 |
0.1.39
|
|
0.1.39
|
| Dependencies | None (CUDA optional for GPU acceleration) |
|---|---|
| Python Import | from pyhelios import EnergyBalanceModel |
| Main Class | EnergyBalanceModel |
| Dependencies | None required (OpenMP recommended for parallel CPU execution) |
|---|---|
| Platforms | Windows, Linux, macOS |
| GPU | Optional - NVIDIA CUDA-capable GPU for GPU acceleration |
| Execution Modes | Three-tier: GPU (CUDA) → OpenMP (parallel CPU) → Serial CPU fallback |
Note: As of helios-core v1.3.61, CUDA is optional. The plugin automatically selects the best available execution mode: GPU acceleration with CUDA if available, parallel CPU with OpenMP if available, or serial CPU as fallback. OpenMP is recommended for most workloads on systems without CUDA.
| Installation Method | Dependency | 
| 
| 
|
|---|---|---|---|---|
| Runtime (pip install) | NVIDIA GPU | Optional - Enables GPU acceleration if available | ||
| CUDA Runtime | Optional - Version 9.0+ for GPU acceleration If not available, uses OpenMP CPU mode or serial fallback | |||
| Runtime (recommended) | OpenMP | Recommended - Enables parallel CPU execution Typically included with GCC/Clang compilers | ||
| Build from Source (optional) | CUDA Toolkit | Optional - Version 9.0+ with nvcc compiler for GPU support Download CUDA Toolkit | ||
Execution Mode Selection (v1.3.61+): The plugin automatically selects the best available execution mode at runtime:
For GPU acceleration setup, see the comprehensive CUDA Setup Guide, which covers:
None.
This model plugin calculates a local energy balance for every primitive, and ultimately predicts sensible, latent, and longwave fluxes as well as surface temperature. The energy balance equation is solved using one of three execution modes:
The execution mode is automatically selected based on available hardware and compiled libraries. GPU acceleration can be explicitly controlled using enableGPUAcceleration() and disableGPUAcceleration().
The model is solving the steady-state budget between absorbed radiation, emitted radiation, sensible heat exchange, and latent heat exchange, which is written as
\[ R-n_s\varepsilon\sigma T_s^4 = c_p g_H \left( T_s-T_a \right) + \lambda g_M \left(\frac{e_s(T_s)f_s-e_s(T_a)h}{p_{atm}}\right)+C_p\frac{dT_s}{dt}+Q_{other}\]
Variables in this equation are listed in this table:
| Variable (units) | Description |
|---|---|
| \(R\) (W/m2) | Absorbed all-wave radiation flux (shortwave+longwave). |
| \(T_s\) (K) | Primitive surface temperature. |
| \(T_a\) (K) | Air temperature just outside primitive boundary-layer. |
| \(g_H\) (mol air/m2-s) | Conductance to heat from primitive surface to outside of boundary-layer. |
| \(g_S\) (mol air/m2-s) | Conductance to moisture between sub-surface air spaces to surface (e.g., for leaves this is stomatal conductance). |
| \(g_M\) (mol air/m2-s) | Conductance to moisture from primitive surface to outside of boundary-layer. |
| \(e_s(T)\) (Pa) | Saturation vapor pressure at temperature T. Calculated from Tetens equation (see en.wikipedia.org/wiki/Tetens_equation) |
| \(f_s\) | Relative humidity of air at primitive surface. |
| \(h\) | Relative humidity of air outside of boundary-layer. |
| \(p_{atm}\) (Pa) | Atmospheric pressure. |
| \(C_p\) (J/m2/oC) | Heat capacity of object. |
| \(Q_{other}\) (W/m2) | Any surface fluxes other than radiation, convection, latent, or storage. |
| \(\varepsilon\) | Surface emissivity for thermal radiation (see below). |
| \(n_s\) | Number of primitive faces with heat transfer, set by the primitive "twosided_flag": \(n_s=1\) for a single-sided primitive (twosided_flag=0, e.g., the ground), and \(n_s=2\) for a two-sided primitive (twosided_flag=1, e.g., leaves). |
The emitted radiation is the same as the emission computed by the RadiationModel in the bands added to the energy balance model, so that the energy a surface emits matches what other surfaces and the sky exchange with it. In the usual case, emission is enabled for a single longwave band without wavelength bounds, which emits \(\varepsilon\sigma T_s^4\), and \(\varepsilon\) is the "emissivity_[band]" of that band. The thermal spectrum may instead be divided among several emitting bands with wavelength bounds, for example to represent an emissivity that varies with wavelength. Each such band emits only its share of the blackbody spectrum, \(F_b(T)\), so the emissivity is the sum over the emitting bands
\[ \varepsilon = \sum_b \varepsilon_b F_b(T_s), \]
with \(F_b=1\) for a band without wavelength bounds. If the scattering depth of an emitting band is 0, the RadiationModel treats surfaces as black in that band and emits \(\sigma T^4\) regardless of the emissivity, so \(\varepsilon_b=1\) is used for it here as well. The blackbody fractions are evaluated at the initial surface temperature, at which the RadiationModel computed the emission. The emitting bands must account for the blackbody spectrum once: run() raises an error if bounded emitting bands overlap, if \(\sum_b F_b\) is less than 0.95 or greater than 1.01, or if more than one band without wavelength bounds emits (or one emits alongside a bounded thermal band). A thermal camera band such as 8-14 \(\mu\)m therefore cannot be the only longwave band of the energy balance, and should not be added to the energy balance model alongside a broadband longwave band.
Constants are given by:
| Constant (units) | Value | Description |
|---|---|---|
| \(c_p\) (J/mol/K) | 29.25 | Heat capacity of air. |
| \(\lambda\) (J/mol) | 44,000 | Latent heat of vaporization of air. |
| Constructors |
|---|
| EnergyBalanceModel |
The EnergyBalanceModel class is initialized by passing the Helios context as an argument to the constructor.
| Primitive Data | Units | Data Type | Description | Available Plug-ins | Default Value |
|---|---|---|---|---|---|
| radiation_flux_[*] | W/m2 | float | Net absorbed radiation flux for band [*] (e.g., radiation_flux_PAR). | Can be computed by RadiationModel plug-in. | N/A (must add at least one band) |
| emissivity_[*] | unitless | float | Surface emissivity for band [*], used for the emitted thermal radiation term. The emissivity of each band for which emission is enabled in the RadiationModel is weighted by the share of the blackbody spectrum that band emits (see Introduction). When no emission information is available (radiation fluxes set manually without running the RadiationModel), the first band that defines an emissivity is used. | Can be set by RadiationModel plug-in. | 1.0 |
| wind_speed | m/s | float | Air wind speed just outside of primitive boundary-layer. | N/A | 1 m/s |
| object_length | m | float | Characteristic dimension of object formed by primitive. | N/A | Square root of primitive surface area |
| boundarylayer_conductance | mol air/m2-s | float | Leaf boundary-layer conductance to heat, summed over both faces of a two-sided primitive. | BoundaryLayerConductanceModel plug-in | Try calculating from model \(0.135\sqrt{\frac{U}{L}}\) |
| air_temperature | Kelvin | float | Ambient air temperature outside of surface boundary layer. | N/A | 300 K |
| moisture_conductance | mol air/m2-s | float | Conductance to moisture between sub-surface air spaces and surface (e.g., for leaves this is stomatal conductance). | Can be computed by StomatalConductanceModel plug-in. | 0 |
| surface_humidity | unitless | float | Relative humidity of air immediately above surface evaporating site. | N/A | 1.0 (saturated) |
| air_humidity | unitless | float | Ambient air relative humidity outside of surface boundary layer. | N/A | 0.5 |
| air_pressure | Pascals | float | Atmospheric pressure. | N/A | 101,000 Pa |
| heat_capacity | J/m2/oC | float | Heat capacity of object. | N/A | 0 |
| other_surface_flux | W/m2 | float | Other surface energy fluxes. | N/A | 0 |
| twosided_flag | N/A | uint | Flag indicating the number of primitive faces with heat transfer (twosided_flag = 0 for one-sided heat transfer; twosided_flag = 1 for two-sided heat transfer). Can also be set via the material system using Context::setMaterialTwosidedFlag(). When a user-assigned material exists, the material's twosided_flag takes precedence over primitive data. | N/A | 1 |
| stomatal_sidedness | \(\zeta\) | float | Ratio of stomatal density on the upper leaf surface to the sum of the stomatal density on upper and lower leaf surfaces. Note: if "twosided_flag" is equal to 0, stomatal_sidedness will be automatically set to 0. | N/A | 0 |
| Primitive Data | Units | Data Type | Description |
|---|---|---|---|
| temperature | Kelvin | float | Primitive surface temperature. |
| sensible_flux | W/m2 | float | Sensible heat flux. |
| latent_flux | W/m2 | float | Latent heat flux. |
| Primitive Data | Units | Data Type | Description |
|---|---|---|---|
| boundarylayer_conductance_out | mol air/m2-s | float | Primitive boundary-layer conductance calculated by this plug-in. |
| vapor_pressure_deficit | mol/mol | float | Surface vapor pressure deficit. |
| storage_flux | W/m2 | float | Storage heat flux. |
| net_radiation_flux | W/m2 | float | Primitive net radiation flux (absorption minus emission). |
Inputs to the model are set by creating primitive variable data in the usual way. If a variable needed for a model input has not been create in the Context, the default value is assumed.
In order for the model to calculate the absorbed all-wave radiation flux, it needs to know the names of all radiation bands that were added to the radiation model. This is done using the function addRadiationBand(). In the following example, assume we have three radiative wavebands "PAR", "NIR", and "LW".
The primitive boundary-layer conductance can either be set using the BoundaryLayerConductanceModel plug-in, or using the default model which is the Polhausen equation. Note also that custom conductance values can also be provided by setting the value of primitive data labeled 'boundarylayer_conductance', which overrides the computed model value.
It is also important to note that, by default, the length scale used to calculate the boundary-layer conductance is taken to be the square root of the primitive surface area. If the size of the object is different from the size of the primitive, then it is important to manually set the length scale to be the size of the object, as this is the relevant scale for boundary-layer development.
The default boundary-layer conductance is calculated as
\(g_H = 0.135 n_s\sqrt{\frac{U}{L}}\),
where \(U\) is the wind speed just outside of the primitive boundary-layer, and \(L\) is the characteristic length/dimension of the object that the primitive belongs to. For a leaf consisting of a single primitive, \(L\) could be assumed to be the length of the primitive. \(n_s\) is the number of primitive faces, which is determined by the value of primitive data "twosided_flag" (twosided_flag=0 is single-sided and \(n_s=1\), twosided_flag=1 is two-sided and \(n_s=2\)).
By default the air temperature and humidity surrounding each leaf are held at whatever the 'air_temperature' and 'air_humidity' primitive data say, so the canopy cannot influence the air that drives its own transpiration. enableCanopyAirspaceModel() resolves the within-canopy air temperature \(T_{ac}\) and vapor pressure \(e_{ac}\) from a vertically layered resistance network instead. The canopy airspace is divided into vertical layers of equal leaf area index, each exchanging sensible heat and water vapor with the leaves it contains, with its neighboring layers, and – for the bottom and top layers – with the soil surface and the above-canopy reference air.
Unlike enableAirEnergyBalance(), which evolves a prognostic atmospheric boundary layer and thereby assumes a horizontally infinite canopy, this model imposes the measured above-canopy air state as a fixed boundary condition. It is therefore appropriate for canopies of limited extent subject to advection, such as an orchard block. The two models determine the same air state and may not both be enabled. Because the canopy airspace model solves for a steady state, it also cannot be combined with the form of run() that takes a timestep.
The above-canopy boundary condition is read from global data 'air_temperature_reference', 'air_humidity_reference' and 'wind_speed_reference' when present, and otherwise falls back to the model's default values.
| Method | Effect |
|---|---|
enableCanopyAirspaceModel(canopy_UUIDs, canopy_height_m, reference_height_m, leaf_area_index, num_layers=1, ground_UUIDs=None) | Enable the model. reference_height_m must exceed canopy_height_m; num_layers=1 gives a single within-canopy node |
disableCanopyAirspaceModel() | Return to a single surface energy balance pass using the current 'air_temperature' and 'air_humidity' |
setCanopyAirspaceConvergence(tolerance_K=0.01, max_iterations=50) | Temperature tolerance (K) and iteration cap of the coupled solution |
enableCanopyAirspaceStabilityCorrection() | Correct turbulent transport above and at the top of the canopy for atmospheric stability (the default) |
disableCanopyAirspaceStabilityCorrection() | Assume neutral stability throughout |
Running the model sets primitive data 'air_temperature', 'air_humidity' and 'wind_speed' on the canopy and ground primitives, and reports global data 'canopy_air_temperature', 'canopy_air_humidity', their per-layer counterparts, 'aerodynamic_resistance', 'canopy_airspace_iterations', 'friction_velocity' (friction velocity above the canopy, m/s) and 'canopy_stability_parameter' (the Monin-Obukhov stability parameter \(\zeta=(z_H-d)/L\) at the reference height; negative is unstable, positive is stable, and zero is neutral – always zero when the stability correction is disabled). It also adds 'boundarylayer_conductance_out' to the optional outputs, which the airspace solution requires.
The aerodynamic resistance above the canopy follows Perrier (1975), which holds for neutral stability. When the canopy is warmer than the air above it, buoyancy enhances turbulent mixing (unstable conditions), and when it is cooler, as on a clear night, buoyancy suppresses it (stable conditions). By default the model corrects for this using Monin-Obukhov similarity theory, so that
\[ u_* = \frac{\kappa u}{\ln\left[\frac{z_m-d}{z_{oM}}\right] - \psi_m\left(\frac{z_m-d}{L}\right) + \psi_m\left(\frac{z_{oM}}{L}\right)}, \qquad r_a = \frac{\ln\left[\frac{z_H-d}{h_c-d}\right] - \psi_h\left(\frac{z_H-d}{L}\right) + \psi_h\left(\frac{h_c-d}{L}\right)}{\kappa u_*}, \]
which reduces to the neutral expression when \(L\rightarrow\infty\). The integrated stability functions \(\psi_m\) and \(\psi_h\) are those of Paulson (1970) with the Businger-Dyer coefficients for unstable conditions and of Beljaars and Holtslag (1991) for stable conditions. The Obukhov length is
\[ L = -\frac{u_*^3\,T}{\kappa g\,\Delta T_v/r_a}, \]
where \(T\) is the reference air temperature and \(\Delta T_v\) is the virtual temperature of the uppermost canopy airspace layer minus that of the reference air. The eddy diffusivity at the canopy top used for the within-canopy layer conductances is likewise corrected, \(K(h_c)=\kappa u_*(h_c-d)/\phi_h\left(\frac{h_c-d}{L}\right)\), so that it is continuous with the resistance above the canopy. The correction is turned off with disableCanopyAirspaceStabilityCorrection(), in which case neutral stability is assumed throughout, and back on with enableCanopyAirspaceStabilityCorrection().
Note: Before helios-core 1.3.93 the canopy airspace model always assumed neutral stability, so daytime canopy air temperatures came out warmer, and nighttime ones closer to the reference air, than they do now.
The ground primitive set may be empty, in which case there is no exchange with the soil surface. Otherwise the soil enters the lowest layer in the same way as a leaf, using the conductances and surface temperatures of the ground primitives as computed by their own surface energy balance, so the heat and water vapor released by the soil are exactly what the airspace receives. Soil fluxes are averaged over the area of the ground primitives, which are taken to represent the soil surface beneath the whole canopy footprint, so ground geometry extending beyond the canopy does not scale the soil contribution. Soil evaporation is governed by the 'moisture_conductance' (zero by default, i.e., a dry soil) and 'surface_humidity' primitive data of the ground primitives. Leaves act as water vapor sources at their 'surface_humidity' (1 by default, i.e., a saturated substomatal airspace) times the saturation vapor pressure at the leaf temperature.
Unless primitive data 'boundarylayer_conductance' is set for a ground primitive, its boundary-layer conductance is the bare-soil relation of Kustas and Norman (1999), \(g_{bl}=41.56\,(0.004+0.012U_s)\) mol/m \(^2\)-s, evaluated at the wind speed at the soil surface \(U_s\), rather than the flat-plate relation used for other primitives. This is the same relation as the "Ground" model of the boundary-layer conductance plug-in. A conductance set by the user or by that plug-in is used unchanged.
Every primitive given to enableCanopyAirspaceModel(), canopy and ground alike, must also be included in the primitives passed to run(). The within-canopy wind speeds are written during run(), so if boundary-layer conductances are instead supplied by the BoundaryLayerConductanceModel run before the energy balance, it reads the 'wind_speed' left by the previous call to run() (or its default value on the first call), and its conductances lag a change in the reference wind speed by one call.
For surfaces that are not completely dry or completely saturated with water, availability of water at the surface determines the rate of moisture transfer from the surface. This is represented by the surface moisture conductance \(g_S\), and corresponds to the conductance between the sub-surface air spaces and the surface.
For a dry surface, \(g_S=0\) and thus there is no moisture transfer at the surface. For a surface that is saturated with water (e.g., a lake, puddle), \(g_S\rightarrow\infty\) as the supply of water at the surface is theoretically infinite.
For surfaces such as leaves and the soil in which surface moisture is present but limited, \(g_S\) will have some intermediate value typically on the order of 0.1 mol air/m2-s.
For leaves, \(g_S\) corresponds to the conductance to water vapor between the intercellular (sub-stomatal) air spaces and the surface of the leaf, which we call the stomatal conductance. When stomata are closed, \(g_S\rightarrow 0\), as water transfer is restricted by the stomata.
For the soil surface, the physical interpretation of \(g_S\) is slightly less intuitive, and corresponds to the rate at which water vapor can diffuse from just below the soil to the surface. This is impacted by the soil texture, and the "tortuosity" of the path water vapor must take through the soil, among other factors.
Once the water vapor has diffused to the surface, the rate of transfer is then determined by the surface boundary-layer conductance. The sub-surface conductance and boundary-layer conductance act in serial, and can be combined to yield an overall moisture conductance between the sub-surface air spaces and the outside of the boundary layer according to
\(g_M = \dfrac{1.08g_H}{n_s}g_S\left[\dfrac{\zeta}{1.08g_H/n_s+g_S\zeta}+\dfrac{(1-\zeta)}{1.08g_H/n_s+g_S(1-\zeta)}\right]\),
where \(g_H\) is the boundary-layer conductance to heat summed over all \(n_s\) faces of the primitive, and \(1.08g_H/n_s\) gives the boundary-layer conductance to moisture of a single face considering the differences in diffusivity between water vapor and heat. Moisture leaving a face crosses only that face's boundary layer, so the stomatal conductance on each face acts in series with the boundary-layer conductance of that face alone.
\(\zeta=\dfrac{D_{upper}}{D_{lower}+D_{upper}}\)
is the stomatal sidedness, which is the ratio of the stomatal density of the upper leaf surface to the sum of the upper and lower leaf surface densities, which is set by the primitive data value "stomatal_sidedness". For leaves, \(\zeta=0\) corresponds to hypostomatous leaves (stomata on one side), and \(\zeta=0.5\) to amphistomatous leaves (stomata equally on both sides). It is important to note that if \(n_s=1\), then the value of \(\zeta\) will be overridden and set to 0.
The model is run assuming steady-state conditions (heat storage term is zero) using the run() function, which will run the model for all primitives in the Context if no argument is given, or will run the model for a subset of primitives if a list of UUIDs is given as the argument. Each UUID may appear only once in that list; a repeated UUID raises an error naming the primitive, as does a UUID that does not exist.
| Model Run Function | Description |
|---|---|
| run() | Run the model for all primitives in the Context. |
| run(uuids) | Run the model for a select set of primitives in the Context, which are specified by a list of their UUIDs. |
Additional functions are available to run the unsteady energy balance model with heat storage. If the timestep argument "dt" is passed with value greater than 0, and the heat capacity \(C_p\) is greater than 0, the unsteady energy balance equation will be applied with the heat storage term included. The following equation is then solved for the surface temperature at time \(t+\Delta t\)
\[ R-n_s\varepsilon\sigma T_s(t+\Delta t)^4 - c_p g_H \left( T_s(t+\Delta t)-T_a \right) - \lambda g_M \left(\frac{e_s(T_s(t+\Delta t))f_s-e_s(T_a)h}{p_{atm}}\right)-Q_{other} - C_p\dfrac{T_s(t+\Delta t) - T_s(t)}{\Delta t} = 0.\]
The timestep should be chosen such that it is significantly smaller than the characteristic time constant of heat storage.
Functions for running the unsteady energy balance model are listed below.
| Model Run Function | Description |
|---|---|
| run(dt) | Run the model for all primitives in the Context and advance in time by dt seconds. |
| run(uuids, dt) | Run the model for a select set of primitives in the Context, which are specified by a list of their UUIDs, and advance in time by dt seconds. |
Note that you can run the steady-state model for some primitives but run the unsteady model for others simply by changing the value of the heat capacity primitive data. This is illustrated in the code example below.
1.13.2