--- title: "Model control options and parameters" output: rmarkdown::html_vignette: toc: true bibliography: refs.bib vignette: > %\VignetteIndexEntry{Model control options and parameters} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ## Introduction This vignette demonstrates the use of different model control option setting (`options_90`-argument of `runLWFB90()`) and parameters (`param_b90`-argument). Aside from the basic technical information of the simulation (`startdate`, `enddate`, `fornetrad`, `prec_interval` and `correct_prec`), the model control options control the annual course of leaf area index (`lai_method`), the phenology models to use (`budburst_method`, `leaffall_method`), the input and interpolations of annual stand properties (`standprop_input`, `standprop_interp`, `use_growthperiod`) and which root density depth distribution function to use (`root_method`). The interplay of options and parameters is shown briefly in the following paragraphs, by describing how options and parameters are passed from the `options_b90` and `param_b90` arguments to the individual functions that are called from within `run_LWFB90()`. For this purpose, we create basic lists of model control options and parameters, as well as soil and climate objects. ```{r, warning = FALSE, message = FALSE} library(LWFBrook90R) library(data.table) options_b90 <- set_optionsLWFB90() param_b90 <- set_paramLWFB90() ``` ## Aboveground vegetation characteristics ### Intra-annual variation of leaf area index The default parameter set (created with `set_paramLWFB90()`) represents a deciduous forest stand, without leafs in winter and maximum leaf area index in summer. The maximum leaf area index is defined by the parameter `param_b90$maxlai`, the minimum value in winter is internally calculated as a fraction (`param_b90$winlaifrac`) of `param_b90$maxlai`. The basic shape of the intra-annual leaf area index dynamics can be selected by the option `options_b90$lai_method`. The default setting `'b90'` makes use of the parameters `budburstdoy`, `leaffalldoy`, `emergedur` and `leaffalldur`, that define the dates of budburst and leaffall, and the durations of leaf unfolding and leaf shedding until `maxlai`, and respectively `winlaifrac` in winter are reached. Within `run_LWFB90()`, the parameters are passed to `make_seasLAI()` that constructs the daily timeseries of leaf area index development for a single year: ```{r} LAI_b90 <- make_seasLAI(method = options_b90$lai_method, year = 2003, maxlai = param_b90$maxlai, winlaifrac = param_b90$winlaifrac, budburst_doy = param_b90$budburstdoy, leaffall_doy = param_b90$leaffalldoy, emerge_dur = param_b90$emergedur, leaffall_dur = param_b90$leaffalldur) ``` `make_seasLAI()` also provides other shape functions that require additional parameters. For example, the option `lai_method = 'linear'` uses value pairs of day-of-year and leaf area index as fraction of `maxlai` passed from parameters `param_b90$lai_doy` and `param_b90$lai_frac`. The doy/value-pairs are then used to interpolate the intra-annual course of leaf area index to a daily time series. ```{r} options_b90$lai_method <- "linear" param_b90$lai_doy <- c(1,110,117,135,175,220,250,290,365) param_b90$lai_frac <- c(0.1,0.1,0.5,0.7,1.2,1.2,1.0,0.1,0.1) LAI_linear <- make_seasLAI(method = options_b90$lai_method, year = 2003, maxlai = param_b90$maxlai, lai_doy = param_b90$lai_doy , lai_frac = param_b90$lai_frac) ``` A third shape-option for the intra-annual variation of leaf area index is called 'Coupmodel' and uses the interpolation method as implemented in the 'Coupmodel' [@jansson_coupled_2004]. With `lai_method ='Coupmodel`, form parameters for leaf unfolding and leaf fall (`shp_budburst`, `shp_leaffall`), and the date when leaf area is at its maximum (`shp_optdoy`) are required, in addition to the parameters required `by lai_method = 'b90'`. ```{r} options_b90$lai_method <- "Coupmodel" param_b90$shp_budburst <- 0.5 param_b90$shp_leaffall <- 5 param_b90$shp_optdoy <- 180 LAI_coupmodel <- make_seasLAI(method = options_b90$lai_method, year = 2003, maxlai = param_b90$maxlai, budburst_doy = param_b90$budburstdoy, leaffall_doy = param_b90$leaffalldoy, shp_budburst = param_b90$shp_budburst, shp_leaffall = param_b90$shp_leaffall, shp_optdoy = param_b90$shp_optdoy) ``` A plot of all three methods shows the roles of the different parameters: ```{r, echo = FALSE, fig.height=5, fig.width=7, fig.cap = "Methods featured by make_seasLAI()" } oldpar <- par(no.readonly = T) par(xpd = TRUE, mar = c(5.1,4.1,2.1,2.1), oma = c(1,1,1,1)) plot(LAI_b90, type = "n", xlab = "doy", ylab = "lai [m²/m²]", ylim = c(0,6)) with(param_b90, abline(v = c(budburstdoy,budburstdoy+emergedur, leaffalldoy, leaffalldoy+leaffalldur), lty = 2, xpd = FALSE)) lines(LAI_b90, col ="green", lwd = 2) lines(LAI_linear, col ="red", lwd = 2) lines(LAI_coupmodel, col ="blue", lwd = 2) with(param_b90, arrows(x0 = c(budburstdoy,leaffalldoy,shp_optdoy,budburstdoy,leaffalldoy), x1 = c(budburstdoy,leaffalldoy, shp_optdoy, budburstdoy+emergedur, leaffalldoy + leaffalldur), y0 = c(-1.6,-1.6,3.5,2.5,2.5),y1 = c(-0.3,-0.3,4.8,2.5,2.5), length = 0.15, code = 2)) with(param_b90, arrows(x0 = c(budburstdoy,leaffalldoy), x1 = c(budburstdoy+emergedur, leaffalldoy + leaffalldur), y0 = c(2.5,2.5),y1 = c(2.5,2.5), length = 0.15, code = 3)) with(param_b90, text(x = c(budburstdoy, leaffalldoy, shp_optdoy), y = c(-1.9,-1.9,3.2), c("budburstdoy", "leaffalldoy", "shp_optdoy"))) with(param_b90, text(x = c(budburstdoy+0.5*emergedur,leaffalldoy+0.5*leaffalldur), y = 2, c("emergedur","leaffalldur"))) with(param_b90, text( x = c(budburstdoy/2, (leaffalldoy - budburstdoy - emergedur)/2 + budburstdoy + emergedur), y = c(winlaifrac*maxlai+0.4,maxlai+0.4), c("winlaifrac * maxlai", "maxlai"))) legend("topleft",c("'b90'","'linear'", "'Coupmodel'"), pch = NULL, lwd =2, col = c("green", "red", "blue"), bty = "n") par(oldpar) ``` ### Inter-annual variation of leaf area index By passing a single value via `param_b90$maxlai` we used the same maximum leaf area index for each year of the simulation period. In order to incorporate between-year variation of the leaf area index, we can simply assign vectors of values for each year of the simulation period to any of the parameters used by function `make_seasLAI()`. In the following example, we pass three values for `maxlai` and `shp_optdoy`, to get different seasonal courses of leaf area index for the three years of the simulation period. Additionally, we add variation to the dates of budburst, by assigning a vector of values to the parameter `budburstdoy`. ```{r,echo = 1:6, results = 'hide',fig.height=5, fig.width=7, fig.cap = "Options and parameters affecting interannual variation of leaf area index."} years <- 2001:2003 param_b90$maxlai <- c(4,6,5) param_b90$shp_optdoy <- c(210,180,240) param_b90$shp_budburst <- c(3,1,0.3) param_b90$budburstdoy <- c(100,135,121) lai_variation <- make_seasLAI(method = options_b90$lai_method, year = years, maxlai = param_b90$maxlai, budburst_doy = param_b90$budburstdoy, leaffall_doy = param_b90$leaffalldoy, shp_budburst = param_b90$shp_budburst, shp_leaffall = param_b90$shp_leaffall, shp_optdoy = param_b90$shp_optdoy) par(mar=c(4.1,4.1,1.1,4.1), oma = c(1,1,1,1)) plot(seq.Date(as.Date("2001-01-01"), as.Date("2003-12-31"), by = "day"), lai_variation, col = "green", ylab = "lai [m²/m²]", type ="l", xlab = "", lwd = 2) arrows( x0 = as.Date(paste(param_b90$budburstdoy, years),format = "%j %Y"), y0 = -1.0, y1 = -0.3, length = 0.15, code = 2, xpd =TRUE) text(x = as.Date(paste(param_b90$budburstdoy, years),format = "%j %Y"), y = -1.3, paste("doy", param_b90$budburstdoy), xpd = TRUE) par(oldpar) ``` Beside the obvious between-year variation of maximum leaf area index, we can also see the effect of the shape parameter for the leaf unfolding phase `shp_budburst`. Values greater 1 result in concave, values below 1 in convex functions, while values of 1 give linear progressions. The budburst day-of-year is varying as specified in the parameters, but can also be estimated using temperature based phenology models. By selecting other settings than the default `options_b90$budburst_method = 'fixed'` and `options_b90$leaffall_method = 'fixed'`, the `vegperiod()` function of the 'vegperiod'-Package is called from within `run_LWFB90`. `budburstdoy` and/or `leaffalldoy` are then calculated for each year from the climate data using the desired methods. See `vegperiod` for a list of available models. The estimated values for `budburstdoy` and/or `leaffalldoy` can be found in the `param_b90` list element of the results object after the simulation. ### Other plant properties (`height`, `sai`, `densef`) Like the leaf area index parameters and budburst/leaffall-dates, it is also possible to provide vectors of values for stand height (`height`), stem area index (`sai`), and stand density (`densef`) to generate between-year variation of stand characteristics. From the yearly values, daily values are interpolated using the function `approx_standprop()`. The `approx.method`- argument of the function defines how to interpolate the yearly values passed by `y`. Within `run_LWFB90()`, the option `options_b90$standprop_interp` is passed to the `approx.method`- argument of `approx_standprop`. The default interpolation method `standprop_interp = 'constant'` results in a yearly changing step function, while `standprop_interp = 'linear'` interpolates the values: ```{r} # constant 'interpolation' options_b90$standprop_interp <- 'constant' param_b90$height <- c(20.2,20.8,21.3) simyears <- 2002:2003 height_c <- approx_standprop(x_yrs=years, y = param_b90$height, approx.method = options_b90$standprop_interp) # linear interpolation options_b90$standprop_interp <- 'linear' param_b90$height_ini <- 19.1 height_l <- approx_standprop(x_yrs=years, y = param_b90$height, y_ini = param_b90$height_ini, approx.method = options_b90$standprop_interp) ``` For linear interpolation, additional parameters `height_ini`, `sai_ini`, `densef_ini` have to be provided to `run_LWFB90()` via the `param_b90`-argument. These parameters define the values at the beginning of the simulation, to which the value of the first year is interpolated to. By default, the yearly values are interpreted to be valid at December 31st of the respective years, so that the interpolated timeseries are linearly increasing or decreasing during the whole year. In order to constrain the interpolation to the growth period only, the option `options_b90$use_growthperiod` was introduced, which requires the arguments `startdoy` and `enddoy`, when set to `TRUE`. Then, values decrease or increase between budburst and leaffall only, and remain constant during winter. ```{r} options_b90$use_growthperiod <- TRUE height_l_gp <- approx_standprop(x_yrs = years, y = param_b90$height, y_ini = param_b90$height_ini, use_growthperiod = options_b90$use_growthperiod, startdoy = param_b90$budburstdoy, enddoy = param_b90$leaffalldoy, approx.method = options_b90$standprop_interp) ``` The following plot explains the differences between the interpolated timeseries of stand height using the different options and parameters: ```{r,echo = FALSE, fig.height=5, fig.width=7, fig.cap="Interpolated stand height derived from parameters using approx_standprop()"} dates <- seq.Date(from = as.Date(paste0(min(years),"-01-01")), to = as.Date(paste0(max(years),"-12-31")), by = "day") par(mar=c(4.1,4.1,1.1,4.1)) plot(dates, height_c, type = "l", lwd = 2, col = "black", ylim = c(19,22), ylab = "height [m]", xlab = "", xpd = TRUE) lines(dates, height_l, col = "blue", lwd = 2) lines(dates, height_l_gp, col = "green", lwd = 2) legend("topleft", legend = c("approx.method = 'constant'", "approx.method = 'linear'", "approx.method = 'linear', use_growthperiod = TRUE"), col = c("black", "blue", "green"), lwd = 2, pch = NULL, bty = "n") arrows( x0 = as.Date(paste(param_b90$budburstdoy, years),format = "%j %Y"), y0 = c(param_b90$height_ini,param_b90$height[1:2])-0.5, y1 = c(param_b90$height_ini,param_b90$height[1:2])-0.1, length = 0.15, code = 2, xpd =TRUE) text(x = as.Date(paste(param_b90$budburstdoy, years), format = "%j %Y"), y = c(param_b90$height_ini,param_b90$height[1:2])-0.7, paste("doy", param_b90$budburstdoy), xpd = TRUE) arrows( x0 = as.Date(paste(param_b90$leaffalldoy, years),format = "%j %Y"), y0 = param_b90$height-0.5, y1 = param_b90$height-0.1, length = 0.15, code = 2, xpd =TRUE) text(x = as.Date(paste(param_b90$leaffalldoy, years), format = "%j %Y"), y = param_b90$height-0.7, paste("doy", param_b90$leaffalldoy), xpd = TRUE) par(oldpar) ``` Another option for incorporating between-year variation of plant properties is to provide a data.frame with yearly values of `height`, `maxlai`, `sai`, `densef` and `age` as list item `standprop_table` in `param_b90`. To take effect, the option `options_b90$standprop_input` has to be set to `'table'`. In this case, the values passed via parameters `height`, `sai`, `densef` and `age_ini` are ignored. As `maxlai` is also provided via the table, the `maxlai` value from parameters is ignored as well, while the other parameters that affect intra-annual leaf area development (e.g., `shp_budburst`) are still active. For demonstration purposes we use the table `slb1_standprop`, that contains observed stand data of the Solling Beech Experimental site from 1966 to 2014, along with estimated leaf and stem area index derived using allometric functions. For creating the daily timeseries of stand properties, we use `run_LWFB90()`, and make use of the option to not run the model (`run = FALSE`), but only return the model input. ```{r, messages = FALSE} #Extend simulation period options_b90$startdate <- as.Date("1980-01-01") options_b90$enddate <- as.Date("1999-12-31") soil <- cbind(slb1_soil, hydpar_wessolek_tab(texture = slb1_soil$texture)) #set up options for table input options_b90$standprop_input <- 'table' param_b90$standprop_table <- slb1_standprop # Set up dynamic budburst and leaf fall options_b90$budburst_method <- "Menzel" options_b90$leaffall_method <- "vonWilpert" param_b90$budburst_species <- "Fagus sylvatica" #run LWF-Brook90 without simulation standprop_daily <- run_LWFB90(options_b90 = options_b90, param_b90 = param_b90, climate = slb1_meteo, soil = soil, output = output, run = FALSE)$standprop_daily ``` ```{r, echo = FALSE, fig.height=5, fig.width=7, fig.cap="Stand properties generated using table input of annual stand characteristics"} par(mar=c(4.1,4.1,1.1,4.1), oma = c(1,1,1,1)) with(standprop_daily,{ plot(dates, lai, type = "l", lwd = 2, col = "green", ylab = "lai, sai m²/m²", xlab = "") lines(dates, sai, lwd = 2, col = "brown") par(new = TRUE) plot(dates, height, type = "l", lwd = 2, col = "blue", ylim = c(27, 32), ylab = "", xlab = "", xaxt = "n", yaxt = "n") }) axis(4, at = seq(27,32, by = 1)) mtext("height [m]", side = 4, line = 3) legend("bottom",horiz = TRUE, bty = "n",inset = -0.25, xpd =TRUE, legend = c("lai", "sai", "height"), col = c("green","brown", "blue"), lwd = 2, pch = NULL) par(oldpar) ``` ## Root density depth distribution The root depth density depth distribution can either be provided in the column `rootden` of the `soil`- argument of `run_LWFB90()`, or can be derived from parameters using the function `make_rootden()`. In order to use root density as specified in the soil data, the `root_method` element of the `options_b90`-list has to be set to `'soilvar'`. Other method names are passed to `make_rootden()`. Currently, the function provides four methods to assign values of relative root density to a vector of soil depths. The default method `'betamodel'` uses the model of Gale & Grigal (-@gale_vertical_1987), which is of the form $y = 1- \beta^d$, where $y$ is the cumulative root fraction down to soil depth $d$ and $\beta$ is the depth coefficient. Larger values of $\beta$ correspond to a greater proportion of roots in deeper soil layers: ```{r, echo = FALSE, fig.height=4, fig.width=4} plot(make_rootden(soilnodes = seq(0,-2, by = -0.01), method = "betamodel", beta = 0.93), seq(-0.01,-2, by = -0.01), ylim = c(-2,0), type="l", lwd = 1.5, xlab = "relative root density", ylab = "soil depth [m]") lines(make_rootden(soilnodes = seq(0,-2, by = -0.01), method = "betamodel", beta = 0.96), seq(-0.01,-2, by = -0.01), lty = 2, lwd = 1.5) lines(make_rootden(soilnodes = seq(0,-2, by = -0.01), method = "betamodel", beta = 0.98), seq(-0.01,-2, by = -0.01), lty = 3, lwd = 1.5) lines(make_rootden(soilnodes = seq(0,-2, by = -0.01), method = "betamodel", beta = 0.99), seq(-0.01,-2, by = -0.01), lty = 4, lwd = 1.5) legend("bottomright",legend = c(expression(paste(beta, "= 0.93")),expression(paste(beta, "= 0.96")), expression(paste(beta, "= 0.98")),expression(paste(beta, "= 0.99"))), lwd = 1.5, pch = NULL, lty = 1:4, bty = "n") ``` For larger values of $\beta$, the root density will reach zero only in very deep soil layers. In order to set the root density to zero at any desired soil depth, the parameter `maxrootdepth` was defined. With this parameter, the root density is set to zero in all soil layers that lie deeper than `maxrootdepth`. Within `run_LWFB90()`, the function is called in the following way: ```{r} param_b90$maxrootdepth <- -1.4 options_b90$root_method <- "betamodel" roots_beta <- make_rootden(soilnodes = c(max(slb1_soil$upper), slb1_soil$lower), maxrootdepth = param_b90$maxrootdepth, beta = param_b90$betaroot, method = options_b90$root_method) ``` A second option to define the root distribution for the soil layers is to provide value pairs of soil depth and root density in a data.frame and assign it to the `rootden_table`-entry in `param_b90`. As an example, we set up a hypothetical root density depth distribution: ```{r} options_b90$root_method <- 'table' param_b90$rootden_table <- data.frame( upper = c(0.03,0,-0.02, -0.15, -0.35, -0.5, -0.65,-0.9,-1.1,-1.3), lower = c(0,-0.02, -0.15, -0.35, -0.5, -0.65,-0.9,-1.1,-1.3,-1.6), rootden = c(10,15, 35, 15, 7.5, 4, 12, 2, 2, 0)) roots_table <- make_rootden(soilnodes = c(max(slb1_soil$upper), slb1_soil$lower), method = options_b90$root_method, rootdat = param_b90$rootden_table) ``` A third option generates a linear root density depth distribution, with the maximum at the uppermost soil layer and a root density of 0 at `maxrootdepth`. If the parameter `relrootden` is provided, the first element of the vector is used as the maximum, otherwise the interpolation is made between 0 and 1. The last option returns a uniform root distribution, with the first vector-element of `relrootden` (if provided) as value for all layers down to `maxrootdepth`. ```{r, echo = 1:4, fig.height=4, fig.width=4} options_b90$root_method <- 'linear' roots_linear <- make_rootden(soilnodes = c(max(slb1_soil$upper), slb1_soil$lower), maxrootdepth = param_b90$maxrootdepth, method = options_b90$root_method) options_b90$root_method <- 'const' roots_constant <- make_rootden(soilnodes = c(max(slb1_soil$upper), slb1_soil$lower), maxrootdepth = param_b90$maxrootdepth, method = options_b90$root_method) plot(roots_constant, slb1_soil$lower, type = 's', lwd = 1.5,ylab = "soil depth [m]",xlab = "relative root density", col = "red") lines(roots_linear, slb1_soil$lower, type = 's', col = "blue", lwd = 1.5) lines(roots_table/100, slb1_soil$lower, type = 's', col = "green", lwd = 1.5) lines(roots_beta*10, slb1_soil$lower, type = 's', col = "brown", lwd = 1.5) legend("bottomright", c("'betamodel'","'table'","'linear'", "'constant'"),seg.len = 1.5, pch = NULL, lwd =1.5, col = c("brown", "green", "blue", "red"), bty = "n") ``` # References