Skip to content

1149 bug fixing abiotic model #1153

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
vgro merged 9 commits into from
Nov 19, 2025
Merged

1149 bug fixing abiotic model #1153

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
1daee3a
flux signes fixes, latvap unit correction
vgro Nov 17, 2025
944ad67
LAI sum only canopy layers
vgro Nov 18, 2025
f623066
note added to flux sign
vgro Nov 18, 2025
e11970c
removed loop, integration over 1h, latvap in J
vgro Nov 19, 2025
282d2b6
comment about integration added
vgro Nov 19, 2025
3a10fb5
docstring unit adjusted
vgro Nov 19, 2025
b5aef66
Merge branch 'develop' into 1149-bug-fixing-abiotic-model
vgro Nov 19, 2025
b3ce676
remove print statement
vgro Nov 19, 2025
8493f92
Merge branch '1149-bug-fixing-abiotic-model' of https://github.com/Im…
vgro Nov 19, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 9 additions & 9 deletions tests/models/abiotic/test_energy_balance.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -207,13 +207,13 @@ def test_energy_balance_residual_only(
absorbed_radiation_canopy=data["shortwave_absorption"][canopy_index].to_numpy(),
specific_heat_air=data["specific_heat_air"][canopy_index].to_numpy(),
density_air=data["density_air"][canopy_index].to_numpy(),
density_water=np.full_like(data["density_air"][canopy_index], 1000),
aerodynamic_resistance=data["aerodynamic_resistance_canopy"][
canopy_index
].to_numpy(),
latent_heat_vapourisation=data["latent_heat_vapourisation"][
canopy_index
].to_numpy(),
].to_numpy()
* 1000,
leaf_emissivity=fixture_abiotic_constants.leaf_emissivity,
stefan_boltzmann_constant=fixture_core_constants.stefan_boltzmann_constant,
zero_Celsius=fixture_core_constants.zero_Celsius,
Expand Down Expand Up @@ -248,13 +248,13 @@ def test_energy_balance_return_fluxes(
absorbed_radiation_canopy=data["shortwave_absorption"][canopy_index].to_numpy(),
specific_heat_air=data["specific_heat_air"][canopy_index].to_numpy(),
density_air=data["density_air"][canopy_index].to_numpy(),
density_water=np.full_like(data["density_air"][canopy_index], 1000),
aerodynamic_resistance=data["aerodynamic_resistance_canopy"][
canopy_index
].to_numpy(),
latent_heat_vapourisation=data["latent_heat_vapourisation"][
canopy_index
].to_numpy(),
].to_numpy()
* 1000,
leaf_emissivity=fixture_abiotic_constants.leaf_emissivity,
stefan_boltzmann_constant=fixture_core_constants.stefan_boltzmann_constant,
zero_Celsius=fixture_core_constants.zero_Celsius,
Expand Down Expand Up @@ -295,13 +295,13 @@ def test_solve_canopy_temperature(
absorbed_radiation_canopy=data["shortwave_absorption"][canopy_index].to_numpy(),
specific_heat_air=data["specific_heat_air"][canopy_index].to_numpy(),
density_air=data["density_air"][canopy_index].to_numpy(),
density_water=1000.0,
aerodynamic_resistance=data["aerodynamic_resistance_canopy"][
canopy_index
].to_numpy(),
latent_heat_vapourisation=data["latent_heat_vapourisation"][
canopy_index
].to_numpy(),
].to_numpy()
* 1000,
emissivity_leaf=0.96,
stefan_boltzmann_constant=fixture_core_constants.stefan_boltzmann_constant,
zero_Celsius=fixture_core_constants.zero_Celsius,
Expand Down Expand Up @@ -352,9 +352,9 @@ def test_update_air_temperature(

exp_result = np.array(
[
[29.883755, 29.883755, 29.857958, np.nan],
[28.902139, 28.886732, np.nan, np.nan],
[27.224235, np.nan, np.nan, np.nan],
[29.806235, 29.806235, 29.832032, np.nan],
[28.840201, 28.855608, np.nan, np.nan],
[27.188575, np.nan, np.nan, np.nan],
]
)
np.testing.assert_allclose(updated_air_temperature, exp_result, rtol=1e-4)
Expand Down
176 changes: 0 additions & 176 deletions tests/models/abiotic/test_microclimate.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,179 +84,3 @@ def test_run_microclimate(
assert (
(valid_values_rel_hum_clean >= 0.0) & (valid_values_rel_hum_clean <= 100.0)
).all()


def test_run_microclimate_subdaily(
dummy_climate_data_varying_canopy,
fixture_core_components,
fixture_core_constants,
fixture_abiotic_constants,
fixture_abiotic_simple_configuration,
fixture_pyrealm_config,
):
"""Test microclimate function iterates over hours - no time index."""

# TODO this test returns different results on windows machines (around 1.5 K),
# likely because of differences in rounding digits.
from virtual_ecosystem.models.abiotic.microclimate import (
run_microclimate,
)

lyr_str = fixture_core_components.layer_structure

# Adjust input data to more realistsic values over this time scale
dummy_climate_data_varying_canopy["canopy_evaporation"] = (
dummy_climate_data_varying_canopy["canopy_evaporation"] / 86400 * 3600
)
dummy_climate_data_varying_canopy["transpiration"] = (
dummy_climate_data_varying_canopy["transpiration"] / 86400 * 3600
)
result = run_microclimate(
data=dummy_climate_data_varying_canopy,
time_index=0,
time_interval=3600 * 4,
cell_area=10000,
layer_structure=lyr_str,
abiotic_constants=fixture_abiotic_constants,
core_constants=fixture_core_constants,
pyrealm_core_constants=fixture_pyrealm_config.core,
abiotic_bounds=fixture_abiotic_simple_configuration.bounds,
)

for var in [
"canopy_temperature",
"air_temperature",
"sensible_heat_flux",
"latent_heat_flux",
]:
assert var in result

# Check that values fall within a reasonable expected range
soil_temps = result["soil_temperature"].isel(layers=lyr_str.index_all_soil)

# To test with varying canopy layers, need to mask
canopy_mask = ~np.isnan(
dummy_climate_data_varying_canopy["canopy_temperature"].isel(
layers=lyr_str.index_filled_canopy
)
)
atm_mask = ~np.isnan(
dummy_climate_data_varying_canopy["air_temperature"].isel(
layers=lyr_str.index_filled_atmosphere
)
)

canopy_temp_result = result["canopy_temperature"].isel(
layers=lyr_str.index_filled_canopy
)
air_temp_result = result["air_temperature"].isel(
layers=lyr_str.index_filled_atmosphere
)
rel_hum_result = result["relative_humidity"].isel(
layers=lyr_str.index_filled_atmosphere
)

# Use the mask as a DataArray for .where()
valid_values_can_temp = canopy_temp_result.where(canopy_mask)
valid_values_air_temp = air_temp_result.where(atm_mask)
valid_values_rel_hum = rel_hum_result.where(atm_mask)

# Now drop the NaNs (i.e., masked values)
valid_values_can_temp_clean = valid_values_can_temp.dropna(dim="layers", how="any")
valid_values_air_temp_clean = valid_values_air_temp.dropna(dim="layers", how="any")
valid_values_rel_hum_clean = valid_values_rel_hum.dropna(dim="layers", how="any")

# Now do the test
assert ((soil_temps >= 18.0) & (soil_temps <= 25.0)).all()
assert (
(valid_values_can_temp_clean >= 15.0) & (valid_values_can_temp_clean <= 40.0)
).all()
assert (
(valid_values_air_temp_clean >= 15.0) & (valid_values_air_temp_clean <= 40.0)
).all()
assert (
(valid_values_rel_hum_clean >= 0.0) & (valid_values_rel_hum_clean <= 100.0)
).all()


def test_run_microclimate_minutes(
dummy_climate_data_varying_canopy,
fixture_core_components,
fixture_core_constants,
fixture_abiotic_constants,
fixture_abiotic_simple_configuration,
fixture_pyrealm_config,
):
"""Test microclimate function iterates once for <1h time interval."""

from virtual_ecosystem.models.abiotic.microclimate import (
run_microclimate,
)

lyr_str = fixture_core_components.layer_structure
# Adjust input data to more realistsic values over this time scale
dummy_climate_data_varying_canopy["canopy_evaporation"] = (
dummy_climate_data_varying_canopy["canopy_evaporation"] / 86400 * 60
)
dummy_climate_data_varying_canopy["transpiration"] = (
dummy_climate_data_varying_canopy["transpiration"] / 86400 * 60
)

result = run_microclimate(
data=dummy_climate_data_varying_canopy,
time_index=0,
time_interval=60,
cell_area=10000,
layer_structure=lyr_str,
abiotic_constants=fixture_abiotic_constants,
core_constants=fixture_core_constants,
pyrealm_core_constants=fixture_pyrealm_config.core,
abiotic_bounds=fixture_abiotic_simple_configuration.bounds,
)

# Check that values fall within a reasonable expected range
soil_temps = result["soil_temperature"].isel(layers=lyr_str.index_all_soil)

# To test with varying canopy layers, need to mask
canopy_mask = ~np.isnan(
dummy_climate_data_varying_canopy["canopy_temperature"].isel(
layers=lyr_str.index_filled_canopy
)
)
atm_mask = ~np.isnan(
dummy_climate_data_varying_canopy["air_temperature"].isel(
layers=lyr_str.index_filled_atmosphere
)
)

canopy_temp_result = result["canopy_temperature"].isel(
layers=lyr_str.index_filled_canopy
)
air_temp_result = result["air_temperature"].isel(
layers=lyr_str.index_filled_atmosphere
)
rel_hum_result = result["relative_humidity"].isel(
layers=lyr_str.index_filled_atmosphere
)

# Use the mask as a DataArray for .where()
valid_values_can_temp = canopy_temp_result.where(canopy_mask)
valid_values_air_temp = air_temp_result.where(atm_mask)
valid_values_rel_hum = rel_hum_result.where(atm_mask)

# Now drop the NaNs (i.e., masked values)
valid_values_can_temp_clean = valid_values_can_temp.dropna(dim="layers", how="any")
valid_values_air_temp_clean = valid_values_air_temp.dropna(dim="layers", how="any")
valid_values_rel_hum_clean = valid_values_rel_hum.dropna(dim="layers", how="any")

# Now do the test
assert ((soil_temps >= 18.0) & (soil_temps <= 25.0)).all()
assert (
(valid_values_can_temp_clean >= 15.0) & (valid_values_can_temp_clean <= 40.0)
).all()
assert (
(valid_values_air_temp_clean >= 15.0) & (valid_values_air_temp_clean <= 40.0)
).all()
assert (
(valid_values_rel_hum_clean >= 0.0) & (valid_values_rel_hum_clean <= 100.0)
).all()
21 changes: 4 additions & 17 deletions virtual_ecosystem/models/abiotic/energy_balance.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -256,7 +256,6 @@ def calculate_energy_balance_residual(
absorbed_radiation_canopy: NDArray[np.floating],
specific_heat_air: NDArray[np.floating],
density_air: NDArray[np.floating],
density_water: float | NDArray[np.floating],
aerodynamic_resistance: NDArray[np.floating],
latent_heat_vapourisation: NDArray[np.floating],
leaf_emissivity: float,
Expand Down Expand Up @@ -288,9 +287,8 @@ def calculate_energy_balance_residual(
[W m-2]
specific_heat_air: Specific heat capacity of air, [J kg-1 K-1]
density_air: Density of air, [kg m-3]
density_water: Density of water, [kg m-3]
aerodynamic_resistance: Aerodynamic resistamce of canopy, [s m-1]
latent_heat_vapourisation: Latent heat of vapourisation, [kJ kg-1]
latent_heat_vapourisation: Latent heat of vapourisation, [J kg-1]
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
Expand Down Expand Up @@ -322,9 +320,8 @@ def calculate_energy_balance_residual(
# Latent heat flux canopy, [W m-2]
# The current implementation converts outputs from plant and hydrology model to
# ensure energy conservation between modules for now.
# TODO cross-check units with plant model, time step currently hour to second
latent_heat_flux_canopy = (
evapotranspiration * density_water * latent_heat_vapourisation
evapotranspiration * latent_heat_vapourisation
) / seconds_to_hour

# Energy balance residual, [W m-2]
Expand Down Expand Up @@ -355,7 +352,6 @@ def solve_canopy_temperature(
absorbed_radiation_canopy: NDArray[np.floating],
specific_heat_air: NDArray[np.floating],
density_air: NDArray[np.floating],
density_water: float | NDArray[np.floating],
aerodynamic_resistance: NDArray[np.floating],
latent_heat_vapourisation: NDArray[np.floating],
emissivity_leaf: float,
Expand Down Expand Up @@ -416,10 +412,9 @@ def solve_canopy_temperature(
[W m-2]
specific_heat_air: Specific heat capacity of air, [J kg-1 K-1]
density_air: Density of air, [kg m-3]
density_water: Density of water, [kg m-3]
aerodynamic_resistance: Aerodynamic resistamce of canopy, [s m-1]
stomatal_resistance: Stomatal resistance, [s m-1]
latent_heat_vapourisation: Latent heat of vapourisation, [kJ kg-1]
latent_heat_vapourisation: Latent heat of vapourisation, [J kg-1]
emissivity_leaf: Leaf emissivity, dimensionless
stefan_boltzmann_constant: Stefan Boltzmann constant, [W m-2 K-4]
zero_Celsius: Factor to convert between Celsius and Kelvin
Expand All @@ -438,11 +433,6 @@ def solve_canopy_temperature(
nrows, ncols = canopy_temperature_initial.shape
solved_temperature = np.empty_like(canopy_temperature_initial, dtype=np.float64)

if isinstance(density_water, float):
density_water_array = np.full(solved_temperature.shape, density_water)
else:
density_water_array = density_water

# TODO this loop might be a potential performance bottleneck.
# The function only takes scalar values
for i in range(nrows):
Expand All @@ -467,9 +457,6 @@ def residual_func(canopy_temp_scalar):
[[specific_heat_air[i, j]]], dtype=np.float64
),
density_air=np.array([[density_air[i, j]]], dtype=np.float64),
density_water=np.array(
[[density_water_array[i, j]]], dtype=np.float64
),
aerodynamic_resistance=np.array(
[[aerodynamic_resistance[i, j]]], dtype=np.float64
),
Expand Down Expand Up @@ -561,7 +548,7 @@ def update_air_temperature(

# Update air temperature over a layer of height z (e.g., canopy height)
new_air_temperature = air_temperature + (
-sensible_heat_flux / (density_air * specific_heat_air * mixing_layer_thickness)
sensible_heat_flux / (density_air * specific_heat_air * mixing_layer_thickness)
)
return new_air_temperature

Expand Down
Loading