Basic water balance

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"

In the basic water balance, only mean temperature, precipitation and potential evapotranspiration is required, but radiation may also be necessary to simulate snow melt.

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 can be controlled using a set of global parameters. The default parameterization is obtained using function defaultControl():

control = defaultControl("Granier")

Some parameters 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 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)

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"             "soil"                "canopy"             
##  [4] "cohorts"             "above"               "below"              
##  [7] "belowLayers"         "paramsPhenology"     "paramsInterception" 
## [10] "paramsTranspiration" "internalPhenology"   "internalWater"

Note that (since ver 2.0) the water balance input object contains an element soil with the soil input.

It is important to remember that, unlike normal objects in R, any water balance simulation will modify the spwbInput object. For example, 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. 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.

Finally, note that users can set cohort-specific parameters for soil water balance (instead of using species-level values) by modifying manually the parameter values in this object. Since some parameters may be coordinated by design, however, it is better to use specific package functions for this purpose.

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, 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"

• cohorts: Table with the species code and species name of each cohort.
• WaterBalance: Contains the soil water balance flows (precipitation, infiltration, transpiration, …)
• Soil: Contains output values by soil layer (i.e. water evaporated from each soil layer, water transpired from each soil layer and the final soil water potential).
• Stand: A list with stand LAI (expanded and dead leaves), canopy water retention capacity and the proportion of light (SWR or PAR) reaching the ground.
• Plants: Contains output values by plant cohort (i.e. LAI values, transpiration, water potential, drought stress index, …).

sd1

## $cohorts
##         SP              Name
## T1_148 148  Pinus halepensis
## T2_168 168      Quercus ilex
## S1_165 165 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.210594        1.210594 
## 
## $Soil
##   SoilEvaporation PlantExtraction         psi
## 1    4.999998e-01       0.2402826 -0.03487639
## 2    1.529512e-07       0.9703118 -0.03442613
## 
## $Stand
##         LAI     LAIlive LAIexpanded     LAIdead          Cm  LgroundPAR 
##   1.9727123   1.9727123   1.9727123   0.0000000   1.4895748   0.3553437 
##  LgroundSWR 
##   0.4646726 
## 
## $Plants
##               LAI AbsorbedSWRFraction Transpiration GrossPhotosynthesis
## T1_148 1.00643723          0.27864285    0.63012543          2.27683706
## T2_168 0.92661573          0.24911558    0.56335249          0.83134791
## S1_165 0.03965932          0.00756893    0.01711646          0.04164607
##           PlantPsi          DDS      StemPLC
## T1_148 -0.03324984 1.449686e-06 4.582245e-07
## T2_168 -0.03324984 9.421210e-07 9.421210e-07
## S1_165 -0.03370783 6.436854e-07 7.904802e-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 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:

x$soil$W

## [1] 0.9891555 0.9916930

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)
x$soil$W

## [1] 1 1

Now we are ready to call function spwb():

S = spwb(x, examplemeteo, latitude = 41.82592, elevation = 100)

## Initial soil water content (mm): 185.069
## Initial snowpack content (mm): 0
## Performing daily simulations 
##  Year 2001:....................................done.
## Final soil water content (mm): 162.249
## Final snowpack content (mm): 0
## Change in soil water content (mm): -22.82
## Soil water balance result (mm): -22.82
## 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) 100 Net rainfall (mm) 362
##   Infiltration (mm) 404 Runoff (mm) 9 Deep drainage (mm) 63
##   Soil evaporation (mm) 38 Transpiration (mm) 326

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"    "WaterBalance" "Soil"        
## [6] "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"             "soil"                "canopy"             
##  [4] "cohorts"             "above"               "below"              
##  [7] "belowLayers"         "paramsPhenology"     "paramsInterception" 
## [10] "paramsTranspiration" "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.2920057        0
## 2001-01-02 2.2185985      2.498292  2.498292    0  0.9563356        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.2310882        0
## 2001-01-05 2.2914449      1.884401  1.884401    0  0.6696101        0
## 2001-01-06 1.7255058     13.359801 13.359801    0 11.4937061        0
##            Infiltration Runoff DeepDrainage Evapotranspiration Interception
## 2001-01-01    3.2920057      0    2.6113540          2.3955225     1.577103
## 2001-01-02    0.9563356      0    0.6493216          2.5766241     1.541956
## 2001-01-03    0.0000000      0    0.0000000          0.6144928     0.000000
## 2001-01-04    4.2310882      0    2.1823438          2.2151420     1.565885
## 2001-01-05    0.6696101      0    0.4673639          2.2670149     1.214791
## 2001-01-06   11.4937061      0    8.3792443          2.7819304     1.866094
##            SoilEvaporation PlantExtraction Transpiration
## 2001-01-01       0.5000000       0.3184193     0.3184193
## 2001-01-02       0.5000000       0.5346680     0.5346680
## 2001-01-03       0.1796157       0.4348772     0.4348772
## 2001-01-04       0.4275280       0.2217292     0.2217292
## 2001-01-05       0.5000000       0.5522235     0.5522235
## 2001-01-06       0.5000000       0.4158359     0.4158359

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_148      T2_168      S1_165
## 2001-01-01 -0.03249324 -0.03249324 -0.03330650
## 2001-01-02 -0.03305830 -0.03305830 -0.03360539
## 2001-01-03 -0.03369993 -0.03369993 -0.03429369
## 2001-01-04 -0.03257571 -0.03257571 -0.03326236
## 2001-01-05 -0.03311174 -0.03311174 -0.03363392
## 2001-01-06 -0.03094718 -0.03094718 -0.03251484

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, gross photosynthesis, water potential or drought stress:

plot(S, type="Transpiration")

plot(S, type="GrossPhotosynthesis")

plot(S, type="PlantPsi")

plot(S, type="PlantStress")

Finally, one can interactively create plots using function shinyplot, e.g.:

shinyplot(S)

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.9913015 0.9980352 67.66950 116.57658 184.2461  997.7262 1.68796802
## 2001-02-01 0.9682675 0.9364257 66.09712 109.38021 175.4773 1000.0000 0.28335308
## 2001-03-01 0.9719232 0.9633419 66.34667 112.52419 178.8709  998.7720 0.01762496
## 2001-04-01 0.9454933 0.8250240 64.54248  96.36782 160.9103 1000.0000 0.59071605
## 2001-05-01 0.9447760 0.7819989 64.49352  91.34223 155.8357 1000.0000 0.00000000
## 2001-06-01 0.8176867 0.5048067 55.81799  58.96449 114.7825 1000.0000 0.00000000
## 2001-07-01 0.9513086 0.5878203 64.93946  68.66098 133.6004 1000.0000 0.00000000
## 2001-08-01 0.9595725 0.6134178 65.50357  71.65092 137.1545 1000.0000 0.00000000
## 2001-09-01 0.9610666 0.6048004 65.60557  70.64436 136.2499 1000.0000 0.00000000
## 2001-10-01 0.9737840 0.5913550 66.47370  69.07386 135.5476 1000.0000 0.00000000
## 2001-11-01 0.9693913 0.7604586 66.17384  88.82619 155.0000 1000.0000 2.68795419
## 2001-12-01 0.9354059 0.9093367 63.85388 106.21605 170.0699 1000.0000 0.00000000
##            PlantExt.1 PlantExt.2       psi.1       psi.2
## 2001-01-01 0.06862091  0.2771057 -0.03462781 -0.03347430
## 2001-02-01 0.14314093  0.5780320 -0.03927984 -0.04836369
## 2001-03-01 0.14873370  0.6006176 -0.03815044 -0.04053473
## 2001-04-01 0.18980069  0.7664371 -0.04422593 -0.09353383
## 2001-05-01 0.22402884  0.9045384 -0.04560445 -0.13481545
## 2001-06-01 0.31938419  1.0115551 -0.09582387 -1.60769546
## 2001-07-01 0.26958872  0.9251424 -0.04413189 -1.03089299
## 2001-08-01 0.27236638  1.0965811 -0.04087960 -0.43312505
## 2001-09-01 0.21376098  0.8604834 -0.04078323 -0.46204112
## 2001-10-01 0.15333136  0.6135933 -0.03801668 -0.56367102
## 2001-11-01 0.11075329  0.4467669 -0.03890918 -0.23349606
## 2001-12-01 0.10358719  0.4183059 -0.04696503 -0.05545887

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_148       T2_168       S1_165
## 2001-01-01 1.468709e-06 9.544836e-07 6.539347e-07
## 2001-02-01 4.246538e-06 2.759740e-06 1.447573e-06
## 2001-03-01 2.653684e-06 1.724576e-06 1.036851e-06
## 2001-04-01 2.762424e-05 1.795256e-05 7.136898e-06
## 2001-05-01 1.171952e-04 7.616851e-05 2.827651e-05
## 2001-06-01 1.774825e-01 1.334371e-01 5.338952e-02
## 2001-07-01 1.297993e-01 9.975341e-02 4.059602e-02
## 2001-08-01 3.031035e-03 1.972700e-03 7.141245e-04
## 2001-09-01 3.486202e-03 2.268459e-03 8.210120e-04
## 2001-10-01 8.690910e-03 5.671831e-03 2.055770e-03
## 2001-11-01 1.003865e-03 6.526789e-04 2.364025e-04
## 2001-12-01 6.345690e-06 4.123938e-06 2.270343e-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     1.352957e-06      6.209666e-07 8.792590e-07
## 2001-01-02     1.424776e-06      6.378347e-07 9.259324e-07
## 2001-01-03     1.509357e-06      6.778349e-07 9.809000e-07
## 2001-01-04     1.363286e-06      6.185011e-07 8.859712e-07
## 2001-01-05     1.431697e-06      6.394604e-07 9.304306e-07
## 2001-01-06     1.168876e-06      5.777313e-07 7.596287e-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     1.468709e-06      6.539347e-07 9.544836e-07
## 2001-02-01     4.246538e-06      1.447573e-06 2.759740e-06
## 2001-03-01     2.653684e-06      1.036851e-06 1.724576e-06
## 2001-04-01     2.762424e-05      7.136898e-06 1.795256e-05
## 2001-05-01     1.171952e-04      2.827651e-05 7.616851e-05
## 2001-06-01     1.774825e-01      5.338952e-02 1.334371e-01
## 2001-07-01     1.297993e-01      4.059602e-02 9.975341e-02
## 2001-08-01     3.031035e-03      7.141245e-04 1.972700e-03
## 2001-09-01     3.486202e-03      8.210120e-04 2.268459e-03
## 2001-10-01     8.690910e-03      2.055770e-03 5.671831e-03
## 2001-11-01     1.003865e-03      2.364025e-04 6.526789e-04
## 2001-12-01     6.345690e-06      2.270343e-06 4.123938e-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_148 T2_168 S1_165
## 2001-01-01      8      3      0

spwb_stress(S, index = "MDS", freq = "years", draw=FALSE)

##               T1_148    T2_168    S1_165
## 2001-01-01 0.7018867 0.5925645 0.2583035

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. In this case it is equal to the value of the input species parameter WUE:

spwb_waterUseEfficiency(S, type = "Stand Ag/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 
##   2.601888   2.601888   2.601888   2.601883   2.601860   2.570343   2.583841 
## 2001-08-01 2001-09-01 2001-10-01 2001-11-01 2001-12-01 
##   2.601350   2.601320   2.600268   2.601696   2.601887