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, canopy evaporation, soil evaporation,
and functions related to surface runoff, bypass flow, and river discharge.
TODO change temperatures to Kelvin
Functions:
|
Calculate preferential bypass flow. |
Calculate evaporation of intercepted water from the canopy, [mm]. |
|
|
Calculate drainage map based on digital elevation model. |
|
Estimate canopy interception. |
|
Calculate soil evaporation based on classical bulk aerodynamic formulation. |
Calculate surface runoff, [mm]. |
|
Convert river discharge from millimeters to m3 s-1. |
|
|
Distribute monthly to daily rainfall using stochastic weather generator. |
|
Find lowest neighbour for each grid cell from digital elevation model. |
|
Find all upstream cell IDs for all grid cells. |
|
Calculate canopy potential evaporation rate using Penman-Monteith equation. |
|
Route horizontal flow for each grid cell (instantaneous channel routing). |
- virtual_ecosystem.models.hydrology.above_ground.calculate_bypass_flow(top_soil_moisture: ndarray[tuple[Any, ...], dtype[floating]], sat_top_soil_moisture: ndarray[tuple[Any, ...], dtype[floating]], available_water: ndarray[tuple[Any, ...], dtype[floating]], bypass_flow_coefficient: float) ndarray[tuple[Any, ...], dtype[floating]][source]#
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]
bypass_flow_coefficient – Bypass flow coefficient, dimensionless
- Returns:
preferential bypass flow, [mm]
- virtual_ecosystem.models.hydrology.above_ground.calculate_canopy_evaporation(leaf_area_index: ndarray[tuple[Any, ...], dtype[floating]], interception: ndarray[tuple[Any, ...], dtype[floating]], net_radiation: ndarray[tuple[Any, ...], dtype[floating]], vapour_pressure_deficit: ndarray[tuple[Any, ...], dtype[floating]], air_temperature: ndarray[tuple[Any, ...], dtype[floating]], density_air_kg: ndarray[tuple[Any, ...], dtype[floating]], specific_heat_air: ndarray[tuple[Any, ...], dtype[floating]], aerodynamic_resistance_canopy: ndarray[tuple[Any, ...], dtype[floating]], stomatal_resistance: ndarray[tuple[Any, ...], dtype[floating]], latent_heat_vapourisation: ndarray[tuple[Any, ...], dtype[floating]], psychrometric_constant: ndarray[tuple[Any, ...], dtype[floating]], saturated_pressure_slope_parameters: tuple[float, float, float, float], time_interval: float, extinction_coefficient_global_radiation: float) ndarray[tuple[Any, ...], dtype[floating]][source]#
Calculate evaporation of intercepted water from the canopy, [mm].
This function calculates evaporation of intercepted water from the canopy following the LISFLOOD model Van Der Knijff et al. (2010). The maximum evaporation per time step \(EW_{max}\) [mm] is proportional to the fraction of vegetated area:
\[EW_{max} = EW_{0} [1 - e^{(-\kappa_{gb} LAI)}] \Delta t\]where \(EW_{0}\) is the potential evaporation rate, the dimensionless constant \(\kappa_{gb}\) is the extinction coefficient for global solar radiation. In LISFLOOD, \(\kappa_{gb}\) is given by the product \(0.75 \cdot \kappa_{df}\), where \(\kappa_{df}\) is the extinction coefficient for diffuse visible light: its value is provided as input to the model and it varies between 0.4 and 1.1.
The actual amount of evaporation \(EW_{int}\) [mm] is limited by the amount of water stored on the leaves \(Int_{cum}\):
\[EW_{int} = min(EW_{max} \Delta t, Int_{cum})\]Leaf drainage is not modelled explicitly given the short residence time of water on the leaves compared to the model time step.
- Parameters:
leaf_area_index – Leaf area index, [m m-1]
interception – Interception of water in canopy, [mm]
net_radiation – Net radiation in canopy, [W m-2]
vapour_pressure_deficit – Vapour pressure deficit, [kPa]
air_temperature – Air temperature in canopy, [C]
density_air_kg – Density of air, [kg m-3]
specific_heat_air – Specific heat of air, [kJ kg-1 K-1]
aerodynamic_resistance_canopy – Aerodynamic resistance in canopy, [s m-1]
stomatal_resistance – Stomatal resistance, [s m-1]
latent_heat_vapourisation – Latent heat of vapourisation, [kJ kg-1]
psychrometric_constant – Psychrometric constant, [kPa K-1]
saturated_pressure_slope_parameters – List of parameters to calculate the slope of the saturated vapour pressure curve
time_interval – Time interval, [s]
extinction_coefficient_global_radiation – Extinction coefficient for global radiation
- Returns:
canopy evaporation [mm per time interval]
- virtual_ecosystem.models.hydrology.above_ground.calculate_drainage_map(grid: Grid, elevation: ndarray) dict[int, list[int]][source]#
Calculate drainage map based on digital elevation model.
This function finds the lowest neighbour for each grid cell, identifies all upstream cell 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[tuple[Any, ...], dtype[floating]], precipitation: ndarray[tuple[Any, ...], dtype[floating]], intercept_parameters: tuple[float, float, float], veg_density_param: float) ndarray[tuple[Any, ...], dtype[floating]][source]#
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 [m1 m-1]. \(k\) is estimated as:
\(k=0.046 \cdot LAI\)
- Parameters:
leaf_area_index – Leaf area index for all canopy layers, [m m-1]
precipitation – Precipitation, [mm]
intercept_parameters – Parameters for equation estimating 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[tuple[Any, ...], dtype[floating]], relative_humidity: ndarray[tuple[Any, ...], dtype[floating]], atmospheric_pressure: ndarray[tuple[Any, ...], dtype[floating]], soil_moisture: ndarray[tuple[Any, ...], dtype[floating]], soil_moisture_residual: float | ndarray[tuple[Any, ...], dtype[floating]], soil_moisture_saturation: float | ndarray[tuple[Any, ...], dtype[floating]], leaf_area_index: ndarray[tuple[Any, ...], dtype[floating]], wind_speed_surface: ndarray[tuple[Any, ...], dtype[floating]], density_air: float | ndarray[tuple[Any, ...], dtype[floating]], latent_heat_vapourisation: float | ndarray[tuple[Any, ...], dtype[floating]], gas_constant_water_vapour: float, drag_coefficient_evaporation: float, extinction_coefficient_global_radiation: float, time_interval: float, pyrealm_core_constants: CoreConst) dict[str, ndarray[tuple[Any, ...], dtype[floating]]][source]#
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}=(C_{E} u_{a})^-1\) is the aerodynamic resistance, with \(C_{E}\) the drag coefficient for evaporation and \(u_{a}\) the wind speed near the surface, \(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} e^{(-\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_saturation – Soil moisture saturation, [unitless]
wind_speed_surface – Wind speed in the bottom air layer, [m s-1]
density_air – Density if air, [kg m-3]
latent_heat_vapourisation – Latent heat of vapourisation, [kJ kg-1]
leaf_area_index – Leaf area index [m m-1]
gas_constant_water_vapour – Gas constant for water vapour, [kJ kg-1 K-1]
drag_coefficient_evaporation – Drag coefficient for evaporation, dimensionless
extinction_coefficient_global_radiation – Extinction coefficient for global radiation, [unitless]
time_interval – Time interval, [s]
pyrealm_core_constants – Core constants from pyrealm package
- Returns:
soil evaporation, [mm per time interval], aerodynamic resistance soil [s m-1]
- virtual_ecosystem.models.hydrology.above_ground.calculate_surface_runoff(precipitation_surface: ndarray[tuple[Any, ...], dtype[floating]], top_soil_moisture: ndarray[tuple[Any, ...], dtype[floating]], top_soil_moisture_saturation: ndarray[tuple[Any, ...], dtype[floating]]) ndarray[tuple[Any, ...], dtype[floating]][source]#
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 saturation , the excess water is added to runoff and top soil moisture is set to soil moisture saturation value; if the top soil is not saturated, precipitation is added to the current soil moisture level and runoff is set to zero.
TODO adjust saturation to account for new set of soil layers #535
- Parameters:
precipitation_surface – Precipitation that reaches surface, [mm]
top_soil_moisture – Water content of top soil layer, [mm]
top_soil_moisture_saturation – Soil mositure saturation of top soil layer, [mm]
- virtual_ecosystem.models.hydrology.above_ground.convert_mm_flow_to_m3_per_second(river_discharge_mm: ndarray[tuple[Any, ...], dtype[floating]], area: int | float, days: int, seconds_to_day: float, meters_to_millimeters: float) ndarray[tuple[Any, ...], dtype[floating]][source]#
Convert river discharge from millimeters to m3 s-1.
- 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, [m3 s-1]
- virtual_ecosystem.models.hydrology.above_ground.distribute_monthly_rainfall(total_monthly_rainfall: ndarray[tuple[Any, ...], dtype[floating]], num_days: int, p_wet_wet: float, p_wet_dry: float, shape_parameter: float, scale_parameter: float, seed: int | None = None) ndarray[tuple[Any, ...], dtype[floating]][source]#
Distribute monthly to daily rainfall using stochastic weather generator.
Daily rainfall occurrence is simulated using a first-order Markov chain , and rainfall intensity on wet days is drawn from a Gamma distribution (Katz, 1977).
Wet/dry occurrence model#
Let \(S_{t}\) be the rainfall state on day \(t\):
\[\begin{split}S_t = \begin{cases} 1 & \text{wet day} \\ 0 & \text{dry day} \end{cases}\end{split}\]The probability of rainfall depends on the previous day’s state:
\[P(S_t = 1 \mid S_{t-1} = 1) = p_{ww}\]\[P(S_t = 1 \mid S_{t-1} = 0) = p_{wd}\]Rainfall intensity model#
Rainfall on wet days is sampled from a Gamma distribution:
\[x_i \sim \text{Gamma}(k, \theta)\]where \(k\) is a shape parameter (dimensionless), and \(\theta\) is a scale parameter, (dimensionless).
The sampled intensities are then scaled so that their sum equals the specified monthly rainfall total:
\[r_{i} = \frac{x_i}{\sum_{j=1}^{n} x_j} P\]where \(r_{i}\) is the rainfall on day \(i\) [mm], \(x_{i}\) is the sampled Gamma intensity, and \(P\) is the total monthly rainfall [mm].
- param total_monthly_rainfall:
Total rainfall per month [mm].
- param num_days:
Number of days in the month.
- param p_wet_wet:
Probability a wet day follows a wet day. Typical values are 0.5-0.7 temperate climates; 0.6-0.8 humid climates; 0.3-0.5 arid climates
- param p_wet_dry:
Probability a wet day follows a dry day. Typical values are 0.1-0.3 arid climates; 0.2-0.4 temperate climates; 0.3-0.5 tropical climates
- param shape_parameter:
Shape parameter of the Gamma distribution controlling rainfall variability. Typical values are 0.7-1.0 intense storms / high variability; 1.0-2.0 moderate variability (common default 1.5); 2.0-4.0 more uniform rainfall
- param scale_parameter:
Scale parameter of the Gamma distribution controlling absolute magnitude of rainfall, typically 1.0.
- param seed:
Seed for random number generator (optional).
Returns: Daily rainfall array (len(total_monthly_rainfall), num_days), [mm].
- raises ValueError:
if any input is invalid (negative rainfall, invalid
- raises probabilities, etc.):
- virtual_ecosystem.models.hydrology.above_ground.find_lowest_neighbour(neighbours: list[ndarray], elevation: ndarray) list[int][source]#
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]][source]#
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
- virtual_ecosystem.models.hydrology.above_ground.potential_evaporation_leaf(net_radiation: ndarray[tuple[Any, ...], dtype[floating]], vapour_pressure_deficit: ndarray[tuple[Any, ...], dtype[floating]], air_temperature: ndarray[tuple[Any, ...], dtype[floating]], density_air_kg: ndarray[tuple[Any, ...], dtype[floating]], specific_heat_air: ndarray[tuple[Any, ...], dtype[floating]], aerodynamic_resistance_canopy: ndarray[tuple[Any, ...], dtype[floating]], stomatal_resistance: ndarray[tuple[Any, ...], dtype[floating]], latent_heat_vapourisation: ndarray[tuple[Any, ...], dtype[floating]], psychrometric_constant: ndarray[tuple[Any, ...], dtype[floating]], saturated_pressure_slope_parameters: tuple[float, float, float, float])[source]#
Calculate canopy potential evaporation rate using Penman-Monteith equation.
The potential evaporation rate \(EW_{0}\) is calculated as follows:
\[EW_{0} = \frac{\Delta R_n + \rho_a c_p \frac{D}{r_a}} {\lambda_v \left(\Delta + \gamma \left(1 + \frac{r_s}{r_a}\right)\right)}\]where \(\Delta\) is the slope of the saturation vapour pressure curve, \(R_n\) is the net radiation, \(\rho_a\) is the density of air, \(c_p\) is the specific heat of air, \(D\) is the vapour pressure deficit, \(r_a\) is the aerodynamic resistance of canopy, \(\lambda_v\) is the latent heat of vapourization, \(\gamma\) is the psychrometric constant, and \(r_s\) is the stomatal resistance.
Note that we do NOT include ground heat flux in the consideration of canopy evaporation; TODO discuss where we use instead the energy flux into NPP
- Parameters:
net_radiation – Net radiation at leaf surface, [W m-2]
vapour_pressure_deficit – Vapour pressure deficit, [kPa]
air_temperature – Air temperature, [C]
density_air_kg – Air density, [kg m-3]
specific_heat_air – Specific heat of air, [kJ kg-1 K-1]
aerodynamic_resistance_canopy – Aerodynamic resistance in canopy, [s m-1]
stomatal_resistance – Stomatal resistance, [s m-1]
latent_heat_vapourisation – Latent heat of vapourisation, [kJ kg-1]
psychrometric_constant – Psychrometric constant, [kPa K-1]
saturated_pressure_slope_parameters – List of parameters to calculate the slope of the saturated vapour pressure curve
- Returns:
potential evaporation rate, [kg m-2 s-1]
- virtual_ecosystem.models.hydrology.above_ground.route_horizontal_flow(drainage_map: dict[int, list[int]], surface_runoff: ndarray, subsurface_runoff: ndarray) ndarray[source]#
Route horizontal flow for each grid cell (instantaneous channel routing).
This function calculates the total river discharge at each grid cell for the current timestep by combining:
Local generation: the water generated in the cell itself during the timestep, including surface runoff and subsurface (lateral + baseflow) runoff.
Inflow from upstream cells: contributions from all cells that drain into the current cell, using their local generation from the same timestep.
No flows from previous timesteps are included, avoiding double-counting, and there is also no time delay, so all the water runs through the whole grid in one time step .
- Parameters:
drainage_map – Dict mapping each cell ID -> list of upstream cell IDs
surface_runoff – Surface runoff for this timestep, [mm]
subsurface_runoff – Subsurface runoff for this timestep, [mm]
- Returns:
Total river discharge at each grid cell, [mm]