GroundwaterDupuitPercolator Component.
@author: G Tucker, D Litwin, K Barnhart
GroundwaterDupuitPercolator
(grid, hydraulic_conductivity=0.001, porosity=0.2, recharge_rate=1e08, regularization_f=0.01, courant_coefficient=0.5, vn_coefficient=0.8, callback_fun=<function GroundwaterDupuitPercolator.<lambda>>)[source]¶Bases: landlab.core.model_component.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 DupuitForcheimer 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 model 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=1E7)
>>> 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:
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:
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:
where \(n\) is the drainable porosity.
An explicit forward in time finite volume method is used to implement a numerical solution. Groundwater flow between neighboring nodes is calculated using the saturated thickness at the upgradient node.
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: 


K
¶hydraulic conductivity at link (m/s)
__init__
(grid, hydraulic_conductivity=0.001, porosity=0.2, recharge_rate=1e08, regularization_f=0.01, courant_coefficient=0.5, vn_coefficient=0.8, callback_fun=<function GroundwaterDupuitPercolator.<lambda>>)[source]¶Parameters: 


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)
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. 

n
¶drainable porosity of the aquifer ()
number_of_substeps
¶The number of substeps used by the run_with_adaptive_time_step_solver method in the latest method call.
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. 

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. 
