858 abiotic set initial constants based on data and add references and units

docs/source/refs.bib

@@ -1,3 +1,66 @@
+@article { molders_plant_2005,
+      author = "Nicole Mölders",
+      title = "Plant- and Soil-Parameter-Caused Uncertainty of Predicted Surface Fluxes",
+      journal = "Monthly Weather Review",
+      year = "2005",
+      publisher = "American Meteorological Society",
+      address = "Boston MA, USA",
+      volume = "133",
+      number = "12",
+      doi = "10.1175/MWR3046.1",
+      pages="3498 - 3516",
+      url = "https://journals.ametsoc.org/view/journals/mwre/133/12/mwr3046.1.xml"
+}
+
+@article{su_aerodynamic_2021,
+	title = {Aerodynamic resistance and Bowen ratio explain the biophysical effects of forest cover on understory air and soil temperatures at the global scale},
+	volume = {308-309},
+	issn = {0168-1923},
+	url = {https://www.sciencedirect.com/science/article/pii/S0168192321003014},
+	doi = {https://doi.org/10.1016/j.agrformet.2021.108615},
+	pages = {108615},
+	journal = {Agricultural and Forest Meteorology},
+	author = {Su, Yongxian and Zhang, Chaoqun and Chen, Xiuzhi and Liu, Liyang and Ciais, Philippe and Peng, Jian and Wu, Shengbiao and Wu, Jianping and Shang, Jiali and Wang, Yingping and Yuan, Wenping and Yang, Yuanzhi and Wu, Zhifeng and Lafortezza, Raffaele},
+	year = {2021}
+}
+
+@article{rasimeng_characterization_2020,
+	title = {Characterization of soil thermal properties for design of underground cable routes at the wind farm power plant area,Panyipatan and Pelaihari sub-district, South Kalimantan},
+	volume = {4},
+	url = {https://journal.unhas.ac.id/index.php/geocelebes/article/view/9599},
+	doi = {10.20956/geocelebes.v4i1.9599},
+	pages = {61--69},
+	number = {1},
+	journal = {{Jurnal} {Geocelebs}},
+	author = {Rasimeng, Syamsurijal},
+	year = {2020},
+}
+
+@article{gupta_soilksatdb_2021,
+	title = {{SoilKsatDB}: global database of soil saturated hydraulic conductivity measurements for geoscience applications},
+	volume = {13},
+	url = {https://essd.copernicus.org/articles/13/1593/2021/},
+	doi = {10.5194/essd-13-1593-2021},
+	pages = {1593--1612},
+	number = {4},
+	journal = {Earth System Science Data},
+	author = {Gupta, S. and Hengl, T. and Lehmann, P. and Bonetti, S. and Or, D.},
+	year = {2021},
+}
+
+@article{ma_an_2019,
+	title = {An improved vegetation emissivity scheme for land surface modeling and its impact on snow cover simulations},
+	volume = {53},
+	issn = {1432-0894},
+	url = {https://doi.org/10.1007/s00382-019-04924-9},
+	doi = {10.1007/s00382-019-04924-9},
+	abstract = {This study developed an improved vegetation emissivity scheme for the Community Land Model ({CLM}) version 4.5 to more accurately simulate the effects of vegetation emissivity on snow processes in the Northern Hemisphere over winter and spring. The original scheme of vegetation emissivity in {CLM} produced an unreasonably low vegetation emissivity with a minimum value of around 0.70 in the cold season. Thus, we developed a new vegetation emissivity scheme based on maximum emissivity and leaf and stem area indices of vegetation, which can simulate vegetation emissivity more realistically than the original scheme. Our simulations were driven by the Climatic Research Unit-National Centers for Environmental Prediction ({CRU}-{NCEP}) reanalysis data. Results show that {CLM} with the new scheme produces stronger longwave radiation to the ground surface and generates more solid water drips off vegetation over winter and spring than with the original scheme. Such changes improve snow cover fraction ({SCF}) simulations for the middle and high latitudes in North America, central Eurasia, and the eastern Tibetan Plateau. About 200 and 350 thousand km2 with the {SCF} changes show a better {SCF} simulation with the new scheme over winter and spring, respectively. However, increased errors were found in {SCF} simulations with the new scheme, and further analysis indicates that such errors may be related to biases in the {CLM} forcing variables from the {CRU}-{NCEP} reanalysis data as compared with those from in situ observations. Moreover, the new emissivity scheme decreases total upward longwave radiation and increases surface net radiation and turbulent fluxes. Overall, the improved vegetation emissivity scheme in this study provides an effective tool to generate better understanding of the effects of vegetation on snow at regional scales and gives strong insight into improved land surface process modeling.},
+	pages = {6215--6226},
+	number = {9},
+	journal = {Climate Dynamics},
+	author = {Ma, Xiaogang and Jin, Jiming and Liu, Jian and Niu, Guo-Yue},
+	year = {2019},
+}
 
 @article{hodnett_marked_2002,
 	title = {Marked differences between van Genuchten soil water-retention parameters for temperate and tropical soils: a new water-retention pedo-transfer functions developed for tropical soils},
@@ -13,17 +76,6 @@ @article{hodnett_marked_2002
 	keywords = {Andosols, Ferralsols, Pedo-transfer function, Tropical soils, van Genuchten parameters, Water-retention curves},
 }
 
-@article{su_aerodynamic_2021,
-	title = {Aerodynamic resistance and Bowen ratio explain the biophysical effects of forest cover on understory air and soil temperatures at the global scale},
-	volume = {308-309},
-	issn = {0168-1923},
-	url = {https://www.sciencedirect.com/science/article/pii/S0168192321003014},
-	doi = {https://doi.org/10.1016/j.agrformet.2021.108615},
-	pages = {108615},
-	journal = {Agricultural and Forest Meteorology},
-	author = {Su, Yongxian and Zhang, Chaoqun and Chen, Xiuzhi and Liu, Liyang and Ciais, Philippe and Peng, Jian and Wu, Shengbiao and Wu, Jianping and Shang, Jiali and Wang, Yingping and Yuan, Wenping and Yang, Yuanzhi and Wu, Zhifeng and Lafortezza, Raffaele},
-	year = {2021}
-}
 
 @article{saldarriaga_solar_1991,
 	title={Solar energy conversion efficiencies during succession of a tropical rain forest in Amazonia},

tests/models/abiotic/test_abiotic_model.py

@@ -372,12 +372,12 @@ def test_setup_abiotic_model(dummy_climate_data, fixture_core_components):
         model.update(time_index=0)
 
     expected_soil_temp1 = lyr_strct.from_template()
-    expected_soil_temp1[lyr_strct.index_all_soil] = np.array([18.730802, 19.989525])[
+    expected_soil_temp1[lyr_strct.index_all_soil] = np.array([17.46605, 19.765925])[
         :, None
     ]
     expected_soil_moist = lyr_strct.from_template()
     expected_soil_moist[lyr_strct.index_all_soil] = np.array([5.0, 500])[:, None]
     xr.testing.assert_allclose(
-        model.data["soil_temperature"], expected_soil_temp1, rtol=0.0001
+        model.data["soil_temperature"], expected_soil_temp1, rtol=1e-3
     )
     xr.testing.assert_allclose(model.data["soil_moisture"], expected_soil_moist)

tests/models/abiotic/test_energy_balance.py

@@ -112,7 +112,7 @@ def test_calculate_longwave_emission():
         emissivity=AbioticConsts.soil_emissivity,
         stefan_boltzmann=CoreConsts.stefan_boltzmann_constant,
     )
-    np.testing.assert_allclose(result, np.repeat(320.84384, 3), rtol=1e-04, atol=1e-04)
+    np.testing.assert_allclose(result, np.repeat(381.002069, 3), rtol=1e-04, atol=1e-04)
 
 
 def test_calculate_sensible_heat_flux(

tests/models/abiotic/test_microclimate.py

@@ -28,7 +28,7 @@ def test_run_microclimate(dummy_climate_data, fixture_core_components):
     )
 
     exp_soiltemp = lyr_str.from_template()
-    exp_soiltemp[lyr_str.index_all_soil] = np.array([17.381122, 19.960125])[:, None]
+    exp_soiltemp[lyr_str.index_all_soil] = np.array([15.882118, 19.524066])[:, None]
     np.testing.assert_allclose(
         result["soil_temperature"][lyr_str.index_all_soil],
         exp_soiltemp[lyr_str.index_all_soil],
@@ -38,7 +38,7 @@ def test_run_microclimate(dummy_climate_data, fixture_core_components):
 
     exp_cantemp = lyr_str.from_template()
     exp_cantemp[lyr_str.index_filled_canopy] = np.array(
-        [21.661984, 20.863471, 19.496992]
+        [21.58994, 20.792164, 19.426936]
     )[:, None]
     np.testing.assert_allclose(
         result["canopy_temperature"][lyr_str.index_filled_canopy],
@@ -49,7 +49,7 @@ def test_run_microclimate(dummy_climate_data, fixture_core_components):
 
     exp_airtemp = lyr_str.from_template()
     exp_airtemp[lyr_str.index_filled_atmosphere] = np.array(
-        [30.0, 25.657778, 24.815103, 23.373881, 19.625]
+        [30.0, 25.621995, 24.779698, 23.373881, 18.890155]
     )[:, None]
     np.testing.assert_allclose(
         result["air_temperature"],
@@ -71,13 +71,10 @@ def test_run_microclimate(dummy_climate_data, fixture_core_components):
 
     exp_vp = lyr_str.from_template()
     exp_vp[lyr_str.index_filled_atmosphere] = np.array(
-        [4.233724, 2.643879, 2.517338, 2.313084, 2.305459]
+        [4.233724, 2.632796, 2.50683, 2.303499, 2.210216]
     )[:, None]
     np.testing.assert_allclose(
-        result["vapour_pressure"],
-        exp_vp,
-        rtol=1e-04,
-        atol=1e-04,
+        result["vapour_pressure"], exp_vp, rtol=1e-04, atol=1e-04
     )
 
     exp_vpd = lyr_str.from_template()
@@ -91,12 +88,12 @@ def test_run_microclimate(dummy_climate_data, fixture_core_components):
 
     # Sensible heat flux, canopy only
     exp_shc = lyr_str.from_template()
-    exp_shc[lyr_str.index_filled_canopy] = np.array(
-        [-370.635892, -366.524782, -359.567395]
+    exp_shc[lyr_str.index_flux_layers] = np.array(
+        [-374.205574, -370.058511, -363.040176, -278.772803]
     )[:, None]
     np.testing.assert_allclose(
-        result["sensible_heat_flux"][1:4],
-        exp_shc[1:4],
+        result["sensible_heat_flux"][-2],
+        exp_shc[-2],
         rtol=1e-04,
         atol=1e-04,
     )

tests/models/abiotic_simple/test_abiotic_simple_model.py

@@ -261,9 +261,9 @@ def test_setup(dummy_climate_data_varying_canopy, fixture_core_components):
 
     exp_netrad = lyr_strct.from_template()
     exp_netrad[lyr_strct.index_flux_layers] = [
-        [74.955331, 74.955171, 74.955011, 74.955011],
-        [74.958261, 74.957146, np.nan, np.nan],
-        [74.962932, np.nan, np.nan, np.nan],
-        [77.491651, 77.490632, 77.489522, 77.489522],
+        [74.954871, 74.954709, 74.954548, 74.954548],
+        [74.95783, 74.956705, np.nan, np.nan],
+        [74.96255, np.nan, np.nan, np.nan],
+        [77.490086, 77.488875, 77.487557, 77.487557],
     ]
     xr.testing.assert_allclose(model.data["net_radiation"], exp_netrad)

virtual_ecosystem/models/abiotic/constants.py

@@ -14,17 +14,32 @@
 class AbioticConsts(ConstantsDataclass):
     """Dataclass to store all constants for the `abiotic` model."""
 
-    leaf_emissivity: float = 0.97
-    """Leaf emissivity."""
+    leaf_emissivity: float = 0.98
+    """Leaf emissivity, unitless.
+
+    Leaf emissivity is a measure of how efficiently a leaf emits thermal radiation
+    compared to a perfect blackbody, typically ranging from 0.95 to 0.99. Value for
+    tropical vegetation is taken from :cite:t:`ma_an_2019`.
+    """
 
     leaf_albedo: float = 0.15
-    """Leaf albedo, dimensionless. """
+    """Leaf albedo, unitless.
+
+    Leaf albedo is the fraction of incoming solar radiation that a leaf reflects,
+    typically ranging from 0.12 to 0.18 in tropical forests due to their dark, broadleaf
+    surfaces. Value here is taken from :cite:t:`su_aerodynamic_2021`.
+    """
 
-    bulk_density_soil: float = 1.5
-    """Bulk density of soil."""
+    bulk_density_soil: float = 1.175 * 1000
+    """Bulk density of soil, [kg m-3].
+
+    Bulk density describes the mass of dry soil per unit volume, including both the
+    solid soil particles and the pore spaces between them. Value for average rainforest
+    soil is taken from :cite:t:`gupta_soilksatdb_2021`.
+    """
 
     wind_reference_height: float = 10.0
-    """Reference height for wind speed above the canopy.
+    """Reference height for wind speed above the canopy, [m].
 
     The reference height for horizontal wind is typically 10m above ground compared to
     2m for other atmospheric variables such as temperature and relative humidity. We
@@ -74,7 +89,11 @@ class AbioticConsts(ConstantsDataclass):
     max_ratio_wind_to_friction_velocity: float = 0.3
     """Maximum ratio of wind velocity to friction velocity, dimensionless.
 
-    Implementation and value from :cite:t:`maclean_microclimc_2021`."""
+    The maximum ratio of wind velocity to friction velocity refers to the highest
+    observed or theoretical value of the ratio between the wind speed at a given height
+    and the surface friction velocity (u*), indicating the efficiency of momentum
+    transfer from the atmosphere to the surface. Implementation and value from
+    :cite:t:`maclean_microclimc_2021`."""
 
     drag_coefficient: float = 0.2
     """Drag coefficient, dimensionless.
@@ -96,18 +115,34 @@ class AbioticConsts(ConstantsDataclass):
     surface. Implementation and value from :cite:t:`maclean_microclimc_2021`."""
 
     canopy_temperature_ini_factor: float = 0.01
-    """Factor used to initialise canopy temperature as a function of air temperature and
-    absorbed shortwave radiation."""
+    """Factor used to initialise canopy temperature, dimensionless.
+
+    This factor is used to initialise canopy temperature in the model setup as a
+    function of air temperature and absorbed shortwave radiation."""
 
     light_extinction_coefficient: float = 0.01
-    """Light extinction coefficient for canopy."""
-
-    soil_thermal_conductivity: float = 0.7
-    """Soil thermal conductivity, [W m-1 K-1], :cite:p:`monteith_principles_1990`.
+    """Light extinction coefficient for canopy, unitless.
+
+    The light extinction coefficient for a canopy quantifies how quickly light
+    diminishes as it passes through vegetation, reflecting the canopy's ability to
+    absorb or scatter incoming radiation. This value is only used in the model setup and
+    later derived in the plant model."""
+
+    soil_thermal_conductivity: float = 1.206
+    """Soil thermal conductivity, [W m-1 K-1].
+
+    Soil thermal conductivity is a measure of the soil's ability to conduct heat,
+    influenced by factors such as moisture content, texture, and density. Value is
+    taken from :cite:t:`rasimeng_characterization_2020`.
     """
 
-    specific_heat_capacity_soil: float = 2.7e6
-    """Specific heat capacity of soil, [J kg-1 K-1], :cite:p:`monteith_principles_1990`.
+    specific_heat_capacity_soil: float = 881
+    """Specific heat capacity of soil, [J kg-1 K-1].
+
+    Specific heat capacity of soil is the amount of heat required to raise the
+    temperature of a unit mass of soil by one degree Celsius (or Kelvin), and depends on
+    soil composition, moisture content, and organic matter. Value taken from
+    :cite:t:`molders_plant_2005`.
     """
 
     surface_albedo: float = 0.125
@@ -116,8 +151,13 @@ class AbioticConsts(ConstantsDataclass):
     The value is takes from a study that compares changes in surface albedo before and
     after deforestation in South East Asia :cite:p:`wilson_role_2020`."""
 
-    soil_emissivity: float = 0.8
-    """Soil emissivity, dimensionless."""
+    soil_emissivity: float = 0.95
+    """Soil emissivity, dimensionless.
+
+    Soil emissivity is a measure of how efficiently the soil surface emits thermal
+    radiation compared to a perfect blackbody, with values typically ranging from 0.90
+    to 0.98 depending on soil texture, moisture, and surface roughness. Value taken
+    from :cite:t:`molders_plant_2005`."""
 
     saturated_pressure_slope_parameters: tuple[float, float, float, float] = (
         4098.0,
@@ -132,7 +172,7 @@ class AbioticConsts(ConstantsDataclass):
 
     This term accounts for the proportion of dry air when computing the partitioning
     of total air pressure. It is the complement of the
-    molecular_weight_ratio_water_to_dry_air in core.constants."""
+    `molecular_weight_ratio_water_to_dry_air` in core.constants."""
 
     initial_flux_value: float = 0.001
     """Initial non-zero fill value for energy fluxes, [W m-2]."""