From 1daee3a46d3bed0fb7aad9f574b500b15802c27b Mon Sep 17 00:00:00 2001 From: vgro Date: Mon, 17 Nov 2025 14:56:59 +0000 Subject: [PATCH 1/7] flux signes fixes, latvap unit correction --- tests/models/abiotic/test_energy_balance.py | 6 +++--- virtual_ecosystem/models/abiotic/energy_balance.py | 2 +- virtual_ecosystem/models/abiotic/microclimate.py | 9 +++++---- 3 files changed, 9 insertions(+), 8 deletions(-) diff --git a/tests/models/abiotic/test_energy_balance.py b/tests/models/abiotic/test_energy_balance.py index 89cdaff96..ae3fd0001 100644 --- a/tests/models/abiotic/test_energy_balance.py +++ b/tests/models/abiotic/test_energy_balance.py @@ -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) diff --git a/virtual_ecosystem/models/abiotic/energy_balance.py b/virtual_ecosystem/models/abiotic/energy_balance.py index 47c3d537a..1b469f4a1 100644 --- a/virtual_ecosystem/models/abiotic/energy_balance.py +++ b/virtual_ecosystem/models/abiotic/energy_balance.py @@ -561,7 +561,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 diff --git a/virtual_ecosystem/models/abiotic/microclimate.py b/virtual_ecosystem/models/abiotic/microclimate.py index ee6779e57..bed428de5 100644 --- a/virtual_ecosystem/models/abiotic/microclimate.py +++ b/virtual_ecosystem/models/abiotic/microclimate.py @@ -269,11 +269,12 @@ def run_microclimate( data["soil_evaporation"].to_numpy() * core_constants.density_water * latent_heat_vapourisation[-1] + / 1000 ) / core_constants.seconds_to_hour # Ground heat flux, [W m-2] ground_heat_flux = ( - net_radiation_soil - latent_heat_flux_soil - sensible_heat_flux_soil + net_radiation_soil + latent_heat_flux_soil + sensible_heat_flux_soil ) # Update soil temperatures, [C] @@ -309,7 +310,7 @@ def run_microclimate( density_air=density_air[1:-1], density_water=core_constants.density_water, aerodynamic_resistance=aerodynamic_resistance_canopy, - latent_heat_vapourisation=latent_heat_vapourisation[1:-1], + latent_heat_vapourisation=latent_heat_vapourisation[1:-1] / 1000, emissivity_leaf=abiotic_constants.leaf_emissivity, stefan_boltzmann_constant=core_constants.stefan_boltzmann_constant, zero_Celsius=core_constants.zero_Celsius, @@ -425,7 +426,7 @@ def run_microclimate( density_air=density_air[1:-1], density_water=core_constants.density_water, aerodynamic_resistance=aerodynamic_resistance_canopy, - latent_heat_vapourisation=latent_heat_vapourisation[1:-1], + latent_heat_vapourisation=latent_heat_vapourisation[1:-1] / 1000, stefan_boltzmann_constant=core_constants.stefan_boltzmann_constant, zero_Celsius=core_constants.zero_Celsius, seconds_to_hour=core_constants.seconds_to_hour, @@ -513,7 +514,7 @@ def run_microclimate( air_temperature_out = layer_structure.from_template() air_temperature_out[layer_structure.index_above] = all_air_temperature[0] air_temperature_out[layer_structure.index_filled_canopy] = air_temperature_canopy - air_temperature_out[layer_structure.index_surface] = surface_air_temperature + air_temperature_out[layer_structure.index_surface_scalar] = surface_air_temperature output["air_temperature"] = air_temperature_out canopy_temperature_out = layer_structure.from_template() From 944ad675e40d7703015353f9de38efbfa6804ce7 Mon Sep 17 00:00:00 2001 From: vgro Date: Tue, 18 Nov 2025 13:03:06 +0000 Subject: [PATCH 2/7] LAI sum only canopy layers --- virtual_ecosystem/models/abiotic/microclimate.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/virtual_ecosystem/models/abiotic/microclimate.py b/virtual_ecosystem/models/abiotic/microclimate.py index bed428de5..38e279424 100644 --- a/virtual_ecosystem/models/abiotic/microclimate.py +++ b/virtual_ecosystem/models/abiotic/microclimate.py @@ -56,7 +56,10 @@ def run_microclimate( # Selection of often used subsets (could be moved to separate function, e.g. tools) # NOTE Canopy height will likely become a separate variable, update as required canopy_height = data["layer_heights"][1].to_numpy() - leaf_area_index_sum = np.nansum(data["leaf_area_index"].to_numpy(), axis=0) + # NOTE currently sums LAI over all canopy layers, not surface grass layer + leaf_area_index_sum = np.nansum( + data["leaf_area_index"][layer_structure.index_filled_canopy].to_numpy(), axis=0 + ) atmospheric_pressure_out = layer_structure.from_template() atmospheric_pressure_out[layer_structure.index_filled_atmosphere] = data[ From f623066cce4bcf96aa12aebbb2f2317dd94737b6 Mon Sep 17 00:00:00 2001 From: vgro Date: Tue, 18 Nov 2025 13:47:21 +0000 Subject: [PATCH 3/7] note added to flux sign --- virtual_ecosystem/models/abiotic/microclimate.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/virtual_ecosystem/models/abiotic/microclimate.py b/virtual_ecosystem/models/abiotic/microclimate.py index 38e279424..b458951c4 100644 --- a/virtual_ecosystem/models/abiotic/microclimate.py +++ b/virtual_ecosystem/models/abiotic/microclimate.py @@ -276,6 +276,8 @@ def run_microclimate( ) / core_constants.seconds_to_hour # Ground heat flux, [W m-2] + # Note the convention is that latent and sensible heat fluxes are negative when + # directed away from the surface, hence added here ground_heat_flux = ( net_radiation_soil + latent_heat_flux_soil + sensible_heat_flux_soil ) From e11970c54d829872a469102d4fb6d61eca2b8045 Mon Sep 17 00:00:00 2001 From: vgro Date: Wed, 19 Nov 2025 11:26:28 +0000 Subject: [PATCH 4/7] removed loop, integration over 1h, latvap in J --- tests/models/abiotic/test_energy_balance.py | 12 +- tests/models/abiotic/test_microclimate.py | 176 -------- .../models/abiotic/energy_balance.py | 15 +- .../models/abiotic/microclimate.py | 380 +++++++++--------- 4 files changed, 188 insertions(+), 395 deletions(-) diff --git a/tests/models/abiotic/test_energy_balance.py b/tests/models/abiotic/test_energy_balance.py index ae3fd0001..c515dc154 100644 --- a/tests/models/abiotic/test_energy_balance.py +++ b/tests/models/abiotic/test_energy_balance.py @@ -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, @@ -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, @@ -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, diff --git a/tests/models/abiotic/test_microclimate.py b/tests/models/abiotic/test_microclimate.py index 7af70fe19..6052c256d 100644 --- a/tests/models/abiotic/test_microclimate.py +++ b/tests/models/abiotic/test_microclimate.py @@ -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() diff --git a/virtual_ecosystem/models/abiotic/energy_balance.py b/virtual_ecosystem/models/abiotic/energy_balance.py index 1b469f4a1..fa1669023 100644 --- a/virtual_ecosystem/models/abiotic/energy_balance.py +++ b/virtual_ecosystem/models/abiotic/energy_balance.py @@ -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, @@ -288,7 +287,6 @@ 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] leaf_emissivity: Leaf emissivity, dimensionless @@ -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] @@ -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, @@ -416,7 +412,6 @@ 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] @@ -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): @@ -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 ), diff --git a/virtual_ecosystem/models/abiotic/microclimate.py b/virtual_ecosystem/models/abiotic/microclimate.py index b458951c4..6b9976649 100644 --- a/virtual_ecosystem/models/abiotic/microclimate.py +++ b/virtual_ecosystem/models/abiotic/microclimate.py @@ -53,7 +53,8 @@ def run_microclimate( output = {} - # Selection of often used subsets (could be moved to separate function, e.g. tools) + # Precompute reusable quantities + # NOTE Canopy height will likely become a separate variable, update as required canopy_height = data["layer_heights"][1].to_numpy() # NOTE currently sums LAI over all canopy layers, not surface grass layer @@ -179,7 +180,7 @@ def run_microclimate( ) # ------------------------------------------------------------------------- - # Initialise variables to iterate energy balance to update temperatures + # Initialise temperatures and humidity profiles # ------------------------------------------------------------------------- all_air_temperature = data["air_temperature"][ @@ -201,220 +202,201 @@ def run_microclimate( layer_structure.index_filled_atmosphere ].to_numpy() - # Get combined evapotranspiration from plant and hydrology model, per hour + # Evapotranspiration from plant and hydrology model, per time interval evapotranspiration = data["canopy_evaporation"] + data["transpiration"] - # If hourly input data is provided, iterate over day, else equilibrium assumption - hourly_time_interval = max(int(time_interval / core_constants.seconds_to_hour), 1) - - # TODO Run diurnal cycle - # enable hourly input in data object, select by time index, iterate over hours, - # and return averages/min/max over the VE time interval, similar to hydrology - if core_constants.seconds_to_hour <= time_interval <= core_constants.seconds_to_day: - iteration = hourly_time_interval - - else: - iteration = 1 - - for _ in range(iteration): - # ------------------------------------------------------------------------- - # Calculate atmospheric background variables - # ------------------------------------------------------------------------- - # Density of air, [kg m-3] - density_air = abiotic_tools.calculate_air_density( - air_temperature=all_air_temperature, - atmospheric_pressure=atmospheric_pressure, - specific_gas_constant_dry_air=core_constants.specific_gas_constant_dry_air, - celsius_to_kelvin=core_constants.zero_Celsius, - ) + # ------------------------------------------------------------------------- + # Calculate atmospheric background variables + # ------------------------------------------------------------------------- + # Density of air, [kg m-3] + density_air = abiotic_tools.calculate_air_density( + air_temperature=all_air_temperature, + atmospheric_pressure=atmospheric_pressure, + specific_gas_constant_dry_air=core_constants.specific_gas_constant_dry_air, + celsius_to_kelvin=core_constants.zero_Celsius, + ) - # Specific heat capacity of air, [J kg-1 K-1] - specific_heat_air = calc_specific_heat( - tc=all_air_temperature, - ) + # Specific heat capacity of air, [J kg-1 K-1] + specific_heat_air = calc_specific_heat( + tc=all_air_temperature, + ) - # Latent heat of vapourisation, [kJ kg-1] - latent_heat_vapourisation = abiotic_tools.calculate_latent_heat_vapourisation( - temperature=all_air_temperature, - celsius_to_kelvin=core_constants.zero_Celsius, - latent_heat_vap_equ_factors=abiotic_constants.latent_heat_vap_equ_factors, - ) + # Latent heat of vapourisation, [kJ kg-1] + latent_heat_vapourisation = abiotic_tools.calculate_latent_heat_vapourisation( + temperature=all_air_temperature, + celsius_to_kelvin=core_constants.zero_Celsius, + latent_heat_vap_equ_factors=abiotic_constants.latent_heat_vap_equ_factors, + ) - # ------------------------------------------------------------------------- - # Soil energy balance - # ------------------------------------------------------------------------- - # Longwave emission from soil, [W m-2] - longwave_emission_soil = energy_balance.calculate_longwave_emission( - temperature=soil_temperature[0] + core_constants.zero_Celsius, - emissivity=abiotic_constants.soil_emissivity, - stefan_boltzmann=core_constants.stefan_boltzmann_constant, - ) + # ------------------------------------------------------------------------- + # Soil energy balance + # ------------------------------------------------------------------------- + # Longwave emission from soil, [W m-2] + longwave_emission_soil = energy_balance.calculate_longwave_emission( + temperature=soil_temperature[0] + core_constants.zero_Celsius, + emissivity=abiotic_constants.soil_emissivity, + stefan_boltzmann=core_constants.stefan_boltzmann_constant, + ) - # Net radiation topsoil, shortwave in - longwave out, [W m-2] - net_radiation_soil = ( - data["shortwave_absorption"][ - layer_structure.index_topsoil_scalar - ].to_numpy() - - longwave_emission_soil - ) + # Net radiation topsoil, shortwave in - longwave out, [W m-2] + net_radiation_soil = ( + data["shortwave_absorption"][layer_structure.index_topsoil_scalar].to_numpy() + - longwave_emission_soil + ) - # Sensible heat flux from topsoil, [W m-2] - sensible_heat_flux_soil = energy_balance.calculate_sensible_heat_flux( - density_air=density_air[-1], - specific_heat_air=specific_heat_air[-1], - air_temperature=surface_air_temperature, - surface_temperature=soil_temperature[0], - aerodynamic_resistance=aerodynamic_resistance_soil, - ) + # Sensible heat flux from topsoil, [W m-2] + sensible_heat_flux_soil = energy_balance.calculate_sensible_heat_flux( + density_air=density_air[-1], + specific_heat_air=specific_heat_air[-1], + air_temperature=surface_air_temperature, + surface_temperature=soil_temperature[0], + aerodynamic_resistance=aerodynamic_resistance_soil, + ) - # Latent heat flux topsoil, [W m-2] - latent_heat_flux_soil = ( - data["soil_evaporation"].to_numpy() - * core_constants.density_water - * latent_heat_vapourisation[-1] - / 1000 - ) / core_constants.seconds_to_hour - - # Ground heat flux, [W m-2] - # Note the convention is that latent and sensible heat fluxes are negative when - # directed away from the surface, hence added here - ground_heat_flux = ( - net_radiation_soil + latent_heat_flux_soil + sensible_heat_flux_soil - ) + # Latent heat flux topsoil, [W m-2] + latent_heat_flux_soil = ( + data["soil_evaporation"].to_numpy() * latent_heat_vapourisation[-1] * 1000 + ) / time_interval - # Update soil temperatures, [C] - # TODO Revisit implementation of soil temperature update, consider Newton - # TODO Soil parameter currently constants, replace with soil maps - # TODO include effect of soil moisture - soil_temperature = energy_balance.update_soil_temperature( - ground_heat_flux=ground_heat_flux, - soil_temperature=soil_temperature, - soil_layer_thickness=layer_structure.soil_layer_thickness, - soil_thermal_conductivity=abiotic_constants.soil_thermal_conductivity, - soil_bulk_density=abiotic_constants.bulk_density_soil, - specific_heat_capacity_soil=abiotic_constants.specific_heat_capacity_soil, - time_interval=core_constants.seconds_to_hour, - ) + # Ground heat flux, [W m-2] + # Note the convention is that latent and sensible heat fluxes are negative when + # directed away from the surface, hence added here + ground_heat_flux = ( + net_radiation_soil + latent_heat_flux_soil + sensible_heat_flux_soil + ) - # ------------------------------------------------------------------------- - # Update canopy and air temperatures using the Newton method - # ------------------------------------------------------------------------- - - # Solve energy balance for canopy temperature, [C] - canopy_temperature = energy_balance.solve_canopy_temperature( - canopy_temperature_initial=canopy_temperature, - air_temperature=air_temperature_canopy, - evapotranspiration=evapotranspiration[ - layer_structure.index_filled_canopy - ].to_numpy() - / hourly_time_interval, - absorbed_radiation_canopy=data["shortwave_absorption"][ - layer_structure.index_filled_canopy - ].to_numpy(), - specific_heat_air=specific_heat_air[1:-1], - density_air=density_air[1:-1], - density_water=core_constants.density_water, - aerodynamic_resistance=aerodynamic_resistance_canopy, - latent_heat_vapourisation=latent_heat_vapourisation[1:-1] / 1000, - emissivity_leaf=abiotic_constants.leaf_emissivity, - stefan_boltzmann_constant=core_constants.stefan_boltzmann_constant, - zero_Celsius=core_constants.zero_Celsius, - seconds_to_hour=core_constants.seconds_to_hour, - return_fluxes=False, - maxiter=10000, - ) + # Update soil temperatures, [C], integration interval 1 hour + # TODO Revisit implementation of soil temperature update, consider Newton + # TODO Soil parameter currently constants, replace with soil maps + # TODO include effect of soil moisture + soil_temperature = energy_balance.update_soil_temperature( + ground_heat_flux=ground_heat_flux, + soil_temperature=soil_temperature, + soil_layer_thickness=layer_structure.soil_layer_thickness, + soil_thermal_conductivity=abiotic_constants.soil_thermal_conductivity, + soil_bulk_density=abiotic_constants.bulk_density_soil, + specific_heat_capacity_soil=abiotic_constants.specific_heat_capacity_soil, + time_interval=core_constants.seconds_to_hour, + ) - # Update air temperature based on new canopy and soil temperatures, [C] - air_temperature_canopy = energy_balance.update_air_temperature( - air_temperature=air_temperature_canopy, - surface_temperature=canopy_temperature, - specific_heat_air=specific_heat_air[1:-1], - density_air=density_air[1:-1], - aerodynamic_resistance=aerodynamic_resistance_canopy, - mixing_layer_thickness=above_ground_layer_thickness[1:-1], - ) + # ------------------------------------------------------------------------- + # Update canopy and air temperatures using the Newton method + # ------------------------------------------------------------------------- - surface_air_temperature = energy_balance.update_air_temperature( - air_temperature=surface_air_temperature, - surface_temperature=soil_temperature[0], - specific_heat_air=specific_heat_air[-1], - density_air=density_air[-1], - aerodynamic_resistance=aerodynamic_resistance_soil, - mixing_layer_thickness=above_ground_layer_thickness[-1], - ) + # Solve energy balance for canopy temperature, [C], integration interval 1 hour + canopy_temperature = energy_balance.solve_canopy_temperature( + canopy_temperature_initial=canopy_temperature, + air_temperature=air_temperature_canopy, + evapotranspiration=evapotranspiration[ + layer_structure.index_filled_canopy + ].to_numpy() + / (time_interval * core_constants.seconds_to_hour), + absorbed_radiation_canopy=data["shortwave_absorption"][ + layer_structure.index_filled_canopy + ].to_numpy(), + specific_heat_air=specific_heat_air[1:-1], + density_air=density_air[1:-1], + aerodynamic_resistance=aerodynamic_resistance_canopy, + latent_heat_vapourisation=latent_heat_vapourisation[1:-1] * 1000, + emissivity_leaf=abiotic_constants.leaf_emissivity, + stefan_boltzmann_constant=core_constants.stefan_boltzmann_constant, + zero_Celsius=core_constants.zero_Celsius, + seconds_to_hour=core_constants.seconds_to_hour, + return_fluxes=False, + maxiter=10000, + ) + print(canopy_temperature) + # Update air temperature based on new canopy and soil temperatures, [C] + air_temperature_canopy = energy_balance.update_air_temperature( + air_temperature=air_temperature_canopy, + surface_temperature=canopy_temperature, + specific_heat_air=specific_heat_air[1:-1], + density_air=density_air[1:-1], + aerodynamic_resistance=aerodynamic_resistance_canopy, + mixing_layer_thickness=above_ground_layer_thickness[1:-1], + ) - all_air_temperature[1 : len(canopy_temperature) + 1] = air_temperature_canopy - all_air_temperature[-1] = surface_air_temperature + surface_air_temperature = energy_balance.update_air_temperature( + air_temperature=surface_air_temperature, + surface_temperature=soil_temperature[0], + specific_heat_air=specific_heat_air[-1], + density_air=density_air[-1], + aerodynamic_resistance=aerodynamic_resistance_soil, + mixing_layer_thickness=above_ground_layer_thickness[-1], + ) - all_air_temperature = wind.mix_and_ventilate( - input_variable=all_air_temperature, - layer_thickness=above_ground_layer_thickness, - ventilation_rate=ventilation_rate, - mixing_coefficient=mixing_coefficient, - limits=abiotic_bounds.air_temperature[:2], - time_interval=core_constants.seconds_to_hour, - ) + all_air_temperature[1 : len(canopy_temperature) + 1] = air_temperature_canopy + all_air_temperature[-1] = surface_air_temperature - # NOTE Advection not implemented as everything is removed with time interval>1h - # and horizontal transfer is not implemented - # advection_rate = ( - # data["wind_speed_ref"].isel(time_index=time_index).to_numpy() - # / np.sqrt(cell_area) - # ) - # advected_fraction = np.clip(advection_rate * time_interval, 0, 1) - # all_air_temperature[0] -=all_air_temperature[0] *advected_fraction - - # Update atmospheric humidity/VPD - # Saturated vapour pressure of air, [kPa] - saturated_vapour_pressure_air = calc_vp_sat( - ta=all_air_temperature, - core_const=pyrealm_core_constants, - ) + all_air_temperature = wind.mix_and_ventilate( + input_variable=all_air_temperature, + layer_thickness=above_ground_layer_thickness, + ventilation_rate=ventilation_rate, + mixing_coefficient=mixing_coefficient, + limits=abiotic_bounds.air_temperature[:2], + time_interval=core_constants.seconds_to_hour, + ) - # Specific humidity of air, [kg kg-1] - specific_humidity_air = abiotic_tools.calculate_specific_humidity( - air_temperature=all_air_temperature, - relative_humidity=relative_humidity, - atmospheric_pressure=atmospheric_pressure, - molecular_weight_ratio_water_to_dry_air=( - core_constants.molecular_weight_ratio_water_to_dry_air - ), - pyrealm_core_constants=pyrealm_core_constants, - ) + # NOTE Advection not implemented as everything is removed with time interval>1h + # and horizontal transfer is not implemented + # advection_rate = ( + # data["wind_speed_ref"].isel(time_index=time_index).to_numpy() + # / np.sqrt(cell_area) + # ) + # advected_fraction = np.clip(advection_rate * time_interval, 0, 1) + # all_air_temperature[0] -=all_air_temperature[0] *advected_fraction + + # Update atmospheric humidity/VPD + # Saturated vapour pressure of air, [kPa] + saturated_vapour_pressure_air = calc_vp_sat( + ta=all_air_temperature, + core_const=pyrealm_core_constants, + ) - # Calculate specific humidity at saturation - mixing_ratio_saturation = ( + # Specific humidity of air, [kg kg-1] + specific_humidity_air = abiotic_tools.calculate_specific_humidity( + air_temperature=all_air_temperature, + relative_humidity=relative_humidity, + atmospheric_pressure=atmospheric_pressure, + molecular_weight_ratio_water_to_dry_air=( core_constants.molecular_weight_ratio_water_to_dry_air - * saturated_vapour_pressure_air - / (atmospheric_pressure - saturated_vapour_pressure_air) - ) - max_specific_humidity = mixing_ratio_saturation / (1 + mixing_ratio_saturation) - - # Update atmospheric humidity variables - new_atmospheric_humidity_vars = energy_balance.update_humidity_vpd( - evapotranspiration=evapotranspiration[ - layer_structure.index_filled_canopy - ].to_numpy(), - soil_evaporation=data["soil_evaporation"].to_numpy(), - saturated_vapour_pressure=saturated_vapour_pressure_air, - specific_humidity=specific_humidity_air, - layer_thickness=above_ground_layer_thickness, - atmospheric_pressure=atmospheric_pressure, - density_air=density_air, - mixing_coefficient=mixing_coefficient, - ventilation_rate=ventilation_rate, - molecular_weight_ratio_water_to_dry_air=( - core_constants.molecular_weight_ratio_water_to_dry_air - ), - dry_air_factor=abiotic_constants.dry_air_factor, - cell_area=cell_area, - limits=(0, max_specific_humidity[0]), # TODO make layer specific - time_interval=core_constants.seconds_to_hour, - ) + ), + pyrealm_core_constants=pyrealm_core_constants, + ) - relative_humidity = new_atmospheric_humidity_vars["relative_humidity"] + # Calculate specific humidity at saturation + mixing_ratio_saturation = ( + core_constants.molecular_weight_ratio_water_to_dry_air + * saturated_vapour_pressure_air + / (atmospheric_pressure - saturated_vapour_pressure_air) + ) + max_specific_humidity = mixing_ratio_saturation / (1 + mixing_ratio_saturation) - # End of loop, write out fluxes and variables + # Update atmospheric humidity variables, integration interval 1 hour + new_atmospheric_humidity_vars = energy_balance.update_humidity_vpd( + evapotranspiration=evapotranspiration[ + layer_structure.index_filled_canopy + ].to_numpy() + / (time_interval * core_constants.seconds_to_hour), + soil_evaporation=data["soil_evaporation"].to_numpy() + / (time_interval * core_constants.seconds_to_hour), + saturated_vapour_pressure=saturated_vapour_pressure_air, + specific_humidity=specific_humidity_air, + layer_thickness=above_ground_layer_thickness, + atmospheric_pressure=atmospheric_pressure, + density_air=density_air, + mixing_coefficient=mixing_coefficient, + ventilation_rate=ventilation_rate, + molecular_weight_ratio_water_to_dry_air=( + core_constants.molecular_weight_ratio_water_to_dry_air + ), + dry_air_factor=abiotic_constants.dry_air_factor, + cell_area=cell_area, + limits=(0, max_specific_humidity[0]), # TODO make layer specific + time_interval=core_constants.seconds_to_hour, + ) + + relative_humidity = new_atmospheric_humidity_vars["relative_humidity"] # Calculate new energy balance and return all fluxes, [W m-2] new_energy_balance_canopy = energy_balance.calculate_energy_balance_residual( @@ -422,16 +404,16 @@ def run_microclimate( air_temperature=air_temperature_canopy, evapotranspiration=evapotranspiration[ layer_structure.index_filled_canopy - ].to_numpy(), + ].to_numpy() + / (time_interval * core_constants.seconds_to_hour), absorbed_radiation_canopy=data["shortwave_absorption"][ layer_structure.index_filled_canopy ].to_numpy(), leaf_emissivity=abiotic_constants.leaf_emissivity, specific_heat_air=specific_heat_air[1:-1], density_air=density_air[1:-1], - density_water=core_constants.density_water, aerodynamic_resistance=aerodynamic_resistance_canopy, - latent_heat_vapourisation=latent_heat_vapourisation[1:-1] / 1000, + latent_heat_vapourisation=latent_heat_vapourisation[1:-1] * 1000, stefan_boltzmann_constant=core_constants.stefan_boltzmann_constant, zero_Celsius=core_constants.zero_Celsius, seconds_to_hour=core_constants.seconds_to_hour, From 282d2b624fde57bc706234b8e75de8f9327cf9bb Mon Sep 17 00:00:00 2001 From: vgro Date: Wed, 19 Nov 2025 11:36:12 +0000 Subject: [PATCH 5/7] comment about integration added --- virtual_ecosystem/models/abiotic/microclimate.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/virtual_ecosystem/models/abiotic/microclimate.py b/virtual_ecosystem/models/abiotic/microclimate.py index 6b9976649..8c4c8e632 100644 --- a/virtual_ecosystem/models/abiotic/microclimate.py +++ b/virtual_ecosystem/models/abiotic/microclimate.py @@ -29,8 +29,10 @@ def run_microclimate( ) -> dict[str, DataArray]: """Run microclimate model. - This function iteratively updates air, soil and canopy temperatures by calculating - the energy balance for each layer. + This function updates air, soil and canopy temperatures by calculating + the energy balance for each layer. We currently make the assumption that over the + time interval of one month, different compartments are in equilibrium. For numerical + stability, the integration interval is 1 hour. ..TODO: Temperatures change between Kelvin and Celsius due to a mix of references, needs to be revisited and converted properly. From 3a10fb520224174e22b0a6ca4bbb301dd9d9bab1 Mon Sep 17 00:00:00 2001 From: vgro Date: Wed, 19 Nov 2025 11:55:14 +0000 Subject: [PATCH 6/7] docstring unit adjusted --- virtual_ecosystem/models/abiotic/energy_balance.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/virtual_ecosystem/models/abiotic/energy_balance.py b/virtual_ecosystem/models/abiotic/energy_balance.py index fa1669023..127aebd36 100644 --- a/virtual_ecosystem/models/abiotic/energy_balance.py +++ b/virtual_ecosystem/models/abiotic/energy_balance.py @@ -288,7 +288,7 @@ def calculate_energy_balance_residual( specific_heat_air: Specific heat capacity of air, [J kg-1 K-1] density_air: Density of air, [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 @@ -414,7 +414,7 @@ def solve_canopy_temperature( density_air: Density of air, [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 From b3ce676631a2c1d7dfc18dc8e568b23af6742ca3 Mon Sep 17 00:00:00 2001 From: vgro Date: Wed, 19 Nov 2025 12:53:34 +0000 Subject: [PATCH 7/7] remove print statement --- virtual_ecosystem/models/abiotic/microclimate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/virtual_ecosystem/models/abiotic/microclimate.py b/virtual_ecosystem/models/abiotic/microclimate.py index 8c4c8e632..e6d15dcbf 100644 --- a/virtual_ecosystem/models/abiotic/microclimate.py +++ b/virtual_ecosystem/models/abiotic/microclimate.py @@ -307,7 +307,7 @@ def run_microclimate( return_fluxes=False, maxiter=10000, ) - print(canopy_temperature) + # Update air temperature based on new canopy and soil temperatures, [C] air_temperature_canopy = energy_balance.update_air_temperature( air_temperature=air_temperature_canopy,