API for the above_ground module

The models.hydrology.above_ground module simulates the above-ground hydrological processes for the Virtual Ecosystem. At the moment, this includes rain water interception by the canopy, soil evaporation, and functions related to surface runoff, bypass flow, and river discharge.

TODO change temperatures to Kelvin

Functions:

accumulate_horizontal_flow(drainage_map, ...)

Calculate accumulated above-/belowground horizontal flow for each grid cell.

calculate_bypass_flow(top_soil_moisture, ...)

Calculate preferential bypass flow.

calculate_drainage_map(grid, elevation)

Calculate drainage map based on digital elevation model.

calculate_interception(leaf_area_index, ...)

Estimate canopy interception.

calculate_soil_evaporation(temperature, ...)

Calculate soil evaporation based on classical bulk aerodynamic formulation.

calculate_surface_runoff(...)

Calculate surface runoff, [mm].

convert_mm_flow_to_m3_per_second(...)

Convert river discharge from millimeters to m3/s.

distribute_monthly_rainfall(...[, seed])

Distributes total monthly rainfall over the specified number of days.

find_lowest_neighbour(neighbours, elevation)

Find lowest neighbour for each grid cell from digital elevation model.

find_upstream_cells(lowest_neighbour)

Find all upstream cell IDs for all grid cells.

virtual_ecosystem.models.hydrology.above_ground.accumulate_horizontal_flow(drainage_map: dict[int, list[int]], current_flow: ndarray, previous_accumulated_flow: ndarray) ndarray

Calculate accumulated above-/belowground horizontal flow for each grid cell.

This function takes the accumulated above-/belowground horizontal flow from the previous timestep and adds all (sub-)surface flow of the current time step from upstream cell IDs.

The function currently raises a ValueError if accumulated flow is negative.

Parameters:
  • drainage_map – dict of all upstream IDs for each grid cell

  • current_flow – (sub-)surface flow of the current time step, [mm]

  • previous_accumulated_flow – accumulated flow from previous time step, [mm]

Returns:

accumulated (sub-)surface flow, [mm]

virtual_ecosystem.models.hydrology.above_ground.calculate_bypass_flow(top_soil_moisture: ndarray[Any, dtype[float32]], sat_top_soil_moisture: ndarray[Any, dtype[float32]], available_water: ndarray[Any, dtype[float32]], infiltration_shape_parameter: float) ndarray[Any, dtype[float32]]

Calculate preferential bypass flow.

Bypass flow is here defined as the flow that bypasses the soil matrix and drains directly to the groundwater. During each time step, a fraction of the water that is available for infiltration is added to the groundwater directly (i.e. without first entering the soil matrix). It is assumed that this fraction is a power function of the relative saturation of the superficial and upper soil layers. This results in the following equation (after Van Der Knijff et al. (2010)):

\(D_{pref, gw} = W_{av} * (\frac{w_{1}}{w_{s1}})^{c_{pref}}\)

where \(D_{pref, gw}\) is the amount of preferential flow per time step [mm], \(W_{av}\) is the amount of water that is available for infiltration, and \(c_{pref}\) is an empirical shape parameter. This parameter affects how much of the water available for infiltration goes directly to groundwater via preferential bypass flow; a value of 0 means all surface water goes directly to groundwater, a value of 1 gives a linear relation between soil moisture and bypass flow. The equation returns a preferential flow component that becomes increasingly important as the soil gets wetter.

Parameters:
  • top_soil_moisture – soil moisture of top soil layer, [mm]

  • sat_top_soil_moisture – soil moisture of top soil layer at saturation, [mm]

  • available_water – amount of water available for infiltration, [mm]

  • infiltration_shape_parameter – shape parameter for infiltration

Returns:

preferential bypass flow, [mm]

virtual_ecosystem.models.hydrology.above_ground.calculate_drainage_map(grid: Grid, elevation: ndarray) dict[int, list[int]]

Calculate drainage map based on digital elevation model.

This function finds the lowest neighbour for each grid cell, identifies all upstream IDs and creates a dictionary that provides all upstream cell IDs for each grid cell. This function currently supports only square grids.

Parameters:
  • grid – grid object

  • elevation – elevation, [m]

Returns:

dictionary of cell IDs and their upstream neighbours

TODO move this to core.grid once we decided on common use

virtual_ecosystem.models.hydrology.above_ground.calculate_interception(leaf_area_index: ndarray[Any, dtype[float32]], precipitation: ndarray[Any, dtype[float32]], intercept_parameters: tuple[float, float, float], veg_density_param: float) ndarray[Any, dtype[float32]]

Estimate canopy interception.

This function estimates canopy interception using the following storage-based equation after Aston (1979) and Merriam (1960) as implemented in Van Der Knijff et al. (2010) :

\(Int = S_{max} * [1 - e \frac{(-k*R*\delta t}{S_{max}})]\)

where \(Int\) [mm] is the interception per time step, \(S_{max}\) [mm] is the maximum interception, \(R\) is the rainfall intensity per time step [mm] and the factor \(k\) accounts for the density of the vegetation.

\(S_{max}\) is calculated using an empirical equation (Von Hoyningen-Huene, 1981):

\[ S_{max} = \begin{cases} 0.935 + 0.498 \cdot \text{LAI} - 0.00575 \cdot \text{LAI}^{2}, & \text{LAI} > 0.1 \\ 0, & \text{LAI} \le 0.1, \end{cases} \]

where LAI is the average Leaf Area Index [m2 m-2]. \(k\) is estimated as:

\(k=0.046 * LAI\)

Parameters:
  • leaf_area_index – leaf area index summed over all canopy layers, [m2 m-2]

  • precipitation – precipitation, [mm]

  • intercept_parameter_1 – Parameter in equation that estimates maximum canopy interception capacity

  • intercept_parameter_2 – Parameter in equation that estimates maximum canopy interception capacity

  • intercept_parameter_3 – Parameter in equation that estimates maximum canopy interception capacity

  • veg_density_param – Parameter used to estimate vegetation density for maximum canopy interception capacity estimate

Returns:

interception, [mm]

virtual_ecosystem.models.hydrology.above_ground.calculate_soil_evaporation(temperature: ndarray[Any, dtype[float32]], relative_humidity: ndarray[Any, dtype[float32]], atmospheric_pressure: ndarray[Any, dtype[float32]], soil_moisture: ndarray[Any, dtype[float32]], soil_moisture_residual: float | ndarray[Any, dtype[float32]], soil_moisture_capacity: float | ndarray[Any, dtype[float32]], leaf_area_index: ndarray[Any, dtype[float32]], wind_speed_surface: ndarray[Any, dtype[float32]], celsius_to_kelvin: float, density_air: float | ndarray[Any, dtype[float32]], latent_heat_vapourisation: float | ndarray[Any, dtype[float32]], gas_constant_water_vapour: float, soil_surface_heat_transfer_coefficient: float, extinction_coefficient_global_radiation: float) dict[str, ndarray[Any, dtype[float32]]]

Calculate soil evaporation based on classical bulk aerodynamic formulation.

This function uses the so-called ‘alpha’ method to estimate the evaporative flux (Mahfouf and Noilhan, 1991). We here use the implementation by Barton (1979):

\(\alpha = \frac{1.8 * \Theta}{\Theta + 0.3}\)

\(E_{g} = \frac{\rho_{air}}{R_{a}} * (\alpha * q_{sat}(T_{s}) - q_{g})\)

where \(\Theta\) is the available top soil moisture (relative volumetric water content), \(E_{g}\) is the evaporation flux (W m-2), \(\rho_{air}\) is the density of air (kg m-3), \(R_{a}\) is the aerodynamic resistance (unitless), \(q_{sat}(T_{s})\) (unitless) is the saturated specific humidity, and \(q_{g}\) is the surface specific humidity (unitless).

In a final step, the bare soil evaporation is adjusted to shaded soil evaporation Supit et al. (1994):

\(E_{act} = E_{g} * exp(-\kappa_{gb}*LAI)\)

where \(\kappa_{gb}\) is the extinction coefficient for global radiation, and \(LAI\) is the total leaf area index.

Parameters:
  • temperature – air temperature at reference height, [C]

  • relative_humidity – relative humidity at reference height, []

  • atmospheric_pressure – atmospheric pressure at reference height, [kPa]

  • soil_moisture – Volumetric relative water content, [unitless]

  • soil_moisture_residual – residual soil moisture, [unitless]

  • soil_moisture_capacity – soil moisture capacity, [unitless]

  • wind_speed_surface – wind speed in the bottom air layer, [m s-1]

  • celsius_to_kelvin – factor to convert teperature from Celsius to Kelvin

  • density_air – density if air, [kg m-3]

  • latent_heat_vapourisation – latent heat of vapourisation, [MJ kg-1]

  • gas_constant_water_vapour – gas constant for water vapour, [J kg-1 K-1]

  • soil_surface_heat_transfer_coefficient – heat transfer coefficient between soil and air, [W m-2 K-1]

  • extinction_coefficient_global_radiation – Extinction coefficient for global radiation, [unitless]

Returns:

soil evaporation, [mm] and aerodynamic resistance near the surface

virtual_ecosystem.models.hydrology.above_ground.calculate_surface_runoff(precipitation_surface: ndarray[Any, dtype[float32]], top_soil_moisture: ndarray[Any, dtype[float32]], top_soil_moisture_capacity: ndarray[Any, dtype[float32]]) ndarray[Any, dtype[float32]]

Calculate surface runoff, [mm].

Surface runoff is calculated with a simple bucket model based on Davis et al. (2017): if precipitation exceeds top soil moisture capacity , the excess water is added to runoff and top soil moisture is set to soil moisture capacity value; if the top soil is not saturated, precipitation is added to the current soil moisture level and runoff is set to zero.

Parameters:
  • precipitation_surface – precipitation that reaches surface, [mm]

  • top_soil_moisture – water content of top soil layer, [mm]

  • top_soil_moisture_capacity – soil mositure capacity of top soil layer, [mm]

virtual_ecosystem.models.hydrology.above_ground.convert_mm_flow_to_m3_per_second(river_discharge_mm: ndarray[Any, dtype[float32]], area: int | float, days: int, seconds_to_day: float, meters_to_millimeters: float) ndarray[Any, dtype[float32]]

Convert river discharge from millimeters to m3/s.

Parameters:
  • river_discharge_mm – total river discharge, [mm]

  • area – area of each grid cell, [m2]

  • days – number of days

  • seconds_to_day – second to day conversion factor

  • meters_to_millimeters – factor to convert between millimeters and meters

Returns:

river discharge rate for each grid cell in m3/s

virtual_ecosystem.models.hydrology.above_ground.distribute_monthly_rainfall(total_monthly_rainfall: ndarray[Any, dtype[float32]], num_days: int, seed: int | None = None) ndarray[Any, dtype[float32]]

Distributes total monthly rainfall over the specified number of days.

At the moment, this function allocates each millimeter of monthly rainfall to a randomly selected day. In the future, this allocation could be based on observed rainfall patterns.

Parameters:
  • total_monthly_rainfall – Total monthly rainfall, [mm]

  • num_days – Number of days to distribute the rainfall over

  • seed – seed for random number generator, optional

Returns:

An array containing the daily rainfall amounts, [mm]

virtual_ecosystem.models.hydrology.above_ground.find_lowest_neighbour(neighbours: list[ndarray], elevation: ndarray) list[int]

Find lowest neighbour for each grid cell from digital elevation model.

This function finds the cell IDs of the lowest neighbour for each grid cell. This can be used to determine in which direction surface runoff flows.

Parameters:
  • neighbours – list of neighbour IDs

  • elevation – elevation, [m]

Returns:

list of lowest neighbour IDs

virtual_ecosystem.models.hydrology.above_ground.find_upstream_cells(lowest_neighbour: list[int]) list[list[int]]

Find all upstream cell IDs for all grid cells.

This function identifies all cell IDs that are upstream of each grid cell. This can be used to calculate the water flow that goes though a grid cell.

Parameters:

lowest_neighbour – list of lowest neighbour cell_ids

Returns:

lists of all upstream IDs for each grid cell