From 023940fa62735d1bd09b1de180de51b77ec2de2a Mon Sep 17 00:00:00 2001 From: Willem van Verseveld Date: Wed, 27 Nov 2024 12:40:57 +0100 Subject: [PATCH] Re-organize routing files --- src/Wflow.jl | 11 +- src/{reservoir_lake.jl => routing/lake.jl} | 243 ----- src/routing/reservoir.jl | 241 +++++ .../routing_process.jl} | 0 src/routing/subsurface.jl | 244 +++++ src/routing/surface_kinwave.jl | 577 +++++++++++ .../surface_local_inertial.jl} | 906 +----------------- src/routing/surface_routing.jl | 66 ++ src/routing/timestepping.jl | 13 + 9 files changed, 1150 insertions(+), 1151 deletions(-) rename src/{reservoir_lake.jl => routing/lake.jl} (64%) create mode 100644 src/routing/reservoir.jl rename src/{horizontal_process.jl => routing/routing_process.jl} (100%) create mode 100644 src/routing/subsurface.jl create mode 100644 src/routing/surface_kinwave.jl rename src/{flow.jl => routing/surface_local_inertial.jl} (62%) create mode 100644 src/routing/surface_routing.jl create mode 100644 src/routing/timestepping.jl diff --git a/src/Wflow.jl b/src/Wflow.jl index 91555a124..fe850971b 100644 --- a/src/Wflow.jl +++ b/src/Wflow.jl @@ -144,8 +144,12 @@ include("parameters.jl") include("groundwater/connectivity.jl") include("groundwater/aquifer.jl") include("groundwater/boundary_conditions.jl") -include("flow.jl") -include("horizontal_process.jl") +include("routing/timestepping.jl") +include("routing/subsurface.jl") +include("routing/surface_kinwave.jl") +include("routing/surface_local_inertial.jl") +include("routing/surface_routing.jl") +include("routing/routing_process.jl") include("vegetation/rainfall_interception.jl") include("vegetation/canopy.jl") include("snow/snow_process.jl") @@ -158,7 +162,8 @@ include("soil/soil_process.jl") include("sbm.jl") include("demand/water_demand.jl") include("sediment.jl") -include("reservoir_lake.jl") +include("routing/reservoir.jl") +include("routing/lake.jl") include("sbm_model.jl") include("sediment_model.jl") include("sbm_gwf_model.jl") diff --git a/src/reservoir_lake.jl b/src/routing/lake.jl similarity index 64% rename from src/reservoir_lake.jl rename to src/routing/lake.jl index 12da86d7c..76d85230c 100644 --- a/src/reservoir_lake.jl +++ b/src/routing/lake.jl @@ -1,246 +1,3 @@ - -@get_units @grid_loc @with_kw struct ReservoirParameters{T} - dt::T # Model time step [s] - maxvolume::Vector{T} | "m3" # maximum storage (above which water is spilled) [m³] - area::Vector{T} | "m2" # reservoir area [m²] - maxrelease::Vector{T} | "m3 s-1" # maximum amount that can be released if below spillway [m³ s⁻¹] - demand::Vector{T} | "m3 s-1" # minimum (environmental) flow requirement downstream of the reservoir [m³ s⁻¹] - targetminfrac::Vector{T} | "-" # target minimum full fraction (of max storage) [-] - targetfullfrac::Vector{T} | "-" # target fraction full (of max storage -end - -function ReservoirParameters(dataset, config, indices_river, n_river_cells, pits, dt) - # read only reservoir data if reservoirs true - # allow reservoirs only in river cells - # note that these locations are only the reservoir outlet pixels - reslocs = ncread( - dataset, - config, - "lateral.river.reservoir.locs"; - optional = false, - sel = indices_river, - type = Int, - fill = 0, - ) - - # this holds the same ids as reslocs, but covers the entire reservoir - rescoverage_2d = ncread( - dataset, - config, - "lateral.river.reservoir.areas"; - optional = false, - allow_missing = true, - ) - # for each reservoir, a list of 2D indices, needed for getting the mean precipitation - inds_res_cov = Vector{CartesianIndex{2}}[] - - rev_inds_reservoir = zeros(Int, size(rescoverage_2d)) - - # construct a map from the rivers to the reservoirs and - # a map of the reservoirs to the 2D model grid - inds_reservoir_map2river = fill(0, n_river_cells) - inds_res = CartesianIndex{2}[] - rescounter = 0 - for (i, ind) in enumerate(indices_river) - res_id = reslocs[i] - if res_id > 0 - push!(inds_res, ind) - rescounter += 1 - inds_reservoir_map2river[i] = rescounter - rev_inds_reservoir[ind] = rescounter - - # get all indices related to this reservoir outlet - # done in this loop to ensure that the order is equal to the order in the - # SimpleReservoir struct - res_cov = findall(isequal(res_id), rescoverage_2d) - push!(inds_res_cov, res_cov) - end - end - - resdemand = ncread( - dataset, - config, - "lateral.river.reservoir.demand"; - optional = false, - sel = inds_res, - type = Float, - fill = 0, - ) - resmaxrelease = ncread( - dataset, - config, - "lateral.river.reservoir.maxrelease"; - optional = false, - sel = inds_res, - type = Float, - fill = 0, - ) - resmaxvolume = ncread( - dataset, - config, - "lateral.river.reservoir.maxvolume"; - optional = false, - sel = inds_res, - type = Float, - fill = 0, - ) - resarea = ncread( - dataset, - config, - "lateral.river.reservoir.area"; - optional = false, - sel = inds_res, - type = Float, - fill = 0, - ) - res_targetfullfrac = ncread( - dataset, - config, - "lateral.river.reservoir.targetfullfrac"; - optional = false, - sel = inds_res, - type = Float, - fill = 0, - ) - res_targetminfrac = ncread( - dataset, - config, - "lateral.river.reservoir.targetminfrac"; - optional = false, - sel = inds_res, - type = Float, - fill = 0, - ) - - # for surface water routing reservoir locations are considered pits in the flow network - # all upstream flow goes to the river and flows into the reservoir - pits[inds_res] .= true - - reservoir_network = ( - indices_outlet = inds_res, - indices_coverage = inds_res_cov, - reverse_indices = rev_inds_reservoir, - river_indices = findall(x -> x ≠ 0, inds_reservoir_map2river), - ) - - parameters = ReservoirParameters{Float}(; - dt = dt, - demand = resdemand, - maxrelease = resmaxrelease, - maxvolume = resmaxvolume, - area = resarea, - targetfullfrac = res_targetfullfrac, - targetminfrac = res_targetminfrac, - ) - - return parameters, reservoir_network, inds_reservoir_map2river, pits -end - -@get_units @grid_loc @with_kw struct ReservoirVariables{T} - volume::Vector{T} | "m3" # reservoir volume [m³] - outflow::Vector{T} | "m3 s-1" # outflow from reservoir [m³ s⁻¹] - totaloutflow::Vector{T} | "m3" # total outflow from reservoir [m³] - percfull::Vector{T} | "-" # fraction full (of max storage) [-] - demandrelease::Vector{T} | "m3 s-1" # minimum (environmental) flow released from reservoir [m³ s⁻¹] - actevap::Vector{T} # average actual evaporation for reservoir area [mm Δt⁻¹] -end - -function ReservoirVariables(n, parameters) - (; targetfullfrac, maxvolume) = parameters - variables = ReservoirVariables{Float}(; - volume = targetfullfrac .* maxvolume, - outflow = fill(mv, n), - totaloutflow = fill(mv, n), - percfull = fill(mv, n), - demandrelease = fill(mv, n), - actevap = fill(mv, n), - ) - return variables -end - -@get_units @grid_loc @with_kw struct ReservoirBC{T} - inflow::Vector{T} | "m3" # total inflow into reservoir [m³] - precipitation::Vector{T} # average precipitation for reservoir area [mm Δt⁻¹] - evaporation::Vector{T} # average potential evaporation for reservoir area [mm Δt⁻¹] -end - -function ReservoirBC(n) - bc = ReservoirBC{Float}(; - inflow = fill(mv, n), - precipitation = fill(mv, n), - evaporation = fill(mv, n), - ) - return bc -end - -@with_kw struct SimpleReservoir{T} - boundary_conditions::ReservoirBC{T} - parameters::ReservoirParameters{T} - variables::ReservoirVariables{T} -end - -function SimpleReservoir(dataset, config, indices_river, n_river_cells, pits, dt) - parameters, reservoir_network, inds_reservoir_map2river, pits = - ReservoirParameters(dataset, config, indices_river, n_river_cells, pits, dt) - - n_reservoirs = length(parameters.area) - @info "Read `$n_reservoirs` reservoir locations." - - variables = ReservoirVariables(n_reservoirs, parameters) - boundary_conditions = ReservoirBC(n_reservoirs) - reservoir = SimpleReservoir{Float}(; boundary_conditions, parameters, variables) - - return reservoir, reservoir_network, inds_reservoir_map2river, pits -end - -""" -Update a single reservoir at position `i`. - -This is called from within the kinematic wave loop, therefore updating only for a single -element rather than all at once. -""" -function update!(model::SimpleReservoir, i, inflow, timestepsecs) - res_bc = model.boundary_conditions - res_p = model.parameters - res_v = model.variables - - # limit lake evaporation based on total available volume [m³] - precipitation = - 0.001 * res_bc.precipitation[i] * (timestepsecs / res_p.dt) * res_p.area[i] - available_volume = res_v.volume[i] + inflow * timestepsecs + precipitation - evap = 0.001 * res_bc.evaporation[i] * (timestepsecs / res_p.dt) * res_p.area[i] - actevap = min(available_volume, evap) # [m³/timestepsecs] - - vol = res_v.volume[i] + (inflow * timestepsecs) + precipitation - actevap - vol = max(vol, 0.0) - - percfull = vol / res_p.maxvolume[i] - # first determine minimum (environmental) flow using a simple sigmoid curve to scale for target level - fac = scurve(percfull, res_p.targetminfrac[i], Float(1.0), Float(30.0)) - demandrelease = min(fac * res_p.demand[i] * timestepsecs, vol) - vol = vol - demandrelease - - wantrel = max(0.0, vol - (res_p.maxvolume[i] * res_p.targetfullfrac[i])) - # Assume extra maximum Q if spilling - overflow_q = max((vol - res_p.maxvolume[i]), 0.0) - torelease = - min(wantrel, overflow_q + res_p.maxrelease[i] * timestepsecs - demandrelease) - vol = vol - torelease - outflow = torelease + demandrelease - percfull = vol / res_p.maxvolume[i] - - # update values in place - res_v.outflow[i] = outflow / timestepsecs - res_bc.inflow[i] += inflow * timestepsecs - res_v.totaloutflow[i] += outflow - res_v.demandrelease[i] = demandrelease / timestepsecs - res_v.percfull[i] = percfull - res_v.volume[i] = vol - res_v.actevap[i] += 1000.0 * (actevap / res_p.area[i]) - - return nothing -end - @get_units @grid_loc @with_kw struct LakeParameters{T} dt::T # Model time step [s] lowerlake_ind::Vector{Int} | "-" # Index of lower lake (linked lakes) diff --git a/src/routing/reservoir.jl b/src/routing/reservoir.jl new file mode 100644 index 000000000..02d998fa9 --- /dev/null +++ b/src/routing/reservoir.jl @@ -0,0 +1,241 @@ +@get_units @grid_loc @with_kw struct ReservoirParameters{T} + dt::T # Model time step [s] + maxvolume::Vector{T} | "m3" # maximum storage (above which water is spilled) [m³] + area::Vector{T} | "m2" # reservoir area [m²] + maxrelease::Vector{T} | "m3 s-1" # maximum amount that can be released if below spillway [m³ s⁻¹] + demand::Vector{T} | "m3 s-1" # minimum (environmental) flow requirement downstream of the reservoir [m³ s⁻¹] + targetminfrac::Vector{T} | "-" # target minimum full fraction (of max storage) [-] + targetfullfrac::Vector{T} | "-" # target fraction full (of max storage +end + +function ReservoirParameters(dataset, config, indices_river, n_river_cells, pits, dt) + # read only reservoir data if reservoirs true + # allow reservoirs only in river cells + # note that these locations are only the reservoir outlet pixels + reslocs = ncread( + dataset, + config, + "lateral.river.reservoir.locs"; + optional = false, + sel = indices_river, + type = Int, + fill = 0, + ) + + # this holds the same ids as reslocs, but covers the entire reservoir + rescoverage_2d = ncread( + dataset, + config, + "lateral.river.reservoir.areas"; + optional = false, + allow_missing = true, + ) + # for each reservoir, a list of 2D indices, needed for getting the mean precipitation + inds_res_cov = Vector{CartesianIndex{2}}[] + + rev_inds_reservoir = zeros(Int, size(rescoverage_2d)) + + # construct a map from the rivers to the reservoirs and + # a map of the reservoirs to the 2D model grid + inds_reservoir_map2river = fill(0, n_river_cells) + inds_res = CartesianIndex{2}[] + rescounter = 0 + for (i, ind) in enumerate(indices_river) + res_id = reslocs[i] + if res_id > 0 + push!(inds_res, ind) + rescounter += 1 + inds_reservoir_map2river[i] = rescounter + rev_inds_reservoir[ind] = rescounter + + # get all indices related to this reservoir outlet + # done in this loop to ensure that the order is equal to the order in the + # SimpleReservoir struct + res_cov = findall(isequal(res_id), rescoverage_2d) + push!(inds_res_cov, res_cov) + end + end + + resdemand = ncread( + dataset, + config, + "lateral.river.reservoir.demand"; + optional = false, + sel = inds_res, + type = Float, + fill = 0, + ) + resmaxrelease = ncread( + dataset, + config, + "lateral.river.reservoir.maxrelease"; + optional = false, + sel = inds_res, + type = Float, + fill = 0, + ) + resmaxvolume = ncread( + dataset, + config, + "lateral.river.reservoir.maxvolume"; + optional = false, + sel = inds_res, + type = Float, + fill = 0, + ) + resarea = ncread( + dataset, + config, + "lateral.river.reservoir.area"; + optional = false, + sel = inds_res, + type = Float, + fill = 0, + ) + res_targetfullfrac = ncread( + dataset, + config, + "lateral.river.reservoir.targetfullfrac"; + optional = false, + sel = inds_res, + type = Float, + fill = 0, + ) + res_targetminfrac = ncread( + dataset, + config, + "lateral.river.reservoir.targetminfrac"; + optional = false, + sel = inds_res, + type = Float, + fill = 0, + ) + + # for surface water routing reservoir locations are considered pits in the flow network + # all upstream flow goes to the river and flows into the reservoir + pits[inds_res] .= true + + reservoir_network = ( + indices_outlet = inds_res, + indices_coverage = inds_res_cov, + reverse_indices = rev_inds_reservoir, + river_indices = findall(x -> x ≠ 0, inds_reservoir_map2river), + ) + + parameters = ReservoirParameters{Float}(; + dt = dt, + demand = resdemand, + maxrelease = resmaxrelease, + maxvolume = resmaxvolume, + area = resarea, + targetfullfrac = res_targetfullfrac, + targetminfrac = res_targetminfrac, + ) + + return parameters, reservoir_network, inds_reservoir_map2river, pits +end + +@get_units @grid_loc @with_kw struct ReservoirVariables{T} + volume::Vector{T} | "m3" # reservoir volume [m³] + outflow::Vector{T} | "m3 s-1" # outflow from reservoir [m³ s⁻¹] + totaloutflow::Vector{T} | "m3" # total outflow from reservoir [m³] + percfull::Vector{T} | "-" # fraction full (of max storage) [-] + demandrelease::Vector{T} | "m3 s-1" # minimum (environmental) flow released from reservoir [m³ s⁻¹] + actevap::Vector{T} # average actual evaporation for reservoir area [mm Δt⁻¹] +end + +function ReservoirVariables(n, parameters) + (; targetfullfrac, maxvolume) = parameters + variables = ReservoirVariables{Float}(; + volume = targetfullfrac .* maxvolume, + outflow = fill(mv, n), + totaloutflow = fill(mv, n), + percfull = fill(mv, n), + demandrelease = fill(mv, n), + actevap = fill(mv, n), + ) + return variables +end + +@get_units @grid_loc @with_kw struct ReservoirBC{T} + inflow::Vector{T} | "m3" # total inflow into reservoir [m³] + precipitation::Vector{T} # average precipitation for reservoir area [mm Δt⁻¹] + evaporation::Vector{T} # average potential evaporation for reservoir area [mm Δt⁻¹] +end + +function ReservoirBC(n) + bc = ReservoirBC{Float}(; + inflow = fill(mv, n), + precipitation = fill(mv, n), + evaporation = fill(mv, n), + ) + return bc +end + +@with_kw struct SimpleReservoir{T} + boundary_conditions::ReservoirBC{T} + parameters::ReservoirParameters{T} + variables::ReservoirVariables{T} +end + +function SimpleReservoir(dataset, config, indices_river, n_river_cells, pits, dt) + parameters, reservoir_network, inds_reservoir_map2river, pits = + ReservoirParameters(dataset, config, indices_river, n_river_cells, pits, dt) + + n_reservoirs = length(parameters.area) + @info "Read `$n_reservoirs` reservoir locations." + + variables = ReservoirVariables(n_reservoirs, parameters) + boundary_conditions = ReservoirBC(n_reservoirs) + reservoir = SimpleReservoir{Float}(; boundary_conditions, parameters, variables) + + return reservoir, reservoir_network, inds_reservoir_map2river, pits +end + +""" +Update a single reservoir at position `i`. + +This is called from within the kinematic wave loop, therefore updating only for a single +element rather than all at once. +""" +function update!(model::SimpleReservoir, i, inflow, timestepsecs) + res_bc = model.boundary_conditions + res_p = model.parameters + res_v = model.variables + + # limit lake evaporation based on total available volume [m³] + precipitation = + 0.001 * res_bc.precipitation[i] * (timestepsecs / res_p.dt) * res_p.area[i] + available_volume = res_v.volume[i] + inflow * timestepsecs + precipitation + evap = 0.001 * res_bc.evaporation[i] * (timestepsecs / res_p.dt) * res_p.area[i] + actevap = min(available_volume, evap) # [m³/timestepsecs] + + vol = res_v.volume[i] + (inflow * timestepsecs) + precipitation - actevap + vol = max(vol, 0.0) + + percfull = vol / res_p.maxvolume[i] + # first determine minimum (environmental) flow using a simple sigmoid curve to scale for target level + fac = scurve(percfull, res_p.targetminfrac[i], Float(1.0), Float(30.0)) + demandrelease = min(fac * res_p.demand[i] * timestepsecs, vol) + vol = vol - demandrelease + + wantrel = max(0.0, vol - (res_p.maxvolume[i] * res_p.targetfullfrac[i])) + # Assume extra maximum Q if spilling + overflow_q = max((vol - res_p.maxvolume[i]), 0.0) + torelease = + min(wantrel, overflow_q + res_p.maxrelease[i] * timestepsecs - demandrelease) + vol = vol - torelease + outflow = torelease + demandrelease + percfull = vol / res_p.maxvolume[i] + + # update values in place + res_v.outflow[i] = outflow / timestepsecs + res_bc.inflow[i] += inflow * timestepsecs + res_v.totaloutflow[i] += outflow + res_v.demandrelease[i] = demandrelease / timestepsecs + res_v.percfull[i] = percfull + res_v.volume[i] = vol + res_v.actevap[i] += 1000.0 * (actevap / res_p.area[i]) + + return nothing +end \ No newline at end of file diff --git a/src/horizontal_process.jl b/src/routing/routing_process.jl similarity index 100% rename from src/horizontal_process.jl rename to src/routing/routing_process.jl diff --git a/src/routing/subsurface.jl b/src/routing/subsurface.jl new file mode 100644 index 000000000..7a4b31dc4 --- /dev/null +++ b/src/routing/subsurface.jl @@ -0,0 +1,244 @@ +@get_units @grid_loc struct KhExponential{T} + # Horizontal hydraulic conductivity at soil surface [m d⁻¹] + kh_0::Vector{T} | "m d-1" + # A scaling parameter [m⁻¹] (controls exponential decline of kh_0) + f::Vector{T} | "m-1" +end + +@get_units @grid_loc struct KhExponentialConstant{T} + # Exponential horizontal hydraulic conductivity profile type + exponential::KhExponential + # Depth [m] from soil surface for which exponential decline of kv_0 is valid + z_exp::Vector{T} | "m" +end + +@get_units @grid_loc struct KhLayered{T} + # Horizontal hydraulic conductivity [m d⁻¹] + kh::Vector{T} | "m d-1" +end + +abstract type SubsurfaceFlow end + +@get_units @grid_loc @with_kw struct LateralSsfParameters{T, Kh} + kh_profile::Kh # Horizontal hydraulic conductivity profile type [-] + khfrac::Vector{T} | "-" # A muliplication factor applied to vertical hydraulic conductivity `kv` [-] + soilthickness::Vector{T} | "m" # Soil thickness [m] + theta_s::Vector{T} | "-" # Saturated water content (porosity) [-] + theta_r::Vector{T} | "-" # Residual water content [-] + dt::T # model time step [d] + slope::Vector{T} | "m m-1" # Slope [m m⁻¹] + flow_length::Vector{T} | "m" # Flow length [m] + flow_width::Vector{T} | "m" # Flow width [m] +end + +function LateralSsfParameters( + dataset, + config, + indices; + soil, + slope, + flow_length, + flow_width, + dt, +) + khfrac = ncread( + dataset, + config, + "lateral.subsurface.ksathorfrac"; + sel = indices, + defaults = 1.0, + type = Float, + ) + n_cells = length(khfrac) + + (; theta_s, theta_r, soilthickness) = soil + soilthickness = soilthickness .* 0.001 + + kh_profile_type = get(config.input.vertical, "ksat_profile", "exponential")::String + if kh_profile_type == "exponential" + (; kv_0, f) = soil.kv_profile + kh_0 = khfrac .* kv_0 .* 0.001 .* dt + kh_profile = KhExponential(kh_0, f .* 1000.0) + elseif kh_profile_type == "exponential_constant" + (; z_exp) = soil.kv_profile + (; kv_0, f) = soil.kv_profile.exponential + kh_0 = khfrac .* kv_0 .* 0.001 .* dt + exp_profile = KhExponential(kh_0, f .* 1000.0) + kh_profile = KhExponentialConstant(exp_profile, z_exp .* 0.001) + elseif kh_profile_type == "layered" || kh_profile_type == "layered_exponential" + kh_profile = KhLayered(fill(mv, n_cells)) + end + parameters = LateralSsfParameters( + kh_profile, + khfrac, + soilthickness, + theta_s, + theta_r, + dt, + slope, + flow_length, + flow_width, + ) + return parameters +end + +@get_units @grid_loc @with_kw struct LateralSsfVariables{T} + zi::Vector{T} | "m" # Pseudo-water table depth [m] (top of the saturated zone) + exfiltwater::Vector{T} | "m dt-1" # Exfiltration [m Δt⁻¹] (groundwater above surface level, saturated excess conditions) + recharge::Vector{T} | "m2 dt-1" # Net recharge to saturated store [m² Δt⁻¹] + ssf::Vector{T} | "m3 d-1" # Subsurface flow [m³ d⁻¹] + ssfin::Vector{T} | "m3 d-1" # Inflow from upstream cells [m³ d⁻¹] + ssfmax::Vector{T} | "m2 d-1" # Maximum subsurface flow [m² d⁻¹] + to_river::Vector{T} | "m3 d-1" # Part of subsurface flow [m³ d⁻¹] that flows to the river + volume::Vector{T} | "m3" # Subsurface volume [m³] +end + +function LateralSsfVariables(ssf, zi, xl, yl) + n = length(zi) + volume = @. (ssf.theta_s - ssf.theta_r) * (ssf.soilthickness - zi) * (xl * yl) + variables = LateralSsfVariables(; + zi, + exfiltwater = fill(mv, n), + recharge = fill(mv, n), + ssf = fill(mv, n), + ssfin = fill(mv, n), + ssfmax = fill(mv, n), + to_river = zeros(n), + volume, + ) + return variables +end + +@get_units @grid_loc @with_kw struct LateralSsfBC{T} + recharge::Vector{T} | "m2 dt-1" # Net recharge to saturated store [m² Δt⁻¹] +end + +@with_kw struct LateralSSF{T, Kh} <: SubsurfaceFlow + boundary_conditions::LateralSsfBC{T} + parameters::LateralSsfParameters{T, Kh} + variables::LateralSsfVariables{T} +end + +function LateralSSF( + dataset, + config, + indices; + soil, + slope, + flow_length, + flow_width, + x_length, + y_length, + dt, +) + parameters = LateralSsfParameters( + dataset, + config, + indices; + soil = soil.parameters, + slope, + flow_length, + flow_width, + dt, + ) + zi = 0.001 * soil.variables.zi + variables = LateralSsfVariables(parameters, zi, x_length, y_length) + boundary_conditions = LateralSsfBC(; recharge = fill(mv, length(zi))) + ssf = LateralSSF(; boundary_conditions, parameters, variables) + return ssf +end + +function update!(model::LateralSSF, network) + (; + order_of_subdomains, + order_subdomain, + subdomain_indices, + upstream_nodes, + area, + frac_to_river, + ) = network + + (; recharge) = model.boundary_conditions + (; ssfin, ssf, ssfmax, to_river, zi, exfiltwater, volume) = model.variables + (; slope, theta_s, theta_r, soilthickness, flow_length, flow_width, dt, kh_profile) = + model.parameters + + ns = length(order_of_subdomains) + for k in 1:ns + threaded_foreach(eachindex(order_of_subdomains[k]); basesize = 1) do i + m = order_of_subdomains[k][i] + for (n, v) in zip(subdomain_indices[m], order_subdomain[m]) + # for a river cell without a reservoir or lake part of the upstream + # subsurface flow goes to the river (frac_to_river) and part goes to the + # subsurface flow reservoir (1.0 - frac_to_river) upstream nodes with a + # reservoir or lake are excluded + ssfin[v] = sum_at( + i -> ssf[i] * (1.0 - frac_to_river[i]), + upstream_nodes[n], + eltype(ssfin), + ) + to_river[v] = sum_at( + i -> ssf[i] * frac_to_river[i], + upstream_nodes[n], + eltype(to_river), + ) + ssf[v], zi[v], exfiltwater[v] = kinematic_wave_ssf( + ssfin[v], + ssf[v], + zi[v], + recharge[v], + slope[v], + theta_s[v] - theta_r[v], + soilthickness[v], + dt, + flow_length[v], + flow_width[v], + ssfmax[v], + kh_profile, + v, + ) + volume[v] = (theta_s[v] - theta_r[v]) * (soilthickness[v] - zi[v]) * area[v] + end + end + end + return nothing +end + +@get_units@grid_loc @with_kw struct GroundwaterExchangeVariables{T} + exfiltwater::Vector{T} | "m dt-1" # Exfiltration [m Δt⁻¹] (groundwater above surface level, saturated excess conditions) + zi::Vector{T} | "m" # Pseudo-water table depth [m] (top of the saturated zone) + to_river::Vector{T} | "m3 d-1" # Part of subsurface flow [m³ d⁻¹] that flows to the river + ssf::Vector{T} | "m3 d-1" # Subsurface flow [m³ d⁻¹] +end + +function GroundwaterExchangeVariables(n) + variables = GroundwaterExchangeVariables{Float}(; + exfiltwater = fill(mv, n), + zi = fill(mv, n), + to_river = fill(mv, n), + ssf = zeros(n), + ) + return variables +end + +@with_kw struct GroundwaterExchangeParameters{T} + dt::T # model time step [d] +end + +@with_kw struct GroundwaterExchange{T} <: SubsurfaceFlow + parameters::GroundwaterExchangeParameters{T} + variables::GroundwaterExchangeVariables{T} +end + +function GroundwaterExchange(n, dt) + parameters = GroundwaterExchangeParameters{Float}(; dt = dt / basetimestep) + variables = GroundwaterExchangeVariables(n) + ssf = GroundwaterExchange{Float}(; parameters, variables) + return ssf +end + +get_water_depth(subsurface::SubsurfaceFlow) = subsurface.variables.zi +get_exfiltwater(subsurface::SubsurfaceFlow) = subsurface.variables.exfiltwater + +get_flux_to_river(subsurface::SubsurfaceFlow) = + subsurface.variables.to_river ./ tosecond(basetimestep) # [m³ s⁻¹] \ No newline at end of file diff --git a/src/routing/surface_kinwave.jl b/src/routing/surface_kinwave.jl new file mode 100644 index 000000000..de94cb14c --- /dev/null +++ b/src/routing/surface_kinwave.jl @@ -0,0 +1,577 @@ +abstract type AbstractRiverFlowModel end + +@get_units @grid_loc @with_kw struct FlowVariables{T} + q::Vector{T} | "m3 s-1" # Discharge [m³ s⁻¹] + qlat::Vector{T} | "m2 s-1" # Lateral inflow per unit length [m² s⁻¹] + qin::Vector{T} | "m3 s-1" # Inflow from upstream cells [m³ s⁻¹] + q_av::Vector{T} | "m3 s-1" # Average discharge [m³ s⁻¹] + volume::Vector{T} | "m3" # Kinematic wave volume [m³] (based on water depth h) + h::Vector{T} | "m" # Water depth [m] + h_av::Vector{T} | "m" # Average water depth [m] +end + +function init_kinematic_wave_timestepping(config, n; domain, dt_fixed) + adaptive = get(config.model, "kin_wave_iteration", false)::Bool + @info "Kinematic wave approach is used for $domain flow, adaptive timestepping = $adaptive." + + if adaptive + stable_timesteps = zeros(n) + timestepping = TimeStepping(; stable_timesteps, adaptive) + else + dt_fixed = get(config.model, "kw_$(domain)_tstep", dt_fixed) + @info "Using a fixed sub-timestep (seconds) $dt_fixed for kinematic wave $domain flow." + timestepping = TimeStepping(; dt_fixed, adaptive) + end + return timestepping +end + +function FlowVariables(n) + variables = FlowVariables(; + q = zeros(Float, n), + qlat = zeros(Float, n), + qin = zeros(Float, n), + q_av = zeros(Float, n), + volume = zeros(Float, n), + h = zeros(Float, n), + h_av = zeros(Float, n), + ) + return variables +end + +@get_units @grid_loc @with_kw struct ManningFlowParameters{T} + beta::T # constant in Manning's equation [-] + slope::Vector{T} | "m m-1" # Slope [m m⁻¹] + mannings_n::Vector{T} | "s m-1/3" # Manning's roughness [s m⁻⅓] + flow_length::Vector{T} | "m" # Flow length [m] + flow_width::Vector{T} | "m" # Flow width [m] + alpha_pow::T # Used in the power part of alpha [-] + alpha_term::Vector{T} | "-" # Term used in computation of alpha [-] + alpha::Vector{T} | "s3/5 m1/5" # Constant in momentum equation A = alpha*Q^beta, based on Manning's equation +end + +function ManningFlowParameters(slope, mannings_n, flow_length, flow_width) + n = length(slope) + parameters = ManningFlowParameters(; + beta = Float(0.6), + slope, + mannings_n, + flow_length, + flow_width, + alpha_pow = Float((2.0 / 3.0) * 0.6), + alpha_term = fill(mv, n), + alpha = fill(mv, n), + ) + return parameters +end + +@get_units @grid_loc @with_kw struct RiverFlowParameters{T} + flow::ManningFlowParameters{T} + bankfull_depth::Vector{T} | "m" # Bankfull water level [m] +end + +function Base.getproperty(v::RiverFlowParameters, s::Symbol) + if s === :bankfull_depth + getfield(v, s) + elseif s === :flow + getfield(v, :flow) + else + getfield(getfield(v, :flow), s) + end +end + +function RiverFlowParameters(dataset, config, indices, river_length, river_width) + mannings_n = ncread( + dataset, + config, + "lateral.river.mannings_n"; + sel = indices, + defaults = 0.036, + type = Float, + ) + bankfull_depth = ncread( + dataset, + config, + "lateral.river.bankfull_depth"; + alias = "lateral.river.h_bankfull", + sel = indices, + defaults = 1.0, + type = Float, + ) + if haskey(config.input.lateral.river, "h_bankfull") + @warn string( + "The `h_bankfull` key in `[input.lateral.river]` is now called ", + "`bankfull_depth`. Please update your TOML file.", + ) + end + slope = ncread( + dataset, + config, + "lateral.river.slope"; + optional = false, + sel = indices, + type = Float, + ) + clamp!(slope, 0.00001, Inf) + + flow_parameter_set = ManningFlowParameters(slope, mannings_n, river_length, river_width) + parameters = + RiverFlowParameters(; flow = flow_parameter_set, bankfull_depth = bankfull_depth) + return parameters +end + +@get_units @grid_loc @with_kw struct RiverFlowBC{T, R, L} + inwater::Vector{T} | "m3 s-1" # Lateral inflow [m³ s⁻¹] + inflow::Vector{T} | "m3 s-1" # External inflow (abstraction/supply/demand) [m³ s⁻¹] + inflow_waterbody::Vector{T} | "m3 s-1" # inflow waterbody (lake or reservoir model) from land part [m³ s⁻¹] + abstraction::Vector{T} | "m3 s-1" # Abstraction (computed as part of water demand and allocation) [m³ s⁻¹] + reservoir::R # Reservoir model struct of arrays + lake::L # Lake model struct of arrays +end + +function RiverFlowBC(n, reservoir, lake) + bc = RiverFlowBC(; + inwater = zeros(Float, n), + inflow = zeros(Float, n), + inflow_waterbody = zeros(Float, n), + abstraction = zeros(Float, n), + reservoir = reservoir, + lake = lake, + ) + return bc +end + +@with_kw struct SurfaceFlowRiver{T, R, L, A} <: AbstractRiverFlowModel + timestepping::TimeStepping{T} + boundary_conditions::RiverFlowBC{T, R, L} + parameters::RiverFlowParameters{T} + variables::FlowVariables{T} + allocation::A # Water allocation +end + +function SurfaceFlowRiver( + dataset, + config, + indices; + river_length, + river_width, + reservoir, + lake, +) + n = length(indices) + + timestepping = + init_kinematic_wave_timestepping(config, n; domain = "river", dt_fixed = 900.0) + + do_water_demand = haskey(config.model, "water_demand") + allocation = do_water_demand ? AllocationRiver(n) : NoAllocationRiver{Float}() + + variables = FlowVariables(n) + parameters = RiverFlowParameters(dataset, config, indices, river_length, river_width) + boundary_conditions = RiverFlowBC(n, reservoir, lake) + + sf_river = SurfaceFlowRiver(; + timestepping, + boundary_conditions, + parameters, + variables, + allocation, + ) + + return sf_river +end + +@get_units @grid_loc @with_kw struct LandFlowVariables{T} + flow::FlowVariables{T} + to_river::Vector{T} | "m3 s-1" # Part of overland flow [m³ s⁻¹] that flows to the river +end + +function Base.getproperty(v::LandFlowVariables, s::Symbol) + if s === :to_river + getfield(v, s) + elseif s === :flow + getfield(v, :flow) + else + getfield(getfield(v, :flow), s) + end +end + +@get_units @grid_loc @with_kw struct LandFlowBC{T} + inwater::Vector{T} | "m3 s-1" # Lateral inflow [m³ s⁻¹] +end + +@with_kw struct SurfaceFlowLand{T} + timestepping::TimeStepping{T} + boundary_conditions::LandFlowBC{T} + parameters::ManningFlowParameters{T} + variables::LandFlowVariables{T} +end + +function SurfaceFlowLand(dataset, config, indices; slope, flow_length, flow_width) + mannings_n = ncread( + dataset, + config, + "lateral.land.mannings_n"; + sel = indices, + defaults = 0.072, + type = Float, + ) + + n = length(indices) + timestepping = + init_kinematic_wave_timestepping(config, n; domain = "land", dt_fixed = 3600.0) + + variables = LandFlowVariables(; flow = FlowVariables(n), to_river = zeros(Float, n)) + parameters = ManningFlowParameters(slope, mannings_n, flow_length, flow_width) + boundary_conditions = LandFlowBC(; inwater = zeros(Float, n)) + sf_land = SurfaceFlowLand(; timestepping, boundary_conditions, variables, parameters) + + return sf_land +end + +function surfaceflow_land_update!(model::SurfaceFlowLand, network, dt) + (; + order_of_subdomains, + order_subdomain, + subdomain_indices, + upstream_nodes, + frac_to_river, + ) = network + + (; beta, alpha, flow_width, flow_length) = model.parameters + (; h, h_av, q, q_av, qin, qlat, to_river) = model.variables + + ns = length(order_of_subdomains) + qin .= 0.0 + for k in 1:ns + threaded_foreach(eachindex(order_of_subdomains[k]); basesize = 1) do i + m = order_of_subdomains[k][i] + for (n, v) in zip(subdomain_indices[m], order_subdomain[m]) + # for a river cell without a reservoir or lake part of the upstream + # surface flow goes to the river (frac_to_river) and part goes to the + # surface flow reservoir (1.0 - frac_to_river), upstream nodes with a + # reservoir or lake are excluded + to_river[v] += + sum_at( + i -> q[i] * frac_to_river[i], + upstream_nodes[n], + eltype(to_river), + ) * dt + if flow_width[v] > 0.0 + qin[v] = sum_at( + i -> q[i] * (1.0 - frac_to_river[i]), + upstream_nodes[n], + eltype(q), + ) + end + + q[v] = kinematic_wave( + qin[v], + q[v], + qlat[v], + alpha[v], + beta, + dt, + flow_length[v], + ) + + # update h, only if flow width > 0.0 + if flow_width[v] > 0.0 + crossarea = alpha[v] * pow(q[v], beta) + h[v] = crossarea / flow_width[v] + end + q_av[v] += q[v] * dt + h_av[v] += h[v] * dt + end + end + end +end + +function update!(model::SurfaceFlowLand, network, dt) + (; inwater) = model.boundary_conditions + (; alpha_term, mannings_n, slope, beta, alpha_pow, alpha, flow_width, flow_length) = + model.parameters + (; h, h_av, q_av, qlat, volume, to_river) = model.variables + (; adaptive) = model.timestepping + + @. alpha_term = pow(mannings_n / sqrt(slope), beta) + # use fixed alpha value based flow width + @. alpha = alpha_term * pow(flow_width, alpha_pow) + @. qlat = inwater / flow_length + + q_av .= 0.0 + h_av .= 0.0 + to_river .= 0.0 + + t = 0.0 + while t < dt + dt_s = adaptive ? stable_timestep(model, 0.02) : model.timestepping.dt_fixed + dt_s = check_timestepsize(dt_s, t, dt) + surfaceflow_land_update!(model, network, dt_s) + t = t + dt_s + end + q_av ./= dt + h_av ./= dt + to_river ./= dt + volume .= flow_length .* flow_width .* h + return nothing +end + +function surfaceflow_river_update!(model::SurfaceFlowRiver, network, doy, dt) + (; + graph, + order_of_subdomains, + order_subdomain, + subdomain_indices, + upstream_nodes, + reservoir_indices, + lake_indices, + ) = network.river + + (; reservoir, lake, inwater, inflow, abstraction, inflow_waterbody) = + model.boundary_conditions + + (; beta, alpha, flow_width, flow_length) = model.parameters + (; h, h_av, q, q_av, qin, qlat, volume) = model.variables + + ns = length(order_of_subdomains) + qin .= 0.0 + for k in 1:ns + threaded_foreach(eachindex(order_of_subdomains[k]); basesize = 1) do i + m = order_of_subdomains[k][i] + for (n, v) in zip(subdomain_indices[m], order_subdomain[m]) + # qin by outflow from upstream reservoir or lake location is added + qin[v] += sum_at(q, upstream_nodes[n]) + # Inflow supply/abstraction is added to qlat (divide by flow length) + # If inflow < 0, abstraction is limited + if inflow[v] < 0.0 + max_abstract = + min((inwater[v] + qin[v] + volume[v] / dt) * 0.80, -inflow[v]) + _inflow = -max_abstract / flow_length[v] + else + _inflow = inflow[v] / flow_length[v] + end + _inflow -= abstraction[v] / flow_length[v] + + q[v] = kinematic_wave( + qin[v], + q[v], + qlat[v] + _inflow, + alpha[v], + beta, + dt, + flow_length[v], + ) + + if !isnothing(reservoir) && reservoir_indices[v] != 0 + # run reservoir model and copy reservoir outflow to inflow (qin) of + # downstream river cell + i = reservoir_indices[v] + update!(reservoir, i, q[v] + inflow_waterbody[v], dt) + + downstream_nodes = outneighbors(graph, v) + n_downstream = length(downstream_nodes) + if n_downstream == 1 + j = only(downstream_nodes) + qin[j] = reservoir.variables.outflow[i] + elseif n_downstream == 0 + error( + """A reservoir without a downstream river node is not supported. + Add a downstream river node or move the reservoir to an upstream node (model schematization). + """, + ) + else + error("bifurcations not supported") + end + + elseif !isnothing(lake) && lake_indices[v] != 0 + # run lake model and copy lake outflow to inflow (qin) of downstream river + # cell + i = lake_indices[v] + update!(lake, i, q[v] + inflow_waterbody[v], doy, dt) + + downstream_nodes = outneighbors(graph, v) + n_downstream = length(downstream_nodes) + if n_downstream == 1 + j = only(downstream_nodes) + qin[j] = lake.variables.outflow[i] + elseif n_downstream == 0 + error( + """A lake without a downstream river node is not supported. + Add a downstream river node or move the lake to an upstream node (model schematization). + """, + ) + else + error("bifurcations not supported") + end + end + # update h + crossarea = alpha[v] * pow(q[v], beta) + h[v] = crossarea / flow_width[v] + volume[v] = flow_length[v] * flow_width[v] * h[v] + q_av[v] += q[v] * dt + h_av[v] += h[v] * dt + end + end + end +end + +function update!(model::SurfaceFlowRiver, network, doy, dt) + (; reservoir, lake, inwater) = model.boundary_conditions + + (; + alpha_term, + mannings_n, + slope, + beta, + alpha_pow, + alpha, + flow_width, + flow_length, + bankfull_depth, + ) = model.parameters + (; h, h_av, q_av, qlat, volume) = model.variables + (; adaptive) = model.timestepping + + @. alpha_term = pow(mannings_n / sqrt(slope), beta) + # use fixed alpha value based on 0.5 * bankfull_depth + @. alpha = alpha_term * pow(flow_width + bankfull_depth, alpha_pow) + @. qlat = inwater / flow_length + + q_av .= 0.0 + h_av .= 0.0 + # because of possible iterations set reservoir and lake inflow and total outflow at + # start to zero, the total sum of inflow and outflow at each sub time step is calculated + if !isnothing(reservoir) + reservoir.boundary_conditions.inflow .= 0.0 + reservoir.variables.totaloutflow .= 0.0 + reservoir.variables.actevap .= 0.0 + end + if !isnothing(lake) + lake.boundary_conditions.inflow .= 0.0 + lake.variables.totaloutflow .= 0.0 + lake.variables.actevap .= 0.0 + end + + t = 0.0 + while t < dt + dt_s = adaptive ? stable_timestep(model, 0.05) : model.timestepping.dt_fixed + dt_s = check_timestepsize(dt_s, t, dt) + surfaceflow_river_update!(model, network, doy, dt_s) + t = t + dt_s + end + q_av ./= dt + h_av ./= dt + volume .= flow_length .* flow_width .* h + return nothing +end + +function stable_timestep(model::S, p) where {S <: Union{SurfaceFlowLand, SurfaceFlowRiver}} + (; q) = model.variables + (; alpha, beta, flow_length) = model.parameters + (; stable_timesteps) = model.timestepping + + n = length(q) + stable_timesteps .= Inf + for i in 1:n + if q[i] > 0.0 + c = 1.0 / (alpha[i] * beta * pow(q[i], (beta - 1.0))) + stable_timesteps[i] = (flow_length[i] / c) + end + end + _stable_timesteps = filter(x -> !isinf(x), stable_timesteps) + + if !isempty(_stable_timesteps) + dt_s = quantile!(_stable_timesteps, p) + else + dt_s = 600.0 + end + + return dt_s +end + +""" +Update boundary condition lateral inflow `inwater` of a `river` flow model for a single +timestep. +""" +function update_lateral_inflow!( + model::AbstractRiverFlowModel, + external_models::NamedTuple, + river_cell_area, + land_area, + river_indices, + dt, +) + (; allocation, runoff, land, subsurface) = external_models + (; inwater) = model.boundary_conditions + (; net_runoff_river) = runoff.variables + + inwater .= ( + get_flux_to_river(subsurface)[river_indices] .+ + land.variables.to_river[river_indices] .+ + (net_runoff_river[river_indices] .* land_area[river_indices] .* 0.001) ./ dt .+ + (get_nonirrigation_returnflow(allocation) .* 0.001 .* river_cell_area) ./ dt + ) + return nothing +end + +""" +Update boundary condition lateral inflow `inwater` of kinematic wave overland flow model +`SurfaceFlowLand` for a single timestep. +""" +function update_lateral_inflow!( + model::SurfaceFlowLand, + external_models::NamedTuple, + area, + config, + dt, +) + (; soil, subsurface, allocation) = external_models + (; net_runoff) = soil.variables + (; inwater) = model.boundary_conditions + + do_drains = get(config.model, "drains", false)::Bool + if do_drains + drainflux = zeros(length(net_runoff)) + drainflux[subsurface.drain.index] = + -subsurface.drain.variables.flux ./ tosecond(basetimestep) + else + drainflux = 0.0 + end + inwater .= + (net_runoff .+ get_nonirrigation_returnflow(allocation)) .* area * 0.001 ./ dt .+ + drainflux + + return nothing +end + +""" +Update boundary condition inflow to a waterbody from land `inflow_waterbody` of a model +`AbstractRiverFlowModel` for a single timestep. +""" +function update_inflow_waterbody!( + model::AbstractRiverFlowModel, + external_models::NamedTuple, + river_indices, +) + (; land, subsurface) = external_models + (; reservoir, lake, inflow_waterbody) = model.boundary_conditions + + if !isnothing(reservoir) || !isnothing(lake) + inflow_land = get_inflow_waterbody(model, land) + inflow_subsurface = get_inflow_waterbody(model, subsurface) + + @. inflow_waterbody = inflow_land[river_indices] + inflow_subsurface[river_indices] + end + return nothing +end + +# For the river kinematic wave, the variable `to_river` can be excluded, because this part +# is added to the river kinematic wave. +get_inflow_waterbody(::SurfaceFlowRiver, model::SurfaceFlowLand) = model.variables.q_av +get_inflow_waterbody(::SurfaceFlowRiver, model::LateralSSF) = + model.variables.ssf ./ tosecond(basetimestep) + +# Exclude subsurface flow for other groundwater components than `LateralSSF`. +get_inflow_waterbody(::AbstractRiverFlowModel, model::GroundwaterFlow) = + model.flow.connectivity.ncell .* 0.0 +get_inflow_waterbody(::SurfaceFlowRiver, model) = model.variables.to_river .* 0.0 \ No newline at end of file diff --git a/src/flow.jl b/src/routing/surface_local_inertial.jl similarity index 62% rename from src/flow.jl rename to src/routing/surface_local_inertial.jl index c442cc6e7..e084b58b5 100644 --- a/src/flow.jl +++ b/src/routing/surface_local_inertial.jl @@ -1,753 +1,3 @@ -abstract type AbstractRiverFlowModel end - -@get_units @grid_loc @with_kw struct FlowVariables{T} - q::Vector{T} | "m3 s-1" # Discharge [m³ s⁻¹] - qlat::Vector{T} | "m2 s-1" # Lateral inflow per unit length [m² s⁻¹] - qin::Vector{T} | "m3 s-1" # Inflow from upstream cells [m³ s⁻¹] - q_av::Vector{T} | "m3 s-1" # Average discharge [m³ s⁻¹] - volume::Vector{T} | "m3" # Kinematic wave volume [m³] (based on water depth h) - h::Vector{T} | "m" # Water depth [m] - h_av::Vector{T} | "m" # Average water depth [m] -end - -@with_kw struct TimeStepping{T} - stable_timesteps::Vector{T} = Float[] - dt_fixed::T = 0.0 - adaptive::Bool = true - cfl::T = 0.70 -end - -function init_kinematic_wave_timestepping(config, n; domain, dt_fixed) - adaptive = get(config.model, "kin_wave_iteration", false)::Bool - @info "Kinematic wave approach is used for $domain flow, adaptive timestepping = $adaptive." - - if adaptive - stable_timesteps = zeros(n) - timestepping = TimeStepping(; stable_timesteps, adaptive) - else - dt_fixed = get(config.model, "kw_$(domain)_tstep", dt_fixed) - @info "Using a fixed sub-timestep (seconds) $dt_fixed for kinematic wave $domain flow." - timestepping = TimeStepping(; dt_fixed, adaptive) - end - return timestepping -end - -function check_timestepsize(timestepsize, currenttime, endtime) - if currenttime + timestepsize > endtime - timestepsize = endtime - currenttime - end - return timestepsize -end - -function FlowVariables(n) - variables = FlowVariables(; - q = zeros(Float, n), - qlat = zeros(Float, n), - qin = zeros(Float, n), - q_av = zeros(Float, n), - volume = zeros(Float, n), - h = zeros(Float, n), - h_av = zeros(Float, n), - ) - return variables -end - -@get_units @grid_loc @with_kw struct ManningFlowParameters{T} - beta::T # constant in Manning's equation [-] - slope::Vector{T} | "m m-1" # Slope [m m⁻¹] - mannings_n::Vector{T} | "s m-1/3" # Manning's roughness [s m⁻⅓] - flow_length::Vector{T} | "m" # Flow length [m] - flow_width::Vector{T} | "m" # Flow width [m] - alpha_pow::T # Used in the power part of alpha [-] - alpha_term::Vector{T} | "-" # Term used in computation of alpha [-] - alpha::Vector{T} | "s3/5 m1/5" # Constant in momentum equation A = alpha*Q^beta, based on Manning's equation -end - -function ManningFlowParameters(slope, mannings_n, flow_length, flow_width) - n = length(slope) - parameters = ManningFlowParameters(; - beta = Float(0.6), - slope, - mannings_n, - flow_length, - flow_width, - alpha_pow = Float((2.0 / 3.0) * 0.6), - alpha_term = fill(mv, n), - alpha = fill(mv, n), - ) - return parameters -end - -@get_units @grid_loc @with_kw struct RiverFlowParameters{T} - flow::ManningFlowParameters{T} - bankfull_depth::Vector{T} | "m" # Bankfull water level [m] -end - -function Base.getproperty(v::RiverFlowParameters, s::Symbol) - if s === :bankfull_depth - getfield(v, s) - elseif s === :flow - getfield(v, :flow) - else - getfield(getfield(v, :flow), s) - end -end - -function RiverFlowParameters(dataset, config, indices, river_length, river_width) - mannings_n = ncread( - dataset, - config, - "lateral.river.mannings_n"; - sel = indices, - defaults = 0.036, - type = Float, - ) - bankfull_depth = ncread( - dataset, - config, - "lateral.river.bankfull_depth"; - alias = "lateral.river.h_bankfull", - sel = indices, - defaults = 1.0, - type = Float, - ) - if haskey(config.input.lateral.river, "h_bankfull") - @warn string( - "The `h_bankfull` key in `[input.lateral.river]` is now called ", - "`bankfull_depth`. Please update your TOML file.", - ) - end - slope = ncread( - dataset, - config, - "lateral.river.slope"; - optional = false, - sel = indices, - type = Float, - ) - clamp!(slope, 0.00001, Inf) - - flow_parameter_set = ManningFlowParameters(slope, mannings_n, river_length, river_width) - parameters = - RiverFlowParameters(; flow = flow_parameter_set, bankfull_depth = bankfull_depth) - return parameters -end - -@get_units @grid_loc @with_kw struct RiverFlowBC{T, R, L} - inwater::Vector{T} | "m3 s-1" # Lateral inflow [m³ s⁻¹] - inflow::Vector{T} | "m3 s-1" # External inflow (abstraction/supply/demand) [m³ s⁻¹] - inflow_waterbody::Vector{T} | "m3 s-1" # inflow waterbody (lake or reservoir model) from land part [m³ s⁻¹] - abstraction::Vector{T} | "m3 s-1" # Abstraction (computed as part of water demand and allocation) [m³ s⁻¹] - reservoir::R # Reservoir model struct of arrays - lake::L # Lake model struct of arrays -end - -function RiverFlowBC(n, reservoir, lake) - bc = RiverFlowBC(; - inwater = zeros(Float, n), - inflow = zeros(Float, n), - inflow_waterbody = zeros(Float, n), - abstraction = zeros(Float, n), - reservoir = reservoir, - lake = lake, - ) - return bc -end - -@with_kw struct SurfaceFlowRiver{T, R, L, A} <: AbstractRiverFlowModel - timestepping::TimeStepping{T} - boundary_conditions::RiverFlowBC{T, R, L} - parameters::RiverFlowParameters{T} - variables::FlowVariables{T} - allocation::A # Water allocation -end - -function SurfaceFlowRiver( - dataset, - config, - indices; - river_length, - river_width, - reservoir, - lake, -) - n = length(indices) - - timestepping = - init_kinematic_wave_timestepping(config, n; domain = "river", dt_fixed = 900.0) - - do_water_demand = haskey(config.model, "water_demand") - allocation = do_water_demand ? AllocationRiver(n) : NoAllocationRiver{Float}() - - variables = FlowVariables(n) - parameters = RiverFlowParameters(dataset, config, indices, river_length, river_width) - boundary_conditions = RiverFlowBC(n, reservoir, lake) - - sf_river = SurfaceFlowRiver(; - timestepping, - boundary_conditions, - parameters, - variables, - allocation, - ) - - return sf_river -end - -@get_units @grid_loc @with_kw struct LandFlowVariables{T} - flow::FlowVariables{T} - to_river::Vector{T} | "m3 s-1" # Part of overland flow [m³ s⁻¹] that flows to the river -end - -function Base.getproperty(v::LandFlowVariables, s::Symbol) - if s === :to_river - getfield(v, s) - elseif s === :flow - getfield(v, :flow) - else - getfield(getfield(v, :flow), s) - end -end - -@get_units @grid_loc @with_kw struct LandFlowBC{T} - inwater::Vector{T} | "m3 s-1" # Lateral inflow [m³ s⁻¹] -end - -@with_kw struct SurfaceFlowLand{T} - timestepping::TimeStepping{T} - boundary_conditions::LandFlowBC{T} - parameters::ManningFlowParameters{T} - variables::LandFlowVariables{T} -end - -function SurfaceFlowLand(dataset, config, indices; slope, flow_length, flow_width) - mannings_n = ncread( - dataset, - config, - "lateral.land.mannings_n"; - sel = indices, - defaults = 0.072, - type = Float, - ) - - n = length(indices) - timestepping = - init_kinematic_wave_timestepping(config, n; domain = "land", dt_fixed = 3600.0) - - variables = LandFlowVariables(; flow = FlowVariables(n), to_river = zeros(Float, n)) - parameters = ManningFlowParameters(slope, mannings_n, flow_length, flow_width) - boundary_conditions = LandFlowBC(; inwater = zeros(Float, n)) - sf_land = SurfaceFlowLand(; timestepping, boundary_conditions, variables, parameters) - - return sf_land -end - -function surfaceflow_land_update!(model::SurfaceFlowLand, network, dt) - (; - order_of_subdomains, - order_subdomain, - subdomain_indices, - upstream_nodes, - frac_to_river, - ) = network - - (; beta, alpha, flow_width, flow_length) = model.parameters - (; h, h_av, q, q_av, qin, qlat, to_river) = model.variables - - ns = length(order_of_subdomains) - qin .= 0.0 - for k in 1:ns - threaded_foreach(eachindex(order_of_subdomains[k]); basesize = 1) do i - m = order_of_subdomains[k][i] - for (n, v) in zip(subdomain_indices[m], order_subdomain[m]) - # for a river cell without a reservoir or lake part of the upstream - # surface flow goes to the river (frac_to_river) and part goes to the - # surface flow reservoir (1.0 - frac_to_river), upstream nodes with a - # reservoir or lake are excluded - to_river[v] += - sum_at( - i -> q[i] * frac_to_river[i], - upstream_nodes[n], - eltype(to_river), - ) * dt - if flow_width[v] > 0.0 - qin[v] = sum_at( - i -> q[i] * (1.0 - frac_to_river[i]), - upstream_nodes[n], - eltype(q), - ) - end - - q[v] = kinematic_wave( - qin[v], - q[v], - qlat[v], - alpha[v], - beta, - dt, - flow_length[v], - ) - - # update h, only if flow width > 0.0 - if flow_width[v] > 0.0 - crossarea = alpha[v] * pow(q[v], beta) - h[v] = crossarea / flow_width[v] - end - q_av[v] += q[v] * dt - h_av[v] += h[v] * dt - end - end - end -end - -function update!(model::SurfaceFlowLand, network, dt) - (; inwater) = model.boundary_conditions - (; alpha_term, mannings_n, slope, beta, alpha_pow, alpha, flow_width, flow_length) = - model.parameters - (; h, h_av, q_av, qlat, volume, to_river) = model.variables - (; adaptive) = model.timestepping - - @. alpha_term = pow(mannings_n / sqrt(slope), beta) - # use fixed alpha value based flow width - @. alpha = alpha_term * pow(flow_width, alpha_pow) - @. qlat = inwater / flow_length - - q_av .= 0.0 - h_av .= 0.0 - to_river .= 0.0 - - t = 0.0 - while t < dt - dt_s = adaptive ? stable_timestep(model, 0.02) : model.timestepping.dt_fixed - dt_s = check_timestepsize(dt_s, t, dt) - surfaceflow_land_update!(model, network, dt_s) - t = t + dt_s - end - q_av ./= dt - h_av ./= dt - to_river ./= dt - volume .= flow_length .* flow_width .* h - return nothing -end - -function surfaceflow_river_update!(model::SurfaceFlowRiver, network, doy, dt) - (; - graph, - order_of_subdomains, - order_subdomain, - subdomain_indices, - upstream_nodes, - reservoir_indices, - lake_indices, - ) = network.river - - (; reservoir, lake, inwater, inflow, abstraction, inflow_waterbody) = - model.boundary_conditions - - (; beta, alpha, flow_width, flow_length) = model.parameters - (; h, h_av, q, q_av, qin, qlat, volume) = model.variables - - ns = length(order_of_subdomains) - qin .= 0.0 - for k in 1:ns - threaded_foreach(eachindex(order_of_subdomains[k]); basesize = 1) do i - m = order_of_subdomains[k][i] - for (n, v) in zip(subdomain_indices[m], order_subdomain[m]) - # qin by outflow from upstream reservoir or lake location is added - qin[v] += sum_at(q, upstream_nodes[n]) - # Inflow supply/abstraction is added to qlat (divide by flow length) - # If inflow < 0, abstraction is limited - if inflow[v] < 0.0 - max_abstract = - min((inwater[v] + qin[v] + volume[v] / dt) * 0.80, -inflow[v]) - _inflow = -max_abstract / flow_length[v] - else - _inflow = inflow[v] / flow_length[v] - end - _inflow -= abstraction[v] / flow_length[v] - - q[v] = kinematic_wave( - qin[v], - q[v], - qlat[v] + _inflow, - alpha[v], - beta, - dt, - flow_length[v], - ) - - if !isnothing(reservoir) && reservoir_indices[v] != 0 - # run reservoir model and copy reservoir outflow to inflow (qin) of - # downstream river cell - i = reservoir_indices[v] - update!(reservoir, i, q[v] + inflow_waterbody[v], dt) - - downstream_nodes = outneighbors(graph, v) - n_downstream = length(downstream_nodes) - if n_downstream == 1 - j = only(downstream_nodes) - qin[j] = reservoir.variables.outflow[i] - elseif n_downstream == 0 - error( - """A reservoir without a downstream river node is not supported. - Add a downstream river node or move the reservoir to an upstream node (model schematization). - """, - ) - else - error("bifurcations not supported") - end - - elseif !isnothing(lake) && lake_indices[v] != 0 - # run lake model and copy lake outflow to inflow (qin) of downstream river - # cell - i = lake_indices[v] - update!(lake, i, q[v] + inflow_waterbody[v], doy, dt) - - downstream_nodes = outneighbors(graph, v) - n_downstream = length(downstream_nodes) - if n_downstream == 1 - j = only(downstream_nodes) - qin[j] = lake.variables.outflow[i] - elseif n_downstream == 0 - error( - """A lake without a downstream river node is not supported. - Add a downstream river node or move the lake to an upstream node (model schematization). - """, - ) - else - error("bifurcations not supported") - end - end - # update h - crossarea = alpha[v] * pow(q[v], beta) - h[v] = crossarea / flow_width[v] - volume[v] = flow_length[v] * flow_width[v] * h[v] - q_av[v] += q[v] * dt - h_av[v] += h[v] * dt - end - end - end -end - -function update!(model::SurfaceFlowRiver, network, doy, dt) - (; reservoir, lake, inwater) = model.boundary_conditions - - (; - alpha_term, - mannings_n, - slope, - beta, - alpha_pow, - alpha, - flow_width, - flow_length, - bankfull_depth, - ) = model.parameters - (; h, h_av, q_av, qlat, volume) = model.variables - (; adaptive) = model.timestepping - - @. alpha_term = pow(mannings_n / sqrt(slope), beta) - # use fixed alpha value based on 0.5 * bankfull_depth - @. alpha = alpha_term * pow(flow_width + bankfull_depth, alpha_pow) - @. qlat = inwater / flow_length - - q_av .= 0.0 - h_av .= 0.0 - # because of possible iterations set reservoir and lake inflow and total outflow at - # start to zero, the total sum of inflow and outflow at each sub time step is calculated - if !isnothing(reservoir) - reservoir.boundary_conditions.inflow .= 0.0 - reservoir.variables.totaloutflow .= 0.0 - reservoir.variables.actevap .= 0.0 - end - if !isnothing(lake) - lake.boundary_conditions.inflow .= 0.0 - lake.variables.totaloutflow .= 0.0 - lake.variables.actevap .= 0.0 - end - - t = 0.0 - while t < dt - dt_s = adaptive ? stable_timestep(model, 0.05) : model.timestepping.dt_fixed - dt_s = check_timestepsize(dt_s, t, dt) - surfaceflow_river_update!(model, network, doy, dt_s) - t = t + dt_s - end - q_av ./= dt - h_av ./= dt - volume .= flow_length .* flow_width .* h - return nothing -end - -function stable_timestep(model::S, p) where {S <: Union{SurfaceFlowLand, SurfaceFlowRiver}} - (; q) = model.variables - (; alpha, beta, flow_length) = model.parameters - (; stable_timesteps) = model.timestepping - - n = length(q) - stable_timesteps .= Inf - for i in 1:n - if q[i] > 0.0 - c = 1.0 / (alpha[i] * beta * pow(q[i], (beta - 1.0))) - stable_timesteps[i] = (flow_length[i] / c) - end - end - _stable_timesteps = filter(x -> !isinf(x), stable_timesteps) - - if !isempty(_stable_timesteps) - dt_s = quantile!(_stable_timesteps, p) - else - dt_s = 600.0 - end - - return dt_s -end - -@get_units @grid_loc struct KhExponential{T} - # Horizontal hydraulic conductivity at soil surface [m d⁻¹] - kh_0::Vector{T} | "m d-1" - # A scaling parameter [m⁻¹] (controls exponential decline of kh_0) - f::Vector{T} | "m-1" -end - -@get_units @grid_loc struct KhExponentialConstant{T} - # Exponential horizontal hydraulic conductivity profile type - exponential::KhExponential - # Depth [m] from soil surface for which exponential decline of kv_0 is valid - z_exp::Vector{T} | "m" -end - -@get_units @grid_loc struct KhLayered{T} - # Horizontal hydraulic conductivity [m d⁻¹] - kh::Vector{T} | "m d-1" -end - -abstract type SubsurfaceFlow end - -@get_units @grid_loc @with_kw struct LateralSsfParameters{T, Kh} - kh_profile::Kh # Horizontal hydraulic conductivity profile type [-] - khfrac::Vector{T} | "-" # A muliplication factor applied to vertical hydraulic conductivity `kv` [-] - soilthickness::Vector{T} | "m" # Soil thickness [m] - theta_s::Vector{T} | "-" # Saturated water content (porosity) [-] - theta_r::Vector{T} | "-" # Residual water content [-] - dt::T # model time step [d] - slope::Vector{T} | "m m-1" # Slope [m m⁻¹] - flow_length::Vector{T} | "m" # Flow length [m] - flow_width::Vector{T} | "m" # Flow width [m] -end - -function LateralSsfParameters( - dataset, - config, - indices; - soil, - slope, - flow_length, - flow_width, - dt, -) - khfrac = ncread( - dataset, - config, - "lateral.subsurface.ksathorfrac"; - sel = indices, - defaults = 1.0, - type = Float, - ) - n_cells = length(khfrac) - - (; theta_s, theta_r, soilthickness) = soil - soilthickness = soilthickness .* 0.001 - - kh_profile_type = get(config.input.vertical, "ksat_profile", "exponential")::String - if kh_profile_type == "exponential" - (; kv_0, f) = soil.kv_profile - kh_0 = khfrac .* kv_0 .* 0.001 .* dt - kh_profile = KhExponential(kh_0, f .* 1000.0) - elseif kh_profile_type == "exponential_constant" - (; z_exp) = soil.kv_profile - (; kv_0, f) = soil.kv_profile.exponential - kh_0 = khfrac .* kv_0 .* 0.001 .* dt - exp_profile = KhExponential(kh_0, f .* 1000.0) - kh_profile = KhExponentialConstant(exp_profile, z_exp .* 0.001) - elseif kh_profile_type == "layered" || kh_profile_type == "layered_exponential" - kh_profile = KhLayered(fill(mv, n_cells)) - end - parameters = LateralSsfParameters( - kh_profile, - khfrac, - soilthickness, - theta_s, - theta_r, - dt, - slope, - flow_length, - flow_width, - ) - return parameters -end - -@get_units @grid_loc @with_kw struct LateralSsfVariables{T} - zi::Vector{T} | "m" # Pseudo-water table depth [m] (top of the saturated zone) - exfiltwater::Vector{T} | "m dt-1" # Exfiltration [m Δt⁻¹] (groundwater above surface level, saturated excess conditions) - recharge::Vector{T} | "m2 dt-1" # Net recharge to saturated store [m² Δt⁻¹] - ssf::Vector{T} | "m3 d-1" # Subsurface flow [m³ d⁻¹] - ssfin::Vector{T} | "m3 d-1" # Inflow from upstream cells [m³ d⁻¹] - ssfmax::Vector{T} | "m2 d-1" # Maximum subsurface flow [m² d⁻¹] - to_river::Vector{T} | "m3 d-1" # Part of subsurface flow [m³ d⁻¹] that flows to the river - volume::Vector{T} | "m3" # Subsurface volume [m³] -end - -function LateralSsfVariables(ssf, zi, xl, yl) - n = length(zi) - volume = @. (ssf.theta_s - ssf.theta_r) * (ssf.soilthickness - zi) * (xl * yl) - variables = LateralSsfVariables(; - zi, - exfiltwater = fill(mv, n), - recharge = fill(mv, n), - ssf = fill(mv, n), - ssfin = fill(mv, n), - ssfmax = fill(mv, n), - to_river = zeros(n), - volume, - ) - return variables -end - -@get_units @grid_loc @with_kw struct LateralSsfBC{T} - recharge::Vector{T} | "m2 dt-1" # Net recharge to saturated store [m² Δt⁻¹] -end - -@with_kw struct LateralSSF{T, Kh} <: SubsurfaceFlow - boundary_conditions::LateralSsfBC{T} - parameters::LateralSsfParameters{T, Kh} - variables::LateralSsfVariables{T} -end - -function LateralSSF( - dataset, - config, - indices; - soil, - slope, - flow_length, - flow_width, - x_length, - y_length, - dt, -) - parameters = LateralSsfParameters( - dataset, - config, - indices; - soil = soil.parameters, - slope, - flow_length, - flow_width, - dt, - ) - zi = 0.001 * soil.variables.zi - variables = LateralSsfVariables(parameters, zi, x_length, y_length) - boundary_conditions = LateralSsfBC(; recharge = fill(mv, length(zi))) - ssf = LateralSSF(; boundary_conditions, parameters, variables) - return ssf -end - -function update!(model::LateralSSF, network) - (; - order_of_subdomains, - order_subdomain, - subdomain_indices, - upstream_nodes, - area, - frac_to_river, - ) = network - - (; recharge) = model.boundary_conditions - (; ssfin, ssf, ssfmax, to_river, zi, exfiltwater, volume) = model.variables - (; slope, theta_s, theta_r, soilthickness, flow_length, flow_width, dt, kh_profile) = - model.parameters - - ns = length(order_of_subdomains) - for k in 1:ns - threaded_foreach(eachindex(order_of_subdomains[k]); basesize = 1) do i - m = order_of_subdomains[k][i] - for (n, v) in zip(subdomain_indices[m], order_subdomain[m]) - # for a river cell without a reservoir or lake part of the upstream - # subsurface flow goes to the river (frac_to_river) and part goes to the - # subsurface flow reservoir (1.0 - frac_to_river) upstream nodes with a - # reservoir or lake are excluded - ssfin[v] = sum_at( - i -> ssf[i] * (1.0 - frac_to_river[i]), - upstream_nodes[n], - eltype(ssfin), - ) - to_river[v] = sum_at( - i -> ssf[i] * frac_to_river[i], - upstream_nodes[n], - eltype(to_river), - ) - ssf[v], zi[v], exfiltwater[v] = kinematic_wave_ssf( - ssfin[v], - ssf[v], - zi[v], - recharge[v], - slope[v], - theta_s[v] - theta_r[v], - soilthickness[v], - dt, - flow_length[v], - flow_width[v], - ssfmax[v], - kh_profile, - v, - ) - volume[v] = (theta_s[v] - theta_r[v]) * (soilthickness[v] - zi[v]) * area[v] - end - end - end - return nothing -end - -@get_units@grid_loc @with_kw struct GroundwaterExchangeVariables{T} - exfiltwater::Vector{T} | "m dt-1" # Exfiltration [m Δt⁻¹] (groundwater above surface level, saturated excess conditions) - zi::Vector{T} | "m" # Pseudo-water table depth [m] (top of the saturated zone) - to_river::Vector{T} | "m3 d-1" # Part of subsurface flow [m³ d⁻¹] that flows to the river - ssf::Vector{T} | "m3 d-1" # Subsurface flow [m³ d⁻¹] -end - -function GroundwaterExchangeVariables(n) - variables = GroundwaterExchangeVariables{Float}(; - exfiltwater = fill(mv, n), - zi = fill(mv, n), - to_river = fill(mv, n), - ssf = zeros(n), - ) - return variables -end - -@with_kw struct GroundwaterExchangeParameters{T} - dt::T # model time step [d] -end - -@with_kw struct GroundwaterExchange{T} <: SubsurfaceFlow - parameters::GroundwaterExchangeParameters{T} - variables::GroundwaterExchangeVariables{T} -end - -function GroundwaterExchange(n, dt) - parameters = GroundwaterExchangeParameters{Float}(; dt = dt / basetimestep) - variables = GroundwaterExchangeVariables(n) - ssf = GroundwaterExchange{Float}(; parameters, variables) - return ssf -end - -get_water_depth(subsurface::SubsurfaceFlow) = subsurface.variables.zi -get_exfiltwater(subsurface::SubsurfaceFlow) = subsurface.variables.exfiltwater - -get_flux_to_river(subsurface::SubsurfaceFlow) = - subsurface.variables.to_river ./ tosecond(basetimestep) # [m³ s⁻¹] - @get_units @grid_loc @with_kw struct ShallowWaterRiverParameters{T} n::Int # number of cells [-] ne::Int # number of edges/links [-] @@ -2074,82 +1324,6 @@ function FloodPlain( return floodplain end -""" -Update boundary condition lateral inflow `inwater` of a `river` flow model for a single -timestep. -""" -function update_lateral_inflow!( - model::AbstractRiverFlowModel, - external_models::NamedTuple, - river_cell_area, - land_area, - river_indices, - dt, -) - (; allocation, runoff, land, subsurface) = external_models - (; inwater) = model.boundary_conditions - (; net_runoff_river) = runoff.variables - - inwater .= ( - get_flux_to_river(subsurface)[river_indices] .+ - land.variables.to_river[river_indices] .+ - (net_runoff_river[river_indices] .* land_area[river_indices] .* 0.001) ./ dt .+ - (get_nonirrigation_returnflow(allocation) .* 0.001 .* river_cell_area) ./ dt - ) - return nothing -end - -""" -Update boundary condition lateral inflow `inwater` of kinematic wave overland flow model -`SurfaceFlowLand` for a single timestep. -""" -function update_lateral_inflow!( - model::SurfaceFlowLand, - external_models::NamedTuple, - area, - config, - dt, -) - (; soil, subsurface, allocation) = external_models - (; net_runoff) = soil.variables - (; inwater) = model.boundary_conditions - - do_drains = get(config.model, "drains", false)::Bool - if do_drains - drainflux = zeros(length(net_runoff)) - drainflux[subsurface.drain.index] = - -subsurface.drain.variables.flux ./ tosecond(basetimestep) - else - drainflux = 0.0 - end - inwater .= - (net_runoff .+ get_nonirrigation_returnflow(allocation)) .* area * 0.001 ./ dt .+ - drainflux - - return nothing -end - -""" -Update boundary condition inflow to a waterbody from land `inflow_waterbody` of a model -`AbstractRiverFlowModel` for a single timestep. -""" -function update_inflow_waterbody!( - model::AbstractRiverFlowModel, - external_models::NamedTuple, - river_indices, -) - (; land, subsurface) = external_models - (; reservoir, lake, inflow_waterbody) = model.boundary_conditions - - if !isnothing(reservoir) || !isnothing(lake) - inflow_land = get_inflow_waterbody(model, land) - inflow_subsurface = get_inflow_waterbody(model, subsurface) - - @. inflow_waterbody = inflow_land[river_indices] + inflow_subsurface[river_indices] - end - return nothing -end - """ Update boundary conditions `runoff` and inflow to a waterbody from land `inflow_waterbody` for overland flow model `ShallowWaterLand` for a single timestep. @@ -2180,87 +1354,9 @@ function update_boundary_conditions!( return nothing end -# For the river kinematic wave, the variable `to_river` can be excluded, because this part -# is added to the river kinematic wave. -get_inflow_waterbody(::SurfaceFlowRiver, model::SurfaceFlowLand) = model.variables.q_av -get_inflow_waterbody(::SurfaceFlowRiver, model::LateralSSF) = - model.variables.ssf ./ tosecond(basetimestep) - -# Exclude subsurface flow for other groundwater components than `LateralSSF`. -get_inflow_waterbody(::Union{SurfaceFlowRiver, ShallowWaterRiver}, model::GroundwaterFlow) = - model.flow.connectivity.ncell .* 0.0 -get_inflow_waterbody(::SurfaceFlowRiver, model) = model.variables.to_river .* 0.0 - # For local inertial river routing, `to_river` is included, as water body cells are excluded # (boundary condition). get_inflow_waterbody(::ShallowWaterRiver, model::SurfaceFlowLand) = model.variables.q_av .+ model.variables.to_river get_inflow_waterbody(::ShallowWaterRiver, model::LateralSSF) = - (model.variables.ssf .+ model.variables.to_river) ./ tosecond(basetimestep) - -""" - surface_routing!(model) - -Run surface routing (land and river) for a single timestep. Kinematic wave for overland flow -and kinematic wave or local inertial model for river flow. -""" -function surface_routing!(model) - (; vertical, lateral, network, config, clock) = model - (; soil, runoff, allocation) = vertical - (; land, river, subsurface) = lateral - - dt = tosecond(clock.dt) - # update lateral inflow for kinematic wave overland flow - update_lateral_inflow!( - land, - (; soil, allocation, subsurface), - network.land.area, - config, - dt, - ) - # run kinematic wave overland flow - update!(land, network.land, dt) - - # update lateral inflow river flow - update_lateral_inflow!( - river, - (; allocation = river.allocation, runoff, land, subsurface), - network.river.cell_area, - network.land.area, - network.river.land_indices, - dt, - ) - update_inflow_waterbody!(river, (; land, subsurface), network.river.land_indices) - update!(river, network, julian_day(clock.time - clock.dt), dt) - return nothing -end - -""" - surface_routing!( - model::Model{N,L,V,R,W,T} - ) where {N,L<:NamedTuple{<:Any,<:Tuple{Any,ShallowWaterLand,ShallowWaterRiver}},V,R,W,T} - -Run surface routing (land and river) for a model type that contains the lateral components -`ShallowWaterLand` and `ShallowWaterRiver` for a single timestep. -""" -function surface_routing!( - model::Model{N, L, V, R, W, T}, -) where { - N, - L <: NamedTuple{<:Any, <:Tuple{Any, ShallowWaterLand, ShallowWaterRiver}}, - V, - R, - W, - T, -} - (; lateral, vertical, network, clock) = model - (; land, river, subsurface) = lateral - (; soil, runoff) = vertical - - dt = tosecond(clock.dt) - update_boundary_conditions!(land, (; river, subsurface, soil, runoff), network, dt) - - update!(land, river, network, julian_day(clock.time - clock.dt), dt) - - return nothing -end + (model.variables.ssf .+ model.variables.to_river) ./ tosecond(basetimestep) \ No newline at end of file diff --git a/src/routing/surface_routing.jl b/src/routing/surface_routing.jl new file mode 100644 index 000000000..325cabb4d --- /dev/null +++ b/src/routing/surface_routing.jl @@ -0,0 +1,66 @@ +""" + surface_routing!(model) + +Run surface routing (land and river) for a single timestep. Kinematic wave for overland flow +and kinematic wave or local inertial model for river flow. +""" +function surface_routing!(model) + (; vertical, lateral, network, config, clock) = model + (; soil, runoff, allocation) = vertical + (; land, river, subsurface) = lateral + + dt = tosecond(clock.dt) + # update lateral inflow for kinematic wave overland flow + update_lateral_inflow!( + land, + (; soil, allocation, subsurface), + network.land.area, + config, + dt, + ) + # run kinematic wave overland flow + update!(land, network.land, dt) + + # update lateral inflow river flow + update_lateral_inflow!( + river, + (; allocation = river.allocation, runoff, land, subsurface), + network.river.cell_area, + network.land.area, + network.river.land_indices, + dt, + ) + update_inflow_waterbody!(river, (; land, subsurface), network.river.land_indices) + update!(river, network, julian_day(clock.time - clock.dt), dt) + return nothing +end + +""" + surface_routing!( + model::Model{N,L,V,R,W,T} + ) where {N,L<:NamedTuple{<:Any,<:Tuple{Any,ShallowWaterLand,ShallowWaterRiver}},V,R,W,T} + +Run surface routing (land and river) for a model type that contains the lateral components +`ShallowWaterLand` and `ShallowWaterRiver` for a single timestep. +""" +function surface_routing!( + model::Model{N, L, V, R, W, T}, +) where { + N, + L <: NamedTuple{<:Any, <:Tuple{Any, ShallowWaterLand, ShallowWaterRiver}}, + V, + R, + W, + T, +} + (; lateral, vertical, network, clock) = model + (; land, river, subsurface) = lateral + (; soil, runoff) = vertical + + dt = tosecond(clock.dt) + update_boundary_conditions!(land, (; river, subsurface, soil, runoff), network, dt) + + update!(land, river, network, julian_day(clock.time - clock.dt), dt) + + return nothing +end \ No newline at end of file diff --git a/src/routing/timestepping.jl b/src/routing/timestepping.jl new file mode 100644 index 000000000..23888d42f --- /dev/null +++ b/src/routing/timestepping.jl @@ -0,0 +1,13 @@ +@with_kw struct TimeStepping{T} + stable_timesteps::Vector{T} = Float[] + dt_fixed::T = 0.0 + adaptive::Bool = true + cfl::T = 0.70 +end + +function check_timestepsize(timestepsize, currenttime, endtime) + if currenttime + timestepsize > endtime + timestepsize = endtime - currenttime + end + return timestepsize +end \ No newline at end of file