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

  1. calculate_water_table_depth(pressure, saturation, dz)

    Calculate water table depth from the land surface.

    param pressure

    An nz by ny by nx ndarray of pressure values (bottom layer to top layer)

    param saturation

    An nz by ny by nx 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 by nx ndarray of water table depth values (measured from the top)

  2. 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:

    1. Incompressible subsurface storage (porosity * saturation * depth of this layer) * dx * dy

    2. Compressible subsurface storage (pressure * saturation * specific storage * depth of this layer) * dx * dy

    param porosity

    An nz by ny by nx ndarray of porosity values (bottom layer to top layer)

    param pressure

    An nz by ny by nx ndarray of pressure values (bottom layer to top layer)

    param saturation

    An nz by ny by nx ndarray of saturation values (bottom layer to top layer)

    param specific_storage

    An nz by ny by nx 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 by ny by nx ndarray of mask values (bottom layer to top layer). If None, assumed to be an nz by ny by nx ndarray of 1s.

    return

    An nz by ny by nx ndarray of subsurface storage values, spanning all layers (bottom to top)

  3. 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 by ny by nx 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 by ny by nx ndarray of mask values (bottom layer to top layer). If None, assumed to be an nz by ny by nx ndarray of 1s.

    return

    An ny by nx ndarray of surface storage values

  4. calculate_evapotranspiration(et, dx, dy, dz, mask=None)

    Calculate gridded evapotranspiration across several layers.

    param et

    An nz by ny by nx 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 by ny by nx ndarray of mask values (bottom layer to top layer). If None, assumed to be an nz by ny by nx ndarray of 1s.

    return

    An nz by ny by nx ndarray of evapotranspiration values (units L^3/T), spanning all layers (bottom to top)

  5. 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 by ny by nx ndarray of pressure values (bottom layer to top layer)

    param slopex

    ny by nx

    param slopey

    ny by nx

    param mannings

    a scalar value, or a ny by nx 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 the Solver.OverlandKinematic.Epsilon key in Parflow.

    param mask

    An nz by ny by nx ndarray of mask values (bottom layer to top layer). If None, assumed to be an nz by ny by nx ndarray of 1s.

    return

    A 2-tuple:

    (qeast: A ny by (nx+1) ndarray of overland flux values,

    qnorth: A (ny+1) by nx 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)
  6. 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 by ny by nx ndarray of pressure values (bottom layer to top layer)

    param slopex

    ny by nx

    param slopey

    ny by nx

    param mannings

    a scalar value, or a ny by nx 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 the Solver.OverlandKinematic.Epsilon key in Parflow.

    param mask

    An nz by ny by nx ndarray of mask values (bottom layer to top layer). If None, assumed to be an nz by ny by nx ndarray of 1s.

    return

    An ny by nx ndarray of overland flow values

  7. calculate_overland_flow(pressure, slopex, slopey, mannings, dx, dy, flow_method='OverlandKinematic', epsilon=1e-5, mask=None)

    param pressure

    An nz by ny by nx ndarray of pressure values (bottom layer to top layer)

    param slopex

    ny by nx

    param slopey

    ny by nx

    param mannings

    a scalar value, or a ny by nx 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 the Solver.OverlandKinematic.Epsilon key in Parflow.

    param mask

    An nz by ny by nx ndarray of mask values (bottom layer to top layer). If None, assumed to be an nz by ny by nx ndarray of 1s.

    return

    A ny by nx ndarray of overland flow values