8.3.8. Hydrology Module
8.3.8.1. Introduction
The Python PFTools Hydrology module provides standalone functions for common hydrologic calculations.
8.3.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
8.3.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