Advanced water and energy balance

Miquel De Caceres

2021-12-16

About this vignette

This document describes how to run a water and energy balance model that uses a more detailed approach for hydraulics and stomatal regulation described in De Cáceres et al. (2021). This document is meant to teach users to run the simulation model within R. All the details of the model design and formulation can be found at https://emf-creaf.github.io/medfatebook/index.html.

Preparing model inputs

Model inputs are explained in greater detail in vignette ‘Simulation inputs’. Here we only review the different steps required to run function spwb().

Soil, vegetation, meteorology and species data

Soil information needs to be entered as a data frame with soil layers in rows and physical attributes in columns. Soil physical attributes can be initialized to default values, for a given number of layers, using function defaultSoilParams():

spar = defaultSoilParams(2)

The soil input for water balance simulation is actually a list of class soil that is created using a function with the same name:

examplesoil = soil(spar)

As explained in the package overview, models included in medfate were primarily designed to be ran on forest inventory plots. Here we use the example object provided with the package:

data(exampleforestMED)
exampleforestMED
## $ID
## [1] "1"
## 
## $patchsize
## [1] 10000
## 
## $treeData
##   Species   N   DBH Height Z50  Z95
## 1     148 168 37.55    800 750 3000
## 2     168 384 14.60    660 750 3000
## 
## $shrubData
##   Species Cover Height Z50  Z95
## 1     165  3.75     80 300 1500
## 
## $herbCover
## [1] 10
## 
## $herbHeight
## [1] 20
## 
## attr(,"class")
## [1] "forest" "list"

Advanced water and energy balance modeling requires daily precipitation, radiation, wind speed, min/max temparatures and relative humitidy as inputs:

data(examplemeteo)
head(examplemeteo)
##            MeanTemperature MinTemperature MaxTemperature Precipitation
## 2001-01-01      3.57668969     -0.5934215       6.287950      4.869109
## 2001-01-02      1.83695972     -2.3662458       4.569737      2.498292
## 2001-01-03      0.09462563     -3.8541036       2.661951      0.000000
## 2001-01-04      1.13866156     -1.8744860       3.097705      5.796973
## 2001-01-05      4.70578690      0.3288287       7.551532      1.884401
## 2001-01-06      4.57036721      0.5461322       7.186784     13.359801
##            MeanRelativeHumidity MinRelativeHumidity MaxRelativeHumidity
## 2001-01-01             78.73709            65.15411           100.00000
## 2001-01-02             69.70800            57.43761            94.71780
## 2001-01-03             70.69610            58.77432            94.66823
## 2001-01-04             76.89156            66.84256            95.80950
## 2001-01-05             76.67424            62.97656           100.00000
## 2001-01-06             89.01940            74.25754           100.00000
##            Radiation WindSpeed WindDirection       PET
## 2001-01-01  12.89251  2.000000           172 1.3212770
## 2001-01-02  13.03079  7.662544           278 2.2185985
## 2001-01-03  16.90722  2.000000           141 1.8045176
## 2001-01-04  11.07275  2.000000           172 0.9200627
## 2001-01-05  13.45205  7.581347           321 2.2914449
## 2001-01-06  12.84841  6.570501           141 1.7255058

Finally, simulations in medfate require a data frame with species parameter values, which we load using defaults for Catalonia (NE Spain):

data("SpParamsMED")

Simulation control

Apart from data inputs, the behaviour of simulation models is controlled using a set of global parameters. The default parameterization is obtained using function defaultControl():

control = defaultControl("Sperry")

To use the complex soil water balance model we must change the values of transpirationMode (to switch from “Granier” to “Sperry”) and soilFunctions (to switch from Saxton’s retention curve, “SX”, to Van Genuchten’s retention curve, “VG”). This is automatically done when calling function defaultControl() with "Sperry" as parameter input.

Water balance input object

A last object is needed before calling simulation functions, called spwbInput. It consists in the compilation of aboveground and belowground parameters and the specification of additional parameter values for each plant cohort. This is done by calling function spwbInput(), but if one has a forest object, the object can be generated more directly using function forest2spwbInput():

x = forest2spwbInput(exampleforestMED, examplesoil, SpParamsMED, control)

The spwbInput object for advanced water and energy balance is similar to that of simple water balance simulations, but contains more elements. Information about the cohort species is found in element cohorts (i.e. code, species and name):

x$cohorts
##         SP              Name
## T1_148 148  Pinus halepensis
## T2_168 168      Quercus ilex
## S1_165 165 Quercus coccifera

Element canopy contains state variables within the canopy:

x$canopy
##    zlow zmid  zup Tair Cair VPair
## 1     0   50  100   NA   NA    NA
## 2   100  150  200   NA   NA    NA
## 3   200  250  300   NA   NA    NA
## 4   300  350  400   NA   NA    NA
## 5   400  450  500   NA   NA    NA
## 6   500  550  600   NA   NA    NA
## 7   600  650  700   NA   NA    NA
## 8   700  750  800   NA   NA    NA
## 9   800  850  900   NA   NA    NA
## 10  900  950 1000   NA   NA    NA
## 11 1000 1050 1100   NA   NA    NA
## 12 1100 1150 1200   NA   NA    NA
## 13 1200 1250 1300   NA   NA    NA
## 14 1300 1350 1400   NA   NA    NA
## 15 1400 1450 1500   NA   NA    NA
## 16 1500 1550 1600   NA   NA    NA
## 17 1600 1650 1700   NA   NA    NA
## 18 1700 1750 1800   NA   NA    NA
## 19 1800 1850 1900   NA   NA    NA
## 20 1900 1950 2000   NA   NA    NA
## 21 2000 2050 2100   NA   NA    NA
## 22 2100 2150 2200   NA   NA    NA
## 23 2200 2250 2300   NA   NA    NA
## 24 2300 2350 2400   NA   NA    NA
## 25 2400 2450 2500   NA   NA    NA
## 26 2500 2550 2600   NA   NA    NA
## 27 2600 2650 2700   NA   NA    NA
## 28 2700 2750 2800   NA   NA    NA

Canopy temperature, water vapour pressure and \(CO_2\) concentration are state variables needed for canopy energy balance. If the canopy energy balance assumes a single canopy layer, the same values will be assumed through the canopy. Variation of within-canopy state variables is modelled if a multi-canopy energy balance is used (see control parameter multiLayerBalance).

As you may already known, element above contains the aboveground structure data that we already know:

x$above
##          H        CR        N   LAI_live LAI_expanded LAI_dead
## T1_148 800 0.6605196 168.0000 1.00643723   1.00643723        0
## T2_168 660 0.6055642 384.0000 0.92661573   0.92661573        0
## S1_165  80 0.8032817 749.4923 0.03965932   0.03965932        0

Belowground parameters can be seen in below:

x$below
##        Z50  Z95 fineRootBiomass coarseRootSoilVolume
## T1_148 750 3000       7041.9790             15.00000
## T2_168 750 3000       3201.6067             19.50000
## S1_165 300 1500        164.5287              0.09375

and in belowLayers:

x$belowLayers
## $V
##                1         2
## T1_148 0.1933638 0.8066362
## T2_168 0.1933638 0.8066362
## S1_165 0.5554389 0.4445611
## 
## $L
##                1        2
## T1_148 2325.6624 4162.854
## T2_168 2609.6105 4606.209
## S1_165  571.4173 1185.629
## 
## $VGrhizo_kmax
##                1        2
## T1_148  42285017 62895956
## T2_168  51808686 79507938
## S1_165 253152670 73383323
## 
## $VCroot_kmax
##                1         2
## T1_148 0.4050231 0.9439254
## T2_168 1.2815989 3.0289117
## S1_165 2.5460170 0.9821121
## 
## $Wpool
##        1 2
## T1_148 1 1
## T2_168 1 1
## S1_165 1 1
## 
## $RhizoPsi
##        1 2
## T1_148 0 0
## T2_168 0 0
## S1_165 0 0

The spwbInputobject also includes cohort parameter values for several kinds of traits. For example, plant anatomy parameters are described in paramsAnatomy:

x$paramsAnatomy
##        Hmax Hmed    Al2As      SLA LeafWidth LeafDensity WoodDensity
## T1_148 1900  850 1716.262 5.348269  0.700000   0.3237865   0.6064556
## T2_168 1300  550 2481.709 6.817833  1.767436   0.5272863   0.8743066
## S1_165  250   80 6180.357 5.027918  1.376108   0.3871262   0.1713280
##        FineRootDensity conduit2sapwood      SRL RLD     r635
## T1_148       0.3237865       0.9236406 3172.572  10 1.964226
## T2_168       0.5272863       0.6238125 4066.803  10 1.805872
## S1_165       0.3871262       0.6238125 4066.803  10 2.289452

Parameters related to plant transpiration and photosynthesis can be seen in paramsTranspiration:

x$paramsTranspiration
##             Gswmin    Gswmax  Vmax298  Jmax298 Kmax_stemxylem Kmax_rootxylem
## T1_148 0.002031250 0.1940833 62.81613 128.8272      0.1500000       0.600000
## T2_168 0.006003126 0.2007222 39.24539 103.0690      0.7735834       3.094334
## S1_165 0.010455247 0.2830167 62.96065 119.4442      0.2900000       1.160000
##        VCleaf_kmax VCleaf_c  VCleaf_d VCstem_kmax VCstem_c  VCstem_d
## T1_148    7.166471 1.497782 -2.142745    1.028343 7.458883 -4.481392
## T2_168    7.354306 1.331076 -1.777291    4.716717 2.558358 -5.103492
## S1_165    9.579077 1.844224 -3.030130    3.117593 3.537784 -4.126512
##        VCroot_kmax VCroot_c  VCroot_d VGrhizo_kmax Plant_kmax
## T1_148    1.348948 1.916254 -3.240714    105180972  0.5395794
## T2_168    4.310511 1.957824 -3.366992    131316623  1.7242042
## S1_165    3.528129 1.760307 -2.797092    326535993  1.4112516

Finally, parameters related to pressure-volume curves and water storage capacity of leaf and stem organs are in paramsWaterStorage:

x$paramsWaterStorage
##          LeafPI0   LeafEPS LeafAF     Vleaf    StemPI0   StemEPS    StemAF
## T1_148 -1.591429  8.918571 0.3525 0.4560549 -2.0028552 13.185187 0.9236406
## T2_168 -1.483333 19.260000 0.1700 0.1829250 -3.1171155 41.465622 0.6238125
## S1_165 -2.370000 17.230000 0.2400 0.3846097 -0.1927245  1.218936 0.6238125
##         Vsapwood
## T1_148 2.8256658
## T2_168 1.1495998
## S1_165 0.1150416

Finally, remember that one can play with plant-specific parameters for soil water balance (instead of using species-level values) by modifying manually the parameter values in this object.

Static analysis of submodels

Before using the advanced water and energy balance model, is important to understand the parameters that influence the different sub-models. Package medfate provides low-level functions corresponding to sub-models (light extinction, hydraulics, transpiration, photosynthesis…). In addition, there are several high-level plotting functions that allow examining several aspects of these processes.

Vulnerability curves

Given a spwbInput object, we can use function hydraulics_vulnerabilityCurvePlot() to inspect vulnerability curves (i.e. how hydraulic conductance of a given segment changes with the water potential) for each plant cohort and each of the different segments of the soil-plant hydraulic network: rhizosphere, roots, stems and leaves:

hydraulics_vulnerabilityCurvePlot(x, type="leaf")

hydraulics_vulnerabilityCurvePlot(x, type="stem")

hydraulics_vulnerabilityCurvePlot(x, type="root")

hydraulics_vulnerabilityCurvePlot(x, examplesoil, type="rhizo")

The maximum values and shape of vulnerability curves for leaves and stems are regulated by parameters in paramsTranspiration. Roots have vulnerability curve parameters in the same data frame, but maximum conductance values need to be specified for each soil layer and are given in belowLayers$VCroot_kmax. Note that the last call to hydraulics_vulnerabilityCurvePlot() includes a soil object. This is because the van Genuchten parameters that define the shape of the vulnerability curve for the rhizosphere are stored in this object. Maximum conductance values in the rhizosphere are given in belowLayers$VGrhizo_kmax.

Supply functions

The vulnerability curves conformng the hydraulic network are used in the model to build the supply function, which relates water flow (i.e. transpiration) with the drop of water potential along the whole hydraulic pathway. The supply function contains not only these two variables, but also the water potential of intermediate nodes in the the hydraulic network. Function hydraulics_supplyFunctionPlot() can be used to inspect any of this variables:

hydraulics_supplyFunctionPlot(x, examplesoil, type="E")

hydraulics_supplyFunctionPlot(x, examplesoil, type="ERhizo")

hydraulics_supplyFunctionPlot(x, examplesoil, type="dEdP")

hydraulics_supplyFunctionPlot(x, examplesoil, type="StemPsi")

Calls to hydraulics_supplyFunctionPlot() always need both a spwbInput object and a soil object. The soil moisture state (i.e. its water potential) is the starting point for the calculation of the supply function, so different curves will be obtained for different values of soil moisture.

Stomatal regulation and photosynthesis

The soil water balance model determines stomatal conductance and transpiration separately for sunlit and shade leaves. Stomatal conductance is determined after building a photosynthesis function corresponding to the supply function and finding the value of stomatal conductance that maximizes carbon revenue while avoiding hydraulic damage (a profit-maximization approach). Given a meteorological and soil inputs and a chosen day and timestep, function transp_stomatalRegulationPlot() allows displaying the supply and photosynthesis curves for sunlit and shade leaves, along with an indication of the values corresponding to the chosen stomatal aperture:

d = 100
transp_stomatalRegulationPlot(x, examplemeteo, day = d, timestep=12,
                              latitude = 41.82592, elevation = 100, type="E")

transp_stomatalRegulationPlot(x, examplemeteo, day = d, timestep=12,
                              latitude = 41.82592, elevation = 100, type="An")

transp_stomatalRegulationPlot(x, examplemeteo, day = d, timestep=12,
                              latitude = 41.82592, elevation = 100, type="Gsw")

transp_stomatalRegulationPlot(x, examplemeteo, day = d, timestep=12,
                              latitude = 41.82592, elevation = 100, type="T")

transp_stomatalRegulationPlot(x, examplemeteo, day = d, timestep=12,
                              latitude = 41.82592, elevation = 100, type="VPD")

Pressure volume curves

moisture_pressureVolumeCurvePlot(x, segment="leaf", fraction="symplastic")

moisture_pressureVolumeCurvePlot(x, segment="leaf", fraction="apoplastic")

moisture_pressureVolumeCurvePlot(x, segment="stem", fraction="symplastic")

moisture_pressureVolumeCurvePlot(x, segment="stem", fraction="apoplastic")

Water balance for a single day

Running the model

Soil water balance simulations will normally span periods of several months or years, but since the model operates at a daily and subdaily temporal scales, it is possible to perform soil water balance for one day only. This is done using function spwb_day(). In the following code we select the same day as before from the meteorological input data and perform soil water balance for that day only:

sd1<-spwb_day(x, rownames(examplemeteo)[d],  
             examplemeteo$MinTemperature[d], examplemeteo$MaxTemperature[d], 
             examplemeteo$MinRelativeHumidity[d], examplemeteo$MaxRelativeHumidity[d], 
             examplemeteo$Radiation[d], examplemeteo$WindSpeed[d], 
             latitude = 41.82592, elevation = 100, 
             slope= 0, aspect = 0, prec = examplemeteo$Precipitation[d])

The output of spwb_day() is a list with several elements:

names(sd1)
##  [1] "cohorts"          "WaterBalance"     "EnergyBalance"    "Soil"            
##  [5] "Stand"            "Plants"           "RhizoPsi"         "SunlitLeaves"    
##  [9] "ShadeLeaves"      "ExtractionInst"   "PlantsInst"       "SunlitLeavesInst"
## [13] "ShadeLeavesInst"  "LightExtinction"  "LWRExtinction"    "CanopyTurbulence"

Water balance output

Element WaterBalance contains the soil water balance flows of the day (precipitation, infiltration, transpiration, …)

sd1$WaterBalance
##                     PET                    Rain                    Snow 
##               5.0233468               0.0000000               0.0000000 
##                 NetRain                Snowmelt                   Runon 
##               0.0000000               0.0000000               0.0000000 
##            Infiltration                  Runoff            DeepDrainage 
##               0.0000000               0.0000000               0.0000000 
##         SoilEvaporation         PlantExtraction           Transpiration 
##               0.5000000               0.5895333               0.5895333 
## HydraulicRedistribution 
##               0.0000000

And Soil contains water evaporated from each soil layer, water transpired from each soil layer and the final soil water potential:

sd1$Soil
##   SoilEvaporation HydraulicInput HydraulicOutput PlantExtraction         psi
## 1    4.999998e-01              0       0.1815242       0.1815242 -0.03478846
## 2    1.529512e-07              0       0.4080091       0.4080091 -0.03355136

Soil and canopy energy balance

Element EnergyBalance contains subdaily variation in atmosphere, canopy and soil temperatures, as well as canopy and soil energy balance components.

names(sd1$EnergyBalance)
## [1] "Temperature"         "CanopyEnergyBalance" "SoilEnergyBalance"  
## [4] "TemperatureLayers"   "VaporPressureLayers"

Package medfate provides a plot function for objects of class spwb_day that can be used to inspect the results of the simulation. We use this function to display subdaily dynamics in plant, soil and canopy variables. For example, we can use it to display temperature variations (only the temperature of the topmost soil layer is drawn):

plot(sd1, type = "Temperature")

plot(sd1, type = "CanopyEnergyBalance")

plot(sd1, type = "SoilEnergyBalance")

Plant output

Element Plants contains output values by plant cohort. Several output variables can be inspected in this element.

sd1$Plants
##               LAI Extraction Transpiration GrossPhotosynthesis
## T1_148 1.00643723 0.25942064    0.25942064          2.07528657
## T2_168 0.92661573 0.31348103    0.31348103          1.38992465
## S1_165 0.03965932 0.01663162    0.01663162          0.06123801
##        NetPhotosynthesis    RootPsi    StemPsi      StemPLC LeafPsiMin
## T1_148        1.93248750 -0.4977087 -0.7991837 4.720222e-07 -1.2299352
## T2_168        1.29836524 -0.2416643 -0.3366812 2.278513e-04 -0.5786835
## S1_165        0.05553232 -0.3176123 -0.4767654 1.095076e-04 -0.7466040
##         LeafPsiMax      dEdP        DDS   StemRWC   LeafRWC StemSympRWC
## T1_148 -0.04848757 0.3518063 0.03466548 0.9983121 0.9237990   0.9779010
## T2_168 -0.05059580 1.1248855 0.03406221 0.9986545 0.9792477   0.9968011
## S1_165 -0.06542421 0.9422049 0.01151417 0.9061952 0.9821248   0.7508251
##        LeafSympRWC  WaterBalance
## T1_148   0.9527228  2.450297e-17
## T2_168   0.9890906 -9.540979e-18
## S1_165   0.9836975  5.556536e-19

While Plants contains one value per cohort and variable that summarizes the whole simulated day, information by disaggregated by time step can be accessed in PlantsInst. Moreover, we can use function plot.spwb_day() to draw plots of sub-daily variation per species of plant transpiration per ground area (L·m\(^{-2}\)), transpiration per leaf area (also in L·m\(^{-2}\)), plant net photosynthesis (in g C·m\(^{-2}\)), and plant water potential (in MPa):

plot(sd1, type = "PlantTranspiration", bySpecies = T)

plot(sd1, type = "TranspirationPerLeaf", bySpecies = T)

plot(sd1, type = "NetPhotosynthesis", bySpecies = T)

plot(sd1, type = "LeafPsiAverage", bySpecies = T)

Output for sunlit and shade leaves

The model distinguishes between sunlit and shade leaves for stomatal regulation. Static properties of sunlit and shade leaves, for each cohort, can be accessed via:

sd1$SunlitLeaves
##                LAI  Vmax298   Jmax298 LeafPsiMin  LeafPsiMax      GSWMin
## T1_148 0.434534247 49.42817 101.37032 -1.6868505 -0.04848757 0.002051063
## T2_168 0.283067367 30.15910  79.20595 -0.9826431 -0.04629099 0.006121761
## S1_165 0.008054994 44.08017  83.62555 -1.6063974 -0.06542421 0.010497379
##            GSWMax  TempMin  TempMax
## T1_148 0.08125225 3.322541 14.80512
## T2_168 0.07887225 3.068057 20.71198
## S1_165 0.12101295 3.492115 20.80481
sd1$ShadeLeaves
##               LAI  Vmax298  Jmax298 LeafPsiMin  LeafPsiMax      GSWMin
## T1_148 0.57190298 41.41563 84.93772 -0.9771766 -0.04848757 0.002037635
## T2_168 0.64354836 27.07762 71.11315 -0.3911507 -0.05059580 0.006109172
## S1_165 0.03160432 44.08017 83.62555 -0.5454198 -0.06542421 0.010497379
##            GSWMax  TempMin  TempMax
## T1_148 0.06385147 3.065392 11.02684
## T2_168 0.06601539 2.975828 10.84967
## S1_165 0.09701438 3.097101 10.77719

Instantaneous values are also stored for sunlit and shade leaves. We can also use the plot function for objects of class spwb_day to draw instantaneous variations in temperature for sunlit and shade leaves:

plot(sd1, type = "LeafTemperature", bySpecies=TRUE)

Note that sunlit leaves of some species reach temperatures higher than the canopy. We can also plot variations in instantaneous gross and net photosynthesis rates:

plot(sd1, type = "LeafGrossPhotosynthesis", bySpecies=TRUE)

plot(sd1, type = "LeafNetPhotosynthesis", bySpecies=TRUE)

Or variations in stomatal conductance:

plot(sd1, type = "LeafStomatalConductance", bySpecies=TRUE)