Solute diffusion from a fracture into a porous rock with an altered zone bordering the fracture is modeled by a system of two diffusion equations (one for the altered zone and another for the intact porous matrix) with different coefficients of effective diffusivity. Since experimental studies of diffusion into rock samples with altered zones indicate that mathematical models of diffusion based on Fick's law do not adequately describe the concentration field in a sample, fractional order diffusion equations are chosen in this study for modeling the anomalous mass transport in the rocks. In the case of significantly higher porosity of the altered zone (e.g., this is typical for carbonates) the effective diffusivity here can be much higher than the effective diffusivity of non-altered rocks. By introducing a small parameter that is the ratio of effective diffusivities in the non-altered and altered regions and applying the technique of perturbations, approximate analytical solutions for concentrations in the altered zone bordering the fracture and in the intact surrounding rocks are obtained. Based on these solutions, different regimes of diffusion into the rocks with different physical properties are modeled and analyzed. It is shown that, using experimentally obtained data, the orders of the fractional derivatives in the differential equations can be readily calibrated for the every specific rock.