Landscape water balance
Miquel De Caceres
2023-03-20
LandscapeWaterBalance.RmdDisplaying properties of spatial objects
We begin by loading an example dataset of 100 forest stands distributed on points in the landscape:
data("examplepointslandscape")We transform the data set to class sf because sp objects are deprecated from ver. 1.0.0 on.
pts_sf <- sp_to_sf(examplepointslandscape)
pts_sf## Simple feature collection with 100 features and 8 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 402000 ymin: 4640000 xmax: 429000 ymax: 4646000
## Projected CRS: ETRS89 / UTM zone 31N
## # A tibble: 100 × 9
## geometry id elevation slope aspect land_cover_type forest soil
## <POINT [m]> <chr> <dbl> <dbl> <dbl> <chr> <list> <list>
## 1 (411000 4646000) 81065 516 9.94 92.7 wildland <forest> <soil>
## 2 (412000 4646000) 81066 537 13.7 115. wildland <forest> <soil>
## 3 (413000 4646000) 81067 578 23.3 152. wildland <forest> <soil>
## 4 (415000 4646000) 81068 615 12.7 360. wildland <forest> <soil>
## 5 (416000 4646000) 81069 644 22.5 186. wildland <forest> <soil>
## 6 (423000 4646000) 81070 520 4.04 225. wildland <forest> <soil>
## 7 (426000 4646000) 81073 608 12.2 337. wildland <forest> <soil>
## 8 (428000 4646000) 81075 745 8.06 332. wildland <forest> <soil>
## 9 (402000 4645000) 81076 527 14.4 245. wildland <forest> <soil>
## 10 (404000 4645000) 81077 483 26.0 36.7 wildland <forest> <soil>
## # … with 90 more rows, and 1 more variable: state <list>Using plot_variable() functions for spatial landscape
objects, we can draw maps of some variables using:
plot_variable(pts_sf, "basal_area")
The set of maps available can be known by inspecting the help of
function extract_variables(). Alternatively, the package
provides function shinyplot_land() to display maps
interactively.
Simulations on a set of forest stands
Package medfate includes functions spwb(),
growth() and fordyn() to simulate soil water
balance, carbon balance and forest dynamics on a single forest stand,
respectively. This section describe how to run simulations on a set of
forest stands in one call. This is done using functions
spwb_spatial(), growth_spatial() and
fordyn_spatial(), respectively.
As an example, we will use function spwb_spatial(),
which simulates soil plant water balance on forests distributed in
particular locations. The function takes an object of class
sf as input.
Weather inputs and local control parameters
As before, we need meteorological data, species parameters and control parameters for local simulations (which will apply to all forest stands):
Calling the simulation function
The call to spwb_spatial() can be done as follows (here
we restrict the dates for simplicity):
dates <- seq(as.Date("2001-01-01"), as.Date("2001-02-28"), by="day")
res <- spwb_spatial(pts_sf, SpParamsMED, examplemeteo,
dates = dates, local_control = local_control, progress = FALSE)Function spwb_spatial() first initializes model inputs
by calling forest2spwbInput() for each forest stand
described in the sf landscape object. Then it calls
function spwb() for each forest stand and stores the
result. The fact that we used examplemeteo as
meteorological input involves that the same weather was applied to all
forest stands, but different weather could have been specified for each
one (see documentation of function spwb_spatial()).
The result of calling spwb_spatial() is an object of
class sf with the following columns:
names(res)## [1] "geometry" "id" "state" "result"Column geometry contains the geometry given as input to
simulations, column id contains the identification label of
each stand, column state contains the
spwbInput corresponding to each forest stand (which can be
used in subsequent simulations) and column result contains
the output of spwb() function for each forest stand
(i.e. its elements are objects of the S3 class spwb).
Temporal summaries, plots and maps
The structure of the output of spwb_spatial() allows
querying information for the simulation of any particular forest stand.
For example, we may use function plot.spwb(), from package
medfate, to display the simulation results on a
particular plot:
plot(res$result[[1]], "Evapotranspiration")
Similarly, if we want a monthly summary of water balance for the
first stand, we can use function summary.spwb() from
package medfate:
summary(res$result[[1]], freq="months",FUN=sum, output="WaterBalance")## PET Precipitation Rain Snow NetRain Snowmelt
## 2001-01-01 45.07945 74.74949 58.098839 16.650650 39.6370782 13.093006
## 2001-02-01 84.71521 4.99943 2.457859 2.541571 0.9699644 5.552842
## Infiltration Runoff DeepDrainage Evapotranspiration Interception
## 2001-01-01 52.730084 0 40.6162567 34.03849 18.461760
## 2001-02-01 6.522806 0 0.1569381 21.18700 1.487895
## SoilEvaporation PlantExtraction Transpiration
## 2001-01-01 5.5785170 9.996655 9.998214
## 2001-02-01 0.9530808 18.740241 18.746027
## HydraulicRedistribution
## 2001-01-01 0.03169547
## 2001-02-01 0.00000000However, a more convenient way of generating summaries is by
calculating them on all forest stands in one step, using function
simulation_summary() on objects issued from
simulations:
res_sum <- simulation_summary(res, summary_function = summary.spwb,
freq="months", output="WaterBalance")where the argument summary_function points to the
function to be used to generate local summaries and the remaining
arguments are those of the local summary function. The result of using
simulation_summary() is again an object of class
sf that contains the spatial geometry and the list of
summaries for all stands:
names(res_sum)## [1] "geometry" "id" "summary"The summary for the first stand can now be accessed through the first
element of column summary:
res_sum$summary[[1]]## PET Precipitation Rain Snow NetRain Snowmelt
## 2001-01-01 45.07945 74.74949 58.098839 16.650650 39.6370782 13.093006
## 2001-02-01 84.71521 4.99943 2.457859 2.541571 0.9699644 5.552842
## Infiltration Runoff DeepDrainage Evapotranspiration Interception
## 2001-01-01 52.730084 0 40.6162567 34.03849 18.461760
## 2001-02-01 6.522806 0 0.1569381 21.18700 1.487895
## SoilEvaporation PlantExtraction Transpiration
## 2001-01-01 5.5785170 9.996655 9.998214
## 2001-02-01 0.9530808 18.740241 18.746027
## HydraulicRedistribution
## 2001-01-01 0.03169547
## 2001-02-01 0.00000000Summary objects are handy because their plot_summary()
function allows us to display maps of summaries for specific dates:
plot_summary(res_sum, "Transpiration", "2001-01-01", limits=c(0,45))
plot_summary(res_sum, "Transpiration", "2001-02-01", limits=c(0,45))
To avoid displaying maps one by one, the package includes function
shinyplot_land() that allows displaying maps of temporal
summaries interactively.
Advanced simulation features
Simulation with integrated temporal summaries
If one needs to save memory, it is possible with
spwb_spatial() to generate temporal summaries automatically
after the simulation of soil water balance of each stand, and storing
those summaries instead of all the output of function
spwb().
For example the following code will keep temporal summaries of water balance components instead of simulation results:
res_2 <- spwb_spatial(pts_sf, SpParamsMED, examplemeteo,
dates = dates, local_control = local_control,
keep_results = FALSE, progress = FALSE,
summary_function = summary.spwb, summary_arguments = list(freq="months"))Parameter keep_results = FALSE tells
spwb_spatial() not to keep the simulation results of forest
stands, whereas summary_function = summary.spwb tells
spwb_spatial() to perform and store summaries before
discarding the results of any stand. The output has slightly different
column names:
names(res_2)## [1] "geometry" "id" "state" "result" "summary"In particular, result is not included. Now the temporal
summaries can be directly accessed through the column
summary:
res_2$summary[[1]]## PET Precipitation Rain Snow NetRain Snowmelt
## 2001-01-01 45.07945 74.74949 58.098839 16.650650 39.6370782 13.093006
## 2001-02-01 84.71521 4.99943 2.457859 2.541571 0.9699644 5.552842
## Infiltration Runoff DeepDrainage Evapotranspiration Interception
## 2001-01-01 52.730084 0 40.6162567 34.03849 18.461760
## 2001-02-01 6.522806 0 0.1569381 21.18700 1.487895
## SoilEvaporation PlantExtraction Transpiration
## 2001-01-01 5.5785170 9.996655 9.998214
## 2001-02-01 0.9530808 18.740241 18.746027
## HydraulicRedistribution
## 2001-01-01 0.03169547
## 2001-02-01 0.00000000And one can produce maps with summary results directly from the output of the simulation function:
plot_summary(res_2, "Transpiration", "2001-02-01", limits=c(0,45))
The possibility of running a summary function after the simulation of
each stand is not limited to summary.spwb(). Users can
define their own summary functions, provided the first argument is
object, which will contain the result of the simulation
(i.e., the result of calling spwb(), growth()
or fordyn()). For example, the following function calls
droughtStress():
f_stress <- function(object, ...) {
return(droughtStress(object, ..., draw = FALSE))
}Now we can call again spwb_spatial:
res_3 <- spwb_spatial(pts_sf, SpParamsMED, examplemeteo,
dates = dates, local_control = local_control,
keep_results = FALSE, progress = FALSE,
summary_function = f_stress,
summary_arguments = list(index = "ADS", freq = "months", bySpecies=TRUE))The drought stress summary of stand #1 is:
res_3$summary[[1]]## Dorycnium spp. Genista spp. Hedera helix Juniperus communis
## 2001-01-01 NaN NaN 0.006165402 0.003861844
## 2001-02-01 NaN NaN 0.007618106 0.004739411
## Lavandula spp. Phillyrea angustifolia Pinus nigra Pistacia lentiscus
## 2001-01-01 0.01720418 0.003735971 0.009425589 0.003528457
## 2001-02-01 0.02141699 0.004658094 0.011569840 0.004327903
## Quercus coccifera Quercus faginea Quercus ilex Rubus spp.
## 2001-01-01 0.003229448 NaN 0.01092188 NaN
## 2001-02-01 0.004017091 NaN 0.01272611 NaN
## Salvia rosmarinus Viburnum spp.
## 2001-01-01 NaN NaN
## 2001-02-01 NaN NaNAnd we can draw a map with the monthly drought stress of Quercus ilex across space using:
plot_summary(res_3, "Quercus ilex", "2001-02-01")
Simulation in several steps
The result of a simulation includes an element state,
which stores the state of soil and stand variables at the end of the
simulation. This information can be used to perform a new simulation
from the point where the first one ended. In order to do so, we need to
update the state variables in spatial object with their values at the
end of the simulation, using function
update_landscape():
pts_sf_mod <- update_landscape(pts_sf, res)The resulting object w is the same as y
except for the state variables. For example we can compare the water
potential in the first soil layer:
plot_variable(pts_sf, "psi1")
plot_variable(pts_sf_mod, "psi1")
By using w as input we can now simulate water balance in
the set of stands for an extra month:
dates <- seq(as.Date("2001-03-01"), as.Date("2001-03-31"), by="day")
res_3 <- spwb_spatial(pts_sf_mod, SpParamsMED, examplemeteo,
dates = dates, local_control = local_control,
summary_function = summary.spwb, summary_arguments = list(freq = "months"),
progress = FALSE)And display a map with the resulting month transpiration:
plot_summary(res_3, "Transpiration", "2001-03-01", limits=c(0,45))
Simulation of watershed hydrology
Simulation of watershed hydrology involves describing a watershed in a gridded mode and simulating soil water balance water while including additional compartments and hydrological processes, such as routing runoff from one cell to the other.
Example data set
To illustrate this kind of simulation, which is carried out using
function spwb_land(), we first load a small example
watershed included with the package, which we transform to
sf:
## Simple feature collection with 66 features and 18 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 401380 ymin: 4671820 xmax: 402880 ymax: 4672620
## Projected CRS: WGS 84 / UTM zone 31N
## # A tibble: 66 × 19
## geometry id elevation slope aspect land_cover_type forest soil
## <POINT [m]> <int> <dbl> <dbl> <dbl> <chr> <list> <list>
## 1 (402630 4672570) 1 1162 11.3 79.2 wildland <forest> <soil>
## 2 (402330 4672470) 2 1214 12.4 98.7 agriculture <forest> <soil>
## 3 (402430 4672470) 3 1197 10.4 102. wildland <forest> <soil>
## 4 (402530 4672470) 4 1180 8.12 83.3 wildland <forest> <soil>
## 5 (402630 4672470) 5 1164 13.9 96.8 wildland <forest> <soil>
## 6 (402730 4672470) 6 1146 11.2 8.47 agriculture <forest> <soil>
## 7 (402830 4672470) 7 1153 9.26 356. agriculture <forest> <soil>
## 8 (402230 4672370) 8 1237 14.5 75.1 wildland <forest> <soil>
## 9 (402330 4672370) 9 1213 13.2 78.7 wildland <forest> <soil>
## 10 (402430 4672370) 10 1198 8.56 75.6 agriculture <forest> <soil>
## # … with 56 more rows, and 11 more variables: state <list>, waterOrder <int>,
## # waterQ <list>, queenNeigh <list>, channel <lgl>, depth_to_bedrock <dbl>,
## # bedrock_conductivity <dbl>, bedrock_porosity <dbl>, snowpack <dbl>,
## # aquifer <dbl>, represented_area <dbl>Here is a map of elevation values:
plot_variable(ws_sf, "elevation")
Since the landscape contains agricultural lands, we need to define crop factors, which will determine transpiration flow as a proportion of potential evapotranspiration:
ws_sf$crop_factor = NA
ws_sf$crop_factor[ws_sf$land_cover_type=="agriculture"] = 0.75Simulation function
As for the call to function spwb_spatial() we will use
the same weather (i.e. examplemeteo) across the watershed.
To speed up calculations we call function spwb_land() for a
single month.
dates <- seq(as.Date("2001-01-01"), as.Date("2001-01-31"), by="day")
res_ws <- spwb_land(ws_sf, SpParamsMED, examplemeteo, dates = dates, summary_frequency = "month")##
## ------------ spwb_land ------------
## Grid cells: 66, patchsize: 100 m2, area: 0.66 ha
## Cell land use wildland: 48 agriculture: 17 artificial: 0 rock: 1 water: 0
## Cells with soil: 65
## Meteorological input class: data.frame
## Number of days to simulate: 31
## Number of summaries: 1
## Number of outlet cells: 1
##
## Preparing spwb input...
## done.
##
## Initial average soil water content (mm): 530.8
## Initial average snowpack water content (mm): 0
## Initial average aquifer water content (mm): 0
## Initial watershed water content (mm): 522.76
##
## Performing daily simulations:
## .+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++.+++++done
##
## Final average soil water content (mm): 530.8
## Final average snowpack water content (mm): 0.05
## Final average aquifer water content (mm): 53.17
## Final watershed water content (mm): 522.76
##
## Change in snowpack water content (mm): 0.05
## Snowpack water balance result (mm): 3.19
## Snowpack water balance components:
## Snow fall (mm) 16.65 Snow melt (mm) 13.46
##
## Change in soil water content (mm): 0
## Soil water balance result (mm): -10.58
## Soil water balance components:
## Net rainfall (mm) 37.92 Snow melt (mm) 13.47
## Runon (mm) 1.1 Runoff (mm) 0
## Subsurface input (mm) 0 Subsurface output (mm) 0
## Deep drainage (mm) 53.99 Aquifer discharge (mm) 0
## Soil evaporation (mm) 5.79 Plant transpiration (mm) 3.29
##
## Change in aquifer water content (mm): 53.17
## Aquifer water balance result (mm): 53.17
## Aquifer water balance components:
## Deep drainage (mm) 53.17 Aquifer discharge (mm) 0
##
## Change in watershed water content (mm): 0
## Watershed water balance result (mm): 45.94
## Watershed water balance components:
## Precipitation (mm) 74.75
## Interception (mm) 19.87 Soil evaporation (mm) 5.7 Plant Transpiration (mm) 3.24
## Export (mm) 0
## Watershed lateral flows:
## Subsurface flow (mm) 0
## Groundwater flow (mm) 7.77
##
## ------------ spwb_land ------------Although simulations are performed using daily temporal steps,
parameter summary_frequency allows storing results at
coarser temporal scales, to reduce the amount of memory in spatial
results.
Summaries and plots
Unlike spwb_spatial() where summaries could be
arbitrarily generated a posteriori from simulation results,
with spwb_land() the summaries are always fixed and
embedded with the simulation result. For example, we can inspect the
summaries for a given landscape cell using:
res_ws$sf$summary[[1]]## Runon Runoff Infiltration Rain NetRain Snow Snowmelt
## 2001-01-01 0 0 66.83459 58.09884 53.1979 16.65065 13.6367
## Interception DeepDrainage AquiferDischarge SaturationExcess
## 2001-01-01 4.900943 57.47486 0 0
## SoilEvaporation Transpiration SWE SoilVol Psi1 WTD DTA
## 2001-01-01 6.995272 3.28633 0 532.3416 -0.03300001 4000 14.78647
## InterflowInput InterflowOutput BaseflowInput BaseflowOutput
## 2001-01-01 0 0 5.935362 4.507225Several plots can be drawn from the result of function
spwb_land() in a similar way as done for
spwb_spatial(). As an example we display a map of soil
water volume for the simulated month:
plot_summary(res_ws$sf, "SoilVol", "2001-01-01")