API for the energy_balance module

The models.abiotic.energy_balance module calculates the energy balance for the Virtual Ecosystem. Given that the time increments of the model are an hour or longer, we can assume that below-canopy heat and vapour exchange attain steady state and heat storage in the canopy does not need to be simulated explicitly. (For application where very fine-temporal resolution data might be needed, heat and vapour exchange must be modelled as transient processes, and heat storage by the canopy, and the exchange of heat between different layers of the canopy, must be considered explicitly, see Maclean and Klinges (2021). This is currently not implemented.)

Under steady-state, the balance equation \(\frac{dQ}{dt}\) for the leaves in each canopy layer is as follows (after Maclean and Klinges (2021)):

\[\frac{dQ}{dt} = R_{abs} - R_{em} - H - \lambda E - PP = R_{abs} - \epsilon_{l} \sigma T_{l}^{4} - \frac{\rho_{a} c_p}{r_a}(T_{l} - T_{a}) - \lambda g_{v} \frac {e_{l} - e_{a}}{p_{a}} - PP = 0\]

where \(R_{abs}\) is absorbed shortwave and longwave radiation, \(R_{em}\) emitted radiation, \(H\) the sensible heat flux, \(\lambda E\) the latent heat flux, \(\epsilon_{l}\) the emissivity of the leaf, \(\sigma\) the Stefan-Boltzmann constant, \(T_{l}\) the absolute temperature of the leaf, \(T_{a}\) the absolute temperature of the air surrounding the leaf, \(\lambda\) the latent heat of vapourisation of water, \(e_{l}\) the effective vapour pressure of the leaf, \(e_{a}\) the vapour pressure of air and \(p_{a}\) atmospheric pressure. \(\rho_a\) is the density of air, \(c_{p}\) is the specific heat capacity of air at constant pressure, \(r_{a}\) is the aerodynamic resistance of the surface (leaf or soil), \(g_{v}\) represents the conductivity for vapour loss from the leaves as a function of the stomatal conductivity, \(PP\) stands for primary productivity.

A challenge in solving this equation is the dependency of latent heat and emitted radiation on leaf temperature. We use a secant method to iteratively solve the energy balance for the canopy temperature. The air temperature is updated based on the sensible heat flux from the canopy and soil in equilibrium, and vertical mixing of air between layers.

Atmospheric humidity is also mixed vertically between atmospheric layers. Advection at the top of the canopy is currently not considered as we don’t have have horizontal exchange between grid cells and air above canopy values would be unrealistic.

Functions:

calculate_absorbed_longwave_radiation(...)	Calculate absorbed longwave radiation per layer using Beer-Lambert attenuation.
calculate_energy_balance_residual(...)	Calculate energy balance residual for canopy.
calculate_latent_heat_flux(...)	Calculate latent heat flux from evapotranspiration.
calculate_longwave_emission(temperature, ...)	Calculate longwave emission using the Stefan Boltzmann law.
calculate_sensible_heat_flux(density_air, ...)	Calculate sensible heat flux.
calculate_total_absorbed_shortwave_radiation(...)	Compute total absorbed shortwave radiation contributing to leaf energy balance.
initialise_canopy_and_soil_fluxes(...)	Initialise canopy temperature and energy fluxes.
make_canopy_residual(state, static, ...)	Creates a residual function for canopy temperature to be used in root finding.
secant_solve_cells_layers(residual_function, ...)	Vectorised secant solver for independent (cell, layer) root problems.
update_canopy_air_temperature(...)	Update air temperature surrounding canopy in steady state.
update_humidity_vpd(...)	Update atmospheric humidity and vapour pressure deficit.
update_soil_temperature(ground_heat_flux, ...)	Update soil temperature using heat diffusion.
update_specific_humidity(evapotranspiration, ...)	Update specific humidity from evapotranspiration and soil evaporation.
update_surface_air_temperature(...)	Update surface air temperature in equilibrium with soil and canopy fluxes.

virtual_ecosystem.models.abiotic.energy_balance.calculate_absorbed_longwave_radiation(downward_longwave: ndarray[tuple[Any, ...], dtype[floating]], leaf_area_index: ndarray[tuple[Any, ...], dtype[floating]], leaf_emissivity: float, soil_emissivity: float, extinction_coefficient_lw: float, surface_index: int, topsoil_index: int) → ndarray[tuple[Any, ...], dtype[floating]]

Calculate absorbed longwave radiation per layer using Beer-Lambert attenuation.

Each canopy layer absorbs downward longwave attenuated from above. The surface layer receives the remainder after full canopy attenuation. The topsoil receives what the surface layer transmitted.

Upward longwave from the soil is NOT included here — it is already accounted for in calculate_soil_fluxes via longwave_emission_soil, which drives the ground heat flux and soil sensible heat flux to the surface air layer. Including it here would double-count the soil longwave emission.

Parameters:

• downward_longwave – Downward longwave at top of canopy [W m-2], shape (cells,)

• leaf_area_index – LAI per layer, shape (layers, cells), NaN for empty layers

• leaf_emissivity – Leaf emissivity, dimensionless

• soil_emissivity – Soil emissivity, dimensionless

• extinction_coefficient_lw – Longwave extinction coefficient, typically ~0.5

• surface_index – Row index of surface layer

• topsoil_index – Row index of topsoil layer

Returns:

Absorbed longwave radiation per layer, shape (layers, cells) [W m-2].

virtual_ecosystem.models.abiotic.energy_balance.calculate_energy_balance_residual(canopy_temperature_initial: ndarray[tuple[Any, ...], dtype[floating]], air_temperature: ndarray[tuple[Any, ...], dtype[floating]], evapotranspiration: ndarray[tuple[Any, ...], dtype[floating]], absorbed_shortwave_radiation: ndarray[tuple[Any, ...], dtype[floating]], absorbed_longwave_radiation: ndarray[tuple[Any, ...], dtype[floating]], specific_heat_air: ndarray[tuple[Any, ...], dtype[floating]], density_air: ndarray[tuple[Any, ...], dtype[floating]], aerodynamic_resistance: ndarray[tuple[Any, ...], dtype[floating]], latent_heat_vapourisation: ndarray[tuple[Any, ...], dtype[floating]], surface_index: int, leaf_emissivity: float, stefan_boltzmann_constant: float, zero_Celsius: float, seconds_to_hour: float, return_fluxes: bool) → ndarray[tuple[Any, ...], dtype[floating]] | dict[str, ndarray[tuple[Any, ...], dtype[floating]]]

Calculate energy balance residual for canopy.

The energy balance residual (\(\frac{dQ}{dt}\)) for the canopy is given by:

\[\frac{dQ}{dt} = R_{abs} - \epsilon_{l} \sigma T_{l}^{4} - H - \lambda E - PP\]

Where \(R_abs\) is the absorbed shortwave and longwave radiation by the canopy, \(\epsilon_{l}\) is the leaf emissivity, \(\sigma\) is the Stefan-Boltzmann constant, \(T_{l}\) is the leaf temperature, \(H\) is the sensible heat flux from the canopy, \(\lambda E\) is the latent heat flux from the canopy, \(PP\) is a fraction of the absorbed light is used in photosynthesis (PAR).

Parameters:

• canopy_temperature_initial – Initial leaf temperature for all canopy layers, [C]

• air_temperature – Initial air temperature in canopy layers, [C]

• evapotranspiration – Evapotranspiration, [mm]

• absorbed_shortwave_radiation – Absorbed shortwave radiation for all canopy layers, [W m-2]

• absorbed_longwave_radiation – Absorbed longwave radiation for all canopy layers, [W m-2]

• specific_heat_air – Specific heat capacity of air, [J kg-1 K-1]

• density_air – Density of air, [kg m-3]

• aerodynamic_resistance – Aerodynamic resistance of canopy, [s m-1]

• latent_heat_vapourisation – Latent heat of vapourisation, [J kg-1]

• surface_index – Row index of surface layer

• leaf_emissivity – Leaf emissivity, dimensionless

• stefan_boltzmann_constant – Stefan Boltzmann constant, [W m-2 K-4]

• zero_Celsius – Factor to convert between Celsius and Kelvin

• seconds_to_hour – Factor to convert between hours and seconds

• return_fluxes – Flag to indicate if all components of the energy balance should be returned. This is false for the newton approach to solve for canopy temperature, but true to create the outputs in a second call afterwards.

Returns:

full energy balance or energy balance residual, [W m-2]

virtual_ecosystem.models.abiotic.energy_balance.calculate_latent_heat_flux(evapotranspiration: ndarray[tuple[Any, ...], dtype[floating]], latent_heat_vapourisation: ndarray[tuple[Any, ...], dtype[floating]], time_interval: float) → ndarray[tuple[Any, ...], dtype[floating]]

Calculate latent heat flux from evapotranspiration.

Parameters:

• evapotranspiration – Evapotranspiration per unit area, [kg m-2] (1 kg m-2 of water = 1 mm of water)

• latent_heat_vapourisation – Latent heat of vaporisation of water, [J kg-1]

• time_interval – Time interval over which flux is computed, [s]

Returns:

Latent heat flux, [W m-2]

virtual_ecosystem.models.abiotic.energy_balance.calculate_longwave_emission(temperature: ndarray[tuple[Any, ...], dtype[floating]], emissivity: float | ndarray[tuple[Any, ...], dtype[floating]], stefan_boltzmann: float) → ndarray[tuple[Any, ...], dtype[floating]]

Calculate longwave emission using the Stefan Boltzmann law.

According to the Stefan Boltzmann law, the amount of radiation emitted per unit time from the area of a black body at absolute temperature is directly proportional to the fourth power of the temperature. Emissivity (which is equal to absorptive power) lies between 0 to 1.

Parameters:

• temperature – Temperature, [K]

• emissivity – Emissivity, dimensionless

• stefan_boltzmann – Stefan Boltzmann constant, [W m-2 K-4]

Returns:

Longwave emission, [W m-2]

virtual_ecosystem.models.abiotic.energy_balance.calculate_sensible_heat_flux(density_air: ndarray[tuple[Any, ...], dtype[floating]], specific_heat_air: ndarray[tuple[Any, ...], dtype[floating]], air_temperature: ndarray[tuple[Any, ...], dtype[floating]], surface_temperature: ndarray[tuple[Any, ...], dtype[floating]], aerodynamic_resistance: float | ndarray[tuple[Any, ...], dtype[floating]]) → ndarray[tuple[Any, ...], dtype[floating]]

Calculate sensible heat flux.

The sensible heat flux \(H\) is calculated using the following equation:

\[H = \frac{\rho_{a} c_{p}}{r_{a}} (T_{s} - T_{a})\]

where \(\rho_{a}\) is the density of air, \(c_{p}\) is the specific heat capacity of air at constant pressure, \(r_{a}\) is the aerodynamic resistance of the surface, \(T_{s}\) is the surface temperature, and \(T_{a}\) is the air temperature.

Parameters:

• density_air – Density of air, [kg m-3]

• specific_heat_air – Specific heat of air, [J kg-1 K-1]

• air_temperature – Air temperature, [C]

• surface_temperature – Surface temperature (canopy or soil), [C]

• aerodynamic_resistance – Aerodynamic resistance, [s m-1]

Returns:

sensible heat flux, [W m-2]

virtual_ecosystem.models.abiotic.energy_balance.calculate_total_absorbed_shortwave_radiation(downward_shortwave_radiation: ndarray[tuple[Any, ...], dtype[floating]], shortwave_absorption_by_canopy: ndarray[tuple[Any, ...], dtype[floating]], fraction_par_used: float, leaf_absorptance_non_par: float, par_fraction: float) → ndarray[tuple[Any, ...], dtype[floating]]

Compute total absorbed shortwave radiation contributing to leaf energy balance.

Shortwave (SW) radiation is partitioned into:

• PAR (photosynthetically active radiation)

• non-PAR (primarily near-infrared, NIR)

The plant model provides absorbed PAR. A fraction of this PAR is used in photosynthesis, while the remainder contributes to heat. Non-PAR radiation is assumed not to be used in photosynthesis and contributes entirely to heat after accounting for leaf absorptance.

This function calculates the total absorbed shortwave radiation that contributes to the leaf energy balance by combining the absorbed PAR (adjusted for photosynthesis) and the absorbed non-PAR radiation (adjusted for leaf absorptance and vertical distribution).

Parameters:

• downward_shortwave_radiation – Incoming shortwave radiation [W m-2]

• shortwave_absorption_by_canopy – Absorbed PAR by canopy [W m-2]

• fraction_par_used – Fraction of absorbed PAR used in photosynthesis (0-1).

• leaf_absorptance_non_par – Fraction of non-PAR radiation absorbed by the leaf (0-1).

• par_fraction – Fraction of total shortwave radiation that is PAR (0-1).

Returns:

Total absorbed shortwave radiation contributing to heat [W m-2]

virtual_ecosystem.models.abiotic.energy_balance.initialise_canopy_and_soil_fluxes(air_temperature: DataArray, layer_structure: LayerStructure, initial_flux_value: float) → dict[str, DataArray]

Initialise canopy temperature and energy fluxes.

This function initializes the following variables to run the first step of the energy balance routine: canopy temperature, [C], sensible and latent heat flux (canopy and soil), and ground heat flux, all in [W m-2].

Parameters:

• air_temperature – Air temperature, [C]

• layer_structure – Instance of LayerStructure

• light_extinction_coefficient – Light extinction coefficient for canopy, unitless

• initial_flux_value – Initial non-zero flux, [W m-2]

Returns:

Dictionary with canopy temperature, [C], sensible and latent heat flux (canopy and soil), [W m-2], and ground heat flux, [W m-2].

virtual_ecosystem.models.abiotic.energy_balance.make_canopy_residual(state: dict[str, Any], static: dict[str, Any], aerodynamic_resistance: ndarray[tuple[Any, ...], dtype[floating]], abiotic_constants: AbioticConstants, core_constants: CoreConstants, surface_index: int) → Callable[[ndarray[tuple[Any, ...], dtype[floating]]], ndarray[tuple[Any, ...], dtype[floating]]]

Creates a residual function for canopy temperature to be used in root finding.

Parameters:

• state – Dictionary containing state variables needed for the energy balance residual.

• static – Dictionary containing static variables needed for the energy balance residual.

• aerodynamic_resistance – Aerodynamic resistance of canopy, [s m-1]

• abiotic_constants – Constants related to abiotic processes.

• core_constants – Core constants.

• surface_index – Row index of surface layer.

Returns:

A function that takes canopy_temperature as input and returns the energy balance residual.

virtual_ecosystem.models.abiotic.energy_balance.secant_solve_cells_layers(residual_function: Callable[[ndarray[tuple[Any, ...], dtype[floating]]], ndarray[tuple[Any, ...], dtype[floating]]], initial_guess: ndarray[tuple[Any, ...], dtype[floating]], maxiter_secant: int, convergence_tolerance: float, small_perturbation_second_guess: float, denominator_tolerance: float) → ndarray[tuple[Any, ...], dtype[floating]]

Vectorised secant solver for independent (cell, layer) root problems.

Parameters:

• residual_function – Function f(T) returning residual with same shape as T.

• initial_guess – Initial guess canopy temperature, shape (n_cells, n_layers).

• maxiter_secant – Maximum secant iterations.

• convergence_tolerance – Convergence tolerance on max absolute update.

• small_perturbation_second_guess – Small perturbation for second initial guess.

• denominator_tolerance – Small value to prevent division by zero in secant update.

Returns:

Root estimate canopy temperature solving f(T)=0 elementwise.

virtual_ecosystem.models.abiotic.energy_balance.update_canopy_air_temperature(air_temperature: ndarray[tuple[Any, ...], dtype[floating]], sensible_heat_flux: ndarray[tuple[Any, ...], dtype[floating]], specific_heat_air: ndarray[tuple[Any, ...], dtype[floating]], density_air: ndarray[tuple[Any, ...], dtype[floating]], mixing_layer_thickness: ndarray[tuple[Any, ...], dtype[floating]]) → ndarray[tuple[Any, ...], dtype[floating]]

Update air temperature surrounding canopy in steady state.

The new air temperature \(T_{a}^{new}\) is updated following Bonan (2019):

\[H = \frac{\rho_a c_p}{r_a}(T_{l} - T_{a})\]

and

\[T_{a}^{new} = T_{a}^{old} + \frac{H}{\rho_a c_p z}\]

where \(\rho_{a}\) is the density of air, \(c_{p}\) is the specific heat capacity of air at constant pressure, \(r_{a}\) is the aerodynamic resistance of the surface, \(T_{s}\) is the surface temperature, \(T_{a}\) is the air temperature, and \(z\) is the thickness of the air layer we are updating.

Parameters:

• air_temperature – Air temperature, [C]

• sensible_heat_flux – Sensible heat flux, [W m-2]

• specific_heat_air – Specific heat capacity of air, [J kg-1 K-1]

• density_air – Density of air, [kg m-3]

• mixing_layer_thickness – thickness of the air layer we are updating, [m]

Returns:

updated air temperatures, [C]

virtual_ecosystem.models.abiotic.energy_balance.update_humidity_vpd(saturated_vapour_pressure: ndarray[tuple[Any, ...], dtype[floating]], specific_humidity_mixed: ndarray[tuple[Any, ...], dtype[floating]], atmospheric_pressure: ndarray[tuple[Any, ...], dtype[floating]], layer_thickness: ndarray[tuple[Any, ...], dtype[floating]], molecular_weight_ratio_water_to_dry_air: float, dry_air_factor: float, density_air: ndarray[tuple[Any, ...], dtype[floating]], cell_area: float, limits_relative_humidity: tuple[float, float, float], limits_vapour_pressure_deficit: tuple[float, float, float], denominator_tolerance: float) → dict[str, ndarray[tuple[Any, ...], dtype[floating]]]

Update atmospheric humidity and vapour pressure deficit.

Parameters:

• saturated_vapour_pressure – Saturated vapour pressure, [kPa]

• specific_humidity_mixed – Specific humidity vertically mixed, [kg kg-1]

• atmospheric_pressure – Atmospheric pressure, [kPa]

• layer_thickness – Layer thickness, [m]

• molecular_weight_ratio_water_to_dry_air – Molecular weight ratio of water to dry air, dimensionless

• dry_air_factor – Complement of water_to_air_mass_ratio, accounting for dry air

• density_air – Density of air, [kg m-3]

• cell_area – Grid cell area, [m2]

• limits_relative_humidity – Realistic bounds of relative humidity, []

• limits_vapour_pressure_deficit – Realistic bounds for vapour pressure deficit, [kPa]

• denominator_tolerance – Small value to prevent division by zero

Returns:

A dictionary containing arrays of updated relative_humidity, specific_humidity, vapour_pressure and vapour_pressure_deficit values.

virtual_ecosystem.models.abiotic.energy_balance.update_soil_temperature(ground_heat_flux: ndarray[tuple[Any, ...], dtype[floating]], soil_temperature: ndarray[tuple[Any, ...], dtype[floating]], soil_layer_thickness: ndarray[tuple[Any, ...], dtype[floating]], soil_thermal_conductivity: float | ndarray[tuple[Any, ...], dtype[floating]], soil_bulk_density: float | ndarray[tuple[Any, ...], dtype[floating]], specific_heat_capacity_soil: float | ndarray[tuple[Any, ...], dtype[floating]], time_interval: float) → ndarray[tuple[Any, ...], dtype[floating]]

Update soil temperature using heat diffusion.

The function applies an explicit finite-difference approach to update soil temperatures based on thermal diffusivity and heat flux.

Governing equations:

Soil thermal diffusivity:

\[\alpha = \frac{\lambda}{\rho_s c_s}\]

where \(\lambda\) is the soil thermal conductivity [W m-1 K-1], \(\rho_s\) is the soil bulk density [kg m-3], \(c_s\) is the specific heat capacity of soil [J kg-1 K-1].

Internal layer update:

\[T_i^{t+\Delta t} = T_i^t + (\Delta t / \Delta z^2) * \alpha * (T_{i+1}^t - 2T_i^t + T_{i-1}^t)\]

Top layer update with ground heat flux:

\[T_0^{t+\Delta t} = T_0^t + (\Delta t / (\rho_s c_s \Delta z)) * G\]

No-heat-flux bottom boundary condition:

\[T_{n-1}^{t+\Delta t} = T_{n-1}^t + (\Delta t / \Delta z^2) * \alpha * (T_{n-2}^t - T_{n-1}^t)\]

Parameters:

• ground_heat_flux – Ground heat flux at top soil, [W m-2]

• soil_temperature – Soil temperature for each soil layer, [C]

• soil_thermal_conductivity – Thermal conductivity of soil, [W m-2 K-1]

• soil_bulk_density – Soil bulk density, [kg m-3]

• specific_heat_capacity_soil – Specific heat capacity of soil, [J kg-1 K-1]

• soil_layer_thickness – Thickness of each soil layer, [m]

• time_interval – Time interval, [s]

Returns:

Updated soil temperatures, [C]

virtual_ecosystem.models.abiotic.energy_balance.update_specific_humidity(evapotranspiration: ndarray[tuple[Any, ...], dtype[floating]], soil_evaporation: ndarray[tuple[Any, ...], dtype[floating]], specific_humidity: ndarray[tuple[Any, ...], dtype[floating]], layer_thickness: ndarray[tuple[Any, ...], dtype[floating]], density_air: ndarray[tuple[Any, ...], dtype[floating]], mm_to_kg: float, cell_area: float, time_interval: float, surface_index: int) → ndarray[tuple[Any, ...], dtype[floating]]

Update specific humidity from evapotranspiration and soil evaporation.

This function adds the water from soil evaporation and canopy evapotranspiration to each atmospheric layer. No limits are applied at this stage and no vertical mixing.

Parameters:

• evapotranspiration – Evapotranspiration, [mm]

• soil_evaporation – Soil evaporation to surface layer, [mm]

• saturated_vapour_pressure – Saturated vapour pressure, [kPa]

• specific_humidity – Specific humidity, [kg kg-1]

• layer_thickness – Layer thickness, [m]

• density_air – Density of air, [kg m-3]

• mm_to_kg – Factor to convert variable unit from millimeters to kilograms of water per square meter

• cell_area – Grid cell area, [m2]

• time_interval – Time interval, [s]

• surface_index – Index of surface layer

Returns:

update specific_humidity, [kg kg-1]

virtual_ecosystem.models.abiotic.energy_balance.update_surface_air_temperature(canopy_air_temperature: ndarray[tuple[Any, ...], dtype[floating]], state: dict[str, ndarray[tuple[Any, ...], dtype[floating]]], idx: SimpleNamespace, denominator_tolerance: float)

Update surface air temperature in equilibrium with soil and canopy fluxes.

The surface air temperature is diagnosed from the soil and canopy bottom conductances and temperatures, assuming equilibrium between the soil and canopy fluxes. This is necessary because the surface layer is too thin to be updated based on fluxes over a 1-hour timestep, and we want to avoid unrealistic surface air temperatures that would arise from a flux-based update.

For cells with fewer canopy layers, the bottom canopy temperature is the last finite value in the canopy temperature array. For cells with no canopy, the above-canopy reference temperature is used instead.

Parameters:

• canopy_air_temperature – Canopy air temperature, [C]

• state – Dictionary of state variables

• idx – Layer structure index

• denominator_tolerance – Small value to prevent division by zero

Returns:

Updated surface air temperature, [C]