GroundwaterDupuitPercolator: model flow in a shallow unconfined aquifer using the Dupuit-Forcheimer approximation#

GroundwaterDupuitPercolator Component.

@author: G Tucker, D Litwin, K Barnhart

class GroundwaterDupuitPercolator(*args, **kwds)[source]#

Bases: Component

Simulate groundwater flow in a shallow unconfined aquifer.

The GroundwaterDupuitPercolator solves the Boussinesq equation for flow in an unconfined aquifer over an impermeable aquifer base and calculates groundwater return flow to the surface. This method uses the Dupuit-Forcheimer approximation. This means that the model assumes the aquifer is laterally extensive in comparison to its thickness, such that the vertical component of flow is negligible. It also assumes that the capillary fringe is small, such that the water table can be modeled as a free surface. Please consider the applicability of these assumptions when using this model. For more details, see component documentation here.

Examples

Import the grid class and component

>>> from landlab import RasterModelGrid
>>> from landlab.components import GroundwaterDupuitPercolator

Initialize the grid and component

>>> grid = RasterModelGrid((10, 10), xy_spacing=10.0)
>>> elev = grid.add_zeros("topographic__elevation", at="node")
>>> abe = grid.add_zeros("aquifer_base__elevation", at="node")
>>> elev[:] = 5.0
>>> gdp = GroundwaterDupuitPercolator(grid)

Run component forward.

>>> dt = 1e4
>>> for i in range(100):
...     gdp.run_one_step(dt)
...

When the model generates groundwater return flow, the surface water flux out of the domain can be calculated only after a FlowAccumulator is run. Below is a more complex example that demonstrates this case.

>>> from landlab.components import FlowAccumulator

Set boundary conditions and initialize grid

>>> grid = RasterModelGrid((5, 41), xy_spacing=10.0)
>>> grid.set_closed_boundaries_at_grid_edges(True, True, False, True)

Make a sloping, 3 m thick aquifer, initially fully saturated

>>> elev = grid.add_zeros("topographic__elevation", at="node")
>>> elev[:] = grid.x_of_node / 100 + 3
>>> base = grid.add_zeros("aquifer_base__elevation", at="node")
>>> base[:] = grid.x_of_node / 100
>>> wt = grid.add_zeros("water_table__elevation", at="node")
>>> wt[:] = grid.x_of_node / 100 + 3

Initialize components

>>> gdp = GroundwaterDupuitPercolator(grid, recharge_rate=1e-7)
>>> fa = FlowAccumulator(grid, runoff_rate="surface_water__specific_discharge")

Advance timestep. Default units are meters and seconds, though the component is unit agnostic.

>>> dt = 1e3
>>> for i in range(1000):
...     gdp.run_one_step(dt)
...

Calculate surface water flux out of domain

>>> fa.run_one_step()
>>> np.testing.assert_almost_equal(gdp.calc_sw_flux_out(), 0.0005077)

Notes

Below is a summary of the theory and numerical implementation of the GroundwaterDupuitPercolator. A complete description can be found here.

Groundwater discharge per unit length, \(q\), is calculated as:

\[q = -K_{sat} h \big( \nabla z \big) \cos^2 (\alpha)\]

where \(K_{sat}\) is the saturated hydraulic conductivity, \(h\) is the aquifer thickness, and \(\alpha\) is the slope angle of the aquifer base.

Surface water discharge per unit area, \(q_s\), is calculated as:

\[q_s = \mathcal{G}_r \bigg( \frac{h}{d} \bigg) \mathcal{R} \big(-\nabla \cdot q + f \big)\]

where \(\mathcal{G}_r\) is a smoothed step function, \(\mathcal{R}\) is the ramp function, \(d\) is the regolith thickness, and \(f\) is the recharge rate.

The evolution of aquifer thickness is then given by:

\[n \frac{\partial h}{\partial t} = f - q_s - \nabla \cdot q\]

where \(n\) is the drainable porosity.

A semi-analytical approach is used to update aquifer thickness \(h\), in which the differential equation above is linearized by assuming \(\nabla \cdot q\) varies minimally over the duration of a timestep. In this approach, \(q\) is reported at the beginning of the timestep, while \(q_s\) is reported at the end of the timestep.

References

Required Software Citation(s) Specific to this Component

Litwin, D. G., Tucker, G.E., Barnhart, K. R., Harman, C. J. (2020). GroundwaterDupuitPercolator: A Landlab component for groundwater flow. Journal of Open Source Software, 5(46), 1935, https://doi.org/10.21105/joss.01935.

Additional References

Marçais, J., de Dreuzy, J. R. & Erhel, J. Dynamic coupling of subsurface and seepage flows solved within a regularized partition formulation. Advances in Water Resources 109, 94–105 (2017).

Childs, E. C. Drainage of Groundwater Resting on a Sloping Bed. Water Resources Research 7, 1256–1263 (1971).

Parameters:

  • grid (ModelGrid) – Landlab ModelGrid object

  • hydraulic_conductivity (float, field name, array of float or function.) – the aquifer saturated hydraulic conductivity, m/s. If function is given, it should take a landlab ModelGrid and return an array of floats at link. This may be used if the lateral hydraulic conductivity is not vertically homogenous and the effective hydraulic conductivity needs to be modified based upon on the position of the water table. See component tests for example. Default = 0.001 m/s

  • porosity (float, field name or array of float) – the drainable porosity of the aquifer [-] Default = 0.2

  • recharge_rate (float, field name, or array of float) – Rate of recharge, m/s Default = 1.0e-8 m/s

  • regularization_f (float) – factor controlling the smoothness of the transition between surface and subsurface flow Default = 0.01

  • courant_coefficient (float (-)) – The muliplying factor on the condition that the timestep is smaller than the minimum link length over groundwater flow velocity. This parameter is only used with run_with_adaptive_time_step_solver and must be greater than zero. Default = 0.5

  • vn_coefficient (float (-)) – The multiplying factor C for the condition \(dt >= C*dx^2/(4D)\), where \(D = Kh/n\) is the diffusivity of the Boussinesq equation. This arises from a von Neumann stability analysis of the Boussinesq equation when the hydraulic gradient is small. This parameter is only used with run_with_adaptive_time_step_solver and must be greater than zero. Default = 0.8

  • callback_fun (function(grid, recharge_rate, substep_dt, **kwargs)) – Optional function that will be executed at the end of each sub-timestep in the run_with_adaptive_time_step_solver method. Intended purpose is to write output not otherwise visible outside of the method call. The function should have three required arguments: grid: the ModelGrid instance used by GroundwaterDupuitPercolator recharge_rate: an array at node that is the specified recharge rate substep_dt: the length of the current substep determined internally by run_with_adaptive_time_step_solver to meet stability criteria.

  • callback_kwds (any additional keyword arguments for the provided callback_fun.) –

property K#

hydraulic conductivity at link (m/s)

__init__(grid, hydraulic_conductivity=0.001, porosity=0.2, recharge_rate=1e-08, regularization_f=0.01, courant_coefficient=0.5, vn_coefficient=0.8, callback_fun=<function GroundwaterDupuitPercolator.<lambda>>, **callback_kwds)[source]#
Parameters:

  • grid (ModelGrid) – Landlab ModelGrid object

  • hydraulic_conductivity (float, field name, array of float or function.) – the aquifer saturated hydraulic conductivity, m/s. If function is given, it should take a landlab ModelGrid and return an array of floats at link. This may be used if the lateral hydraulic conductivity is not vertically homogenous and the effective hydraulic conductivity needs to be modified based upon on the position of the water table. See component tests for example. Default = 0.001 m/s

  • porosity (float, field name or array of float) – the drainable porosity of the aquifer [-] Default = 0.2

  • recharge_rate (float, field name, or array of float) – Rate of recharge, m/s Default = 1.0e-8 m/s

  • regularization_f (float) – factor controlling the smoothness of the transition between surface and subsurface flow Default = 0.01

  • courant_coefficient (float (-)) – The muliplying factor on the condition that the timestep is smaller than the minimum link length over groundwater flow velocity. This parameter is only used with run_with_adaptive_time_step_solver and must be greater than zero. Default = 0.5

  • vn_coefficient (float (-)) – The multiplying factor C for the condition \(dt >= C*dx^2/(4D)\), where \(D = Kh/n\) is the diffusivity of the Boussinesq equation. This arises from a von Neumann stability analysis of the Boussinesq equation when the hydraulic gradient is small. This parameter is only used with run_with_adaptive_time_step_solver and must be greater than zero. Default = 0.8

  • callback_fun (function(grid, recharge_rate, substep_dt, **kwargs)) – Optional function that will be executed at the end of each sub-timestep in the run_with_adaptive_time_step_solver method. Intended purpose is to write output not otherwise visible outside of the method call. The function should have three required arguments: grid: the ModelGrid instance used by GroundwaterDupuitPercolator recharge_rate: an array at node that is the specified recharge rate substep_dt: the length of the current substep determined internally by run_with_adaptive_time_step_solver to meet stability criteria.

  • callback_kwds (any additional keyword arguments for the provided callback_fun.) –

calc_gw_flux_at_node()[source]#

Calculate the sum of the groundwater flux leaving a node.

(m2/s)

calc_gw_flux_out()[source]#

Groundwater flux through open boundaries may be positive (out of the domain) or negative (into the domain).

This function determines the correct sign for specific discharge based upon this convention, and sums the flux across the boundary faces. (m3/s)

calc_recharge_flux_in()[source]#

Calculate flux into the domain from recharge.

Includes recharge that may immediately become saturation excess overland flow. (m3/s)

calc_sw_flux_out()[source]#

Surface water flux out of the domain through seepage and saturation excess.

Note that model does not allow for reinfiltration. (m3/s)

calc_total_storage()[source]#

calculate the current water storage in the aquifer (m3)

property callback_fun#

callback function for adaptive timestep solver

Parameters:

callback_fun (function(grid, recharge_rate, substep_dt, **callback_kwds)) – Optional function that will be executed at the end of each sub-timestep in the run_with_adaptive_time_step_solver method. Intended purpose is to write output not otherwise visible outside of the method call. The function should have three required arguments: grid: the ModelGrid instance used by GroundwaterDupuitPercolator recharge_rate: an array at node that is the specified recharge rate substep_dt: the length of the current substep determined internally by run_with_adaptive_time_step_solver to meet stability criteria.

property courant_coefficient#

Courant coefficient for adaptive time step.

Parameters:

courant_coefficient (float (-)) – The muliplying factor on the condition that the timestep is smaller than the minimum link length over groundwater flow velocity. This parameter is only used with run_with_adaptive_time_step_solver and must be greater than zero.

property n#

drainable porosity of the aquifer (-)

property number_of_substeps#

The number of substeps used by the run_with_adaptive_time_step_solver method in the latest method call.

property recharge#

recharge rate (m/s)

run_one_step(dt)[source]#

Advance component by one time step of size dt.

Parameters:

dt (float) – The imposed timestep.

run_with_adaptive_time_step_solver(dt)[source]#

Advance component by one time step of size dt, subdividing the timestep into substeps as necessary to meet stability conditions. Note this method returns the fluxes at the last substep, but also returns a new field, average_surface_water__specific_discharge, that is averaged over all subtimesteps. To return state during substeps, provide a callback_fun.

Parameters:

dt (float) – The imposed timestep.

property vn_coefficient#

Coefficient for the diffusive timestep condition in the adaptive timestep solver.

Parameters:

vn_coefficient (float (-)) – The multiplying factor C for the condition dt >= C*dx^2/(4D), where D = Kh/n is the diffusivity of the Boussinesq equation. This arises from a von Neumann stability analysis of the Boussinesq equation when the hydraulic gradient is small. This parameter is only used with run_with_adaptive_time_step_solver and must be greater than zero.

Returns array of hydraulic conductivity on links, allowing for aquifers with laterally anisotropic hydraulic conductivity.

Parameters:

K ((2x2) array of floats (m/s)) – The hydraulic conductivity tensor: [[Kxx, Kxy],[Kyx,Kyy]]