When water from the surface, e.g. from a lake, fl ows down through porous carbonate rocks, through a region with high hydraulic conductivity and encounters the water table of a phreatic aquifer, both waters mix by diffusion along their boundary. In a carbonate aquifer, where both surface and phreatic waters are saturated with respect to calcite, mixing corrosion causes renewed dissolution capacity Δceq of the carbonate rock in the diffusion-mixing zone, extending from the boundary separating the phreatic water from the surface water encountering it. A numerical model is presented from which the initial change of porosity in such a diffusion-mixing zone is obtained. The initial change of porosity can be calculated from the local distribution of the mixing ratio, and the second derivative of Δceq with respect to m. m(x,y)is the spatial distribution of the mixing ratio m= Vsurf /(Vsuf+ Vprh), where the V’s assign the corresponding volumes of surface and phreatic water. The second derivative has been calculated for three geochemical scenarios with differing CO2-concentrations of surface and phreatic water by the use of PHREEQC-2. The spatial distribution m(x,y) is obtained by using MODFLOW and MT3DMS in a modelling domain with constant hydraulic conductivity for various fl ow velocities of the phreatic aquifer. The time scale of cave evolution is estimated from the results. Passages of dimensions of about one metre in width and several 10 cm in height, extending in length along the boundary line, where surface and phreatic waters meet, can be created in time scales of 10 000 years. These caves are horizontal with blind ending passages and closely resemble the isolated caves observed in Central West Florida.