Section 10 Model Parameterization

These sections are not detailed enough to warrant their own dedicated chapter.

10.1 Parameter Changes

These sections mostly involve changing a few parameter values and are therefore in small sub sections without R-code. Some of this might change later.

10.1.1 Channel parameters revised

This an ongoing unresolved issue #49.

So far nothing has been changed.

10.1.2 Crop parameters verified

The crop database has been updated and stored in the following file:

"model_data/cs10_setup/swat_input_mod/plants.plt"

Discussions to this topic can be found in issue #27.

From the whole database, we only use crops defined in our management files:

Crop management: oats, barl, wwht, swht, pota, fesc

Generic: frst, fesc and others which are minor

The parameters for these crops have been updated to reflect the colder growing conditions.

10.1.3 Soil physical parameters in final form

This has been mentioned in issue #28, and has been closed without discussion. Seems like the parameters are in final form, but no further documentation exists at this point.

10.1.4 Soil chemical parameters in final form

This step will be verified in issue #46. It has been discussed in #19.

The following changes have been implemented, but are subject to change:

Changes made to soilnut1 of nutrients.sol

Parameter	Default	CS10
lab_p	5	20
nitrate	7	8.5
inorgp	3.5	35.1
watersol_p	0.15	0.4

nutrients_sol <- readLines("model_data/cs10_setup/swat_input/nutrients.sol")

header <- nutrients_sol[1]

nutrients_sol_df <- read.table("model_data/cs10_setup/swat_input/nutrients.sol", header = T, skip = 1, sep = "", fill = T, as.is = T, colClasses = "character")

nutrients_sol_df$lab_p <- 20
nutrients_sol_df$nitrate <- 8.5
nutrients_sol_df$inorgp <- 35.1
nutrients_sol_df$watersol_p <- 0.4

write.table(
  nutrients_sol_df,
  "model_data/cs10_setup/swat_input_mod/nutrients.sol",
  quote = F,
  sep = "   ",
  row.names = F, 
)

header_new <- paste("modified by the CS10 workflow in section 9.1")

nutrients_sol <- readLines("model_data/cs10_setup/swat_input_mod/nutrients.sol")

writeLines(c(header_new, nutrients_sol), "model_data/cs10_setup/swat_input_mod/nutrients.sol")

Testing if the setup still works:

simout <- test_swat_mod()
cat(tail(simout), sep = "\n")

10.1.5 Impoundment parameters defined

Has been postponed by #20, will be dealt with in #54

10.1.6 Water diversions defined

Deemed as not relevant for this catchment

10.1.7 Point sources parameters added

Is an ongoing issue in #21

10.1.8 Tile drainage parameters defined

This has been discussed in Issue #7 and is an ongoing issue in #53.

old_path <- "model_data/cs10_setup/swat_input/tiledrain.str"
new_path <-  "model_data/cs10_setup/swat_input_mod/tiledrain.str"

tiledrain_str <- readLines(old_path)

header <- paste("modified by CS10 workflow in section 9.1")

tiledrain_df <-
  read.table(
    old_path,
    sep = "",
    header = T,
    skip = 1,
    fill = T,
    colClasses = "character"
  )

tiledrain_df$dp <- 800 # from 1000
tiledrain_df$t_fc <- 12 # from 24
tiledrain_df$lag <- 30 # from 96
tiledrain_df$rad <- 200 # from 100
tiledrain_df$dist <- 8000 # from 30
tiledrain_df$drain <- 40 # from 10
tiledrain_df$pump <- 0 # from 1
tiledrain_df$lat_ksat <- 2 # from 2 (no change)

write.table(header, file = new_path, quote = F, col.names = F, row.names = F)
write.table(
  x = tiledrain_df,
  file = new_path,
  sep = "   ",
  quote = F,
  append = T,
  row.names = F
)

## Warning in write.table(x = tiledrain_df, file = new_path, sep = " ", quote = F,
## : appending column names to file

Testing if the setup still works:

simout <- test_swat_mod()
cat(tail(simout), sep = "\n")

10.1.9 Atmospheric deposition defined

Has been done in Section 3.4

10.1.10 Additional settings verified

Refers to chapter 3.11 in The Protocol and files parameters.bsn codes.bsn Has been discussed in #22

The following changes have been made:

Changes made to codes.bsn

Parameter	Default	CS10
pet	1	2
rte_cha	0	1
cn	0	1
tiledrain	0	1
soil_p	0	1
atmo_dep	a	m

old_path <- "model_data/cs10_setup/swat_input/codes.bsn"
new_path <- "model_data/cs10_setup/swat_input_mod/codes.bsn"

codes_bsn <- readLines(old_path)

header <- codes_bsn[1]

header <- paste("modified by CS10 workflow in section 9.1")

codes_df <- read.table(old_path, skip = 1, header = T, sep = "", 
                       colClasses = "character")

PET method (pet) The OPTAIN project recommends the PET calculation of Pennman-Monteith, which is code 2.

codes_df$pet <- 2 

Channel routing network (rte_cha) Recommended is to start with Muskingum (code 1) and apply variable storage method if it causes problems.

codes_df$rte_cha <- 1

Stream water quality (wq_cha) Recommended to test both for OPTAIN (1/0)

Daily curve number calculation (cn) Code 1 is recommended. (Plant ET)

codes_df$cn <- 1

OPTAIN recommends new version

codes_df$soil_p <- 1

Lapse rate is being tested in 12.1.1

Plant growth stress is being testing in 12.3

Tile drainage

This should be set to 1, but that causes a crash. We will fix this bug in a dedicated section

codes_df$tiledrain <- 0

Atmospheric Deposition

atmodep was done in Section 3.4 and is annual

codes_df$atmo_dep <- "y"

Writing the changes

write.table(header, file = new_path, col.names = F, row.names = F, quote = F) 


write.table(
  codes_df,
  file = new_path,
  col.names = T,
  append = T,
  quote = F,
  sep = "   ",
  row.names = F
)

## Warning in write.table(codes_df, file = new_path, col.names = T, append = T, :
## appending column names to file

Testing if the setup still works:

simout <- test_swat_mod()
cat(tail(simout), sep = "\n")

10.2 Fixing tiledrain

run_this_chapter = FALSE

Problem: tiledrain=1 in codes.bsn causes crash.

Reason: Some HRUs have drains, but are on a shallow soil type. This causes the drains to be below the soil, leading to a (nondescript) error.

Likely source of error: Misalignment of the soil and land use map.

Discussion of this issue can be found in #53

require(mapview)
require(sf)
require(stringr)
require(dplyr)
require(purrr)
require(DT)

10.2.1 Finding the problem HRUs

Our drains are set to 80cm depth. Which of our soils are less than that?

usersoil <- readr::read_csv("model_data/input/soil/usersoil_26_v2.csv", show_col_types = F)

shallow_soils <- usersoil$SNAM[which(usersoil$SOL_ZMX < 800)]

print(shallow_soils)

## [1] "ANT" "END" "HUM" "Mon" "RG"  "SDs"

Which of our HRUs use these soils, and are drained?

hru_data <- read.table( "model_data/cs10_setup/swat_input/hru-data.hru",
                        skip = 1, header = T, sep = "")

shallow_hrus <- which(hru_data$soil %in% shallow_soils)

drained_hrus <- grepl(x = hru_data$lu_mgt, pattern = "_drn") %>% which()

# Intersection of HRUs which are both shallow and drained
problem_hrus <- base::intersect(shallow_hrus, drained_hrus)

length(problem_hrus)

## [1] 6

Just 6 HRUs have an issue. Lets get some info on these 6:

problem_hru_data <- hru_data[problem_hrus,]

datatable(problem_hru_data)

We need to link the soil type to the hru vector file

hru_shp <- read_sf("model_data/cs10_setup/optain-cs10/data/vector/hru.shp")
hru_soil <- hru_data %>% dplyr::select(soil, name)
hru_shp <- left_join(hru_shp, hru_soil, by = "name")

And separate out our problematic HRUs and affected LUs

problem_shps <- hru_shp %>% filter(name %in% problem_hru_data$name)
problem_hru_ids <- problem_hru_data$lu_mgt %>% str_remove("_drn_lum")
problem_lus <- hru_shp %>% filter(type %in% problem_hru_ids)

Next we remove the problem HRUs and affect LUs from our hru_shp (this is for the upcoming map, to remove overlaps).

hru_shp <- hru_shp %>% filter((hru_shp$name %in% problem_lus$name) == FALSE)
hru_shp <- hru_shp %>% filter((hru_shp$name %in% problem_shps$name) == FALSE)

We also remove the problem HRUs from the affected LUs

problem_lus <- problem_lus %>% filter((problem_lus$name %in% problem_shps$name) == FALSE)

Data processing done, lets have a look:

hru_map <-mapview(
  hru_shp,
  map.types = "Esri.WorldImagery",
  color = "green",
  lwd = 1,
  legend = FALSE,
  zcol = "soil",
  alpha.region = 0
)

problem_map <-
  mapview(
    problem_shps,
    color = "red",
    lwd = 3,
    legend = FALSE,
    zcol = "name",
    layer.name = "Problem HRUs",
    alpha.region = 0
  )

problem_lus_map <- 
  mapview(
    problem_lus,
    color = "orange",
    lwd = 1,
    legend = FALSE,
    zcol = "name",
    layer.name = "Problem LUs",
    alpha.region = 0
  )

full_map <- hru_map + problem_lus_map+ problem_map

full_map

Resolution to this problem has been discussed in issue #53

We will extract the hru and landuse file to modify further.

lum_fixed <- read.table("model_data/cs10_setup/swat_input_mod/landuse.lum", skip = 1, sep = "", header = T) %>% tibble()

hru_fixed <- read.table("model_data/cs10_setup/swat_input/hru-data.hru", skip = 1, sep = "", header = T) %>% tibble()

10.2.2 Fixing the problem HRUs

hru0704 has the wrong soil type for its landuse a_037f_2_drn_lum, which should be STV instead of SDs

hru_fixed$soil[which(hru_fixed$name=="hru0704")] = "STV"

hru1157 has the wrong soil type for its landuse a_067f_1, which should be STR instead of SDs

hru_fixed$soil[which(hru_fixed$name=="hru1157")] = "STR"

hru1233 has the wrong soil type for its landuse a_074f_1, which should be STV instead of Mon

hru_fixed$soil[which(hru_fixed$name=="hru1233")] = "STV"

hru1755 has the wrong land use, while it was classified as a field, it should actually be a road. This is due to misaligned LU maps and should be fixed at the source

hru_fixed$lu_mgt[which(hru_fixed$name=="hru1755")] = "utrn_lum"

hru3523 has the wrong soil type for its landuse a_053f_1, which should be STR instead of RG

hru_fixed$soil[which(hru_fixed$name=="hru3523")] = "STR"

hru5861 of landuse a_210_f2 does not share the same soil (end moraine) as the rest of its field.

hru_fixed$soil[which(hru_fixed$name=="hru5861")] = "STV"

10.2.2.1 Legacy problems

The following HRUs were problematic before we updated our soil map (see issue #68), now they are no longer problematic, therefore the code has been commented out.

hru0218 with landuse a_008f_2_drn_lum – Remove the drains since it is an isolated LU

#lum_fixed$tile[which(lum_fixed$name=="a_008f_2_drn_lum")] = "null"

hru3637 with landuse a_148f_1_drn_lum – Remove the drains since it is an isolated LU.

#lum_fixed$tile[which(lum_fixed$name=="a_148f_1_drn_lum")] = "null"

hru3686 of landuse a_115f_1_drn_lum has the wrong soil type, changing it to "STlv-sl” of its connected field.

#hru_fixed$soil[which(hru_fixed$name=="hru3686")] = "STlv-sl"

hru4880 of landuse a_182f_5_drn_lum has the wrong soil type for the land use. Changing to the “ARdy” of neighboring fields

#hru_fixed$soil[which(hru_fixed$name=="hru4880")] = "ARdy"

hru5033 with isolated landuse a_112f_1_drn_lum can have its drains removed.

#lum_fixed$tile[which(lum_fixed$name=="a_112f_1_drn_lum")] = "null"

hru0692 of isolated land use a_037f_1_drn_lum can safely has its drains removed

#lum_fixed$tile[which(lum_fixed$name=="a_037f_1_drn_lum")] = "null"

hru3314 of landuse a_060f_2_drn_lum has the wrong soil type. It should be STlv-sl, not STlen-rt-sl.

The same is true for hru3376.

#hru_fixed$soil[which(hru_fixed$name=="hru3314")] = "STlv-sl"
#hru_fixed$soil[which(hru_fixed$name=="hru3376")] = "STlv-sl"

hru2784 (and the within hru3367) is a sizable field with landuse a_072f_1_drn_lum and soil RGlen-dy-ar. Tough call here, but best bet is probably to assign soil of the other major field it shares a landuse with (STlv-sl)

#hru_fixed$soil[which(hru_fixed$name=="hru2784")] = "STlv-sl"
#hru_fixed$soil[which(hru_fixed$name=="hru3367")] = "STlv-sl"

hru3465 of LU a_146f_2_drn_lum has the wrong soil for the greater field. Changing it from RGlep-dy to STrt-sl

#hru_fixed$soil[which(hru_fixed$name=="hru3465")] = "STrt-sl"

hru5214 of isolated land use a_020f_4_drn_lum can safely have its drains removed.

#lum_fixed$tile[which(lum_fixed$name=="a_020f_4_drn_lum")] = "null"

hru5926 of landuse a_211f_1_drn_lum has the wrong soil type (Sea-dep-thin), changing to that of the greater field (STrt-sl)

#hru_fixed$soil[which(hru_fixed$name=="hru5926")] = "STlv-sl"

hru2633 of isolated LU a_131f_1_drn_lum can safely have its drains removed.

#lum_fixed$tile[which(lum_fixed$name=="a_131f_1_drn_lum")] = "null"

10.2.3 Checking changes

Do we now have 0 problem HRUs?

shallow_hrus <- which(hru_fixed$soil %in% shallow_soils)

drained_hrus <- grepl(x = hru_fixed$lu_mgt, pattern = "_drn") %>% which()

# Intersection of HRUs which are both shallow and drained
problem_hrus <- base::intersect(shallow_hrus, drained_hrus)

length(problem_hrus)

## [1] 0

Yes

10.2.4 Writing modified changes

Removing the appended columns, and duplicates from the left-join, adding header and writing to modified input files directory.

lu_header <- paste("lu table after fixing up by OPTAIN-CS10 workflow in section 9.2")

write.table(lu_header, "model_data/cs10_setup/swat_input_mod/landuse.lum", quote = F, row.names = F, col.names = F )

write.table(lum_fixed, file = "model_data/cs10_setup/swat_input_mod/landuse.lum", sep = "\t", quote = F, append = T, row.names = F, col.names = T)

## Warning in write.table(lum_fixed, file =
## "model_data/cs10_setup/swat_input_mod/landuse.lum", : appending column names to
## file

Six drains have been removed.

hru_header <- paste("header header header, delete me and regret it forever! fixed soil types by cs10 workflow in section 9.2")

write.table(hru_header, "model_data/cs10_setup/swat_input_mod/hru-data.hru", quote = F, row.names = F, col.names = F )

write.table(hru_fixed, file = "model_data/cs10_setup/swat_input_mod/hru-data.hru", sep = "\t", quote = F, append = T, row.names = F, col.names = T)

## Warning in write.table(hru_fixed, file =
## "model_data/cs10_setup/swat_input_mod/hru-data.hru", : appending column names
## to file

Ten soil types have been changed.

10.2.5 Re-enabling tiledrain

Load from modified input files

old_path <- "model_data/cs10_setup/swat_input_mod/codes.bsn"
new_path <- "model_data/cs10_setup/swat_input_mod/codes.bsn"

codes_bsn <- readLines(old_path)

header <- codes_bsn[1]

header <- paste("modified by CS10 workflow in section 9.2")

codes_df <- read.table(old_path, skip = 1, header = T, sep = "", 
                       colClasses = "character")

Re-enable by setting to 1

codes_df$tiledrain <- 1

Writing to modified input files directory.

write.table(header, file = new_path, col.names = F, row.names = F, quote = F) 

write.table(
  codes_df,
  file = new_path,
  col.names = T,
  append = T,
  quote = F,
  sep = "   ",
  row.names = F
)

## Warning in write.table(codes_df, file = new_path, col.names = T, append = T, :
## appending column names to file

10.2.6 Testing SWAT+

Now we can see if the model runs

source("model_data/code/test_swat_mod.R")

simout <- test_swat_mod()
cat(tail(simout), sep = "\n")

Does it work? No. Of course not because tiledrain=1 is current issue in SWAT+.

10.3 Setting Soil Phosphorous Initial Conditions

# require(dplyr)
# require(sf)
# require(mapview)
# require(ggplot2)
# require(stringr)
# require(DT)
require(magrittr)

sf_theme <- ggplot2::theme(axis.text.x=ggplot2::element_blank(), #remove x axis labels
        axis.ticks.x=ggplot2::element_blank(), #remove x axis ticks
        axis.text.y=ggplot2::element_blank(),  #remove y axis labels
        axis.ticks.y=ggplot2::element_blank()  #remove y axis ticks
        )

In this section we will set the initial value for labile phosphorous in the soil.

SWAT+ reads in inital P values from lab_p which is in a file called nutrients.sol which is read by the nutrients column in the file soil_plant.ini which is read by the soil_plant_init column of the hru-data.hru.

We have a map containing the initial values, created by Dominika using P-AL method, here:

pal <- sf::read_sf("model_data/input/nutrients/cs10_p_al.shp")
ggplot2::ggplot()+ggplot2::geom_sf(pal, mapping = ggplot2::aes(fill = pal))+sf_theme

It does not have full coverage of the basin:

ws <- sf::read_sf("model_data/input/shape/cs10_basin.shp")
ggplot2::ggplot() + ggplot2::geom_sf(ws, mapping = ggplot2::aes(fill = "Watershed"))+
  ggplot2::geom_sf(pal, mapping = ggplot2::aes(fill = "P map")) + sf_theme

Now we need to calculate which HRUs have which amount of lab_p .

hru_shp <- sf::read_sf("model_data/cs10_setup/optain-cs10/data/vector/hru.shp")
# set correct CRS
pal <- sf::st_transform(pal, sf::st_crs(hru_shp))
# join
hru_pal <- sf::st_intersection(hru_shp, pal)

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

# average per hru
hru_pal_avg <- hru_pal %>% dplyr::group_by(id) %>% dplyr::summarize(p_al_mean = mean(pal))
mapview::mapview(hru_pal_avg, zcol = "p_al_mean")

We are finished with the spatial work. We now just need to extract the ID and the lab_p.

# drop geometry and grab only the ID and lab_p
hru_p_df_raw <- hru_pal_avg %>% sf::st_drop_geometry() %>% dplyr::select(id, p_al_mean)
# left join to a complete dataframe
hru_full <- tibble::tibble(id = hru_shp$id)
hru_p_df <- dplyr::left_join(hru_full, hru_p_df_raw , by = "id")
# sort by ID
hru_p_df <- dplyr::arrange(hru_p_df, id)

We have 84 hrus with no data.

 which(hru_p_df$p_al_mean %>% is.na()) %>% length()

## [1] 84

We will fix this by giving the average value of the land use. Calculating average values of land use:

lu_shp <- sf::read_sf("model_data/input/land/CS10_LU.shp")
hru_lu <- sf::st_intersection(hru_pal_avg, lu_shp)

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

generics <- c("frst", "past", "rngb", "urml", "utrn", "watr", "wetf")
hru_lu_generic <- hru_lu %>% dplyr::filter(type %in% generics)
hru_lu_agri <- hru_lu %>% dplyr::filter((type %in% generics)==FALSE) %>% dplyr::summarize (mean = mean(p_al_mean)) %>% sf::st_drop_geometry() %>% dplyr::pull() %>% round(2)
generic_p_avg <- hru_lu_generic %>% dplyr::group_by(type) %>% dplyr::summarise(type_p = mean(p_al_mean)) %>% sf::st_drop_geometry() %>% dplyr::mutate(type_p = round(type_p, 2))
lu_p_avg_df <- rbind(generic_p_avg, c("agri", hru_lu_agri))
DT::datatable(lu_p_avg_df)

Now we need to find the types of the hrus with the missing data and give them the average value, then add those values to the full dataframe.

missing <- hru_p_df$id[which(hru_p_df$p_al_mean %>% is.na())]
missing_df <- hru_p_df %>% dplyr::filter(id %in% missing)
missing_df <- missing_df %>% dplyr::arrange(id)
types <- hru_shp %>% sf::st_drop_geometry() %>% dplyr::filter(id %in% missing) %>% dplyr::arrange(id)
types$type[which((types$type %in% generics) == FALSE)] = "agri"
missing_df$type <- types$type
missing_df <- dplyr::left_join(missing_df, lu_p_avg_df, by = "type") %>% dplyr::mutate(p_al_mean = type_p) %>% dplyr::select(id, p_al_mean)
hru_p_df_no_na <- hru_p_df %>% dplyr::filter(p_al_mean %>% is.na() == FALSE)
hru_p_df <- rbind(hru_p_df_no_na, missing_df) %>% dplyr::arrange(id)

Now we need a soil_plant_ini class for each hru, as well as a nutrient class for each hru.

The soil_plant.ini file requires the following columns:

hru_names <- paste0("spi_hru",hru_p_df$id)
sw_frac <- "0.00000"
nutrients <- paste0("sn_hru", hru_p_df$id)
pest <- "null" 
path <- "null" 
hmet <- "null" 
salt <- "null"
soil_plant_ini <- tibble::tibble(name = hru_names, sw_frac = sw_frac, nutrients = nutrients, pest = pest, path = path, hmet = hmet, salt = salt)
write.table(x = soil_plant_ini, file = "model_data/cs10_setup/swat_input_mod/soil_plant.ini", quote = F, row.names = F)
DT::datatable(soil_plant_ini)

## Warning in instance$preRenderHook(instance): It seems your data is too big for
## client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html

Now we need to add a header and reformat the column names.

header = "soil_plant.ini: last updated by CS10 workflow section XX"
txt <- readLines("model_data/cs10_setup/swat_input_mod/soil_plant.ini")
txt <- c(header, txt)
writeLines(txt, "model_data/cs10_setup/swat_input_mod/soil_plant.ini")

And now we need to modify the nutrients.sol file

name <- paste0("sn_hru", hru_p_df$id)
lab_p = hru_p_df$p_al_mean

cpath <- "model_data/cs10_setup/swat_input_mod/nutrients.sol"

nutrients_df <- read.table(cpath, skip = 1, fill = T, header = T, colClasses = rep("character", 13)) %>% tibble::tibble()

nutrients_df <- nutrients_df[1,]

nutrients_df_add <- tibble::tibble(
  name = name,
  exp_co = nutrients_df$exp_co,
  lab_p = lab_p,
  nitrate = nutrients_df$nitrate,
  fr_hum_act = nutrients_df$fr_hum_act,
  hum_c_n = nutrients_df$hum_c_n,
  hum_c_p = nutrients_df$hum_c_p,
  inorgp = nutrients_df$inorgp,
  watersol_p = nutrients_df$watersol_p,
  h3a_p = nutrients_df$h3a_p,
  mehlich_p = nutrients_df$mehlich_p,
  bray_strong_p = nutrients_df$bray_strong_p,
  description = ""
)
  
nutrients_df_full <- rbind(nutrients_df, nutrients_df_add)

write.table(nutrients_df_full, file = cpath, quote = F, row.names = F)

header = "nutrients.sol: last updated by CS10 workflow section XX"

txt <- readLines(cpath)
txt <- c(header, txt)
writeLines(txt, cpath)

DT::datatable(nutrients_df_full)

Now we need to assign the HRU table in hru-data.hru the correct soil_plant_init column from.

cpath <- "model_data/cs10_setup/swat_input_mod/hru-data.hru"
hru_data <- read.table(cpath, fill = T, header = T, skip = 1) %>% tibble::tibble()

hru_data$soil_plant_init <- hru_names

write.table(hru_data, file = cpath, quote = F, row.names = F)

txt <- readLines(cpath)
txt <- c("hru-data last modified by CS10 workflow @ section XX", txt)
writeLines(txt, con = cpath)

DT::datatable(head(hru_data))

Related issues:

#88, #46, #19