Skip to content

Commit

Permalink
Process review comments by Willem
Browse files Browse the repository at this point in the history
  • Loading branch information
@CFBaptista
CFBaptista committed Apr 17, 2024
1 parent 239c7a0 commit a44dace
Show file tree
Hide file tree
Showing 14 changed files with 119 additions and 120 deletions.
10 changes: 5 additions & 5 deletions docs/src/model_docs/params_lateral.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -158,7 +158,7 @@ in bold represent model parameters that can be set through static input data (ne
soil related parameters `f`, `soilthickness`, `z_exp`, `theta_s` and `theta_r` are derived from the
vertical `SBM` concept (including unit conversion for `f`, `z_exp` and `soilthickness`), and
can be listed in the TOML configuration file under `[input.vertical]`, to map the internal
model parameter to the external netCDF variable. The internal slope model parameter `beta_l` is
model parameter to the external netCDF variable. The internal slope model parameter `slope` is
set through the TOML file as follows:

```toml
Expand Down Expand Up @@ -199,7 +199,7 @@ the `layered_exponential` profile `kv` is used and `z_exp` is required as part o
| **`theta_s`** | saturated water content (porosity) | - | 0.6 |
| **`theta_r`** | residual water content | - | 0.01 |
| `dt` | model time step | d | - |
| **`beta_l`** (`slope`) | slope | m m``^{-1}`` | - |
| **`slope`** | slope | m m``^{-1}`` | - |
| `dl` | drain length | m | - |
| `dw` | drain width | m | - |
| `zi` | pseudo-water table depth (top of the saturated zone) | m | - |
Expand Down Expand Up @@ -256,9 +256,9 @@ model parameter `mannings_n`, and is listed in the Table below between parenthes
| `zb_max` | maximum channel bed elevation | m | - |
| `mannings_n_sq` | Manning's roughness squared at edge/link | (s m``^{-\frac{1}{3}}``)``^2`` | - |
| `h` | water depth | m | - |
| `eta_max` | maximum water elevation | m | - |
| `eta_src` | water elevation of source node of edge | m | - |
| `eta_dst` | water elevation of downstream node of edge | m | - |
| `zs_max` | maximum water elevation | m | - |
| `zs_src` | water elevation of source node of edge | m | - |
| `zs_dst` | water elevation of downstream node of edge | m | - |
| `hf` | water depth at edge/link | m | - |
| `h_av` | average water depth | m | - |
| `dl` | river length | m | - |
Expand Down
8 changes: 4 additions & 4 deletions src/flextopo_model.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -517,17 +517,17 @@ function initialize_flextopo_model(config::Config)
ldd = set_pit_ldd(pits_2d, ldd, inds)
end

beta_l =
landslope =
ncread(nc, config, "lateral.land.slope"; optional = false, sel = inds, type = Float)
clamp!(beta_l, 0.00001, Inf)
clamp!(landslope, 0.00001, Inf)

dl = map(detdrainlength, ldd, xl, yl)
dw = (xl .* yl) ./ dl
olf = initialize_surfaceflow_land(
nc,
config,
inds;
sl = beta_l,
sl = landslope,
dl = dl,
width = map(det_surfacewidth, dw, riverwidth, river),
iterate = kinwave_it,
Expand All @@ -548,7 +548,7 @@ function initialize_flextopo_model(config::Config)

# the indices of the river cells in the land(+river) cell vector
index_river = filter(i -> !isequal(river[i], 0), 1:n)
frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, beta_l, n)
frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, landslope, n)

rf = initialize_surfaceflow_river(
nc,
Expand Down
60 changes: 30 additions & 30 deletions src/flow.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -397,7 +397,7 @@ end
theta_s::Vector{T} | "-" # Saturated water content (porosity) [-]
theta_r::Vector{T} | "-" # Residual water content [-]
dt::T | "d" | 0 | "none" | "none" # model time step [d]
beta_l::Vector{T} | "m m-1" # Slope [m m⁻¹]
slope::Vector{T} | "m m-1" # Slope [m m⁻¹]
dl::Vector{T} | "m" # Drain length [m]
dw::Vector{T} | "m" # Flow width [m]
zi::Vector{T} | "m" # Pseudo-water table depth [m] (top of the saturated zone)
Expand Down Expand Up @@ -446,7 +446,7 @@ function update(ssf::LateralSSF, network, frac_toriver, ksat_profile)
ssf.zi[v],
ssf.recharge[v],
ssf.kh_0[v],
ssf.beta_l[v],
ssf.slope[v],
ssf.theta_s[v] - ssf.theta_r[v],
ssf.f[v],
ssf.soilthickness[v],
Expand All @@ -465,7 +465,7 @@ function update(ssf::LateralSSF, network, frac_toriver, ksat_profile)
ssf.zi[v],
ssf.recharge[v],
ssf.kh[v],
ssf.beta_l[v],
ssf.slope[v],
ssf.theta_s[v] - ssf.theta_r[v],
ssf.soilthickness[v],
ssf.dt,
Expand Down Expand Up @@ -504,9 +504,9 @@ end
mannings_n_sq::Vector{T} | "(s m-1/3)2" | _ | "edge" # Manning's roughness squared at edge/link
mannings_n::Vector{T} | "s m-1/3" # Manning's roughness at node
h::Vector{T} | "m" # water depth
eta_max::Vector{T} | "m" | _ | "edge" # maximum water elevation at edge
eta_src::Vector{T} | "m" # water elevation of source node of edge
eta_dst::Vector{T} | "m" # water elevation of downstream node of edge
zs_max::Vector{T} | "m" | _ | "edge" # maximum water elevation at edge
zs_src::Vector{T} | "m" # water elevation of source node of edge
zs_dst::Vector{T} | "m" # water elevation of downstream node of edge
hf::Vector{T} | "m" | _ | "edge" # water depth at edge/link
h_av::Vector{T} | "m" # average water depth
dl::Vector{T} | "m" # river length
Expand Down Expand Up @@ -685,9 +685,9 @@ function initialize_shallowwater_river(
mannings_n_sq = mannings_n_sq,
mannings_n = n_river,
h = h,
eta_max = zeros(_ne),
eta_src = zeros(_ne),
eta_dst = zeros(_ne),
zs_max = zeros(_ne),
zs_src = zeros(_ne),
zs_dst = zeros(_ne),
hf = zeros(_ne),
h_av = zeros(n),
width = width,
Expand Down Expand Up @@ -736,11 +736,11 @@ function shallowwater_river_update(sw::ShallowWaterRiver, network, dt, doy, upda
i = sw.active_e[j]
i_src = nodes_at_link.src[i]
i_dst = nodes_at_link.dst[i]
sw.eta_src[i] = sw.zb[i_src] + sw.h[i_src]
sw.eta_dst[i] = sw.zb[i_dst] + sw.h[i_dst]
sw.zs_src[i] = sw.zb[i_src] + sw.h[i_src]
sw.zs_dst[i] = sw.zb[i_dst] + sw.h[i_dst]

sw.eta_max[i] = max(sw.eta_src[i], sw.eta_dst[i])
sw.hf[i] = (sw.eta_max[i] - sw.zb_max[i])
sw.zs_max[i] = max(sw.zs_src[i], sw.zs_dst[i])
sw.hf[i] = (sw.zs_max[i] - sw.zb_max[i])

sw.a[i] = sw.width_at_link[i] * sw.hf[i] # flow area (rectangular channel)
sw.r[i] = sw.a[i] / (sw.width_at_link[i] + 2.0 * sw.hf[i]) # hydraulic radius (rectangular channel)
Expand All @@ -749,8 +749,8 @@ function shallowwater_river_update(sw::ShallowWaterRiver, network, dt, doy, upda
sw.hf[i] > sw.h_thresh,
local_inertial_flow(
sw.q0[i],
sw.eta_src[i],
sw.eta_dst[i],
sw.zs_src[i],
sw.zs_dst[i],
sw.hf[i],
sw.a[i],
sw.r[i],
Expand All @@ -770,7 +770,7 @@ function shallowwater_river_update(sw::ShallowWaterRiver, network, dt, doy, upda
sw.q_av[i] += sw.q[i] * dt
end
if !isnothing(sw.floodplain)
@tturbo @. sw.floodplain.hf = max(sw.eta_max - sw.floodplain.zb_max, 0.0)
@tturbo @. sw.floodplain.hf = max(sw.zs_max - sw.floodplain.zb_max, 0.0)

n = 0
@inbounds for i in sw.active_e
Expand Down Expand Up @@ -830,8 +830,8 @@ function shallowwater_river_update(sw::ShallowWaterRiver, network, dt, doy, upda
sw.floodplain.a[i] > 1.0e-05,
local_inertial_flow(
sw.floodplain.q0[i],
sw.eta_src[i],
sw.eta_dst[i],
sw.zs_src[i],
sw.zs_dst[i],
sw.floodplain.hf[i],
sw.floodplain.a[i],
sw.floodplain.r[i],
Expand Down Expand Up @@ -1230,10 +1230,10 @@ function shallowwater_update(
# for flow in x dir) and floodplain flow is not calculated.
if xu <= sw.n && sw.ywidth[i] != T(0.0)

eta_x = sw.z[i] + sw.h[i]
eta_xu = sw.z[xu] + sw.h[xu]
eta_max = max(eta_x, eta_xu)
hf = (eta_max - sw.zx_max[i])
zs_x = sw.z[i] + sw.h[i]
zs_xu = sw.z[xu] + sw.h[xu]
zs_max = max(zs_x, zs_xu)
hf = (zs_max - sw.zx_max[i])

if hf > sw.h_thresh
length = T(0.5) * (sw.xl[i] + sw.xl[xu]) # can be precalculated
Expand All @@ -1242,8 +1242,8 @@ function shallowwater_update(
sw.qx0[i],
sw.qx0[xd],
sw.qx0[xu],
eta_x,
eta_xu,
zs_x,
zs_xu,
hf,
sw.ywidth[i],
length,
Expand All @@ -1270,10 +1270,10 @@ function shallowwater_update(
# for flow in y dir) and floodplain flow is not calculated.
if yu <= sw.n && sw.xwidth[i] != T(0.0)

eta_y = sw.z[i] + sw.h[i]
eta_yu = sw.z[yu] + sw.h[yu]
eta_max = max(eta_y, eta_yu)
hf = (eta_max - sw.zy_max[i])
zs_y = sw.z[i] + sw.h[i]
zs_yu = sw.z[yu] + sw.h[yu]
zs_max = max(zs_y, zs_yu)
hf = (zs_max - sw.zy_max[i])

if hf > sw.h_thresh
length = T(0.5) * (sw.yl[i] + sw.yl[yu]) # can be precalculated
Expand All @@ -1282,8 +1282,8 @@ function shallowwater_update(
sw.qy0[i],
sw.qy0[yd],
sw.qy0[yu],
eta_y,
eta_yu,
zs_y,
zs_yu,
hf,
sw.xwidth[i],
length,
Expand Down
10 changes: 5 additions & 5 deletions src/groundwater/aquifer.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -273,9 +273,9 @@ function conductance(
connectivity.width[nzi],
)
elseif conductivity_profile == "uniform"
phi_i = aquifer.head[i]
phi_j = aquifer.head[j]
if phi_i >= phi_j
head_i = aquifer.head[i]
head_j = aquifer.head[j]
if head_i >= head_j
saturation =
saturated_thickness(aquifer, i) / (aquifer.top[i] - aquifer.bottom[i])
else
Expand All @@ -298,9 +298,9 @@ function flux!(Q, aquifer, connectivity, conductivity_profile)
for nzi in connections(connectivity, i)
# connection from i -> j
j = connectivity.rowval[nzi]
delta_phi = aquifer.head[i] - aquifer.head[j]
delta_head = aquifer.head[i] - aquifer.head[j]
cond = conductance(aquifer, i, j, nzi, conductivity_profile, connectivity)
Q[i] -= cond * delta_phi
Q[i] -= cond * delta_head
end
end
return Q
Expand Down
18 changes: 9 additions & 9 deletions src/groundwater/boundary_conditions.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,16 +25,16 @@ end

function flux!(Q, river::River, aquifer)
for (i, index) in enumerate(river.index)
phi= aquifer.head[index]
head = aquifer.head[index]
stage = river.stage[i]
if stage > phi
if stage > head
cond = river.infiltration_conductance[i]
delta_phi = min(stage - river.bottom[i], stage - phi)
delta_head = min(stage - river.bottom[i], stage - head)
else
cond = river.exfiltration_conductance[i]
delta_phi = stage - phi
delta_head = stage - head
end
river.flux[i] = check_flux(cond * delta_phi, aquifer, index)
river.flux[i] = check_flux(cond * delta_head, aquifer, index)
Q[index] += river.flux[i]
end
return Q
Expand All @@ -53,8 +53,8 @@ end
function flux!(Q, drainage::Drainage, aquifer)
for (i, index) in enumerate(drainage.index)
cond = drainage.conductance[i]
delta_phi = min(0, drainage.elevation[i] - aquifer.head[index])
drainage.flux[i] = check_flux(cond * delta_phi, aquifer, index)
delta_head = min(0, drainage.elevation[i] - aquifer.head[index])
drainage.flux[i] = check_flux(cond * delta_head, aquifer, index)
Q[index] += drainage.flux[i]
end
return Q
Expand All @@ -72,8 +72,8 @@ end
function flux!(Q, headboundary::HeadBoundary, aquifer)
for (i, index) in enumerate(headboundary.index)
cond = headboundary.conductance[i]
delta_phi = headboundary.head[i] - aquifer.head[index]
headboundary.flux[i] = check_flux(cond * delta_phi, aquifer, index)
delta_head = headboundary.head[i] - aquifer.head[index]
headboundary.flux[i] = check_flux(cond * delta_head, aquifer, index)
Q[index] += headboundary.flux[i]
end
return Q
Expand Down
8 changes: 4 additions & 4 deletions src/hbv_model.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -249,17 +249,17 @@ function initialize_hbv_model(config::Config)
ldd = set_pit_ldd(pits_2d, ldd, inds)
end

beta_l =
landslope =
ncread(nc, config, "lateral.land.slope"; optional = false, sel = inds, type = Float)
clamp!(beta_l, 0.00001, Inf)
clamp!(landslope, 0.00001, Inf)

dl = map(detdrainlength, ldd, xl, yl)
dw = (xl .* yl) ./ dl
olf = initialize_surfaceflow_land(
nc,
config,
inds;
sl = beta_l,
sl = landslope,
dl = dl,
width = map(det_surfacewidth, dw, riverwidth, river),
iterate = kinwave_it,
Expand All @@ -282,7 +282,7 @@ function initialize_hbv_model(config::Config)

# the indices of the river cells in the land(+river) cell vector
index_river = filter(i -> !isequal(river[i], 0), 1:n)
frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, beta_l, n)
frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, landslope, n)

rf = initialize_surfaceflow_river(
nc,
Expand Down
Loading

0 comments on commit a44dace

Please sign in to comment.