-
Notifications
You must be signed in to change notification settings - Fork 2
/
groundwater_recharge.R
144 lines (114 loc) · 5.61 KB
/
groundwater_recharge.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
#' Calculate groundwater recharge of a soil
#'
#' This function calculates an index score for groundwater storage based on precipitation surplus, infiltration at saturation, sealing risk, drainage and subsoil compaction
#'
#'
#' @param B_LU_BRP (numeric) The crop code from the BRP
#' @param D_PSP (numeric) The precipitation surplus per crop calculated by \code{\link{calc_psp}}
#' @param D_WRI_K (numeric) The value for top soil permeability (cm/d) as calculated by \code{\link{calc_permeability}}
#' @param I_P_SE (numeric) The indicator value for soil sealing
#' @param I_P_CO (numeric) The indicator value for occurrence of subsoil compaction
#' @param B_DRAIN (boolean) Are drains installed to drain the field (options: yes or no)
#' @param B_GWL_CLASS (character) The groundwater table class
#'
#' @examples
#' ind_gw_recharge(B_LU_BRP = 265,D_PSP = 200, D_WRI_K = 10, I_P_SE = 0.6, I_P_CO = 0.9,
#' B_DRAIN = FALSE, B_GWL_CLASS = 'GtV')
#' ind_gw_recharge(B_LU_BRP = 233, D_PSP = 400, D_WRI_K = 10, I_P_SE = 0.4, I_P_CO = 0.2,
#' B_DRAIN = TRUE, B_GWL_CLASS = 'GtII')
#'
#' @return
#' The evaluated score for the soil function to improve groundwater recharge. A numeric value between 0 and 1.
#'
#' @export
ind_gw_recharge <- function(B_LU_BRP, D_PSP, D_WRI_K, I_P_SE, I_P_CO, B_DRAIN, B_GWL_CLASS){
I_H_GWR = D_I_WRI_K = cf_compaction = cf_drain = D_I_PSP = NULL
# Check inputs
arg.length <- max(length(B_LU_BRP),length(D_WRI_K),length(D_PSP),length(I_P_SE),length(I_P_CO),length(B_DRAIN))
checkmate::assert_numeric(B_LU_BRP, any.missing = FALSE, min.len = 1, len = arg.length)
checkmate::assert_subset(B_LU_BRP, choices = unique(OBIC::crops.obic$crop_code), empty.ok = FALSE)
checkmate::assert_numeric(D_WRI_K, any.missing = FALSE, len = arg.length)
checkmate::assert_numeric(D_PSP, any.missing = FALSE, len = arg.length)
checkmate::assert_numeric(I_P_SE, any.missing = FALSE, len = arg.length)
checkmate::assert_numeric(I_P_CO, any.missing = FALSE, len = arg.length)
checkmate::assert_logical(B_DRAIN, any.missing = FALSE, len = arg.length)
checkmate::assert_subset(B_GWL_CLASS, choices = c('GtI','GtII','GtIII','GtIV','GtV','GtVI','GtVII','GtVIII',
'GtIIb','GtIIIb','GtVb'), empty.ok = FALSE)
# import data into table
dt <- data.table(B_LU_BRP = B_LU_BRP,
D_PSP = D_PSP,
D_WRI_K = D_WRI_K,
I_P_SE = I_P_SE,
I_P_CO = I_P_CO,
B_DRAIN = B_DRAIN,
B_GWL_CLASS = B_GWL_CLASS)
# Calculate the indicator value for water holding capacity
dt[,D_I_WRI_K := ind_permeability(D_WRI_K)]
# Calculate indicator value for precipitation surplus
dt[,D_I_PSP := ind_psp(D_PSP,B_LU_BRP)]
# Correct for subsoil compaction or drainage
dt[B_DRAIN == TRUE & B_GWL_CLASS %in% c('GtIIIb','GtIV'), c('cf_drain','cf_compaction') := list(0.6,1)]
dt[B_DRAIN == FALSE | B_GWL_CLASS %in% c('GtI','GtII','GtIII','GtV','GtVI','GtVII','GtVIII','GtIIb','GtVb'),
c('cf_drain','cf_compaction') := list(1,fifelse(I_P_CO < 0.5,0.8,1))]
# Calculate aggregated score
dt[,I_H_GWR := (0.7 * D_I_PSP + 0.15 * D_I_WRI_K + 0.15 * I_P_SE) * cf_compaction * cf_drain]
# Calculate mean
out <- dt[,I_H_GWR]
return(out)
}
#' Calculate the permeability of the top soil
#'
#' This function calculates the permeability of the top soil
#'
#'
#' @param A_CLAY_MI (numeric) The clay content of the soil (\%)
#' @param A_SAND_MI (numeric) The sand content of the soil (\%)
#' @param A_SILT_MI (numeric) The silt content of the soil (\%)
#' @param A_SOM_LOI (numeric) The organic matter content of the soil (\%)
#'
#'
#' @export
calc_permeability <- function(A_CLAY_MI,A_SAND_MI,A_SILT_MI,A_SOM_LOI){
id = thetaS = thetaR = alfa = n = ksat = l = Pleem = mineral = D_WRI_K = NULL
# Check inputs
arg.length <- max(length(A_CLAY_MI), length(A_SAND_MI),length(A_SILT_MI), length(A_SOM_LOI))
checkmate::assert_numeric(A_CLAY_MI, lower = 0, upper = 100, any.missing = FALSE)
checkmate::assert_numeric(A_SAND_MI, lower = 0, upper = 100, any.missing = FALSE)
checkmate::assert_numeric(A_SILT_MI, lower = 0, upper = 100, any.missing = FALSE)
checkmate::assert_numeric(A_SOM_LOI, lower = 0, upper = 100, any.missing = FALSE)
# Collect data in a table
dt <- data.table(
id = 1:arg.length,
A_CLAY_MI = A_CLAY_MI,
A_SAND_MI = A_SAND_MI,
A_SILT_MI = A_SILT_MI,
A_SOM_LOI = A_SOM_LOI,
value = NA_real_
)
# Calculate mineral fractions
dt[,mineral := A_CLAY_MI + A_SAND_MI + A_SILT_MI]
dt[,A_CLAY_MI := A_CLAY_MI * 100 / mineral]
dt[,A_SILT_MI := A_SILT_MI * 100 / mineral]
dt[,Pleem := A_CLAY_MI + A_SILT_MI]
# Calculate unsaturated permeability
dt[, c("Dichtheid", "thetaR", "thetaS", "alfa", "n", "ksat", "l") := pFpara_ptf_Wosten2001(A_CLAY_MI, Pleem, A_SOM_LOI, 150, 1)]
# Calculate permeabilty, based on Wosten et al., 2001
dt[,D_WRI_K := ksat * (((1 + alfa * thetaR^n)^(1 - 1/n) - alfa * thetaR^(n - 1))^2) / ((1 + alfa * thetaR^n)^((1 - 1/n) * (l + 2)))]
# output
return(dt[,D_WRI_K])
}
#' Calculate the indicator score for the permeability of the top soil
#'
#' This function calculates the indicator score for the permeability of the top soil
#'
#'
#' @param D_WRI_K (numeric) The value for top soil permeability (cm/d) as calculated by \code{\link{calc_permeability}}
#'
#'
#' @export
ind_permeability <- function(D_WRI_K){
# Check inputs
checkmate::assert_numeric(D_WRI_K, lower = 0, any.missing = FALSE)
D_I_WRI_K <- evaluate_logistic(D_WRI_K,0.08,50,0.4)
return(D_I_WRI_K)
}