set air temperature as a local leaf attribute

src/hydroshoot/energy.py

@@ -144,6 +144,22 @@ def set_wind_speed(g, meteo, leaf_lbl_prefix='L'):
     return {vid: u for vid in leaves}
 
 
+def set_local_air_temperature(g, meteo, leaf_lbl_prefix='L') -> dict:
+    """Basic model for air temperature at the leaf level, considered equal to air temperature at reference height for
+    all leaves.
+
+    Args:
+        g: a multiscale tree graph object
+        meteo (DataFrame): forcing meteorological variables
+        leaf_lbl_prefix (str): the prefix of the leaf label
+
+    Returns:
+        (dict): keys are leaves vertices ids and their values are all equal to air temperature at reference height
+
+    """
+    return {vid: meteo['Tac'].iloc[0] for vid in get_leaves(g=g, leaf_lbl_prefix=leaf_lbl_prefix)}
+
+
 def _gbH(leaf_length, wind_speed):
     """Computes boundary layer conductance to heat
 
@@ -187,12 +203,11 @@ def calc_heat_boundary_layer_conductance(g, leaf_ids, unit_scene_length):
     return g
 
 
-def calc_leaf_temperature(g, meteo, t_soil, t_sky_eff, leaf_ids, max_iter=100, t_error_crit=0.01, t_step=0.5):
+def calc_leaf_temperature(g, t_soil, t_sky_eff, leaf_ids, max_iter=100, t_error_crit=0.01, t_step=0.5):
     """Computes the temperature of each individual leaf and soil elements.
 
     Args:
         g: a multiscale tree graph object
-        meteo (DataFrame): forcing meteorological variables
         t_soil (float): [°C] soil surface temperature
         t_sky_eff (float): [°C] effective sky temperature
         solo (bool):
@@ -211,7 +226,6 @@ def calc_leaf_temperature(g, meteo, t_soil, t_sky_eff, leaf_ids, max_iter=100, t
     """
 
     temp_sky = utils.celsius_to_kelvin(t_sky_eff)
-    temp_air = utils.celsius_to_kelvin(meteo['Tac'])
     temp_soil = utils.celsius_to_kelvin(t_soil)
 
     t_prev = g.property('Tlc')
@@ -231,6 +245,7 @@ def calc_leaf_temperature(g, meteo, t_soil, t_sky_eff, leaf_ids, max_iter=100, t
             gb_h = mtg_node.properties()['gbH']
             evap = mtg_node.properties()['E']
             t_leaf = utils.celsius_to_kelvin(mtg_node.properties()['Tlc'])
+            temp_air = utils.celsius_to_kelvin(mtg_node.properties()['Tac'])
 
             longwave_gain_from_leaves = ff_leaves * cst.sigma * t_leaf ** 4
 
@@ -264,13 +279,12 @@ def _VineEnergyX(t_leaf):
     return t_new, it
 
 
-def calc_leaf_temperature2(g, meteo, t_soil, t_sky_eff, leaf_ids):
+def calc_leaf_temperature2(g, t_soil, t_sky_eff, leaf_ids):
     """Computes the temperature of each individual leaf and soil elements whilst considering explicitly each other
     leaf temperature energy in a matrix solution using symbolic solver.
 
     Args:
         g: a multiscale tree graph object
-        meteo (DataFrame): forcing meteorological variables
         t_soil (float): [°C] soil surface temperature
         t_sky_eff (float): [°C] effective sky temperature
         leaf_ids (list of int): leaf ids
@@ -282,7 +296,6 @@ def calc_leaf_temperature2(g, meteo, t_soil, t_sky_eff, leaf_ids):
     """
 
     temp_sky = utils.celsius_to_kelvin(t_sky_eff)
-    temp_air = utils.celsius_to_kelvin(meteo['Tac'])
     temp_soil = utils.celsius_to_kelvin(t_soil)
 
     t_prev = g.property('Tlc')
@@ -299,6 +312,7 @@ def calc_leaf_temperature2(g, meteo, t_soil, t_sky_eff, leaf_ids):
         ff_sky = mtg_node.properties()['ff_sky']
         ff_leaves = mtg_node.properties()['ff_leaves']
         ff_soil = mtg_node.properties()['ff_soil']
+        temp_air = utils.celsius_to_kelvin(mtg_node.properties()['Tac'])
 
         gb_h = mtg_node.properties()['gbH']
         evap = mtg_node.properties()['E']

src/hydroshoot/exchange.py

@@ -533,7 +533,7 @@ def calc_gas_exchange_rates(leaf_water_potential, leaf_temperature, carboxylatio
         leaf_ids (list): leaf ids (default None)
         photo_params (dict): values at 25 °C of Farquhar's model
         gs_params (dict): parameters of the stomatal conductance model (model, g0, m0, psi0, D0, n)
-        air_temperature (float): [°C] air temperature
+        air_temperature (dict): [°C] air temperature
         relative_humidity (float): (%) air relative humidity (between 0 and 100 for dry and saturated air, respectively)
         air_co2 (float): [ppm] air CO2 concentration
         atmospheric_pressure (float): [kPa] atmospheric pressure
@@ -581,7 +581,7 @@ def calc_gas_exchange_rates(leaf_water_potential, leaf_temperature, carboxylatio
             dHd=entalpy_deactivation[vid]))
 
         a_n, c_c, c_i, gs = an_gs_ci(
-            air_temperature=air_temperature,
+            air_temperature=air_temperature[vid],
             absorbed_ppfd=absorbed_ppfd[vid],
             relative_humidity=relative_humidity,
             leaf_temperature=temperature_leaf,
@@ -595,13 +595,13 @@ def calc_gas_exchange_rates(leaf_water_potential, leaf_temperature, carboxylatio
             leaf_length=leaf_length[vid],
             wind_speed=wind_speed[vid],
             atm_pressure=atmospheric_pressure,
-            air_temp=air_temperature,
+            air_temp=air_temperature[vid],
             ideal_gas_cst=r)
 
         # Transpiration
         e = calc_transpiration_rate(
             leaf_temperature=temperature_leaf,
-            ea=utils.calc_air_vapor_pressure(air_temperature=air_temperature, relative_humidity=relative_humidity),
+            ea=utils.calc_air_vapor_pressure(air_temperature=air_temperature[vid], relative_humidity=relative_humidity),
             gs=gs,
             gb=gb,
             atm_pressure=atmospheric_pressure)
@@ -616,15 +616,14 @@ def calc_gas_exchange_rates(leaf_water_potential, leaf_temperature, carboxylatio
             boundary_layer_conductance_rate, transpiration_rate)
 
 
-def set_gas_exchange_rates(g, photo_params, gs_params, air_temperature, relative_humidity, air_co2,
+def set_gas_exchange_rates(g, photo_params, gs_params, relative_humidity, air_co2,
                            atmospheric_pressure, E_type2, leaf_ids, rbt=2. / 3.):
     """Sets gas exchange fluxes at the leaf scale analytically.
 
     Args:
         g: a multiscale tree graph object
         photo_params (dict): values at 25 °C of Farquhar's model
         gs_params (dict): parameters of the stomatal conductance model (model, g0, m0, psi0, D0, n)
-        air_temperature (float): [°C] air temperature
         relative_humidity (float): (%) air relative humidity (between 0 and 100 for dry and saturated air, respectively)
         air_co2 (float): [ppm] air CO2 concentration
         atmospheric_pressure (float): [kPa] atmospheric pressure
@@ -656,6 +655,7 @@ def set_gas_exchange_rates(g, photo_params, gs_params, air_temperature, relative
     all_absorbed_ppfd = g.property(E_type2)
     all_leaf_length = g.property('Length')
     all_wind_speed = g.property('u')
+    all_air_temperature = g.property('Tac')
 
     (g.properties()['An'], g.properties()['Ci'], g.properties()['gs'], g.properties()['gb'],
      g.properties()['E']) = calc_gas_exchange_rates(
@@ -672,7 +672,7 @@ def set_gas_exchange_rates(g, photo_params, gs_params, air_temperature, relative
         leaf_ids=leaf_ids,
         photo_params=photo_params,
         gs_params=gs_params,
-        air_temperature=air_temperature,
+        air_temperature=all_air_temperature,
         relative_humidity=relative_humidity,
         air_co2=air_co2,
         atmospheric_pressure=atmospheric_pressure,

src/hydroshoot/initialisation.py

@@ -8,7 +8,7 @@
 
 from hydroshoot import soil
 from hydroshoot.architecture import get_mtg_base, add_soil_surface_mesh, get_leaves
-from hydroshoot.energy import set_form_factors_simplified, set_wind_speed
+from hydroshoot.energy import set_form_factors_simplified, set_wind_speed, set_local_air_temperature
 from hydroshoot.exchange import leaf_Na
 from hydroshoot.io import HydroShootInputs, HydroShootHourlyInputs
 from hydroshoot.irradiance import irradiance_distribution, hsCaribu, set_optical_properties
@@ -173,6 +173,10 @@ def init_hourly(g: MTG, inputs_hourly: HydroShootHourlyInputs, leaf_ppfd: dict,
     g.properties()['u'] = set_wind_speed(
         g=g, meteo=inputs_hourly.weather, leaf_lbl_prefix=params.mtg_api.leaf_lbl_prefix)
 
+    # initiate local air temperature
+    g.properties()['Tac'] = set_local_air_temperature(
+        g=g, meteo=inputs_hourly.weather, leaf_lbl_prefix=params.mtg_api.leaf_lbl_prefix)
+
     if leaf_ppfd is not None:
         diffuse_to_total_irradiance_ratio = leaf_ppfd[g.date]['diffuse_to_total_irradiance_ratio']
         g.properties()['Ei'] = leaf_ppfd[g.date]['Ei']

src/hydroshoot/solver.py

@@ -80,7 +80,7 @@ def solve_interactions(g, meteo, psi_soil, t_soil, t_sky_eff, params, calc_colla
                 # Compute gas-exchange fluxes. Leaf T and Psi are from prev calc loop
                 exchange.set_gas_exchange_rates(
                     g=g, photo_params=photosynthesis_params, gs_params=stomatal_conductance_params,
-                    air_temperature=meteo['Tac'], relative_humidity=meteo['hs'], air_co2=meteo['Ca'],
+                    relative_humidity=meteo['hs'], air_co2=meteo['Ca'],
                     atmospheric_pressure=meteo['Pa'], E_type2=irradiance_type2, leaf_ids=leaf_ids,
                     rbt=turbulence_resistance)
 
@@ -135,7 +135,7 @@ def solve_interactions(g, meteo, psi_soil, t_soil, t_sky_eff, params, calc_colla
             # Compute gas-exchange fluxes. Leaf T and Psi are from prev calc loop
             exchange.set_gas_exchange_rates(
                 g=g, photo_params=photosynthesis_params, gs_params=stomatal_conductance_params,
-                air_temperature=meteo['Tac'], relative_humidity=meteo['hs'], air_co2=meteo['Ca'],
+                relative_humidity=meteo['hs'], air_co2=meteo['Ca'],
                 atmospheric_pressure=meteo['Pa'], E_type2=irradiance_type2, leaf_ids=leaf_ids,
                 rbt=turbulence_resistance)
 
@@ -150,11 +150,11 @@ def solve_interactions(g, meteo, psi_soil, t_soil, t_sky_eff, params, calc_colla
             it += 1
             if params.energy.solo:
                 g.properties()['Tlc'], t_iter = energy.calc_leaf_temperature(
-                    g=g, meteo=meteo, t_soil=t_soil, t_sky_eff=t_sky_eff, leaf_ids=leaf_ids,
+                    g=g, t_soil=t_soil, t_sky_eff=t_sky_eff, leaf_ids=leaf_ids,
                     max_iter=max_iter, t_error_crit=temp_error_threshold, t_step=temp_step)
             else:
                 g.properties()['Tlc'], t_iter = energy.calc_leaf_temperature2(
-                    g=g, meteo=meteo, t_soil=t_soil, t_sky_eff=t_sky_eff, leaf_ids=leaf_ids)
+                    g=g, t_soil=t_soil, t_sky_eff=t_sky_eff, leaf_ids=leaf_ids)
 
             # t_iter_list.append(t_iter)
             t_new = deepcopy(g.property('Tlc'))

test/test_energy.py

@@ -65,9 +65,10 @@ def test_leaf_temperature():
         node.ff_soil = 0.4
         node.gbH = 1
         node.E = 0.
-        node.Tlc = met.Tac.values[0]
+        node.Tac = met.Tac.values[0]
+        node.Tlc = node.Tac
 
-    tleaf, it = calc_leaf_temperature(g, met, tsoil, tsky, get_leaves(g=g, leaf_lbl_prefix='L'))
+    tleaf, it = calc_leaf_temperature(g, tsoil, tsky, get_leaves(g=g, leaf_lbl_prefix='L'))
     assert len(tleaf) == 46
     first = list(tleaf.keys())[0]
     for vid in tleaf: