FLAC3D Theory and Background • Fluid-Mechanical Interaction

Spreading of a Groundwater Mound

Note

To view this project in FLAC3D, use the menu command Help ► Examples…. The project’s main data files are shown at the end of this example.

This problem studies the transient evolution of a groundwater mound within a porous medium. The mound spreads out and flows along an impervious base under the influence of gravity. It is assumed that the fluid is incompressible, and that the water is contained initially in a cylindrical region with radius \(r_{0}\) and height \(h_{0}\). The water saturation within the mound is equal to one, and Darcy’s law is applicable. The mound elevation is compared to an analytic solution as the elevation evolves with time.

Kochina et al. (1983) have derived the solution for the height, \(h\), of the mound. The solution, as given by Barenblatt (1987), assumes a hydrostatic pore-pressure distribution within the mound. In the case of zero residual saturation, it may be expressed in the form

(1)\[\hat{h}={{1}\over{8\sqrt{\hat{t}}}} \ [4-{{\hat{r}^{2}}\over{\sqrt{\hat{t}}}}]; \ \ \ \hat{r}\leq 2 \root 4 \of {\hat{t}}\]

where \(\hat{h} = h/h_{0}\), \(\hat{r} = r/r_{0}\), \(\hat{t} = t/t_{c}\), and the characteristic time is given as \(t_{c} = r_{0}^{2}/(2\kappa h_{0})\) with \(\kappa = k\rho _{w}g/(2n)\) (\(k\) is the mobility coefficient, \(\rho _{w}\) is water density, \(g\) is gravity and \(n\) is porosity).

This solution applies to long time scales, when the influence of the details of the initial mound geometry have disappeared.

The results are presented in dimensionless form: the scaled geometrical parameters \(r_{0}\) = 1, \(h_{0}\) = 1 are used, and the scaled water properties \(k\) = 0.5 × 10^-4, \(\rho _{w}\) = 10³, and \(n\) = 0.5 are prescribed in the numerical simulation. To model incompressible flow, the bulk modulus of the fluid, \(K_{f}\), is given a value that is large compared to the pore-pressure variations in the simulation (\(K_{f}\) = 2 × 10⁵); the value \(g\) = 10 is used for gravity. The FLAC3D grid corresponds to a quarter cylinder and contains 2000 zones (see Figure 1). The radius of the model is 2, and its height is 1 unit. The initial saturation is 1 within the mound (radius = 1 unit, height = 1 unit), and zero outside. The initial pore-pressure distribution within the mound is hydrostatic. All boundaries are impermeable by default. As time goes on, the mound spreads out under its own weight. The simulation is conducted for a total dimensionless time value of \(\hat{t}\) = 0.85, with intermediate results at \(\hat{t}\) = 0.35, 0.45 and 0.65.

Figure 1: FLAC3D grid and initial state of saturated column.

Saturation contours are sketched in Figure 2 to Figure 5. The analytic prediction for the mound height is calculated by a set of FISH functions and plotted for comparison in the figures (bold line).

Figure 2 corresponds to \(\hat{t}\) = 0.35; at that time, the initial shape of the mound still persists, and a comparison with the analytic solution is probably not yet appropriate. For larger times (see Figure 3 to Figure 5), the spreading of the groundwater mound described by Equation (1) is captured by the numerical solution with reasonable accuracy. The numerical estimate lags behind the analytical prediction; the discrepancy may be explained by the occurrence of residual saturation in the numerical solution and by the coarse discretization used in the simulation.

Figure 2: Saturation contours and analytical mound elevation at \(\hat{t}\) = 0.35.

Figure 3: Saturation contours and analytical mound elevation at \(\hat{t}\) = 0.45.

Figure 4: Saturation contours and analytical mound elevation at \(\hat{t}\) = 0.65.

Figure 5: Saturation contours and analytical mound elevation at \(\hat{t}\) = 0.85.

The sketch of flow vectors and head contours in Figure 6 corresponds to \(\hat{t}\) = 0.85. (The lack of smoothness in the contour plot is caused by the head jump across the phreatic surface.) Two regions with water flowing predominantly downward in the core of the mound and outward in its periphery can be seen in the figure. (The vertical bold line is drawn at the location where the time derivative of the analytical mound elevation vanishes: \(\partial \hat{h}/\partial t\) = 0.)

Figure 6: Head contours and analytical mound elevation at \(\hat{t}\) = 0.85

References

Barenblatt, G. Dimensional Analysis. Gordon and Breach Science Publishers (1987).

Kochina, I., N. Mikhailov and M. Filinov. “Groundwater Mound Damping,” Int. J. Engng. Sci., 21, 413-421 (1983).

Data File

GroundwaterMound.dat

model new
; --- geometry ---
zone create cylinder size 30 15 15 point 1 (0,2,0) point 2 (0,0,1) ...
                                   point 3 (2,0,0)
zone group 'mound' range cylinder ...
    end-1 (0,0,-1) end-2 (0,0,2) radius 1.01
; --- fluid flow model ---
model configure fluid-flow
zone fluid unsaturated cutoff 0
zone fluid property mobility-coefficient 0.5e-4 biot-modulus 1e5 ...
                    fluid-density 1e3 
zone gridpoint initialize saturation 0.0
zone gridpoint initialize saturation 1.0 range group 'mound'
zone gridpoint initialize pore-pressure 1e4 gradient 0 0 -1e4 range group 'mound'
; --- settings ---
model gravity 10
model save 'mound0' 
; --- test ---
zone fluid implicit servo on 
model solve-fluid time 0.35 
model title "Spreading of a groundwater mound (t = 0.35)"
model save 'mound1'
model solve-fluid time 0.10
model title "Spreading of a groundwater mound (t = 0.45)"
model save 'mound2'
model solve-fluid time 0.20
model title "Spreading of a groundwater mound (t = 0.65)"
model save 'mound3'
model solve-fluid time 0.20
model title "Spreading of a groundwater mound (t = 0.85)"
model save 'mound4'