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.
Model inputs are explained in greater detail in vignette ‘Simulation inputs’. Here we only review the different steps required to run function spwb()
.
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()
:
= defaultSoilParams(2) spar
The soil input for water balance simulation is actually a list of class soil
that is created using a function with the same name:
= soil(spar) examplesoil
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")
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()
:
= defaultControl("Sperry") control
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.
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()
:
= forest2spwbInput(exampleforestMED, examplesoil, SpParamsMED, control) x
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):
$cohorts x
## 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:
$canopy x
## 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:
$above x
## 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
:
$below x
## 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
:
$belowLayers x
## $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 spwbInput
object also includes cohort parameter values for several kinds of traits. For example, plant anatomy parameters are described in paramsAnatomy
:
$paramsAnatomy x
## 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
:
$paramsTranspiration x
## 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
:
$paramsWaterStorage x
## 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.
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.
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
.
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.
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:
= 100
d 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")
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")
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:
<-spwb_day(x, rownames(examplemeteo)[d],
sd1$MinTemperature[d], examplemeteo$MaxTemperature[d],
examplemeteo$MinRelativeHumidity[d], examplemeteo$MaxRelativeHumidity[d],
examplemeteo$Radiation[d], examplemeteo$WindSpeed[d],
examplemeteolatitude = 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"
Element WaterBalance
contains the soil water balance flows of the day (precipitation, infiltration, transpiration, …)
$WaterBalance sd1
## 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:
$Soil sd1
## 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
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")
Element Plants
contains output values by plant cohort. Several output variables can be inspected in this element.
$Plants sd1
## 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)
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:
$SunlitLeaves sd1
## 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
$ShadeLeaves sd1
## 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)
Or variations in vapour pressure deficit:
plot(sd1, type = "LeafVPD", bySpecies=TRUE)
Or variations in leaf water potential:
plot(sd1, type = "LeafPsi", bySpecies=TRUE)
plot(sd1, type = "LeafCi", bySpecies=TRUE)
plot(sd1, type = "LeafIntrinsicWUE", bySpecies=TRUE)
Users will often use function spwb()
to run the soil water balance model for several days. This function requires the spwbInput
object, the soil
object and the meteorological data frame. However, running spwb_day()
modified the input objects. In particular, the soil moisture at the end of the simulation was:
$soil$W x
## [1] 0.9893779 0.9965024
And the temperature of soil layers:
$soil$Temp x
## [1] 9.148856 7.392471
We can also see the current state of canopy variables:
$canopy x
## zlow zmid zup Tair Cair VPair
## 1 0 50 100 6.573796 386 0.5170718
## 2 100 150 200 6.573796 386 0.5170718
## 3 200 250 300 6.573796 386 0.5170718
## 4 300 350 400 6.573796 386 0.5170718
## 5 400 450 500 6.573796 386 0.5170718
## 6 500 550 600 6.573796 386 0.5170718
## 7 600 650 700 6.573796 386 0.5170718
## 8 700 750 800 6.573796 386 0.5170718
## 9 800 850 900 6.573796 386 0.5170718
## 10 900 950 1000 6.573796 386 0.5170718
## 11 1000 1050 1100 6.573796 386 0.5170718
## 12 1100 1150 1200 6.573796 386 0.5170718
## 13 1200 1250 1300 6.573796 386 0.5170718
## 14 1300 1350 1400 6.573796 386 0.5170718
## 15 1400 1450 1500 6.573796 386 0.5170718
## 16 1500 1550 1600 6.573796 386 0.5170718
## 17 1600 1650 1700 6.573796 386 0.5170718
## 18 1700 1750 1800 6.573796 386 0.5170718
## 19 1800 1850 1900 6.573796 386 0.5170718
## 20 1900 1950 2000 6.573796 386 0.5170718
## 21 2000 2050 2100 6.573796 386 0.5170718
## 22 2100 2150 2200 6.573796 386 0.5170718
## 23 2200 2250 2300 6.573796 386 0.5170718
## 24 2300 2350 2400 6.573796 386 0.5170718
## 25 2400 2450 2500 6.573796 386 0.5170718
## 26 2500 2550 2600 6.573796 386 0.5170718
## 27 2600 2650 2700 6.573796 386 0.5170718
## 28 2700 2750 2800 6.573796 386 0.5170718
We simply use function resetInputs()
to reset state variables to their default values, so that the new simulation is not affected by the end state of the previous simulation:
resetInputs(x)
$soil$W x
## [1] 1 1
$soil$Temp x
## [1] NA NA
$canopy x
## zlow zmid zup Tair Cair VPair
## 1 0 50 100 NA 386 NA
## 2 100 150 200 NA 386 NA
## 3 200 250 300 NA 386 NA
## 4 300 350 400 NA 386 NA
## 5 400 450 500 NA 386 NA
## 6 500 550 600 NA 386 NA
## 7 600 650 700 NA 386 NA
## 8 700 750 800 NA 386 NA
## 9 800 850 900 NA 386 NA
## 10 900 950 1000 NA 386 NA
## 11 1000 1050 1100 NA 386 NA
## 12 1100 1150 1200 NA 386 NA
## 13 1200 1250 1300 NA 386 NA
## 14 1300 1350 1400 NA 386 NA
## 15 1400 1450 1500 NA 386 NA
## 16 1500 1550 1600 NA 386 NA
## 17 1600 1650 1700 NA 386 NA
## 18 1700 1750 1800 NA 386 NA
## 19 1800 1850 1900 NA 386 NA
## 20 1900 1950 2000 NA 386 NA
## 21 2000 2050 2100 NA 386 NA
## 22 2100 2150 2200 NA 386 NA
## 23 2200 2250 2300 NA 386 NA
## 24 2300 2350 2400 NA 386 NA
## 25 2400 2450 2500 NA 386 NA
## 26 2500 2550 2600 NA 386 NA
## 27 2600 2650 2700 NA 386 NA
## 28 2700 2750 2800 NA 386 NA
Now we are ready to call function spwb()
. In this example, we only simulate 61 days to save computational time:
= spwb(x, examplemeteo[110:170,], latitude = 41.82592, elevation = 100) S
## Initial soil water content (mm): 180.814
## Initial snowpack content (mm): 0
## Performing daily simulations
## Year 2001:......done.
## Final soil water content (mm): 138.812
## Final snowpack content (mm): 0
## Change in soil water content (mm): -42.0022
## Soil water balance result (mm): -42.0022
## Change in snowpack water content (mm): 0
## Snowpack water balance result (mm): 0
## Water balance components:
## Precipitation (mm) 36
## Rain (mm) 26 Snow (mm) 10
## Interception (mm) 7 Net rainfall (mm) 19
## Infiltration (mm) 29 Runoff (mm) 0 Deep drainage (mm) 20
## Soil evaporation (mm) 4 Transpiration (mm) 47
## Plant extraction from soil (mm) 47 Plant water balance (mm) 0 Hydraulic redistribution (mm) 0
Function spwb()
returns an object of class spwb. If we inspect its elements, we realize that the output is arranged differently than in spwb_day()
:
names(S)
## [1] "latitude" "topography" "spwbInput"
## [4] "WaterBalance" "EnergyBalance" "Temperature"
## [7] "TemperatureLayers" "Soil" "Stand"
## [10] "Plants" "SunlitLeaves" "ShadeLeaves"
## [13] "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" "paramsAnatomy"
## [10] "paramsInterception" "paramsTranspiration" "paramsWaterStorage"
## [13] "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-04-20 2.641801 2.1625135 0.0000000 2.162513 0.00000000 0.0000000
## 2001-04-21 1.875251 3.7992356 0.0000000 3.799236 0.00000000 0.0000000
## 2001-04-22 2.903129 4.1962782 0.0000000 4.196278 0.00000000 0.5608083
## 2001-04-23 3.633982 1.4698434 1.4698434 0.000000 0.52229955 9.5972190
## 2001-04-24 3.891957 0.1538991 0.1538991 0.000000 0.05468707 0.0000000
## 2001-04-25 4.171116 0.1317847 0.1317847 0.000000 0.04682885 0.0000000
## Infiltration Runoff DeepDrainage Evapotranspiration Interception
## 2001-04-20 0.00000000 0 0.00000000 0.2810470 0.00000000
## 2001-04-21 0.00000000 0 0.00000000 0.2641887 0.00000000
## 2001-04-22 0.56080827 0 0.01235285 0.3361506 0.00000000
## 2001-04-23 10.11951859 0 7.76312312 1.8292369 0.94754388
## 2001-04-24 0.05468707 0 1.40460498 0.8771818 0.09921203
## 2001-04-25 0.04682885 0 0.00000000 0.9526656 0.08495583
## SoilEvaporation PlantExtraction Transpiration
## 2001-04-20 0.0000000 0.2810470 0.2810470
## 2001-04-21 0.0000000 0.2641887 0.2641887
## 2001-04-22 0.0000000 0.3361506 0.3361506
## 2001-04-23 0.5000000 0.3816930 0.3816930
## 2001-04-24 0.1871698 0.5907999 0.5907999
## 2001-04-25 0.1286249 0.7390849 0.7390849
## HydraulicRedistribution
## 2001-04-20 0.0000000000
## 2001-04-21 0.0000000000
## 2001-04-22 0.0000000000
## 2001-04-23 0.0008133194
## 2001-04-24 0.0002850573
## 2001-04-25 0.0000000000
Elements Plants
is itself a list with several elements that contain daily output results by plant cohorts, for example leaf minimum (midday) water potentials are:
head(S$Plants$LeafPsiMin)
## T1_148 T2_168 S1_165
## 2001-04-20 -0.6064335 -0.2505370 -0.4252339
## 2001-04-21 -0.5855836 -0.2093489 -0.4364805
## 2001-04-22 -0.7981248 -0.2982040 -0.5219960
## 2001-04-23 -0.9299837 -0.3493486 -0.5466284
## 2001-04-24 -1.1504263 -0.4794885 -0.6609777
## 2001-04-25 -1.3720054 -0.5586618 -0.7941075
Package medfate
also provides a plot
function for objects of class spwb
. It can be used to show the meteorological input. Additionally, it can also be used to draw soil and plant variables. In the code below we draw water fluxes, soil water potentials, plant transpiration and plant (mid-day) water potential:
plot(S, type="Evapotranspiration")
plot(S, type="SoilPsi", bySpecies = TRUE)
plot(S, type="PlantTranspiration", bySpecies = TRUE)
plot(S, type="LeafPsiMin", bySpecies = TRUE)
Alternatively, one can interactively create plots using function shinyplot
, e.g.:
shinyplot(S)
While the simulation model uses daily steps, users may 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-04-01 0.9844236 0.9951190 63.16140 116.08359 179.2450 997.5961 1.611044
## 2001-05-01 0.9417348 0.9523744 60.42245 111.09730 171.5198 998.9143 0.000000
## 2001-06-01 0.8246673 0.8165651 52.91131 95.25474 148.1661 1000.0000 0.000000
## PlantExt.1 PlantExt.2 HydraulicInput.1 HydraulicInput.2 psi.1
## 2001-04-01 0.1521013 0.3457103 1.003289e-04 0 -0.03638390
## 2001-05-01 0.2275819 0.5178126 1.568071e-05 0 -0.04591984
## 2001-06-01 0.2907275 0.6556277 0.000000e+00 0 -0.08657902
## psi.2
## 2001-04-01 -0.03412939
## 2001-05-01 -0.04297083
## 2001-06-01 -0.08457711
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-04-01 0.02901917 0.02751211 0.01001912
## 2001-05-01 0.04277728 0.04040150 0.01515990
## 2001-06-01 0.05870397 0.05643002 0.02227212
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-04-20 0.01531745 0.006168423 0.01578876
## 2001-04-21 0.01596393 0.007008712 0.01549581
## 2001-04-22 0.01947257 0.007226023 0.01958509
## 2001-04-23 0.01982170 0.006062780 0.01894226
## 2001-04-24 0.03111000 0.009776324 0.02844583
## 2001-04-25 0.04569488 0.015521868 0.04124936
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-04-01 0.02901917 0.01001912 0.02751211
## 2001-05-01 0.04277728 0.01515990 0.04040150
## 2001-06-01 0.05870397 0.02227212 0.05643002