knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) is_windows <- .Platform$OS.type == "windows"
This workflow is based on using the R package kwb.hydrus1d, which is compatible with
the last official release of HYDRUS1D v4.17.0140. The
development of this new R package
was necessary, due to the fact that already
available model wrappers were developed for outdated HYDRUS1D
versions (e.g.
python package phydrus requires HYDRUS1D v4.08)
and/or are not well maintained (e.g. open issues/pull requests for R package
hydrusR).
The workflow below describes all the steps required for:
Input data preparation
Modifying (i.e. atmospheric
input data defined in file ATMOSPHERE.in
)
an already pre-prepared HYDRUS1D
model (using the HYDRUS1D GUI
)
Running it from within R
by using the command line (cmd
) and
Importing and analysing the HYDRUS1D
results within R.
# Enable this universe options(repos = c( kwbr = 'https://kwb-r.r-universe.dev', CRAN = 'https://cloud.r-project.org')) # Install R package install.packages('flextreat.hydrus1d')
paths_list <- list( extdata = system.file("extdata", package = "flextreat.hydrus1d"), exe_dir = "<extdata>/model", model_name = "test", model_dir = "<exe_dir>/<model_name>", scenario = "status-quo_no-irrigation", atmosphere = "<model_dir>/ATMOSPH.IN", a_level = "<model_dir>/A_LEVEL.out", t_level = "<model_dir>/T_LEVEL.out", runinf = "<model_dir>/Run_Inf.out", solute_id = 1, solute = "<model_dir>/solute<solute_id>.out", soil_data = "<extdata>/input-data/soil/soil_geolog.csv" ) paths <- kwb.utils::resolve(paths_list)
Soil data is derived from depth-dependent grain-size analysis of soil samples taken
in Braunschweig. The following required input parameters for the van Genuchten
model used in HYDRUS1D
were derived by using the pedotransfer function ROSETTA-API version 3, which was developed
by Zhang and Schaap, 2017.
library(flextreat.hydrus1d) soil_dat <- readr::read_csv(paths$soil_data, show_col_types = FALSE) %>% dplyr::mutate(layer_id = dplyr::if_else(.data$tiefe_cm_ende < 60, 1, 2), thickness_cm = .data$tiefe_cm_ende - .data$tiefe_cm_start, sand = .data$S_prozent + .data$G_prozent, silt = .data$U_prozent, clay = round(100 - .data$sand - .data$silt, 1)) %>% dplyr::mutate(clay = dplyr::if_else(.data$clay < 0, 0, .data$clay)) %>% soilDB:::ROSETTA(vars = c("sand", "silt", "clay"), v = "3") %>% dplyr::mutate(alpha = 10^.data$alpha, npar = 10^.data$npar, ksat = 10^.data$ksat) knitr::kable(soil_dat)
Due to similar soil characteristics, only two different layers (column layer_id
)
were defined:
- layer 1
: ranging from 0 cm
- 55 cm
depth
- layer 2
: ranging from 55 cm
- 210 cm
depth
The spatially aggregation for each layer (geometric mean) of the input parameters
for the van Genuchten
model is performed with the code below and shown in the
subsequent table.
soil_dat_per_layer <- soil_dat %>% dplyr::group_by(.data$layer_id) %>% dplyr::summarise(layer_thickness = max(.data$tiefe_cm_ende) - min(.data$tiefe_cm_start), theta_r = sum(.data$theta_r*.data$thickness_cm)/.data$layer_thickness, theta_s = sum(.data$theta_s*.data$thickness_cm)/.data$layer_thickness, alpha = sum(.data$alpha*.data$thickness_cm)/.data$layer_thickness, npar = sum(.data$npar*.data$thickness_cm)/.data$layer_thickness, ksat = sum(.data$ksat*.data$thickness_cm)/.data$layer_thickness) knitr::kable(soil_dat_per_layer)
In total three different atmospheric input boundary conditions (precipitation
, potential evaporation
and irrigation
) are used within the HYDRUS1D
model and described in the following subchapters in more detail.
Rainfall is based on historical hourly raw data downloaded from DWD open-data ftp server for station_id
= 662
(i.e. Braunschweig). Data is aggregated within R
to daily sums by using the code below:
install.packages("rdwd") library(dplyr) rdwd::updateRdwd() rdwd::findID("Braunschweig") rdwd::selectDWD(name = "Braunschweig", res = "daily") url_bs_rain <- rdwd::selectDWD(name = "Braunschweig", res = "hourly", var = "precipitation", per = "historical" ) bs_rain <- rdwd::dataDWD(url_bs_rain) precipitation_hourly <- rdwd::dataDWD(url_bs_rain) %>% dplyr::select(.data$MESS_DATUM, .data$R1) %>% dplyr::rename("datetime" = "MESS_DATUM", "precipitation_mm" = "R1") usethis::use_data(precipitation_hourly) precipitation_daily <- precipitation_hourly %>% dplyr::mutate("date" = as.Date(datetime)) %>% dplyr::group_by(date) %>% dplyr::summarise(rain_mm = sum(precipitation_mm)) usethis::use_data(precipitation_daily)
DT::datatable(head(flextreat.hydrus1d::precipitation_daily))
Is based on daily potential evaporation grids (1km x 1km) downloaded from DWD open-data server, where all grids (i.e. 46) within the Abwasserverregnungsgebiet.shp
were selected and spatially aggregated (mean
,
sd
, min
, max
) by using the code below:
shape_file <- system.file("extdata/input-data/gis/Abwasserverregnungsgebiet.shp", package = "flextreat.hydrus1d") # Only data of full months can currently be read! evapo_p <- kwb.dwd::read_daily_data_over_shape( file = shape_file, variable = "evapo_p", from = "201701", to = "202012" ) usethis::use_data(evapo_p)
DT::datatable(head(flextreat.hydrus1d::evapo_p))
Monthly irrigation volumes were provided by Abwasserverband Braunschweig for
the time period r sprintf("%s - %s", min(flextreat.hydrus1d::irrigation$date_start), max(flextreat.hydrus1d::irrigation$date_end))
as csv
file and separates between
two sources (groundwater
and clearwater
). Data preparation was carried
out with the code below:
irrigation_file <- system.file("extdata/input-data/Beregnungsmengen_AVB.csv", package = "flextreat.hydrus1d") # irrigation_area <- rgdal::readOGR(dsn = shape_file) # irrigation_area_sqm <- irrigation_area$area # 44111068m2 ## 2700ha (https://www.abwasserverband-bs.de/de/was-wir-machen/verregnung/) irrigation_area_sqm <- 27000000 irrigation <- read.csv2(irrigation_file) %>% dplyr::select(- .data$Monat) %>% dplyr::rename(irrigation_m3 = .data$Menge_m3, source = .data$Typ, month = .data$Monat_num, year = .data$Jahr) %>% dplyr::mutate(date_start = as.Date(sprintf("%d-%02d-01", .data$year, .data$month)), days_in_month = as.numeric(lubridate::days_in_month(.data$date_start)), date_end = as.Date(sprintf("%d-%02d-%02d", .data$year, .data$month, .data$days_in_month)), source = kwb.utils::multiSubstitute(.data$source, replacements = list("Grundwasser" = "groundwater.mmPerDay", "Klarwasser" = "clearwater.mmPerDay")), irrigation_cbmPerDay = .data$irrigation_m3/.data$days_in_month, irrigation_area_sqm = irrigation_area_sqm, irrigation_mmPerDay = 1000*irrigation_cbmPerDay/irrigation_area_sqm) %>% dplyr::select(.data$year, .data$month, .data$days_in_month, .data$date_start, .data$date_end, .data$source, .data$irrigation_mmPerDay, .data$irrigation_area_sqm) %>% tidyr::pivot_wider(names_from = .data$source, values_from = .data$irrigation_mmPerDay) usethis::use_data(irrigation)
DT::datatable(head(flextreat.hydrus1d::irrigation))
All model input files were initially setup by using the HYDRUS1D-GUI
but the
following below were modified manually with the
Soil input data were entered manually in SELECTOR.in
for the two layers defined
the second table in Chapter: Input Data - Soil.
Based on the atmospheric input data (see Chapter: Atmospheric Boundary Conditions) the ATMOSPHERE.in
file for HYDRUS1D
is prepared with the code below and starts with the hydrological summer half year
(assumption: soil is fully wetted at end of winter half year):
atm <- flextreat.hydrus1d::prepare_atmosphere_data() ### Set irrigation to 0 mm/day atm$groundwater.mmPerDay <- 0 atm$clearwater.mmPerDay <- 0 atm_selected <- flextreat.hydrus1d::select_hydrologic_years(atm) atm_selected_hydro_wide <- flextreat.hydrus1d::aggregate_atmosphere(atm_selected, "wide") DT::datatable(atm_selected_hydro_wide) atm_selected_hydro_long <- flextreat.hydrus1d::aggregate_atmosphere(atm_selected, "long") atm_hydro_plot <- flextreat.hydrus1d::plot_atmosphere(atm_selected_hydro_long) atm_hydro_plot kwb.utils::preparePdf(pdfFile = sprintf("boundaries-temporal_%s.pdf", paths$scenario), width.cm = 18, height.cm = 10) atm_hydro_plot dev.off() atm_prep <- flextreat.hydrus1d::prepare_atmosphere( atm_selected, conc_irrig_clearwater = 0, # use as "clearwater" tracer conc_irrig_groundwater = 0, conc_rain = 1) atmos <- kwb.hydrus1d::write_atmosphere( atm = atm_prep, MaxAL = nrow(atm_prep), round_digits = 6 # increase precision for "clearwater" tracer! ) writeLines(atmos, paths$atmosphere)
The selected time period covers r as.integer(diff(range(atm_selected$date)))
days
(r sprintf("%s - %s", min(atm_selected$date), max(atm_selected$date))
), i.e. it covers r round(as.integer(diff(range(atm_selected$date)))/(365/2), 0)
hydrological half years.
DT::datatable(atm_selected)
These time-series were converted with the function flextreat.hydrus1d::prepare_atmosphere()
into the data format required
by HYDRUS1D
. Due to the fact, that irrigation rates
(i.e. sum of clearwater.mmPerDay
and groundwater.mmPerDay
) cannot be entered separately as input column within HYDRUS1D
, these were simply added to th prec
(i.e. precipitation
) column. The whole time series defined in ATMOSPHERE.in
is shown below:
DT::datatable(atm_prep)
Finally the model is run automatically by using the following code:
exe_path <- kwb.hydrus1d::check_hydrus_exe(dir = paths$exe_dir, skip_preinstalled = TRUE) kwb.hydrus1d:::run_model(exe_path = exe_path, model_path = paths$model_dir)
runinf <- kwb.hydrus1d::read_runinf(paths$runinf) summary(runinf)
t_level <- kwb.hydrus1d::read_tlevel(paths$t_level) t_level ## t_level aggregate tlevel_aggr_date <- flextreat.hydrus1d::aggregate_tlevel(t_level) tlevel_aggr_yearmonth <- flextreat.hydrus1d::aggregate_tlevel(t_level, col_aggr = "yearmonth") tlevel_aggr_year_hydrologic <- flextreat.hydrus1d::aggregate_tlevel(t_level, col_aggr = "year_hydrologic") %>% dplyr::filter(.data$diff_time >= 364) ### filter out as only may-october DT::datatable(tlevel_aggr_year_hydrologic)
wb_date_plot <- flextreat.hydrus1d::plot_waterbalance(tlevel_aggr_date) wb_yearmonth_plot <- flextreat.hydrus1d::plot_waterbalance(tlevel_aggr_yearmonth) wb_yearhydrologic_plot <- flextreat.hydrus1d::plot_waterbalance(tlevel_aggr_year_hydrologic) wb_date_plot wb_yearmonth_plot wb_yearhydrologic_plot kwb.utils::preparePdf(pdfFile = sprintf("water-balance_yearmonth_%s.pdf", paths$scenario), width.cm = 19, height.cm = 10) wb_yearmonth_plot dev.off() saveRDS(wb_yearmonth_plot, file = sprintf("wb_yearmonth_%s.Rds", paths$scenario))
plotly::ggplotly(wb_date_plot) plotly::ggplotly(wb_yearmonth_plot)
a_level <- kwb.hydrus1d::read_alevel(paths$a_level) a_level
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.