landlab.components.groundwater.dupuit_percolator¶

GroundwaterDupuitPercolator Component.

@author: G Tucker, D Litwin, K Barnhart

class GroundwaterDupuitPercolator[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.)

static __new__(cls, *args, **kwds)¶

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.

cite_as = '@article{litwin2020groundwater,\n doi = {10.21105/joss.01935},\n url = {https://doi.org/10.21105/joss.01935},\n year = {2020},\n publisher = {The Open Journal},\n volume = {5},\n number = {46},\n pages = {1935},\n author = {David Litwin and Gregory Tucker and Katherine Barnhart and Ciaran Harman},\n title = {GroundwaterDupuitPercolator: A Landlab component for groundwater flow},\n journal = {Journal of Open Source Software}\n }'¶

property coords¶: Return the coordinates of nodes on grid attached to the component.

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 current_time¶

Current time.

Some components may keep track of the current time. In this case, the current_time attribute is incremented. Otherwise it is set to None.

Return type:: current_time

definitions = (('aquifer__thickness', 'thickness of saturated zone'), ('aquifer_base__elevation', 'elevation of impervious layer'), ('aquifer_base__gradient', 'gradient of the aquifer base in the link direction'), ('average_surface_water__specific_discharge', 'average surface water specific discharge over variable timesteps'), ('groundwater__specific_discharge', 'discharge per width in link dir'), ('groundwater__velocity', 'velocity of groundwater in link direction'), ('hydraulic__gradient', 'gradient of water table in link direction'), ('surface_water__specific_discharge', 'rate of seepage to surface'), ('topographic__elevation', 'Land surface topographic elevation'), ('water_table__elevation', 'elevation of water table'))¶

classmethod from_path(grid, path)¶

Create a component from an input file.

Parameters:

grid (ModelGrid) – A landlab grid.
path (str or file_like) – Path to a parameter file, contents of a parameter file, or a file-like object.

Returns:

A newly-created component.

Return type:

Component

property grid¶: Return the grid attached to the component.

initialize_optional_output_fields()¶

Create fields for a component based on its optional field outputs, if declared in _optional_var_names.

This method will create new fields (without overwrite) for any fields output by the component as optional. New fields are initialized to zero. New fields are created as arrays of floats, unless the component also contains the specifying property _var_type.

initialize_output_fields(values_per_element=None)¶

Create fields for a component based on its input and output var names.

This method will create new fields (without overwrite) for any fields output by, but not supplied to, the component. New fields are initialized to zero. Ignores optional fields. New fields are created as arrays of floats, unless the component specifies the variable type.

Parameters:: values_per_element (int (optional)) – On occasion, it is necessary to create a field that is of size (n_grid_elements, values_per_element) instead of the default size (n_grid_elements,). Use this keyword argument to acomplish this task.

input_var_names = ('aquifer_base__elevation', 'topographic__elevation')¶

property n¶: drainable porosity of the aquifer (-)

name = 'GroundwaterDupuitPercolator'¶

property number_of_substeps¶: The number of substeps used by the run_with_adaptive_time_step_solver method in the latest method call.

optional_var_names = ()¶

output_var_names = ('aquifer__thickness', 'aquifer_base__gradient', 'average_surface_water__specific_discharge', 'groundwater__specific_discharge', 'groundwater__velocity', 'hydraulic__gradient', 'surface_water__specific_discharge', 'water_table__elevation')¶

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 shape¶: Return the grid shape attached to the component, if defined.

unit_agnostic = True¶

units = (('aquifer__thickness', 'm'), ('aquifer_base__elevation', 'm'), ('aquifer_base__gradient', 'm/m'), ('average_surface_water__specific_discharge', 'm/s'), ('groundwater__specific_discharge', 'm2/s'), ('groundwater__velocity', 'm/s'), ('hydraulic__gradient', 'm/m'), ('surface_water__specific_discharge', 'm/s'), ('topographic__elevation', 'm'), ('water_table__elevation', 'm'))¶

classmethod var_definition(name)¶

Get a description of a particular field.

Parameters:: name (str) – A field name.
Returns:: A description of each field.
Return type:: tuple of (name, *description*)

classmethod var_help(name)¶

Print a help message for a particular field.

Parameters:: name (str) – A field name.

classmethod var_loc(name)¶

Location where a particular variable is defined.

Parameters:: name (str) – A field name.
Returns:: The location (‘node’, ‘link’, etc.) where a variable is defined.
Return type:: str

var_mapping = (('aquifer__thickness', 'node'), ('aquifer_base__elevation', 'node'), ('aquifer_base__gradient', 'link'), ('average_surface_water__specific_discharge', 'node'), ('groundwater__specific_discharge', 'link'), ('groundwater__velocity', 'link'), ('hydraulic__gradient', 'link'), ('surface_water__specific_discharge', 'node'), ('topographic__elevation', 'node'), ('water_table__elevation', 'node'))¶

classmethod var_type(name)¶

Returns the dtype of a field (float, int, bool, str…).

Parameters:: name (str) – A field name.
Returns:: The dtype of the field.
Return type:: dtype

classmethod var_units(name)¶

Get the units of a particular field.

Parameters:: name (str) – A field name.
Returns:: Units for the given field.
Return type:: str

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.

get_link_hydraulic_conductivity(grid, K)[source]¶

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]]