9.8. Hydrology Module¶
9.8.1. Introduction¶
The Python PFTools Hydrology module provides standalone functions for common hydrologic calculations.
9.8.2. Usage of Hydrology
¶
First, we’ll show some examples of using the Hydrology class within a ParFlow Python script:
import numpy as np
from parflow import Run
from parflow.tools.hydrology import calculate_surface_storage, calculate_subsurface_storage, \
calculate_water_table_depth, calculate_evapotranspiration, calculate_overland_flow_grid
# Create a Run object from the .pfidb file
run = Run.from_definition('/path/to/pfidb/file')
# Get the DataAccessor object corresponding to the Run object
data = run.data_accessor
# ----------------------------------------------
# Get relevant information from the DataAccessor
# ----------------------------------------------
# Resolution
dx = data.dx
dy = data.dy
# Thickness of each layer, bottom to top
dz = data.dz
# Extent
nx = data.shape[2]
ny = data.shape[1]
nz = data.shape[0]
# ------------------------------------------
# Time-invariant values
# ------------------------------------------
porosity = data.computed_porosity
specific_storage = data.specific_storage
mask = data.mask
et = data.et # shape (nz, ny, nx) - units 1/T.
slopex = data.slope_x # shape (ny, nx)
slopey = data.slope_y # shape (ny, nx)
mannings = data.mannings # scalar value
# ------------------------------------------
# Time-variant values
# ------------------------------------------
# no. of time steps
nt = len(data.times)
# ------------------------------------------
# Initialization
# ------------------------------------------
# Arrays for total values (across all layers), with time as the first axis
subsurface_storage = np.zeros(nt)
surface_storage = np.zeros(nt)
wtd = np.zeros((nt, ny, nx))
et = np.zeros(nt)
overland_flow = np.zeros((nt, ny, nx))
# ------------------------------------------
# Loop through time steps
# i goes from 0 to n_timesteps - 1
# ------------------------------------------
for i in data.times:
pressure = data.pressure
saturation = data.saturation
# Total subsurface storage for this time step is the summation of substorage surface across all x/y/z slices
subsurface_storage[i, ...] = np.sum(
calculate_subsurface_storage(porosity, pressure, saturation, specific_storage, dx, dy, dz, mask=mask),
axis=(0, 1, 2)
)
# Total surface storage for this time step is the summation of substorage surface across all x/y slices
surface_storage[i, ...] = np.sum(
calculate_surface_storage(pressure, dx, dy, mask=mask),
axis=(0, 1)
)
wtd[i, ...] = calculate_water_table_depth(pressure, saturation, dz)
if et is not None:
# Total ET for this time step is the summation of ET values across all x/y/z slices
et[i, ...] = np.sum(
calculate_evapotranspiration(et_flux_values, dx, dy, dz, mask=mask),
axis=(0, 1, 2)
)
overland_flow[i, ...] = calculate_overland_flow_grid(pressure, slopex, slopey, mannings, dx, dy, mask=mask)
data.time += 1
9.8.3. Full API¶
calculate_water_table_depth(pressure, saturation, dz)
Calculate water table depth from the land surface.
- param
pressure
An
nz
byny
bynx
ndarray
of pressure values (bottom layer to top layer)- param
saturation
An
nz
byny
bynx
ndarray
ndarray of saturation values (bottom layer to top layer)- param
dz
An
ndarray
of shape(nz,)
of thickness values (bottom layer to top layer)- return
A
ny
bynx
ndarray
of water table depth values (measured from the top)
- param
calculate_subsurface_storage(porosity, pressure, saturation, specific_storage, dx, dy, dz, mask=None)
Calculate gridded subsurface storage across several layers. For each layer in the subsurface, storage consists of two parts:
Incompressible subsurface storage (
porosity
*saturation
* depth of this layer) *dx
*dy
Compressible subsurface storage (
pressure
*saturation
*specific storage
* depth of this layer) *dx
*dy
- param
porosity
An
nz
byny
bynx
ndarray
of porosity values (bottom layer to top layer)- param
pressure
An
nz
byny
bynx
ndarray
of pressure values (bottom layer to top layer)- param
saturation
An
nz
byny
bynx
ndarray
of saturation values (bottom layer to top layer)- param
specific_storage
An
nz
byny
bynx
ndarray
of specific storage values (bottom layer to top layer)- param
dx
Length of a grid element in the
x
direction- param
dy
Length of a grid element in the
y
direction- param
dz
Thickness of a grid element in the
z
direction (bottom layer to top layer)- param
mask
An
nz
byny
bynx
ndarray
of mask values (bottom layer to top layer). IfNone
, assumed to be annz
byny
bynx
ndarray
of 1s.- return
An
nz
byny
bynx
ndarray
of subsurface storage values, spanning all layers (bottom to top)
calculate_surface_storage(pressure, dx, dy, mask=None)
Calculate gridded surface storage on the top layer. Surface storage is given by: Pressure at the top layer *
dx
*dy
(for pressure values > 0)- param
pressure
An
nz
byny
bynx
ndarray
of pressure values (bottom layer to top layer)- param
dx
Length of a grid element in the x direction
- param
dy
Length of a grid element in the y direction
- param
mask
An
nz
byny
bynx
ndarray
of mask values (bottom layer to top layer). IfNone
, assumed to be annz
byny
bynx
ndarray
of 1s.- return
An
ny
bynx
ndarray
of surface storage values
- param
calculate_evapotranspiration(et, dx, dy, dz, mask=None)
Calculate gridded evapotranspiration across several layers.
- param
et
An
nz
byny
bynx
ndarray
of evapotranspiration flux values with units 1/T (bottom layer to top layer)- param
dx
Length of a grid element in the
x
direction- param
dy
Length of a grid element in the
y
direction- param
dz
Thickness of a grid element in the
z
direction (bottom layer to top layer)- param
mask
An
nz
byny
bynx
ndarray
of mask values (bottom layer to top layer). IfNone
, assumed to be annz
byny
bynx
ndarray
of 1s.- return
An
nz
byny
bynx
ndarray
of evapotranspiration values (units L^3/T), spanning all layers (bottom to top)
- param
calculate_overland_fluxes(pressure, slopex, slopey, mannings, dx, dy, flow_method='OverlandKinematic', epsilon=1e-5, mask=None)
Calculate overland fluxes across grid faces.
- param
pressure
An
nz
byny
bynx
ndarray
of pressure values (bottom layer to top layer)- param
slopex
ny
bynx
- param
slopey
ny
bynx
- param
mannings
a scalar value, or a
ny
bynx
ndarray
- param
dx
Length of a grid element in the
x
direction- param
dy
Length of a grid element in the
y
direction- param
flow_method
Either ‘OverlandFlow’ or ‘OverlandKinematic’. ‘OverlandKinematic’ by default.
- param
epsilon
Minimum slope magnitude for solver. Only applicable if
flow_method='OverlandKinematic'
. This is set using theSolver.OverlandKinematic.Epsilon
key in Parflow.- param
mask
An
nz
byny
bynx
ndarray
of mask values (bottom layer to top layer). IfNone
, assumed to be annz
byny
bynx
ndarray
of 1s.- return
A 2-tuple:
(
qeast
: Any
by(nx+1)
ndarray
of overland flux values,qnorth
: A(ny+1)
bynx
ndarray
of overland flux values)
Numpy array origin is at the top left. The cardinal direction along axis 0 (rows) is North (going down!!). The cardinal direction along axis 1 (columns) is East (going right). qnorth ``(ny+1,nx)`` and qeast ``(ny,nx+1)`` values are to be interpreted as follows. +-------------------------------------> (East) | | qnorth_i,j (outflow if negative) | +-----+------+ | | | | | | | | | qeast_i,j (outflow if negative) |--> v |---> qeast_i,j+1 (outflow if positive) | | | | | Cell i,j | | +-----+------+ | | | | | v | qnorth_i+1,j (outflow if positive) v (North)
- param
calculate_overland_flow_grid(pressure, slopex, slopey, mannings, dx, dy, flow_method='OverlandKinematic', epsilon=1e-5, mask=None)
Calculate overland outflow per grid cell of a domain.
- param
pressure
An
nz
byny
bynx
ndarray
of pressure values (bottom layer to top layer)- param
slopex
ny
bynx
- param
slopey
ny
bynx
- param
mannings
a scalar value, or a
ny
bynx
ndarray
- param
dx
Length of a grid element in the
x
direction- param
dy
Length of a grid element in the
y
direction- param
flow_method
Either ‘OverlandFlow’ or ‘OverlandKinematic’. ‘OverlandKinematic’ by default.
- param
epsilon
Minimum slope magnitude for solver. Only applicable if
kinematic=True
. This is set using theSolver.OverlandKinematic.Epsilon
key in Parflow.- param
mask
An
nz
byny
bynx
ndarray
of mask values (bottom layer to top layer). IfNone
, assumed to be annz
byny
bynx
ndarray
of 1s.- return
An
ny
bynx
ndarray
of overland flow values
- param
calculate_overland_flow(pressure, slopex, slopey, mannings, dx, dy, flow_method='OverlandKinematic', epsilon=1e-5, mask=None)
- param
pressure
An
nz
byny
bynx
ndarray
of pressure values (bottom layer to top layer)- param
slopex
ny
bynx
- param
slopey
ny
bynx
- param
mannings
a scalar value, or a
ny
bynx
ndarray
- param
dx
Length of a grid element in the
x
direction- param
dy
Length of a grid element in the
y
direction- param
flow_method
Either ‘OverlandFlow’ or ‘OverlandKinematic’. ‘OverlandKinematic’ by default.
- param
epsilon
Minimum slope magnitude for solver. Only applicable if
flow_method='OverlandKinematic'
. This is set using theSolver.OverlandKinematic.Epsilon
key in Parflow.- param
mask
An
nz
byny
bynx
ndarray
of mask values (bottom layer to top layer). If None, assumed to be annz
byny
bynx
ndarray
of 1s.- return
A
ny
bynx
ndarray
of overland flow values
- param