493 hydrology add canopy evaporation

docs/source/virtual_ecosystem/implementation/abiotic_simple_implementation.md

@@ -128,7 +128,7 @@ mystnb:
 display_markdown(
     generate_variable_table(
         'AbioticSimpleModel', 
-        ['vars_populated_by_init', 'vars_populated_by_first_update']
+        ['vars_populated_by_init']
     ), 
     raw=True
 )

tests/conftest.py

@@ -342,7 +342,8 @@ def dummy_climate_data(fixture_core_components):
     spatially_constant = {
         "sensible_heat_flux_soil": 1,
         "latent_heat_flux_soil": 1,
-        "zero_displacement_height": 20.0,
+        "zero_plane_displacement": 20.0,
+        "roughness_length_momentum": 0.5,
         "diabatic_correction_heat_above": 0.1,
         "diabatic_correction_heat_canopy": 1.0,
         "diabatic_correction_momentum_above": 0.1,

tests/models/abiotic/test_abiotic_model.py

@@ -366,15 +366,14 @@ def test_update_abiotic_model(dummy_climate_data, fixture_core_components):
     # Wind speed
     exp_wind_speed = lyr_strct.from_template()
     exp_wind_speed[lyr_strct.index_filled_atmosphere] = np.array(
-        # [0.727122, 0.615474, 0.587838, 0.537028, 0.50198]
         [0.72712164, 0.61547404, 0.57491436, 0.47258967, 0.41466282]
     )[:, None]
     xr.testing.assert_allclose(model.data["wind_speed"], exp_wind_speed)
 
     # Soil temperature
     exp_new_soiltemp = lyr_strct.from_template()
     exp_new_soiltemp[lyr_strct.index_all_soil] = np.array(
-        [  # [20.713167, 20.708367, 20.707833, 20.707833],
+        [
             [20.712458, 20.712457, 20.712456, 20.712456],
             [20.0, 20.0, 20.0, 20.0],
         ]
@@ -384,24 +383,21 @@ def test_update_abiotic_model(dummy_climate_data, fixture_core_components):
     # Leaf vapour conductivity
     exp_gv = lyr_strct.from_template()
     exp_gv[lyr_strct.index_filled_canopy] = np.array(
-        # [0.496563, 0.485763, 0.465142]
         [0.4965627, 0.48056564, 0.43718369]
     )[:, None]
     xr.testing.assert_allclose(model.data["leaf_vapour_conductivity"], exp_gv)
 
     # Air temperature
     exp_air_temp = lyr_strct.from_template()
     exp_air_temp[lyr_strct.index_filled_atmosphere] = np.array(
-        # [30.0, 29.999943, 29.992298, 29.623399, 20.802228]
-        [30.0, 29.99994326, 29.99237944, 29.6604941, 20.80193877]
+        [30.0, 29.999943, 29.992379, 29.660494, 20.801939]
     )[:, None]
     xr.testing.assert_allclose(model.data["air_temperature"], exp_air_temp)
 
     # Canopy temperature
     exp_leaf_temp = lyr_strct.from_template()
     exp_leaf_temp[lyr_strct.index_filled_canopy] = np.array(
-        # [28.787061, 28.290299, 28.15982]
-        [28.78850297, 28.29326228, 28.19789174]
+        [28.788503, 28.293262, 28.197892]
     )[:, None]
     xr.testing.assert_allclose(model.data["canopy_temperature"], exp_leaf_temp)
 
@@ -410,16 +406,14 @@ def test_update_abiotic_model(dummy_climate_data, fixture_core_components):
     # Latent heat flux
     exp_latent_heat = lyr_strct.from_template()
     exp_latent_heat[lyr_strct.index_filled_canopy] = np.array(
-        # [28.07077, 27.568715, 16.006325]
-        [28.07077012, 27.35735709, 14.97729136]
+        [28.07077, 27.357357, 14.977291]
     )[:, None]
     exp_latent_heat[lyr_strct.index_topsoil] = np.array([2.254, 22.54, 225.4, 225.4])
     xr.testing.assert_allclose(model.data["latent_heat_flux"], exp_latent_heat)
 
     # Sensible heat flux
     exp_sens_heat = lyr_strct.from_template()
     exp_sens_heat[lyr_strct.index_flux_layers] = np.array(
-        # [-16.970825, -16.47644, -5.637233, -192.074608]
-        [-16.9708248, -16.26697999, -4.65665595, -192.07460835]
+        [-16.970825, -16.26698, -4.656656, -192.074608]
     )[:, None]
     xr.testing.assert_allclose(model.data["sensible_heat_flux"], exp_sens_heat)

tests/models/abiotic/test_conductivities.py

@@ -127,8 +127,8 @@ def test_calculate_air_heat_conductivity_above(dummy_climate_data):
 
     result = calculate_air_heat_conductivity_above(
         height_above_canopy=dummy_climate_data["layer_heights"][0],
-        zero_displacement_height=(
-            dummy_climate_data["zero_displacement_height"].to_numpy()
+        zero_plane_displacement=(
+            dummy_climate_data["zero_plane_displacement"].to_numpy()
         ),
         canopy_height=dummy_climate_data["layer_heights"][1],
         friction_velocity=dummy_climate_data["friction_velocity"].to_numpy(),

tests/models/abiotic/test_wind.py

@@ -2,6 +2,8 @@
 
 import numpy as np
 import pytest
+import xarray as xr
+from xarray import DataArray
 
 from virtual_ecosystem.core.constants import CoreConsts
 from virtual_ecosystem.models.abiotic.constants import AbioticConsts
@@ -125,13 +127,15 @@ def test_canopy_correction_conditions(
 def test_calculate_mean_mixing_length(dummy_climate_data):
     """Test mixing length with and without vegetation."""
 
+    from virtual_ecosystem.models.abiotic.constants import AbioticConsts
     from virtual_ecosystem.models.abiotic.wind import calculate_mean_mixing_length
 
+    abiotic_consts = AbioticConsts()
     result = calculate_mean_mixing_length(
         canopy_height=dummy_climate_data["layer_heights"][1].to_numpy(),
         zero_plane_displacement=np.array([0.0, 25.312559, 27.58673, 27.58673]),
         roughness_length_momentum=np.array([0.017, 1.4533, 0.9591, 0.9591]),
-        mixing_length_factor=AbioticConsts.mixing_length_factor,
+        mixing_length_factor=abiotic_consts.mixing_length_factor,
     )
 
     np.testing.assert_allclose(
@@ -426,3 +430,88 @@ def test_calculate_wind_profile(
     np.testing.assert_allclose(
         wind_update["wind_speed"], wind_speed_exp, rtol=1e-3, atol=1e-3
     )
+
+
+def test_calculate_aerodynamic_resistance_canopy(dummy_climate_data):
+    """Test calculation of aerodynamic resistance between air and canopy."""
+
+    from virtual_ecosystem.core.core_components import CoreConsts
+    from virtual_ecosystem.models.abiotic.wind import (
+        calculate_aerodynamic_resistance_canopy,
+    )
+
+    result = calculate_aerodynamic_resistance_canopy(
+        canopy_height=dummy_climate_data["layer_heights"][1].to_numpy(),
+        wind_speed_above=dummy_climate_data["wind_speed"][0].to_numpy(),
+        zero_plane_displacement=dummy_climate_data["zero_plane_displacement"],
+        roughness_length_momentum=dummy_climate_data["roughness_length_momentum"],
+        above_canopy_offset=2.0,
+        von_karmans_constant=CoreConsts.von_karmans_constant,
+    )
+
+    exp_result = np.array([1088.610345, 1088.610345, 1088.610345, 1088.610345])
+    np.testing.assert_allclose(result, exp_result, rtol=1e-4, atol=1e-4)
+
+
+def test_update_wind(dummy_climate_data_varying_canopy, fixture_core_components):
+    """Test wind update for abiotic simple model."""
+
+    from virtual_ecosystem.models.abiotic.constants import AbioticConsts
+    from virtual_ecosystem.models.abiotic.wind import update_wind
+
+    data = dummy_climate_data_varying_canopy
+    lyr_strct = fixture_core_components.layer_structure
+
+    microclimate_data = {}
+    microclimate_data["air_temperature"] = data["air_temperature"]
+    microclimate_data["atmospheric_pressure"] = data["atmospheric_pressure"]
+    microclimate_data["sensible_heat_flux"] = data["sensible_heat_flux"]
+
+    result = update_wind(
+        data=data,
+        microclimate_data=microclimate_data,
+        layer_structure=lyr_strct,
+        time_index=0,
+        abiotic_constants=AbioticConsts(),
+        core_constants=fixture_core_components.core_constants,
+    )
+
+    exp_wind_speed = lyr_strct.from_template()
+    exp_wind_speed[lyr_strct.index_filled_atmosphere] = [
+        [0.727122, 0.743643, 0.772241, 0.772241],
+        [0.615474, 0.64478, 0.691463, 0.691463],
+        [0.574914, 0.609452, np.nan, np.nan],
+        [0.47259, np.nan, np.nan, np.nan],
+        [0.414663, 0.544804, 0.635719, 0.635719],
+    ]
+    xr.testing.assert_allclose(result["wind_speed"], exp_wind_speed)
+
+    exp_molar_density = lyr_strct.from_template()
+    exp_molar_density[lyr_strct.index_filled_atmosphere] = [
+        [38.110259, 38.110259, 38.110259, 38.110259],
+        [38.129755, 38.129755, 38.129755, 38.129755],
+        [38.252699, 38.252699, np.nan, np.nan],
+        [38.46472, np.nan, np.nan, np.nan],
+        [39.935316, 39.935316, 39.935316, 39.935316],
+    ]
+    xr.testing.assert_allclose(result["molar_density_air"], exp_molar_density)
+
+    exp_spec_heat = lyr_strct.from_template()
+    exp_spec_heat[lyr_strct.index_filled_atmosphere] = [
+        [29.214, 29.214, 29.214, 29.214],
+        [29.213783, 29.213783, 29.213783, 29.213783],
+        [29.212445, 29.212445, np.nan, np.nan],
+        [29.210245, np.nan, np.nan, np.nan],
+        [29.198443, 29.198443, 29.198443, 29.198443],
+    ]
+    xr.testing.assert_allclose(result["specific_heat_air"], exp_spec_heat)
+
+    exp_d = DataArray(
+        [23.730524, 22.41512, 19.753863, 19.753863], coords={"cell_id": data["cell_id"]}
+    )
+    xr.testing.assert_allclose(result["zero_plane_displacement"], exp_d)
+
+    exp_zm = DataArray(
+        [1.362552, 1.64843, 2.226804, 2.226804], coords={"cell_id": data["cell_id"]}
+    )
+    xr.testing.assert_allclose(result["roughness_length_momentum"], exp_zm)

tests/models/abiotic_simple/test_abiotic_simple_model.py

@@ -15,6 +15,12 @@
 MODEL_VAR_CHECK_LOG = [
     (DEBUG, "abiotic_simple model: required var 'air_temperature_ref' checked"),
     (DEBUG, "abiotic_simple model: required var 'relative_humidity_ref' checked"),
+    (DEBUG, "abiotic_simple model: required var 'atmospheric_pressure_ref' checked"),
+    (DEBUG, "abiotic_simple model: required var 'atmospheric_co2_ref' checked"),
+    (DEBUG, "abiotic_simple model: required var 'leaf_area_index' checked"),
+    (DEBUG, "abiotic_simple model: required var 'layer_heights' checked"),
+    (DEBUG, "abiotic_simple model: required var 'wind_speed_ref' checked"),
+    (DEBUG, "abiotic_simple model: required var 'mean_annual_temperature' checked"),
 ]
 
 
@@ -33,6 +39,7 @@ def test_abiotic_simple_model_initialization(
 ):
     """Test `AbioticSimpleModel` initialization."""
     from virtual_ecosystem.core.base_model import BaseModel
+    from virtual_ecosystem.models.abiotic.constants import AbioticConsts
     from virtual_ecosystem.models.abiotic_simple.abiotic_simple_model import (
         AbioticSimpleModel,
     )
@@ -54,6 +61,7 @@ def test_abiotic_simple_model_initialization(
         assert model.model_name == "abiotic_simple"
         assert repr(model) == "AbioticSimpleModel(update_interval=1209600 seconds)"
         assert model.bounds == AbioticSimpleBounds()
+        assert model.abiotic_constants == AbioticConsts()
 
     # Final check that expected logging entries are produced
     log_check(caplog, expected_log_entries)
@@ -178,6 +186,14 @@ def test_setup(dummy_climate_data_varying_canopy, fixture_core_components):
     model.setup()
 
     exp_soil_temp = lyr_strct.from_template()
+    exp_soil_temp[lyr_strct.index_all_soil] = [
+        [20.712458, 21.317566, 21.922674, 21.922674],
+        [20.0, 20.0, 20.0, 20.0],
+    ]
+    exp_soil_temp[lyr_strct.index_all_soil] = [
+        [20.712458, 21.317566, 21.922674, 21.922674],
+        [20.0, 20.0, 20.0, 20.0],
+    ]
     xr.testing.assert_allclose(model.data["soil_temperature"], exp_soil_temp)
 
     xr.testing.assert_allclose(
@@ -188,20 +204,49 @@ def test_setup(dummy_climate_data_varying_canopy, fixture_core_components):
             coords={"cell_id": [0, 1, 2, 3]},
         ),
     )
+    exp_wind_speed = lyr_strct.from_template()
+    exp_wind_speed[lyr_strct.index_filled_atmosphere] = [
+        [0.727122, 0.743643, 0.772241, 0.772241],
+        [0.615474, 0.64478, 0.691463, 0.691463],
+        [0.574914, 0.609452, np.nan, np.nan],
+        [0.47259, np.nan, np.nan, np.nan],
+        [0.414663, 0.544804, 0.635719, 0.635719],
+    ]
+    xr.testing.assert_allclose(model.data["wind_speed"], exp_wind_speed)
 
-    # Run the update step
-    model.update(time_index=0)
+    exp_sens_heat_flux = lyr_strct.from_template()
+    exp_sens_heat_flux[1] = [0, 0, 0, 0]
+    xr.testing.assert_allclose(model.data["sensible_heat_flux"], exp_sens_heat_flux)
 
     for var in [
+        "soil_temperature",
+        "vapour_pressure_ref",
+        "vapour_pressure_deficit_ref",
+        "vapour_pressure_deficit",
         "air_temperature",
         "relative_humidity",
-        "vapour_pressure_deficit",
-        "soil_temperature",
         "atmospheric_pressure",
         "atmospheric_co2",
+        "wind_speed",
+        "sensible_heat_flux",
+        "molar_density_air",
+        "specific_heat_air",
+        "relative_turbulence_intensity",
+        "attenuation_coefficient",
+        "zero_plane_displacement",
+        "roughness_length_momentum",
+        "mean_mixing_length",
+        "friction_velocity",
+        "diabatic_correction_heat_above",
+        "diabatic_correction_momentum_above",
+        "diabatic_correction_heat_canopy",
+        "diabatic_correction_momentum_canopy",
     ]:
         assert var in model.data
 
+    # Run the update step
+    model.update(time_index=0)
+
     exp_air_temp = lyr_strct.from_template()
     exp_air_temp[lyr_strct.index_filled_atmosphere] = [
         [30.0, 30.0, 30.0, 30.0],
@@ -218,3 +263,13 @@ def test_setup(dummy_climate_data_varying_canopy, fixture_core_components):
         [20.0, 20.0, 20.0, 20.0],
     ]
     xr.testing.assert_allclose(model.data["soil_temperature"], exp_soil_temp)
+
+    exp_wind_speed = lyr_strct.from_template()
+    exp_wind_speed[lyr_strct.index_filled_atmosphere] = [
+        [0.727122, 0.743643, 0.772241, 0.772241],
+        [0.615474, 0.64478, 0.691463, 0.691463],
+        [0.574914, 0.609452, np.nan, np.nan],
+        [0.47259, np.nan, np.nan, np.nan],
+        [0.414663, 0.544804, 0.635719, 0.635719],
+    ]
+    xr.testing.assert_allclose(model.data["wind_speed"], exp_wind_speed)

tests/models/hydrology/test_hydrology_model.py

@@ -277,12 +277,12 @@ def test_setup(
         # Test 2d variables
         expected_2d = {
             "soil_moisture": [
-                [67.0621, 67.0829, 67.05435, 67.04017],
-                [209.8470, 209.8500, 209.8491, 209.8467],
+                [261.477192, 261.597878, 261.469709, 261.296173],
+                [228.479304, 228.557065, 228.459026, 228.394845],
             ],
             "matric_potential": [
-                [-1.532961e07, -1.536408e07, -1.528976e07, -1.53231e07],
-                [-1.250262e03, -1.250131e03, -1.250172e03, -1.250276e3],
+                [-194.444459, -193.68851, -194.543984, -195.407705],
+                [-537.245569, -535.773569, -537.57492, -539.027757],
             ],
         }
 
@@ -299,11 +299,11 @@ def test_setup(
 
         # Test one dimensional variables
         expected_1d = {
-            "vertical_flow": [0.69471, 0.695691, 0.695682, 0.694436],
-            "total_river_discharge": [0, 0, 63361, 20925],
+            "vertical_flow": [1.113488, 1.116363, 1.116338, 1.112685],
+            "total_river_discharge": [0, 0, 64876, 21514],
             "surface_runoff": [0, 0, 0, 0],
             "surface_runoff_accumulated": [0, 0, 0, 0],
-            "soil_evaporation": [345.1148, 344.759928, 345.15422, 344.90802],
+            "soil_evaporation": [12.038902, 12.038902, 12.038902, 12.038902],
         }
 
         for var_name, expected_vals in expected_1d.items():

tests/models/hydrology/test_hydrology_tools.py

@@ -33,7 +33,6 @@ def test_setup_hydrology_input_current_timestep(
     # Check if all variables were created TODO switch back to subcanopy
     var_list = [
         "latent_heat_vapourisation",
-        "molar_density_air",
         "current_precipitation",
         "surface_temperature",
         "surface_humidity",