Basic water balance simulation

Miquel De Caceres

2020-11-05

About this vignette

The present document describes how to run the soil plant water balance model described in De Cáceres et al. (2015) using package medfate. The document indicates how to prepare the inputs, use the simulation functions and inspect the outputs. All the details of the model design and formulation can be found at https://vegmod.ctfc.cat/software/medfate. Because it introduces many basic features of package medfate, this document should be read before addressing advanced topics of water balance simulations or growth simulations.

Preparing model inputs

Any forest water balance model needs information on climate, vegetation and soils of the forest stand to be simulated. Moreover, since the soil water balance in medfate differentiates between species, species-specific parameters are also needed. In this subsection we explain the different steps (including optional pathways) to prepare all the data needed to run function spwb().

Soils

Required soil data

Simulation models in medfate require information on the physical attributes of soil, namely soil depth, texture, bulk density and rock fragment content. Soil information needs to be entered as a data frame with soil layers in rows and physical attributes in columns. The model accepts one to five soil layers with arbitrary widths. Soil physical attributes can be initialized to default values, for a given number of layers, using function defaultSoilParams():

spar = defaultSoilParams(2)
print(spar)
##   widths clay sand om  bd rfc
## 1    300   25   25 NA 1.5  25
## 2    700   25   25 NA 1.5  45

where widths are soil layer widths in mm; clay and sand are the percentage of clay and sand, in percent of dry weight, om stands for organic matter, bd is bulk density (in g·cm\(^{-3}\)) and rfc the percentage of rock fragments. Because soil properties vary strongly at fine spatial scales, ideally soil physical attributes should be measured on samples taken at the forest stand to be simulated. For those users lacking such data, soil properties modelled at larger scales are available via soilgrids.org (see function soilgridsParams()).

Soil input object

The soil input for function spwb() is actually a list of class soil that is created using a function with the same name:

examplesoil = soil(spar)
class(examplesoil)
## [1] "soil" "list"

The list contains the following elements:

names(examplesoil)
##  [1] "SoilDepth"    "W"            "SWE"          "Temp"         "Ksoil"       
##  [6] "Gsoil"        "dVec"         "sand"         "clay"         "om"          
## [11] "VG_alpha"     "VG_n"         "VG_theta_res" "VG_theta_sat" "Ksat"        
## [16] "Kdrain"       "macro"        "bd"           "rfc"

In addition to the physical soil description, this object contains soil parameters needed for soil water balance simulations. For example, macro specifies the macroporosity of each layer; Gsoil and Ksoil are parameters needed to model the process of water infiltration into the soil. The meaning of all elements in the soil object can be found in the help page for function soil().

Importantly, the soil object is used to store the degree of moisture of each soil layer. In particular, element W contains the state variable that represents moisture content - the proportion of moisture relative to field capacity - which is normally initialized to 1 for each layer:

examplesoil$W
## [1] 1 1

It is important to remember that, unlike normal objects in R, any water balance simulation will modify the moisture values in the input soil object. That is, the state of the soil at the end of the simulated process (i.e. W) will be stored. Hence, one can use the same object to simulate water balance sequentially and the final state of one simulation is the initial state of the next.

Water retention curves

At any time, one can show the characteristics and status of the soil object using its print function:

print(examplesoil, model = "SX")
## Soil depth (mm): 1000 
## 
## Layer  1  [ 0  to  300 mm ] 
##     clay (%): 25 silt (%): 50 sand (%): 25 organic matter (%): NA [ Silt loam ]
##     Rock fragment content (%): 25 Macroporosity (%): 5 
##     Theta WP (%): 14 Theta FC (%): 30 Theta SAT (%): 49 Theta current (%) 30 
##     Vol. WP (mm): 32 Vol. FC (mm): 68 Vol. SAT (mm): 111 Vol. current (mm): 68 
##     Temperature (Celsius): NA 
## 
## Layer  2  [ 300  to  1000 mm ] 
##     clay (%): 25 silt (%): 50 sand (%): 25 organic matter (%): NA [ Silt loam ]
##     Rock fragment content (%): 45 Macroporosity (%): 5 
##     Theta WP (%): 14 Theta FC (%): 30 Theta SAT (%): 49 Theta current (%) 30 
##     Vol. WP (mm): 55 Vol. FC (mm): 117 Vol. SAT (mm): 190 Vol. current (mm): 117 
##     Temperature (Celsius): NA 
## 
## Total soil saturated capacity (mm): 300 
## Total soil water holding capacity (mm): 185 
## Total soil extractable water (mm): 116 
## Total soil current Volume (mm): 185 
## 
## Snow pack water equivalent (mm): 0 
## Soil water table depth (mm): 1000

The modelled moisture content of the soil depends on the water retention curve used to represent the relationship between soil volumetric water content (\(\theta\); %) and soil water potential (\(\Psi\); MPa). By default the Saxton (model = "SX") equations are used to model the water retention curve, but the user may choose to follow Van Genuchten - Mualem equations, which will give slightly different values for the same texture:

print(examplesoil, model="VG")
## Soil depth (mm): 1000 
## 
## Layer  1  [ 0  to  300 mm ] 
##     clay (%): 25 silt (%): 50 sand (%): 25 organic matter (%): NA [ Silt loam ]
##     Rock fragment content (%): 25 Macroporosity (%): 5 
##     Theta WP (%): 13 Theta FC (%): 29 Theta SAT (%): 42 Theta current (%) 29 
##     Vol. WP (mm): 30 Vol. FC (mm): 64 Vol. SAT (mm): 95 Vol. current (mm): 64 
##     Temperature (Celsius): NA 
## 
## Layer  2  [ 300  to  1000 mm ] 
##     clay (%): 25 silt (%): 50 sand (%): 25 organic matter (%): NA [ Silt loam ]
##     Rock fragment content (%): 45 Macroporosity (%): 5 
##     Theta WP (%): 13 Theta FC (%): 30 Theta SAT (%): 42 Theta current (%) 30 
##     Vol. WP (mm): 49 Vol. FC (mm): 117 Vol. SAT (mm): 163 Vol. current (mm): 117 
##     Temperature (Celsius): NA 
## 
## Total soil saturated capacity (mm): 258 
## Total soil water holding capacity (mm): 181 
## Total soil extractable water (mm): 118 
## Total soil current Volume (mm): 181 
## 
## Snow pack water equivalent (mm): 0 
## Soil water table depth (mm): 1000

While Saxton equations use texture and organic matter as inputs, the Van Genuchten-Mualem equations need other parameters, which are estimated using pedotransfer functions and their names start with VG_ (two alternative options are provided in function soil to estimate Van Genuchten parameters). The following code calls function soil_retentionCurvePlot() to illustrate the difference between the two water retention models in this soil:

soil_retentionCurvePlot(examplesoil, model="both")

Low-level functions, such as soil_psi2thetaSX() and soil_psi2thetaVG() (and their counterparts soil_theta2psiSX() and soil_theta2psiVG()), can be used to calculate volumetric soil moisture from the water potential (and viceversa) using the two models. When simulating soil water balance, the user can choose among the two models (see control parameters below).

Species-specific parameters

Simulation models in medfate require a data frame with species parameter values. The package provides a default data set of parameter values for 89 Mediterranean species (rows), resulting from bibliographic search, fit to empirical data or expert-based guesses:

data("SpParamsMED")

These species commonly occur in the Spanish forest inventory of Catalonia, but may not be sufficient for other areas. A large number of parameters (columns) can be found in SpParamsMED:

names(SpParamsMED)
##  [1] "Name"               "IFNcodes"           "SpIndex"           
##  [4] "Group"              "Order"              "Family"            
##  [7] "GrowthForm"         "Hmed"               "Hmax"              
## [10] "Z50"                "Z95"                "a_ash"             
## [13] "a_bsh"              "b_bsh"              "cr"                
## [16] "a_fbt"              "b_fbt"              "c_fbt"             
## [19] "d_fbt"              "a_cr"               "b_1cr"             
## [22] "b_2cr"              "b_3cr"              "c_1cr"             
## [25] "c_2cr"              "a_cw"               "b_cw"              
## [28] "PhenologyType"      "LeafDuration"       "Sgdd"              
## [31] "SLA"                "LeafDensity"        "WoodDensity"       
## [34] "FineRootDensity"    "r635"               "pDead"             
## [37] "Al2As"              "LeafWidth"          "SRL"               
## [40] "RLD"                "minFMC"             "maxFMC"            
## [43] "LeafPI0"            "LeafEPS"            "LeafAF"            
## [46] "StemPI0"            "StemEPS"            "StemAF"            
## [49] "LigninPercent"      "ParticleDensity"    "LeafLitterFuelType"
## [52] "Flammability"       "SAV"                "HeatContent"       
## [55] "gammaSWR"           "alphaSWR"           "kPAR"              
## [58] "g"                  "Psi_Extract"        "Psi_Critic"        
## [61] "WUE"                "pRootDisc"          "Gwmin"             
## [64] "Gwmax"              "VCleaf_kmax"        "VCleaf_c"          
## [67] "VCleaf_d"           "Kmax_stemxylem"     "VCstem_c"          
## [70] "VCstem_d"           "Kmax_rootxylem"     "VCroot_c"          
## [73] "VCroot_d"           "Narea"              "Vmax298"           
## [76] "Jmax298"            "WoodC"              "RGRsapwoodmax"     
## [79] "fHDmin"             "fHDmax"

Not all parameters are needed for all models. The user can find parameter definitions in the help page of this data set. However, to fully understand the role of parameters in the model, the user should read the details of model design and formulation (https://vegmod.ctfc.cat/software/medfate).

Vegetation

Forest plot data

Models included in medfate were primarily designed to be ran on forest inventory plots. In this kind of data, the vegetation of a sampled area is described in terms of woody plants (trees and shrubs) along with their size and species identity. Forest plots in medfate are assumed to be in a format that follows closely the Spanish forest inventory. Each forest plot is represented in an object of class forest, a list that contains several elements. Among them, the most important items are two data frames, treeData (for trees) and shrubData for shrubs:

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

Trees are expected to be primarily described in terms of species, diameter (DBH) and height, whereas shrubs are described in terms of species, percent cover and mean height.

Aboveground and belowground data

Because the forest plot format is rather specific, simulation functions in medfate allow starting in a more general way using two data frames, one with aboveground information (i.e. the leave area and size of plants) and the other with belowground information (i.e. root distribution). The aboveground data frame does not distinguish between trees and shrubs. It includes, for each plant cohort to be considered in rows, its species identity, height and leaf area index (LAI). While users can build their input data themselves, we use function forest2aboveground() on the object exampleforestMED to show how should the data look like:

above = forest2aboveground(exampleforestMED, SpParamsMED)
above
##       SP        N   DBH Cover   H        CR  LAI_live LAI_expanded LAI_dead
## T1_54 54 168.0000 37.55    NA 800 0.7150421 0.8167012    0.8167012        0
## T2_68 68 384.0000 14.60    NA 660 0.6055642 0.7977952    0.7977952        0
## S1_65 65 773.9235    NA  3.75  80 0.9738889 0.0653033    0.0653033        0

Note that the call to forest2aboveground() included species parameters, because species-specific values are needed to calculate leaf area from tree diameters or shrub cover. Columns N, DBH and Cover are required for simulating growth, but not for soil water balance, which only requires columns SP, H (in cm), CR (i.e. the crown ratio), LAI_live, LAI_expanded and LAI_dead. Here plant cohorts are given unique codes that tell us whether they correspond to trees or shrubs, but the user can use other row identifiers as long as they are unique. In practice, the user only needs to worry to calculate the values for LAI_live. LAI_live and LAI_expanded can contain the same LAI values, and LAI_dead is normally zero. This is so because models update LAI_expanded and LAI_dead according to the leaf phenology of species.

Aboveground leaf area distribution (with or without distinguishing among cohorts) can be examined by calling function vprofile_leafAreaDensity():

vprofile_leafAreaDensity(above, byCohorts = F)

vprofile_leafAreaDensity(above, byCohorts = T)

Regarding belowground information, we need vectors with depths corresponding to 50% and 95% of fine roots, which we simply concatenate from our forest data:

Z50 = c(exampleforestMED$treeData$Z50, exampleforestMED$shrubData$Z50)
Z95 = c(exampleforestMED$treeData$Z95, exampleforestMED$shrubData$Z95)

These parameters specify a continuous distribution of fine roots. Users can visually inspect the distribution of fine roots of forest objects by calling function vprofile_rootDistribution():

vprofile_rootDistribution(exampleforestMED, SpParamsMED)

Meteorological forcing

Soil water simulations require daily weather inputs. The weather variables that are required depend on the complexity of the soil water balance model we are using. In the simplest case, only mean temperature, precipitation and potential evapotranspiration is required, but the more complex simulation model also requires radiation, wind speed, min/max temparature and relative humitidy. Here we show an example of meteorological forcing data.

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

Simulation models in medfate have been designed to work along with data generated from package meteoland. The user is strongly recommended to resort to this package to obtain suitable weather input for soil water balance simulations.

Simulation control

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

control = defaultControl()
control
## $verbose
## [1] TRUE
## 
## $subdailyResults
## [1] FALSE
## 
## $defaultWindSpeed
## [1] 2.5
## 
## $soilFunctions
## [1] "SX"
## 
## $snowpack
## [1] TRUE
## 
## $rockyLayerDrainage
## [1] TRUE
## 
## $unlimitedSoilWater
## [1] FALSE
## 
## $plantWaterPools
## [1] FALSE
## 
## $leafPhenology
## [1] TRUE
## 
## $unfoldingDD
## [1] 300
## 
## $transpirationMode
## [1] "Granier"
## 
## $cavitationRefill
## [1] "total"
## 
## $refillMaximumRate
## [1] 0.05
## 
## $verticalLayerSize
## [1] 100
## 
## $gainModifier
## [1] 1
## 
## $costModifier
## [1] 1
## 
## $costWater
## [1] "dEdP"
## 
## $cochard
## [1] FALSE
## 
## $capacitance
## [1] FALSE
## 
## $klatstem
## [1] 0.01
## 
## $klatleaf
## [1] 0.01
## 
## $taper
## [1] TRUE
## 
## $numericParams
## $numericParams$maxNsteps
## [1] 400
## 
## $numericParams$ntrial
## [1] 200
## 
## $numericParams$psiTol
## [1] 1e-04
## 
## $numericParams$ETol
## [1] 1e-07
## 
## 
## $fracLeafResistance
## [1] NA
## 
## $fracRootResistance
## [1] 0.4
## 
## $averageFracRhizosphereResistance
## [1] 0.15
## 
## $Catm
## [1] 386
## 
## $ndailysteps
## [1] 24
## 
## $thermalCapacityLAI
## [1] 1e+06
## 
## $allocationStrategy
## [1] "Plant_kmax"
## 
## $nonStomatalPhotosynthesisLimitation
## [1] TRUE
## 
## $floemConductanceFactor
## [1] 0.2
## 
## $nonSugarConcentration
## [1] 0.25
## 
## $equilibriumOsmoticConcentration
##    leaf sapwood 
##     0.8     0.6 
## 
## $minimumSugarForGrowth
##     leaf  sapwood fineroot 
##     0.25     0.25     0.25 
## 
## $respirationRates
##         leaf      sapwood     fineroot 
## 2.602740e-03 6.849315e-05 2.054795e-03 
## 
## $turnoverRates
##      sapwood     fineroot 
## 0.0001261398 0.0018972310 
## 
## $constructionCosts
##     leaf  sapwood fineroot 
##     1.50     1.47     1.30 
## 
## $maximumRelativeGrowthRates
##     leaf  sapwood fineroot 
##    0.010    0.001    0.100

Most of these parameters should normally be left to their default value. However, there are three that deserve explanation here:

  1. Console output can be turned off by setting verbose = FALSE.
  2. The soil water retention curves can be switched between Saxton’s and Van Genuchten’s using parameter soilFunctions.
  3. The complexity of the soil water balance calculations will be very different if we set transpirationMode = "Sperry". Most of the other options apply in the case of advanced soil water balance only.

Water balance input object

A last step is needed before calling simulation functions. It consists in the compilation of aboveground and belowground parameters and the specification of additional parameter values for each plant cohort, such as their light extinction coefficient or their response to drought. This is done by calling function spwbInput() and taking species parameter values from species parameter data:

x = spwbInput(above, Z50, Z95, examplesoil, SpParamsMED, control)

Different parameter variables will be drawn depending on the value of transpirationMode. For the simple water balance model (transpirationMode = "Granier"), relatively few parameters are needed. All the input information for forest data and species parameter values can be inspected by accessing the different elements of this object, whose names are.

names(x)
##  [1] "control"             "canopy"              "cohorts"            
##  [4] "above"               "below"               "belowLayers"        
##  [7] "paramsPhenology"     "paramsInterception"  "paramsTranspiration"
## [10] "internalPhenology"   "internalWater"

As with the soil input object, the spwbInput object will be modified during simulations. In the case of soil water balance, these modifications are small, for example concerning LAI_expanded, but in the case of growth simulations the object is used to store the status of vegetation during and at the end of simulations. If one has a forest object, the spwbInput object can be generated in directly from it, avoiding the need to explicitly define aboveground and belowground input data:

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

Finally, note 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.

Executing the soil water balance model

Water balance for a single day

Soil water balance simulations will normally span periods of several months or years, but since the model operates at a daily temporal scale, 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 day 100 from the meteorological input data and perform soil water balance for that day only:

d = 100
sd1<-spwb_day(x, examplesoil, 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])

Function spwb_day() is most useful when working with the complex transpiration model. This is why so many meteorological variables are required. The output of spwb_day() is a list with five elements:

names(sd1)
## [1] "cohorts"      "WaterBalance" "Soil"         "Stand"        "Plants"
sd1
## $cohorts
##       SP              Name
## T1_54 54  Pinus halepensis
## T2_68 68      Quercus ilex
## S1_65 65 Quercus coccifera
## 
## $WaterBalance
##             PET            Rain            Snow         NetRain        Snowmelt 
##        5.023347        0.000000        0.000000        0.000000        0.000000 
##           Runon    Infiltration          Runoff    DeepDrainage SoilEvaporation 
##        0.000000        0.000000        0.000000        0.000000        0.500000 
## PlantExtraction   Transpiration 
##        1.045672        1.045672 
## 
## $Soil
##   SoilEvaporation PlantExtraction         psi
## 1    4.999998e-01       0.2132496 -0.03480566
## 2    1.529512e-07       0.8324222 -0.03421903
## 
## $Stand
##         LAI     LAIlive LAIexpanded     LAIdead          Cm  LgroundPAR 
##   1.6797997   1.6797997   1.6797997   0.0000000   1.2319246   0.4135178 
##  LgroundSWR 
##   0.5199022 
## 
## $Plants
##             LAI AbsorbedSWRFraction Transpiration Photosynthesis    PlantPsi
## T1_54 0.8167012          0.23326063    0.50805015      0.9499074 -0.03324984
## T2_68 0.7977952          0.23281954    0.50709058      0.9481133 -0.03324984
## S1_65 0.0653033          0.01401766    0.03053107      0.0570843 -0.03370783
##                DDS      StemPLC
## T1_54 3.184957e-06 1.876311e-07
## T2_68 9.436921e-07 7.029520e-08
## S1_65 4.147992e-07 5.311238e-08
## 
## attr(,"class")
## [1] "spwb_day" "list"

Water balance for multiple days

Most often, users will use function spwb() to run the soil water balance model. This function requires the spwbInput object, the soil object and the meteorological data frame. However, function spwb_day() modified the state variables of the input objects. In particular, the values of soil moisture are now:

examplesoil$W
## [1] 0.9895515 0.9928735

We simply reset state variables to their default values so that new simulations are not affected by the end state of the previous simulation:

resetInputs(x, examplesoil)
examplesoil$W
## [1] 1 1

Now we are ready to call function spwb():

S = spwb(x, examplesoil, examplemeteo, latitude = 41.82592, elevation = 100)
## Initial soil water content (mm): 185.069
## Initial snowpack content (mm): 0
## Performing daily simulations .....................................done.
## Final soil water content (mm): 165.599
## Final snowpack content (mm): 0
## Change in soil water content (mm): -19.47
## Soil water balance result (mm): -19.47
## Change in snowpack water content (mm): 0
## Snowpack water balance result (mm): 0
## Water balance components:
##   Precipitation (mm) 513
##   Rain (mm) 462 Snow (mm) 51
##   Interception (mm) 86 Net rainfall (mm) 376
##   Infiltration (mm) 418 Runoff (mm) 10 Deep drainage (mm) 104
##   Soil evaporation (mm) 41 Transpiration (mm) 292

Function spwb() returns an object of class with the same name, actually a list:

class(S)
## [1] "spwb" "list"

If we inspect its elements, we realize that the output is arranged differently than in spwb_day():

names(S)
## [1] "latitude"     "topography"   "spwbInput"    "soilInput"    "WaterBalance"
## [6] "Soil"         "Stand"        "Plants"       "subdaily"

In particular, element spwbInput contains a copy of the input parameters that were used to run the model:

names(S$spwbInput)
##  [1] "control"             "canopy"              "cohorts"            
##  [4] "above"               "below"               "belowLayers"        
##  [7] "paramsPhenology"     "paramsInterception"  "paramsTranspiration"
## [10] "internalPhenology"   "internalWater"

As before, WaterBalance contains water balance components, but in this case in form of a data frame with days in rows:

head(S$WaterBalance)
##                  PET Precipitation      Rain Snow    NetRain Snowmelt
## 2001-01-01 1.3212770      4.869109  4.869109    0  3.5552648        0
## 2001-01-02 2.2185985      2.498292  2.498292    0  1.2148349        0
## 2001-01-03 1.8045176      0.000000  0.000000    0  0.0000000        0
## 2001-01-04 0.9200627      5.796973  5.796973    0  4.4940363        0
## 2001-01-05 2.2914449      1.884401  1.884401    0  0.7792335        0
## 2001-01-06 1.7255058     13.359801 13.359801    0 11.7823320        0
##            Infiltration Runoff DeepDrainage Evapotranspiration Interception
## 2001-01-01    3.5552648      0    2.8201820          2.0888842     1.313844
## 2001-01-02    1.2148349      0    0.9319611          2.2452854     1.283457
## 2001-01-03    0.0000000      0    0.0000000          0.5580148     0.000000
## 2001-01-04    4.4940363      0    2.5519431          1.9728017     1.302937
## 2001-01-05    0.7792335      0    0.6143926          2.0821606     1.105168
## 2001-01-06   11.7823320      0    8.6982699          2.4366541     1.577469
##            SoilEvaporation PlantExtraction Transpiration
## 2001-01-01       0.5000000       0.2750402     0.2750402
## 2001-01-02       0.5000000       0.4618288     0.4618288
## 2001-01-03       0.1823821       0.3756327     0.3756327
## 2001-01-04       0.4783426       0.1915225     0.1915225
## 2001-01-05       0.5000000       0.4769926     0.4769926
## 2001-01-06       0.5000000       0.3591855     0.3591855

Element Plants is in turn a list with several dataframes with plant output variables, for example plant water potentials are in:

head(S$Plants$PlantPsi)
##                  T1_54       T2_68       S1_65
## 2001-01-01 -0.03243430 -0.03243430 -0.03327561
## 2001-01-02 -0.03297563 -0.03297563 -0.03356136
## 2001-01-03 -0.03354097 -0.03354097 -0.03419796
## 2001-01-04 -0.03249816 -0.03249816 -0.03328277
## 2001-01-05 -0.03306854 -0.03306854 -0.03361086
## 2001-01-06 -0.03086650 -0.03086650 -0.03247461

Inspecting model outputs

Plots

Package medfate provides a simple plot function for objects of class spwb. It can be used to show meteorological inputs, snow dynamics, and different components of the water balance:

plot(S, type = "PET_Precipitation")

plot(S, type = "Snow")

plot(S, type = "Export")

plot(S, type = "Evapotranspiration")

Function plot is also allows displaying soil moisture dynamics by layer, which can be done in four different ways (the first two only imply a change in axis units):

plot(S, type="SoilTheta")

plot(S, type="SoilRWC")

plot(S, type="SoilPsi")

plot(S, type="SoilVol")

Finally, the same function can also be used to draw the dynamics of plant variables by cohorts, such as transpiration, photosynthesis, water potential or drought stress:

plot(S, type="PlantTranspiration")

plot(S, type="PlantPhotosynthesis")

plot(S, type="PlantPsi")

plot(S, type="PlantStress")

General summaries

While the simulation model uses daily steps, users will normally be interested in outputs at larger time scales. The package provides a summary for objects of class spwb. This function can be used to summarize the model’s output at different temporal steps (i.e. weekly, annual, …). For example, to obtain the average soil moisture and water potentials by months one can use:

summary(S, freq="months",FUN=mean, output="Soil")
##                  W.1       W.2     ML.1      ML.2    MLTot       WTD        SWE
## 2001-01-01 0.9919235 0.9993122 67.71196 116.72574 184.4377  997.4094 1.64682703
## 2001-02-01 0.9708988 0.9508814 66.27675 111.06872 177.3455 1000.0000 0.27394810
## 2001-03-01 0.9744181 0.9709104 66.51699 113.40823 179.9252  998.4364 0.01762496
## 2001-04-01 0.9552660 0.8586657 65.20960 100.29737 165.5070 1000.0000 0.57778227
## 2001-05-01 0.9489861 0.8599514 64.78092 100.44755 165.2285 1000.0000 0.00000000
## 2001-06-01 0.8354257 0.6144991 57.02890  71.77723 128.8061 1000.0000 0.00000000
## 2001-07-01 0.9606114 0.6763598 65.57450  79.00293 144.5774 1000.0000 0.00000000
## 2001-08-01 0.9635080 0.7619345 65.77222  88.99858 154.7708 1000.0000 0.00000000
## 2001-09-01 0.9640173 0.8075583 65.80700  94.32772 160.1347 1000.0000 0.00000000
## 2001-10-01 0.9765642 0.8350745 66.66349  97.54177 164.2053 1000.0000 0.00000000
## 2001-11-01 0.9711153 0.9329374 66.29152 108.97275 175.2643  997.8439 2.58447175
## 2001-12-01 0.9390468 0.9260861 64.10242 108.17248 172.2749 1000.0000 0.00000000
##            PlantExt.1 PlantExt.2       psi.1       psi.2
## 2001-01-01 0.06090072  0.2377266 -0.03450821 -0.03322981
## 2001-02-01 0.12703683  0.4958886 -0.03868661 -0.04417256
## 2001-03-01 0.13200040  0.5152646 -0.03760862 -0.03874578
## 2001-04-01 0.16844714  0.6575233 -0.04188925 -0.07445462
## 2001-05-01 0.19882431  0.7760849 -0.04435245 -0.07799632
## 2001-06-01 0.28345114  1.0859160 -0.08501592 -0.54603530
## 2001-07-01 0.23925873  0.9026254 -0.04132495 -0.49000152
## 2001-08-01 0.24172374  0.9433890 -0.03996076 -0.14287855
## 2001-09-01 0.18971173  0.7405052 -0.04009309 -0.10035748
## 2001-10-01 0.13608080  0.5311632 -0.03744645 -0.09127960
## 2001-11-01 0.09829296  0.3836857 -0.03851898 -0.04907084
## 2001-12-01 0.09193304  0.3588609 -0.04596838 -0.05013918

Parameter output is used to indicate the element of the spwb object for which we desire summaries. Similarly, it is possible to calculate the average stress of plant cohorts by months:

summary(S, freq="months",FUN=mean, output="PlantStress")
##                   T1_54        T2_68        S1_65
## 2001-01-01 3.176310e-06 9.411300e-07 4.162065e-07
## 2001-02-01 7.122282e-06 2.110313e-06 7.679433e-07
## 2001-03-01 5.028273e-06 1.489862e-06 6.016818e-07
## 2001-04-01 2.996742e-05 8.879371e-06 2.416503e-06
## 2001-05-01 4.322411e-05 1.280765e-05 3.416084e-06
## 2001-06-01 2.130526e-02 6.540441e-03 1.537135e-03
## 2001-07-01 3.997226e-02 1.277906e-02 3.028427e-03
## 2001-08-01 2.535987e-04 7.515681e-05 1.777967e-05
## 2001-09-01 7.215664e-05 2.138061e-05 5.285293e-06
## 2001-10-01 7.299964e-05 2.163069e-05 5.278236e-06
## 2001-11-01 1.081028e-05 3.203065e-06 1.026886e-06
## 2001-12-01 1.051643e-05 3.115992e-06 1.195373e-06

The summary function can be also used to aggregate the output by species. In this case, the values of plant cohorts belonging to the same species will be averaged using LAI values as weights. For example, we may average the daily drought stress across cohorts of the same species (here there is only one cohort by species, so this does not modify the output):

head(summary(S, freq="day", output="PlantStress", bySpecies = TRUE))
##            Pinus halepensis Quercus coccifera Quercus ilex
## 2001-01-01     2.956298e-06      3.990466e-07 8.759411e-07
## 2001-01-02     3.106806e-06      4.094152e-07 9.205360e-07
## 2001-01-03     3.269352e-06      4.331577e-07 9.686979e-07
## 2001-01-04     2.973795e-06      3.993042e-07 8.811252e-07
## 2001-01-05     3.133141e-06      4.112295e-07 9.283391e-07
## 2001-01-06     2.547987e-06      3.709174e-07 7.549598e-07

Or we can combine the aggregation by species with a temporal aggregation (here monthly averages):

summary(S, freq="month", FUN = mean, output="PlantStress", bySpecies = TRUE)
##            Pinus halepensis Quercus coccifera Quercus ilex
## 2001-01-01     3.176310e-06      4.162065e-07 9.411300e-07
## 2001-02-01     7.122282e-06      7.679433e-07 2.110313e-06
## 2001-03-01     5.028273e-06      6.016818e-07 1.489862e-06
## 2001-04-01     2.996742e-05      2.416503e-06 8.879371e-06
## 2001-05-01     4.322411e-05      3.416084e-06 1.280765e-05
## 2001-06-01     2.130526e-02      1.537135e-03 6.540441e-03
## 2001-07-01     3.997226e-02      3.028427e-03 1.277906e-02
## 2001-08-01     2.535987e-04      1.777967e-05 7.515681e-05
## 2001-09-01     7.215664e-05      5.285293e-06 2.138061e-05
## 2001-10-01     7.299964e-05      5.278236e-06 2.163069e-05
## 2001-11-01     1.081028e-05      1.026886e-06 3.203065e-06
## 2001-12-01     1.051643e-05      1.195373e-06 3.115992e-06

Specific output functions

The package provides some functions to extract or transform specific outputs from soil plant water balance simulations. In particular, function spwb_stress() allows calculating several plant stress indices, such as the number of days with drought stress > 0.5 or the maximum drought stress:

spwb_stress(S, index = "NDD", freq = "years", draw=FALSE)
##            T1_54 T2_68 S1_65
## 2001-01-01     0     0     0
spwb_stress(S, index = "MDS", freq = "years", draw=FALSE)
##                T1_54      T2_68      S1_65
## 2001-01-01 0.2254972 0.07467767 0.01785612

As the general summary function, spwb_stress() allows calculating stress indices at several temporal scales. For example the water stress index (integral of water potential values) can be calculated and drawn for every month:

spwb_stress(S, index = "WSI", freq = "months", draw=TRUE)

Another specific summary function is spwb_waterUseEfficiency(). This is most useful with advanced water and energy balance modeling, but for simple water balance it calculates the ratio between photosynthesis and transpiration at the desired scale:

spwb_waterUseEfficiency(S, type = "Stand An/E", freq = "months", draw=FALSE)
## 2001-01-01 2001-02-01 2001-03-01 2001-04-01 2001-05-01 2001-06-01 2001-07-01 
##  0.4831519  0.9739391  1.6897987  1.4888231  3.0347171  4.3398732  4.5336767 
## 2001-08-01 2001-09-01 2001-10-01 2001-11-01 2001-12-01 
##  5.0114108  3.2328872  3.2968819  1.0235083  0.8007528