Feature/energy balance init

docs/source/api/abiotic/energy_balance.md

@@ -0,0 +1,24 @@
+---
+jupytext:
+  cell_metadata_filter: -all
+  formats: md:myst
+  main_language: python
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.13.8
+kernelspec:
+  display_name: vr_python3
+  language: python
+  name: vr_python3
+---
+
+#  API for the {mod}`~virtual_rainforest.models.abiotic.energy_balance` module
+
+```{eval-rst}
+.. automodule:: virtual_rainforest.models.abiotic.energy_balance
+    :autosummary:
+    :members:
+    :special-members: __init__
+```

docs/source/index.md

@@ -110,6 +110,7 @@ team.
   Abiotic Mechanistic Constants <api/abiotic/abiotic_constants.md>
   Abiotic Mechanistic Tools <api/abiotic/abiotic_tools.md>
   Abiotic Mechanistic Wind <api/abiotic/wind.md>
+  Abiotic Mechanistic Energy Balance <api/abiotic/energy_balance.md>
   Hydrology Overview <api/hydrology.md>
   Hydrology Model <api/hydrology/hydrology_model.md>
   Hydrology Above-ground <api/hydrology/above_ground.md>

tests/models/abiotic/test_abiotic_model.py

@@ -44,19 +44,8 @@ def test_abiotic_model_initialization(
     log_check(
         caplog,
         expected_log=(
-            (DEBUG, "abiotic model: required var 'air_temperature' checked"),
-            (DEBUG, "abiotic model: required var 'canopy_height' checked"),
-            (DEBUG, "abiotic model: required var 'layer_heights' checked"),
-            (DEBUG, "abiotic model: required var 'leaf_area_index' checked"),
-            (DEBUG, "abiotic model: required var 'atmospheric_pressure' checked"),
-            (
-                DEBUG,
-                (
-                    "abiotic model: required var 'sensible_heat_flux_topofcanopy'"
-                    " checked"
-                ),
-            ),
-            (DEBUG, "abiotic model: required var 'wind_speed_ref' checked"),
+            (DEBUG, "abiotic model: required var 'air_temperature_ref' checked"),
+            (DEBUG, "abiotic model: required var 'relative_humidity_ref' checked"),
         ),
     )
 
@@ -89,25 +78,14 @@ def test_abiotic_model_initialization_no_data(caplog, dummy_climate_data):
     log_check(
         caplog,
         expected_log=(
-            (ERROR, "abiotic model: init data missing required var 'air_temperature'"),
-            (ERROR, "abiotic model: init data missing required var 'canopy_height'"),
-            (ERROR, "abiotic model: init data missing required var 'layer_heights'"),
-            (ERROR, "abiotic model: init data missing required var 'leaf_area_index'"),
             (
                 ERROR,
-                (
-                    "abiotic model: init data missing required var"
-                    " 'atmospheric_pressure'"
-                ),
+                "abiotic model: init data missing required var 'air_temperature_ref'",
             ),
             (
                 ERROR,
-                (
-                    "abiotic model: init data missing required var"
-                    " 'sensible_heat_flux_topofcanopy'"
-                ),
+                "abiotic model: init data missing required var 'relative_humidity_ref'",
             ),
-            (ERROR, "abiotic model: init data missing required var 'wind_speed_ref'"),
             (ERROR, "abiotic model: error checking required_init_vars, see log."),
         ),
     )
@@ -129,22 +107,8 @@ def test_abiotic_model_initialization_no_data(caplog, dummy_climate_data):
                     "Information required to initialise the abiotic model successfully "
                     "extracted.",
                 ),
-                (DEBUG, "abiotic model: required var 'air_temperature' checked"),
-                (DEBUG, "abiotic model: required var 'canopy_height' checked"),
-                (DEBUG, "abiotic model: required var 'layer_heights' checked"),
-                (DEBUG, "abiotic model: required var 'leaf_area_index' checked"),
-                (
-                    DEBUG,
-                    "abiotic model: required var 'atmospheric_pressure' checked",
-                ),
-                (
-                    DEBUG,
-                    (
-                        "abiotic model: required var 'sensible_heat_flux_topofcanopy'"
-                        " checked"
-                    ),
-                ),
-                (DEBUG, "abiotic model: required var 'wind_speed_ref' checked"),
+                (DEBUG, "abiotic model: required var 'air_temperature_ref' checked"),
+                (DEBUG, "abiotic model: required var 'relative_humidity_ref' checked"),
             ),
             id="default_config",
         ),
@@ -162,22 +126,8 @@ def test_abiotic_model_initialization_no_data(caplog, dummy_climate_data):
                     "Information required to initialise the abiotic model successfully "
                     "extracted.",
                 ),
-                (DEBUG, "abiotic model: required var 'air_temperature' checked"),
-                (DEBUG, "abiotic model: required var 'canopy_height' checked"),
-                (DEBUG, "abiotic model: required var 'layer_heights' checked"),
-                (DEBUG, "abiotic model: required var 'leaf_area_index' checked"),
-                (
-                    DEBUG,
-                    "abiotic model: required var 'atmospheric_pressure' checked",
-                ),
-                (
-                    DEBUG,
-                    (
-                        "abiotic model: required var 'sensible_heat_flux_topofcanopy'"
-                        " checked"
-                    ),
-                ),
-                (DEBUG, "abiotic model: required var 'wind_speed_ref' checked"),
+                (DEBUG, "abiotic model: required var 'air_temperature_ref' checked"),
+                (DEBUG, "abiotic model: required var 'relative_humidity_ref' checked"),
             ),
             id="modified_config_correct",
         ),
@@ -322,6 +272,33 @@ def test_update_abiotic_model(
     )
 
     model.setup()
+
+    np.testing.assert_allclose(
+        model.data["soil_temperature"],
+        DataArray(
+            np.full((15, 3), np.nan),
+            dims=["layers", "cell_id"],
+            coords={
+                "layers": np.arange(0, 15),
+                "layer_roles": (
+                    "layers",
+                    model.layer_roles,
+                ),
+                "cell_id": [0, 1, 2],
+            },
+        ),
+    )
+    np.testing.assert_allclose(
+        model.data["vapour_pressure_deficit_ref"],
+        DataArray(
+            np.full((3, 3), 0.141727),
+            dims=["cell_id", "time_index"],
+            coords={
+                "cell_id": [0, 1, 2],
+            },
+        ),
+    )
+
     model.update(time_index=0)
 
     friction_velocity_exp = np.array(

virtual_rainforest/models/abiotic/__init__.py

@@ -20,6 +20,10 @@
   above- and within-canopy wind profiles for the Virtual Rainforest. These profiles will
   determine the exchange of heat, water, and :math:`CO_{2}` between soil and atmosphere
   below the canopy as well as the exchange with the atmsophere above the canopy.
+
+* The :mod:`~virtual_rainforest.models.abiotic.energy_balance` submodule calculates the
+  energy balance of the Virtual Rainforest. The module returns vertical profiles of air
+  temperature, relative humidity, and soil temperature.
 """  # noqa: D205, D415
 
 from virtual_rainforest.models.abiotic.abiotic_model import AbioticModel  # noqa: F401

virtual_rainforest/models/abiotic/abiotic_model.py

@@ -33,6 +33,8 @@ class as a child of the :class:`~virtual_rainforest.core.base_model.BaseModel` c
 from virtual_rainforest.core.utils import set_layer_roles
 from virtual_rainforest.models.abiotic import wind
 from virtual_rainforest.models.abiotic.constants import AbioticConsts
+from virtual_rainforest.models.abiotic_simple import microclimate
+from virtual_rainforest.models.abiotic_simple.constants import AbioticSimpleConsts
 
 
 class AbioticModel(BaseModel):
@@ -53,17 +55,12 @@ class AbioticModel(BaseModel):
     upper_bound_on_time_scale = "1 day"
     """Longest time scale that abiotic model can sensibly capture."""
     required_init_vars = (
-        ("air_temperature", ("spatial",)),
-        ("canopy_height", ("spatial",)),
-        ("layer_heights", ("spatial",)),
-        ("leaf_area_index", ("spatial",)),
-        ("atmospheric_pressure", ("spatial",)),
-        ("sensible_heat_flux_topofcanopy", ("spatial",)),
-        ("wind_speed_ref", ("spatial",)),
+        ("air_temperature_ref", ("spatial",)),
+        ("relative_humidity_ref", ("spatial",)),
     )
-    """The required variables and axes for the abiotic model"""
+    """The required variables and axes for the abiotic model."""
     vars_updated = ()
-    """Variables updated by the abiotic model"""
+    """Variables updated by the abiotic model."""
 
     def __init__(
         self,
@@ -80,12 +77,8 @@ def __init__(
         # create a list of layer roles
         layer_roles = set_layer_roles(canopy_layers, soil_layers)
 
-        self.data
-        """A Data instance providing access to the shared simulation data."""
         self.layer_roles = layer_roles
         """A list of vertical layer roles."""
-        self.update_interval
-        """The time interval between model updates."""
         self.constants = constants
         """Set of constants for the abiotic model."""
         self.core_constants = core_constants
@@ -129,7 +122,37 @@ def from_config(
         )
 
     def setup(self) -> None:
-        """Function to set up the abiotic model."""
+        """Function to set up the abiotic model.
+
+        This function initializes soil temperature and canopy temperature for all
+        corresponding layers and calculates the reference vapour pressure deficit for
+        all time steps of the simulation. All variables are added directly to the
+        self.data object.
+        """
+
+        # Initialise soil temperature
+        self.data["soil_temperature"] = DataArray(
+            np.full((len(self.layer_roles), len(self.data.grid.cell_id)), np.nan),
+            dims=["layers", "cell_id"],
+            coords={
+                "layers": np.arange(0, len(self.layer_roles)),
+                "layer_roles": ("layers", self.layer_roles),
+                "cell_id": self.data.grid.cell_id,
+            },
+            name="soil_temperature",
+        )
+
+        # Calculate vapour pressure deficit at reference height for all time steps
+        # TODO sort out constants argument in simple abiotic model
+        self.data[
+            "vapour_pressure_deficit_ref"
+        ] = microclimate.calculate_vapour_pressure_deficit(
+            temperature=self.data["air_temperature_ref"],
+            relative_humidity=self.data["relative_humidity_ref"],
+            constants=AbioticSimpleConsts(),
+        ).rename(
+            "vapour_pressure_deficit_ref"
+        )
 
     def spinup(self) -> None:
         """Placeholder function to spin up the abiotic model."""

virtual_rainforest/models/abiotic/energy_balance.py

@@ -0,0 +1,45 @@
+r"""The ``models.abiotic.energy_balance`` module calculates the energy balance for the
+Virtual Rainforest. Given that the time increments of the model are an hour or longer,
+we can assume that below-canopy heat and vapor exchange attain steady state and heat
+storage in the canopy does not need to be simulated explicitly. Under steady-state,
+the balance equation for the leaves in each canopy layer is as follows:
+
+.. math::
+    R_{abs} - R_{em} - H - \lambda E
+    = R_{abs} - \epsilon_{s} \sigma T_{L}^{4} - c_{P}g_{Ha}(T_{L} - T_{A})
+    - \lambda g_{v} \frac {e_{L} - e_{A}}{p_{A}} = 0
+
+where :math:`R_{abs}` is absorbed radiation, :math:`R_{em}` emitted radiation, :math:`H`
+the sensible heat flux, :math:`\lambda E` the latent heat flux, :math:`\epsilon_{s}` the
+emissivity of the leaf, :math:`\sigma` the Stefan-Boltzmann constant, :math:`T_{L}` the
+absolute temperature of the leaf, :math:`T_{A}` the absolute temperature of the air
+surrounding the leaf, :math:`\lambda` the latent heat of vaporisation of water,
+:math:`e_{L}` the effective vapour pressure of the leaf, :math:`e_{A}` the vapor
+pressure of air and :math:`p_{A}` atmospheric pressure. :math:`g_{Ha}` is the heat
+conductance between leaf and atmosphere, :math:`g_{v}` represents the conductance
+for vapor loss from the leaves as a function of the stomatal conductance.
+
+A challenge in solving this equation is the dependency of latent heat and emitted
+radiation on leaf temperature. We use a linearisation approach to solve the equation for
+leaf temperature and air temperature simultaneously, see
+:cite:t:`maclean_microclimc_2021`.
+
+In the soil, heat storage is almost always significant. Thus, Fourier's Law is combined
+with the continuity equation to obtain a time dependant differential equation that
+describes soil temperature as a function of depth and time. A numerical solution can be
+achieved by dividing the soil into discrete layers. Each layer is assigned a node,
+:math:`i`, at depth, :math:`zi`, and with heat storage, :matha;`C_{h_{i}}`, and nodes
+are numbered sequentially downward such that node :math:`i+1` represents the
+node for the soil layer immediately below. Conductivity, :math:`k_{i}`, represents
+conductivity between nodes :math:`i` and :math:`i+1`. The energy balance equation for
+node :math:`i` is then given by
+
+.. math::
+    \kappa_{i}(T_{i+1} - T_{i})- \kappa_{i-1}(T_{i} - T_{i-1})
+    = \frac{C_{h_{i}}(T_{i}^{j+1} - T_{i}^{j})(z_{i+1} - z_{i-1})}{2 \Delta t}
+
+where \Delta t is the time increment, conductance,
+:math:`\kappa_{i}=k_{i}/(z_{i+1} - z_{i})`,  and superscript :math:`j` indicates the
+time at which temperature is determined. This equation can be re-arranged and solved for
+:math:`T_{j+1}` by Gaussian elimination using the Thomas algorithm.
+"""  # noqa: D205, D415

virtual_rainforest/models/abiotic_simple/abiotic_simple_model.py

@@ -86,12 +86,8 @@ def __init__(
         # create a list of layer roles
         layer_roles = set_layer_roles(canopy_layers, soil_layers)
 
-        self.data
-        """A Data instance providing access to the shared simulation data."""
         self.layer_roles = layer_roles
         """A list of vertical layer roles."""
-        self.update_interval
-        """The time interval between model updates."""
         self.constants = constants
         """Set of constants for the abiotic simple model"""