From 40dd0fac864516c35f4def7bf09322b7c8164524 Mon Sep 17 00:00:00 2001 From: Willem van Verseveld Date: Thu, 14 Nov 2024 21:28:09 +0100 Subject: [PATCH] Refactor groundwater flow structs --- src/demand/water_demand.jl | 2 +- src/flow.jl | 2 +- src/groundwater/aquifer.jl | 183 +++++++++++++++------- src/groundwater/boundary_conditions.jl | 155 +++++++++++++++---- src/sbm_gwf_model.jl | 122 +++------------ src/states.jl | 2 +- test/groundwater.jl | 206 +++++++++++++++---------- test/io.jl | 2 +- test/run_sbm_gwf.jl | 36 ++--- test/sbm_gwf_config.toml | 12 +- test/sbm_gwf_piave_demand_config.toml | 4 +- 11 files changed, 428 insertions(+), 298 deletions(-) diff --git a/src/demand/water_demand.jl b/src/demand/water_demand.jl index f7737ceab..74bf74880 100644 --- a/src/demand/water_demand.jl +++ b/src/demand/water_demand.jl @@ -784,7 +784,7 @@ return_flow(non_irri::NoNonIrrigationDemand, nonirri_demand_gross, nonirri_alloc # wrapper methods groundwater_volume(model::LateralSSF) = model.variables.volume -groundwater_volume(model) = model.flow.aquifer.volume +groundwater_volume(model) = model.flow.aquifer.variables.volume """ update_water_allocation!(land_allocation, demand::Demand, lateral, network, dt) diff --git a/src/flow.jl b/src/flow.jl index 45309eb85..22be1813f 100644 --- a/src/flow.jl +++ b/src/flow.jl @@ -2026,7 +2026,7 @@ function set_land_inwater!( do_water_demand = haskey(config.model, "water_demand") if do_drains drainflux[lateral.subsurface.drain.index] = - -lateral.subsurface.drain.flux ./ tosecond(basetimestep) + -lateral.subsurface.drain.variables.flux ./ tosecond(basetimestep) end if do_water_demand @. inwater = diff --git a/src/groundwater/aquifer.jl b/src/groundwater/aquifer.jl index 3a8619bf6..6def5f358 100644 --- a/src/groundwater/aquifer.jl +++ b/src/groundwater/aquifer.jl @@ -84,18 +84,26 @@ NOTA BENE: **specific** storage is per m of aquifer (conf. specific weight). **Storativity** or (**storage coefficient**) is for the entire aquifer (conf. transmissivity). """ -@get_units @grid_loc struct ConfinedAquifer{T} <: Aquifer - head::Vector{T} | "m" # hydraulic head [m] +@get_units @grid_loc @with_kw struct ConfinedAquiferParameters{T} k::Vector{T} | "m d-1" # horizontal conductivity [m d⁻¹] top::Vector{T} | "m" # top of groundwater layer [m] bottom::Vector{T} | "m" # bottom of groundwater layer area::Vector{T} | "m2" # area of cell specific_storage::Vector{T} | "m m-1 m-1" # [m m⁻¹ m⁻¹] storativity::Vector{T} | "m m-1" # [m m⁻¹] +end + +@get_units @grid_loc @with_kw struct ConfinedAquiferVariables{T} + head::Vector{T} | "m" # hydraulic head [m] conductance::Vector{T} | "m2 d-1" # Confined aquifer conductance is constant volume::Vector{T} | "m3" # total volume of water that can be released end +@with_kw struct ConfinedAquifer{T} <: Aquifer + parameters::ConfinedAquiferParameters{T} + variables::ConfinedAquiferVariables{T} +end + """ UnconfinedAquifer{T} <: Aquifer @@ -107,22 +115,57 @@ aquifer will yield when all water drains and the pore volume is filled by air instead. Specific yield will vary roughly between 0.05 (clay) and 0.45 (peat) (Johnson, 1967). """ -@get_units @grid_loc struct UnconfinedAquifer{T} <: Aquifer - head::Vector{T} | "m" # hydraulic head [m] +@get_units @grid_loc @with_kw struct UnconfinedAquiferParameters{T} k::Vector{T} | "m d-1" # reference horizontal conductivity [m d⁻¹] top::Vector{T} | "m" # top of groundwater layer [m] bottom::Vector{T} | "m" # bottom of groundwater layer area::Vector{T} | "m2" specific_yield::Vector{T} | "m m-1" # [m m⁻¹] - conductance::Vector{T} | "m2 d-1" # - volume::Vector{T} | "m3" # total volume of water that can be released f::Vector{T} | "-" # factor controlling the reduction of reference horizontal conductivity # Unconfined aquifer conductance is computed with degree of saturation (only when # conductivity_profile is set to "exponential") end -storativity(A::UnconfinedAquifer) = A.specific_yield -storativity(A::ConfinedAquifer) = A.storativity +function UnconfinedAquiferParameters(nc, config, inds, top, bottom, area) + k = ncread(nc, config, "lateral.subsurface.conductivity"; sel = inds, type = Float) + specific_yield = + ncread(nc, config, "lateral.subsurface.specific_yield"; sel = inds, type = Float) + f = ncread( + nc, + config, + "lateral.subsurface.gwf_f"; + sel = inds, + type = Float, + defaults = 3.0, + ) + + parameters = + UnconfinedAquiferParameters{Float}(; k, top, bottom, area, specific_yield, f) + return parameters +end + +@get_units @grid_loc @with_kw struct UnconfinedAquiferVariables{T} + head::Vector{T} | "m" # hydraulic head [m] + conductance::Vector{T} | "m2 d-1" # conductance + volume::Vector{T} | "m3" # total volume of water that can be released m +end + +@with_kw struct UnconfinedAquifer{T} <: Aquifer + parameters::UnconfinedAquiferParameters{T} + variables::UnconfinedAquiferVariables{T} +end + +function UnconfinedAquifer(nc, config, inds, top, bottom, area, conductance, head) + parameters = UnconfinedAquiferParameters(nc, config, inds, top, bottom, area) + + volume = @. (min(top, head) - bottom) * area * parameters.specific_yield + variables = UnconfinedAquiferVariables{Float}(head, conductance, volume) + aquifer = UnconfinedAquifer{Float}(parameters, variables) + return aquifer +end + +storativity(A::UnconfinedAquifer) = A.parameters.specific_yield +storativity(A::ConfinedAquifer) = A.parameters.storativity """ harmonicmean_conductance(kH1, kH2, l1, l2, width) @@ -147,19 +190,20 @@ function harmonicmean_conductance(kH1, kH2, l1, l2, width) end function saturated_thickness(aquifer::UnconfinedAquifer, index::Int) - return min(aquifer.top[index], aquifer.head[index]) - aquifer.bottom[index] + return min(aquifer.parameters.top[index], aquifer.variables.head[index]) - + aquifer.parameters.bottom[index] end function saturated_thickness(aquifer::ConfinedAquifer, index::Int) - return aquifer.top[index] - aquifer.bottom[index] + return aquifer.parameters.top[index] - aquifer.parameters.bottom[index] end function saturated_thickness(aquifer::UnconfinedAquifer) - @. min(aquifer.top, aquifer.head) - aquifer.bottom + @. min(aquifer.parameters.top, aquifer.variables.head) - aquifer.parameters.bottom end function saturated_thickness(aquifer::ConfinedAquifer) - @. aquifer.top - aquifer.bottom + @. aquifer.parameters.top - aquifer.parameters.bottom end """ @@ -177,10 +221,10 @@ function horizontal_conductance( aquifer::A, connectivity::Connectivity, ) where {A <: Aquifer} - k1 = aquifer.k[i] - k2 = aquifer.k[j] - H1 = aquifer.top[i] - aquifer.bottom[i] - H2 = aquifer.top[j] - aquifer.bottom[j] + k1 = aquifer.parameters.k[i] + k2 = aquifer.parameters.k[j] + H1 = aquifer.parameters.top[i] - aquifer.parameters.bottom[i] + H2 = aquifer.parameters.top[j] - aquifer.parameters.bottom[j] length1 = connectivity.length1[nzi] length2 = connectivity.length2[nzi] width = connectivity.width[nzi] @@ -204,7 +248,7 @@ function initialize_conductance!( # Loop over connections for cell j for nzi in connections(connectivity, i) j = connectivity.rowval[nzi] - aquifer.conductance[nzi] = + aquifer.variables.conductance[nzi] = horizontal_conductance(i, j, nzi, aquifer, connectivity) end end @@ -218,7 +262,7 @@ function conductance( conductivity_profile::String, connectivity::Connectivity, ) - return aquifer.conductance[nzi] + return aquifer.variables.conductance[nzi] end """ @@ -257,17 +301,21 @@ function conductance( ) if conductivity_profile == "exponential" # Extract required variables - zi1 = aquifer.top[i] - aquifer.head[i] - zi2 = aquifer.top[j] - aquifer.head[j] - thickness1 = aquifer.top[i] - aquifer.bottom[i] - thickness2 = aquifer.top[j] - aquifer.bottom[j] + zi1 = aquifer.parameters.top[i] - aquifer.variables.head[i] + zi2 = aquifer.parameters.top[j] - aquifer.variables.head[j] + thickness1 = aquifer.parameters.top[i] - aquifer.parameters.bottom[i] + thickness2 = aquifer.parameters.top[j] - aquifer.parameters.bottom[j] # calculate conductivity values corrected for depth of water table k1 = - (aquifer.k[i] / aquifer.f[i]) * - (exp(-aquifer.f[i] * zi1) - exp(-aquifer.f[i] * thickness1)) + (aquifer.parameters.k[i] / aquifer.parameters.f[i]) * ( + exp(-aquifer.parameters.f[i] * zi1) - + exp(-aquifer.parameters.f[i] * thickness1) + ) k2 = - (aquifer.k[j] / aquifer.f[j]) * - (exp(-aquifer.f[j] * zi2) - exp(-aquifer.f[j] * thickness2)) + (aquifer.parameters.k[j] / aquifer.parameters.f[j]) * ( + exp(-aquifer.parameters.f[j] * zi2) - + exp(-aquifer.parameters.f[j] * thickness2) + ) return harmonicmean_conductance( k1, k2, @@ -276,16 +324,18 @@ function conductance( connectivity.width[nzi], ) elseif conductivity_profile == "uniform" - head_i = aquifer.head[i] - head_j = aquifer.head[j] + head_i = aquifer.variables.head[i] + head_j = aquifer.variables.head[j] if head_i >= head_j saturation = - saturated_thickness(aquifer, i) / (aquifer.top[i] - aquifer.bottom[i]) + saturated_thickness(aquifer, i) / + (aquifer.parameters.top[i] - aquifer.parameters.bottom[i]) else saturation = - saturated_thickness(aquifer, j) / (aquifer.top[j] - aquifer.bottom[j]) + saturated_thickness(aquifer, j) / + (aquifer.parameters.top[j] - aquifer.parameters.bottom[j]) end - return saturation * aquifer.conductance[nzi] + return saturation * aquifer.variables.conductance[nzi] else error( """An unknown "conductivity_profile" is specified in the TOML file ($conductivity_profile). @@ -301,7 +351,7 @@ function flux!(Q, aquifer, connectivity, conductivity_profile) for nzi in connections(connectivity, i) # connection from i -> j j = connectivity.rowval[nzi] - delta_head = aquifer.head[i] - aquifer.head[j] + delta_head = aquifer.variables.head[i] - aquifer.variables.head[j] cond = conductance(aquifer, i, j, nzi, conductivity_profile, connectivity) Q[i] -= cond * delta_head end @@ -309,11 +359,32 @@ function flux!(Q, aquifer, connectivity, conductivity_profile) return Q end -@get_units @grid_loc struct ConstantHead{T} +@get_units @grid_loc @with_kw struct ConstantHeadVariables{T} head::Vector{T} | "m" +end + +@get_units @grid_loc @with_kw struct ConstantHead{T} + variables::ConstantHeadVariables{T} index::Vector{Int} | "-" end +function ConstantHead(nc, config, inds) + constanthead = ncread( + nc, + config, + "lateral.subsurface.constant_head"; + sel = inds, + type = Float, + fill = mv, + ) + n = length(inds) + index_constanthead = filter(i -> !isequal(constanthead[i], mv), 1:n) + head = constanthead[index_constanthead] + variables = ConstantHeadVariables{Float}(head) + constant_head = ConstantHead{Float}(; variables, index = index_constanthead) + return constant_head +end + """ stable_timestep(aquifer) @@ -324,25 +395,28 @@ The following criterion can be found in Chu & Willis (1984) """ function stable_timestep(aquifer, conductivity_profile::String) dt_min = Inf - for i in eachindex(aquifer.head) + for i in eachindex(aquifer.variables.head) if conductivity_profile == "exponential" - zi = aquifer.top[i] - aquifer.head[i] - thickness = aquifer.top[i] - aquifer.bottom[i] + zi = aquifer.parameters.top[i] - aquifer.variables.head[i] + thickness = aquifer.parameters.top[i] - aquifer.parameters.bottom[i] value = - (aquifer.k[i] / aquifer.f[i]) * - (exp(-aquifer.f[i] * zi) - exp(-aquifer.f[i] * thickness)) + (aquifer.parameters.k[i] / aquifer.parameters.f[i]) * ( + exp(-aquifer.parameters.f[i] * zi) - + exp(-aquifer.parameters.f[i] * thickness) + ) elseif conductivity_profile == "uniform" - value = aquifer.k[i] * saturated_thickness(aquifer, i) + value = aquifer.parameters.k[i] * saturated_thickness(aquifer, i) end - dt = aquifer.area[i] * storativity(aquifer)[i] / value + dt = aquifer.parameters.area[i] * storativity(aquifer)[i] / value dt_min = dt < dt_min ? dt : dt_min end return 0.25 * dt_min end -minimum_head(aquifer::ConfinedAquifer) = aquifer.head -minimum_head(aquifer::UnconfinedAquifer) = max.(aquifer.head, aquifer.bottom) +minimum_head(aquifer::ConfinedAquifer) = aquifer.variables.head +minimum_head(aquifer::UnconfinedAquifer) = + max.(aquifer.variables.head, aquifer.parameters.bottom) function update!(gwf, Q, dt, conductivity_profile) Q .= 0.0 # TODO: Probably remove this when linking with other components @@ -350,13 +424,15 @@ function update!(gwf, Q, dt, conductivity_profile) for boundary in gwf.boundaries flux!(Q, boundary, gwf.aquifer) end - gwf.aquifer.head .+= (Q ./ gwf.aquifer.area .* dt ./ storativity(gwf.aquifer)) + gwf.aquifer.variables.head .+= + (Q ./ gwf.aquifer.parameters.area .* dt ./ storativity(gwf.aquifer)) # Set constant head (dirichlet) boundaries - gwf.aquifer.head[gwf.constanthead.index] .= gwf.constanthead.head + gwf.aquifer.variables.head[gwf.constanthead.index] .= gwf.constanthead.variables.head # Make sure no heads ends up below an unconfined aquifer bottom - gwf.aquifer.head .= minimum_head(gwf.aquifer) - gwf.aquifer.volume .= - saturated_thickness(gwf.aquifer) .* gwf.aquifer.area .* storativity(gwf.aquifer) + gwf.aquifer.variables.head .= minimum_head(gwf.aquifer) + gwf.aquifer.variables.volume .= + saturated_thickness(gwf.aquifer) .* gwf.aquifer.parameters.area .* + storativity(gwf.aquifer) return nothing end @@ -381,9 +457,10 @@ end function get_water_depth( gwf::GroundwaterFlow{A, C, CH, B}, ) where {A <: UnconfinedAquifer, C, CH, B} - gwf.aquifer.head .= min.(gwf.aquifer.head, gwf.aquifer.top) - gwf.aquifer.head[gwf.constanthead.index] .= gwf.constanthead.head - wtd = gwf.aquifer.top .- gwf.aquifer.head + gwf.aquifer.variables.head .= + min.(gwf.aquifer.variables.head, gwf.aquifer.parameters.top) + gwf.aquifer.variables.head[gwf.constanthead.index] .= gwf.constanthead.variables.head + wtd = gwf.aquifer.parameters.top .- gwf.aquifer.variables.head return wtd end @@ -391,8 +468,10 @@ function get_exfiltwater( gwf::GroundwaterFlow{A, C, CH, B}, ) where {A <: UnconfinedAquifer, C, CH, B} exfiltwater = - (gwf.aquifer.head .- min.(gwf.aquifer.head, gwf.aquifer.top)) .* - storativity(gwf.aquifer) + ( + gwf.aquifer.variables.head .- + min.(gwf.aquifer.variables.head, gwf.aquifer.parameters.top) + ) .* storativity(gwf.aquifer) exfiltwater[gwf.constanthead.index] .= 0 return exfiltwater end \ No newline at end of file diff --git a/src/groundwater/boundary_conditions.jl b/src/groundwater/boundary_conditions.jl index 48e883527..5e3a45740 100644 --- a/src/groundwater/boundary_conditions.jl +++ b/src/groundwater/boundary_conditions.jl @@ -1,6 +1,6 @@ function check_flux(flux, aquifer::UnconfinedAquifer, index::Int) # Check if cell is dry - if aquifer.head[index] <= aquifer.bottom[index] + if aquifer.variables.head[index] <= aquifer.parameters.bottom[index] # If cell is dry, no negative flux is allowed return max(0, flux) else @@ -11,91 +11,190 @@ end # Do nothing for a confined aquifer: aquifer can always provide flux check_flux(flux, aquifer::ConfinedAquifer, index::Int) = flux -@get_units @grid_loc struct River{T} <: AquiferBoundaryCondition - stage::Vector{T} | "m" +@get_units @grid_loc @with_kw struct RiverParameters{T} infiltration_conductance::Vector{T} | "m2 d-1" exfiltration_conductance::Vector{T} | "m2 d-1" bottom::Vector{T} | "m" +end + +@get_units @grid_loc @with_kw struct RiverVariables{T} + stage::Vector{T} | "m" flux::Vector{T} | "m3 d-1" +end + +function RiverVariables(n) + variables = RiverVariables{Float}(; stage = fill(mv, n), flux = fill(mv, n)) + return variables +end + +@get_units @grid_loc @with_kw struct River{T} <: AquiferBoundaryCondition + parameters::RiverParameters{T} + variables::RiverVariables{T} index::Vector{Int} | "-" end +function River(nc, config, inds, index) + infiltration_conductance = ncread( + nc, + config, + "lateral.subsurface.infiltration_conductance"; + sel = inds, + type = Float, + ) + exfiltration_conductance = ncread( + nc, + config, + "lateral.subsurface.exfiltration_conductance"; + sel = inds, + type = Float, + ) + bottom = ncread(nc, config, "lateral.subsurface.river_bottom"; sel = inds, type = Float) + + parameters = + RiverParameters{Float}(infiltration_conductance, exfiltration_conductance, bottom) + n = length(inds) + variables = RiverVariables(n) + river = River(parameters, variables, index) + return river +end + function flux!(Q, river::River, aquifer) for (i, index) in enumerate(river.index) - head = aquifer.head[index] - stage = river.stage[i] + head = aquifer.variables.head[index] + stage = river.variables.stage[i] if stage > head - cond = river.infiltration_conductance[i] - delta_head = min(stage - river.bottom[i], stage - head) + cond = river.parameters.infiltration_conductance[i] + delta_head = min(stage - river.parameters.bottom[i], stage - head) else - cond = river.exfiltration_conductance[i] + cond = river.parameters.exfiltration_conductance[i] delta_head = stage - head end - river.flux[i] = check_flux(cond * delta_head, aquifer, index) - Q[index] += river.flux[i] + river.variables.flux[i] = check_flux(cond * delta_head, aquifer, index) + Q[index] += river.variables.flux[i] end return Q end -@get_units @grid_loc struct Drainage{T} <: AquiferBoundaryCondition +@get_units @grid_loc @with_kw struct DrainageParameters{T} elevation::Vector{T} | "m" conductance::Vector{T} | "m2 d-1" +end + +@get_units @grid_loc @with_kw struct DrainageVariables{T} flux::Vector{T} | "m3 d-1" +end + +@get_units @grid_loc @with_kw struct Drainage{T} <: AquiferBoundaryCondition + parameters::DrainageParameters{T} + variables::DrainageVariables{T} index::Vector{Int} | "-" end +function Drainage(nc, config, inds, index) + drain_elevation = ncread( + nc, + config, + "lateral.subsurface.drain_elevation"; + sel = inds, + type = Float, + fill = mv, + ) + drain_conductance = ncread( + nc, + config, + "lateral.subsurface.drain_conductance"; + sel = inds, + type = Float, + fill = mv, + ) + elevation = drain_elevation[index] + conductance = drain_conductance[index] + parameters = DrainageParameters{Float}(; elevation, conductance) + variables = DrainageVariables{Float}(; flux = fill(mv, length(index))) + + drains = Drainage{Float}(parameters, variables, index) + return drains +end + function flux!(Q, drainage::Drainage, aquifer) for (i, index) in enumerate(drainage.index) - cond = drainage.conductance[i] - 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] + cond = drainage.parameters.conductance[i] + delta_head = + min(0, drainage.parameters.elevation[i] - aquifer.variables.head[index]) + drainage.variables.flux[i] = check_flux(cond * delta_head, aquifer, index) + Q[index] += drainage.variables.flux[i] end return Q end -@get_units @grid_loc struct HeadBoundary{T} <: AquiferBoundaryCondition - head::Vector{T} | "m" +@get_units @grid_loc @with_kw struct HeadBoundaryParameters{T} conductance::Vector{T} | "m2 d-1" +end + +@get_units @grid_loc @with_kw struct HeadBoundaryVariables{T} + head::Vector{T} | "m" flux::Vector{T} | "m3 d-1" +end + +@get_units @grid_loc @with_kw struct HeadBoundary{T} <: AquiferBoundaryCondition + parameters::HeadBoundaryParameters{T} + variables::HeadBoundaryVariables{T} index::Vector{Int} | "-" end function flux!(Q, headboundary::HeadBoundary, aquifer) for (i, index) in enumerate(headboundary.index) - cond = headboundary.conductance[i] - delta_head = headboundary.head[i] - aquifer.head[index] - headboundary.flux[i] = check_flux(cond * delta_head, aquifer, index) - Q[index] += headboundary.flux[i] + cond = headboundary.parameters.conductance[i] + delta_head = headboundary.variables.head[i] - aquifer.variables.head[index] + headboundary.variables.flux[i] = check_flux(cond * delta_head, aquifer, index) + Q[index] += headboundary.variables.flux[i] end return Q end -@get_units @grid_loc struct Recharge{T} <: AquiferBoundaryCondition +@get_units @grid_loc @with_kw struct RechargeVariables{T} rate::Vector{T} | "m d-1" flux::Vector{T} | "m3 d-1" +end + +@get_units @grid_loc @with_kw struct Recharge{T} <: AquiferBoundaryCondition + variables::RechargeVariables{T} index::Vector{Int} | "-" end +function Recharge(rate, flux, index) + variables = RechargeVariables{Float}(rate, flux) + recharge = Recharge{Float}(variables, index) + return recharge +end + function flux!(Q, recharge::Recharge, aquifer) for (i, index) in enumerate(recharge.index) - recharge.flux[i] = - check_flux(recharge.rate[i] * aquifer.area[index], aquifer, index) - Q[index] += recharge.flux[i] + recharge.variables.flux[i] = check_flux( + recharge.variables.rate[i] * aquifer.parameters.area[index], + aquifer, + index, + ) + Q[index] += recharge.variables.flux[i] end return Q end -@get_units @grid_loc struct Well{T} <: AquiferBoundaryCondition +@get_units @grid_loc @with_kw struct WellVariables{T} volumetric_rate::Vector{T} | "m3 d-1" flux::Vector{T} | "m3 d-1" +end + +@get_units @grid_loc @with_kw struct Well{T} <: AquiferBoundaryCondition + variables::WellVariables{T} index::Vector{Int} | "-" end function flux!(Q, well::Well, aquifer) for (i, index) in enumerate(well.index) - well.flux[i] = check_flux(well.volumetric_rate[i], aquifer, index) - Q[index] += well.flux[i] + well.variables.flux[i] = + check_flux(well.variables.volumetric_rate[i], aquifer, index) + Q[index] += well.variables.flux[i] end return Q end diff --git a/src/sbm_gwf_model.jl b/src/sbm_gwf_model.jl index b6b05a8ff..4b3a518bd 100644 --- a/src/sbm_gwf_model.jl +++ b/src/sbm_gwf_model.jl @@ -82,7 +82,7 @@ function initialize_sbm_gwf_model(config::Config) pits = zeros(Bool, modelsize_2d) if do_reservoirs reservoirs, resindex, reservoir, pits = - initialize_simple_reservoir(config, nc, inds_riv, nriv, pits, tosecond(dt)) + SimpleReservoir(config, nc, inds_riv, nriv, pits, tosecond(dt)) else reservoir = () reservoirs = nothing @@ -91,8 +91,7 @@ function initialize_sbm_gwf_model(config::Config) # lakes if do_lakes - lakes, lakeindex, lake, pits = - initialize_lake(config, nc, inds_riv, nriv, pits, tosecond(dt)) + lakes, lakeindex, lake, pits = Lake(config, nc, inds_riv, nriv, pits, tosecond(dt)) else lake = () lakes = nothing @@ -209,91 +208,39 @@ function initialize_sbm_gwf_model(config::Config) # unconfined aquifer if do_constanthead - constanthead = ncread( - nc, - config, - "lateral.subsurface.constant_head"; - sel = inds, - type = Float, - fill = mv, - ) - index_constanthead = filter(i -> !isequal(constanthead[i], mv), 1:n) - constant_head = ConstantHead(constanthead[index_constanthead], index_constanthead) + constant_head = ConstantHead(nc, config, inds) else - constant_head = ConstantHead{Float}(Float[], Int64[]) + variables = ConstantHeadVariables{Float}(; head = Float[]) + constant_head = ConstantHead{Float}(; variables, index = Int64[]) end - conductivity = - ncread(nc, config, "lateral.subsurface.conductivity"; sel = inds, type = Float) - specific_yield = - ncread(nc, config, "lateral.subsurface.specific_yield"; sel = inds, type = Float) - gwf_f = ncread( - nc, - config, - "lateral.subsurface.gwf_f"; - sel = inds, - type = Float, - defaults = 3.0, - ) - conductivity_profile = - get(config.input.lateral.subsurface, "conductivity_profile", "uniform")::String - connectivity = Connectivity(inds, rev_inds, xl, yl) + initial_head = altitude .- lhm.soil.variables.zi / 1000.0 # cold state for groundwater head based on SBM zi initial_head[index_river] = altitude[index_river] - if do_constanthead - initial_head[constant_head.index] = constant_head.head + initial_head[constant_head.index] = constant_head.variables.head end bottom = altitude .- lhm.soil.parameters.soilthickness ./ Float(1000.0) area = xl .* yl - volume = @. (min(altitude, initial_head) - bottom) * area * specific_yield # total volume than can be released - + conductance = zeros(Float, connectivity.nconnection) aquifer = UnconfinedAquifer( - initial_head, - conductivity, + nc, + config, + inds, altitude, bottom, area, - specific_yield, - zeros(Float, connectivity.nconnection), # conductance - volume, - gwf_f, + conductance, + initial_head, ) # river boundary of unconfined aquifer - infiltration_conductance = ncread( - nc, - config, - "lateral.subsurface.infiltration_conductance"; - sel = inds_riv, - type = Float, - ) - exfiltration_conductance = ncread( - nc, - config, - "lateral.subsurface.exfiltration_conductance"; - sel = inds_riv, - type = Float, - ) - river_bottom = - ncread(nc, config, "lateral.subsurface.river_bottom"; sel = inds_riv, type = Float) - - river_flux = fill(mv, nriv) - river_stage = fill(mv, nriv) - river = River( - river_stage, - infiltration_conductance, - exfiltration_conductance, - river_bottom, - river_flux, - index_river, - ) + river = River(nc, config, inds_riv, index_river) # recharge boundary of unconfined aquifer - r = fill(mv, n) - recharge = Recharge(r, zeros(Float, n), collect(1:n)) + recharge = Recharge(fill(mv, n), zeros(Float, n), collect(1:n)) # drain boundary of unconfined aquifer (optional) if do_drains @@ -310,32 +257,11 @@ function initialize_sbm_gwf_model(config::Config) @info "$n_false_drain drain locations are removed that occur where overland flow is not possible (overland flow width is zero)" end - inds_drain, rev_inds_drain = active_indices(drain_2d, 0) - drain_elevation = ncread( - nc, - config, - "lateral.subsurface.drain_elevation"; - sel = inds, - type = Float, - fill = mv, - ) - drain_conductance = ncread( - nc, - config, - "lateral.subsurface.drain_conductance"; - sel = inds, - type = Float, - fill = mv, - ) + inds_drain, rev_inds_drain = active_indices(drain_2d, 0) index_drain = filter(i -> !isequal(drain[i], 0), 1:n) - drain_flux = fill(mv, length(index_drain)) - drains = Drainage( - drain_elevation[index_drain], - drain_conductance[index_drain], - drain_flux, - index_drain, - ) + + drains = Drainage(nc, config, inds, index_drain) drain = (indices = inds_drain, reverse_indices = rev_inds_drain) aquifer_boundaries = AquiferBoundaryCondition[recharge, river, drains] else @@ -539,8 +465,8 @@ function update!(model::Model{N, L, V, R, W, T}) where {N, L, V, R, W, T <: SbmG update!(vertical, lateral, network, config) # set river stage (groundwater) to average h from kinematic wave - lateral.subsurface.river.stage .= - lateral.river.variables.h_av .+ lateral.subsurface.river.bottom + lateral.subsurface.river.variables.stage .= + lateral.river.variables.h_av .+ lateral.subsurface.river.parameters.bottom # determine stable time step for groundwater flow conductivity_profile = @@ -556,9 +482,10 @@ function update!(model::Model{N, L, V, R, W, T}) where {N, L, V, R, W, T <: SbmG Q = zeros(lateral.subsurface.flow.connectivity.ncell) # exchange of recharge between SBM soil model and groundwater flow domain # recharge rate groundwater is required in units [m d⁻¹] - @. lateral.subsurface.recharge.rate = soil.variables.recharge / 1000.0 * (1.0 / dt_sbm) + @. lateral.subsurface.recharge.variables.rate = + soil.variables.recharge / 1000.0 * (1.0 / dt_sbm) if do_water_demand - @. lateral.subsurface.recharge.rate -= + @. lateral.subsurface.recharge.variables.rate -= vertical.allocation.variables.act_groundwater_abst / 1000.0 * (1.0 / dt_sbm) end # update groundwater domain @@ -568,7 +495,8 @@ function update!(model::Model{N, L, V, R, W, T}) where {N, L, V, R, W, T <: SbmG update!(soil, (; runoff, demand, subsurface = lateral.subsurface.flow)) ssf_toriver = zeros(length(soil.variables.zi)) - ssf_toriver[inds_riv] = -lateral.subsurface.river.flux ./ tosecond(basetimestep) + ssf_toriver[inds_riv] = + -lateral.subsurface.river.variables.flux ./ tosecond(basetimestep) surface_routing!(model; ssf_toriver = ssf_toriver) return nothing diff --git a/src/states.jl b/src/states.jl index 13e598e90..b0515c3a6 100644 --- a/src/states.jl +++ b/src/states.jl @@ -217,7 +217,7 @@ function extract_required_states(config::Config) add_to_required_states(required_states, (:vertical, :soil, :variables), soil_states) # Add subsurface states to dict if model_type == "sbm_gwf" - key_entry = (:lateral, :subsurface, :flow, :aquifer) + key_entry = (:lateral, :subsurface, :flow, :aquifer, :variables) else key_entry = (:lateral, :subsurface, :variables) end diff --git a/test/groundwater.jl b/test/groundwater.jl index def5ce7e5..427c4c96b 100644 --- a/test/groundwater.jl +++ b/test/groundwater.jl @@ -35,28 +35,35 @@ function homogenous_aquifer(nrow, ncol) connectivity = Wflow.Connectivity(indices, reverse_indices, dx, dy) ncell = connectivity.ncell - conf_aqf = Wflow.ConfinedAquifer( - [0.0, 7.5, 20.0], # head - fill(10.0, ncell), # k - fill(10.0, ncell), # top - fill(0.0, ncell), # bottom - fill(100.0, ncell), # area - fill(0.1, ncell), # specific storage - fill(1.0, ncell), # storativity - fill(0.0, connectivity.nconnection), # conductance - fill(0.0, ncell), # total volume that can be released + parameters = Wflow.ConfinedAquiferParameters(; + k = fill(10.0, ncell), + top = fill(10.0, ncell), + bottom = fill(0.0, ncell), + area = fill(100.0, ncell), + specific_storage = fill(0.1, ncell), + storativity = fill(1.0, ncell), ) - unconf_aqf = Wflow.UnconfinedAquifer( - [0.0, 7.5, 20.0], # head - fill(10.0, ncell), # k - fill(10.0, ncell), # top - fill(0.0, ncell), # bottom - fill(100.0, ncell), # area - fill(0.15, ncell), # specific yield - fill(0.0, connectivity.nconnection), # conductance - fill(0.0, ncell), # total volume that can be released - fill(3.0, ncell), # conductance reduction factor + variables = Wflow.ConfinedAquiferVariables(; + head = [0.0, 7.5, 20.0], + conductance = fill(0.0, connectivity.nconnection), + volume = fill(0.0, ncell), ) + conf_aqf = Wflow.ConfinedAquifer(; parameters, variables) + + parameters = Wflow.UnconfinedAquiferParameters(; + k = fill(10.0, ncell), + top = fill(10.0, ncell), + bottom = fill(0.0, ncell), + area = fill(100.0, ncell), + specific_yield = fill(0.15, ncell), + f = fill(3.0, ncell), + ) + variables = Wflow.UnconfinedAquiferVariables(; + head = [0.0, 7.5, 20.0], + conductance = fill(0.0, connectivity.nconnection), + volume = fill(0.0, ncell), + ) + unconf_aqf = Wflow.UnconfinedAquifer(; parameters, variables) return (connectivity, conf_aqf, unconf_aqf) end @@ -244,19 +251,19 @@ end end @testset "minimum_head-confined" begin - original_head = copy(conf_aqf.head) - conf_aqf.head[1] = -10.0 + original_head = copy(conf_aqf.variables.head) + conf_aqf.variables.head[1] = -10.0 @test Wflow.check_flux(-1.0, conf_aqf, 1) == -1.0 @test Wflow.minimum_head(conf_aqf)[1] == -10.0 - conf_aqf.head .= original_head + conf_aqf.variables.head .= original_head end @testset "minimum_head-unconfined" begin - original_head = copy(unconf_aqf.head) - unconf_aqf.head[1] = -10.0 + original_head = copy(unconf_aqf.variables.head) + unconf_aqf.variables.head[1] = -10.0 @test Wflow.check_flux(-1.0, unconf_aqf, 1) == 0.0 @test Wflow.minimum_head(conf_aqf)[1] == 0.0 - unconf_aqf.head .= original_head + unconf_aqf.variables.head .= original_head end @testset "stable_timestep" begin @@ -288,14 +295,13 @@ end end @testset "river" begin - river = Wflow.River( - [2.0, 2.0], - [100.0, 100.0], - [200.0, 200.0], - [1.0, 1.0], - [0.0, 0.0], - [1, 3], + parameters = Wflow.RiverParameters(; + infiltration_conductance = [100.0, 100.0], + exfiltration_conductance = [200.0, 200.0], + bottom = [1.0, 1.0], ) + variables = Wflow.RiverVariables(; stage = [2.0, 2.0], flux = [0.0, 0.0]) + river = Wflow.River(; parameters, variables, index = [1, 3]) Q = zeros(3) Wflow.flux!(Q, river, conf_aqf) # infiltration, below bottom, flux is (stage - bottom) * inf_cond @@ -305,7 +311,12 @@ end end @testset "drainage" begin - drainage = Wflow.Drainage([2.0, 2.0], [100.0, 100.0], [0.0, 0.0], [1, 2]) + parameters = Wflow.DrainageParameters(; + elevation = [2.0, 2.0], + conductance = [100.0, 100.0], + ) + variables = Wflow.DrainageVariables(; flux = [0.0, 0.0]) + drainage = Wflow.Drainage(; parameters, variables, index = [1, 2]) Q = zeros(3) Wflow.flux!(Q, drainage, conf_aqf) @test Q[1] == 0.0 @@ -313,8 +324,10 @@ end end @testset "headboundary" begin - headboundary = - Wflow.HeadBoundary([2.0, 2.0], [100.0, 100.0], [0.0, 0.0], [1, 2]) + parameters = Wflow.HeadBoundaryParameters(; conductance = [100.0, 100.0]) + variables = Wflow.HeadBoundaryVariables(; head = [2.0, 2.0], flux = [0.0, 0.0]) + + headboundary = Wflow.HeadBoundary(; parameters, variables, index = [1, 2]) Q = zeros(3) Wflow.flux!(Q, headboundary, conf_aqf) @test Q[1] == 100.0 * (2.0 - 0.0) @@ -322,14 +335,19 @@ end end @testset "recharge" begin - recharge = Wflow.Recharge([1.0e-3, 1.0e-3, 1.0e-3], [0.0, 0.0, 0.0], [1, 2, 3]) + variables = Wflow.RechargeVariables(; + rate = [1.0e-3, 1.0e-3, 1.0e-3], + flux = [0.0, 0.0, 0.0], + ) + recharge = Wflow.Recharge(; variables, index = [1, 2, 3]) Q = zeros(3) Wflow.flux!(Q, recharge, conf_aqf) @test all(Q .== 1.0e-3 * 100.0) end @testset "well" begin - well = Wflow.Well([-1000.0], [0.0], [1]) + variables = Wflow.WellVariables(; volumetric_rate = [-1000.0], flux = [0.0]) + well = Wflow.Well(; variables, index = [1]) Q = zeros(3) Wflow.flux!(Q, well, conf_aqf) @test Q[1] == -1000.0 @@ -338,7 +356,8 @@ end @testset "integration: steady 1D" begin connectivity, aquifer, _ = homogenous_aquifer(3, 1) - constanthead = Wflow.ConstantHead([2.0, 4.0], [1, 3]) + variables = Wflow.ConstantHeadVariables(; head = [2.0, 4.0]) + constanthead = Wflow.ConstantHead(; variables, index = [1, 3]) conductivity_profile = "uniform" gwf = Wflow.GroundwaterFlow{Wflow.Float}(; aquifer = aquifer, @@ -347,7 +366,8 @@ end boundaries = Wflow.AquiferBoundaryCondition[], ) # Set constant head (dirichlet) boundaries - gwf.aquifer.head[gwf.constanthead.index] .= gwf.constanthead.head + gwf.aquifer.variables.head[gwf.constanthead.index] .= + gwf.constanthead.variables.head Q = zeros(3) dt = 0.25 # days @@ -355,12 +375,13 @@ end Wflow.update!(gwf, Q, dt, conductivity_profile) end - @test gwf.aquifer.head ≈ [2.0, 3.0, 4.0] + @test gwf.aquifer.variables.head ≈ [2.0, 3.0, 4.0] end @testset "integration: steady 1D, exponential conductivity" begin connectivity, aquifer, _ = homogenous_aquifer(3, 1) - constanthead = Wflow.ConstantHead([2.0, 4.0], [1, 3]) + variables = Wflow.ConstantHeadVariables(; head = [2.0, 4.0]) + constanthead = Wflow.ConstantHead(; variables, index = [1, 3]) conductivity_profile = "exponential" gwf = Wflow.GroundwaterFlow{Wflow.Float}(; aquifer = aquifer, @@ -369,7 +390,8 @@ end boundaries = Wflow.AquiferBoundaryCondition[], ) # Set constant head (dirichlet) boundaries - gwf.aquifer.head[gwf.constanthead.index] .= gwf.constanthead.head + gwf.aquifer.variables.head[gwf.constanthead.index] .= + gwf.constanthead.variables.head Q = zeros(3) dt = 0.25 # days @@ -377,7 +399,7 @@ end Wflow.update!(gwf, Q, dt, conductivity_profile) end - @test gwf.aquifer.head ≈ [2.0, 3.0, 4.0] + @test gwf.aquifer.variables.head ≈ [2.0, 3.0, 4.0] end @testset "integration: unconfined transient 1D" begin @@ -402,19 +424,25 @@ end connectivity = Wflow.Connectivity(indices, reverse_indices, dx, dy) ncell = connectivity.ncell xc = collect(range(0.0; stop = aquifer_length - cellsize, step = cellsize)) - aquifer = Wflow.UnconfinedAquifer( - initial_head.(xc), - fill(conductivity, ncell), - fill(top, ncell), - fill(bottom, ncell), - fill(cellsize * cellsize, ncell), - fill(specific_yield, ncell), - fill(0.0, connectivity.nconnection), - fill(0.0, ncell), - fill(gwf_f, ncell), + + variables = Wflow.UnconfinedAquiferVariables(; + head = initial_head.(xc), + conductance = fill(0.0, connectivity.nconnection), + volume = fill(0.0, ncell), + ) + parameters = Wflow.UnconfinedAquiferParameters(; + k = fill(conductivity, ncell), + top = fill(top, ncell), + bottom = fill(bottom, ncell), + area = fill(cellsize * cellsize, ncell), + specific_yield = fill(specific_yield, ncell), + f = fill(gwf_f, ncell), ) + + aquifer = Wflow.UnconfinedAquifer(; parameters, variables) # constant head on left boundary, 0 at 0 - constanthead = Wflow.ConstantHead([0.0], [1]) + variables = Wflow.ConstantHeadVariables(; head = [0.0]) + constanthead = Wflow.ConstantHead(; variables, index = [1]) gwf = Wflow.GroundwaterFlow{Wflow.Float}(; aquifer = aquifer, connectivity = connectivity, @@ -431,7 +459,7 @@ end for i in 1:nstep Wflow.update!(gwf, Q, dt, conductivity_profile) # Gradient dh/dx is positive, all flow to the left - @test all(diff(gwf.aquifer.head) .> 0.0) + @test all(diff(gwf.aquifer.variables.head) .> 0.0) end head_analytical = [ @@ -444,7 +472,7 @@ end beta, ) for x in xc ] - difference = gwf.aquifer.head .- head_analytical + difference = gwf.aquifer.variables.head .- head_analytical # @test all(difference .< ?) #TODO end @@ -470,19 +498,25 @@ end connectivity = Wflow.Connectivity(indices, reverse_indices, dx, dy) ncell = connectivity.ncell xc = collect(range(0.0; stop = aquifer_length - cellsize, step = cellsize)) - aquifer = Wflow.UnconfinedAquifer( - initial_head.(xc), - fill(conductivity, ncell), - fill(top, ncell), - fill(bottom, ncell), - fill(cellsize * cellsize, ncell), - fill(specific_yield, ncell), - fill(0.0, connectivity.nconnection), - fill(0.0, ncell), - fill(gwf_f, ncell), + + variables = Wflow.UnconfinedAquiferVariables(; + head = initial_head.(xc), + conductance = fill(0.0, connectivity.nconnection), + volume = fill(0.0, ncell), + ) + parameters = Wflow.UnconfinedAquiferParameters(; + k = fill(conductivity, ncell), + top = fill(top, ncell), + bottom = fill(bottom, ncell), + area = fill(cellsize * cellsize, ncell), + specific_yield = fill(specific_yield, ncell), + f = fill(gwf_f, ncell), ) + + aquifer = Wflow.UnconfinedAquifer(; parameters, variables) # constant head on left boundary, 0 at 0 - constanthead = Wflow.ConstantHead([0.0], [1]) + variables = Wflow.ConstantHeadVariables(; head = [0.0]) + constanthead = Wflow.ConstantHead(; variables, index = [1]) gwf = Wflow.GroundwaterFlow{Wflow.Float}(; aquifer = aquifer, connectivity = connectivity, @@ -499,7 +533,7 @@ end for i in 1:nstep Wflow.update!(gwf, Q, dt, conductivity_profile) # Gradient dh/dx is positive, all flow to the left - @test all(diff(gwf.aquifer.head) .> 0.0) + @test all(diff(gwf.aquifer.variables.head) .> 0.0) end head_analytical = [ @@ -512,7 +546,7 @@ end beta, ) for x in xc ] - difference = gwf.aquifer.head .- head_analytical + difference = gwf.aquifer.variables.head .- head_analytical # @test all(difference .< ?) #TODO end @@ -541,23 +575,29 @@ end indices, reverse_indices = Wflow.active_indices(domain, false) connectivity = Wflow.Connectivity(indices, reverse_indices, dx, dy) ncell = connectivity.ncell - aquifer = Wflow.ConfinedAquifer( - fill(startinghead, ncell), - fill(conductivity, ncell), - fill(top, ncell), - fill(bottom, ncell), - fill(cellsize * cellsize, ncell), - fill(specific_storage, ncell), - fill(storativity, ncell), - fill(0.0, connectivity.nconnection), # conductance, to be set - fill(0.0, ncell), # total volume that can be released, to be set + + parameters = Wflow.ConfinedAquiferParameters(; + k = fill(conductivity, ncell), + top = fill(top, ncell), + bottom = fill(bottom, ncell), + area = fill(cellsize * cellsize, ncell), + specific_storage = fill(specific_storage, ncell), + storativity = fill(storativity, ncell), + ) + variables = Wflow.ConfinedAquiferVariables(; + head = fill(startinghead, ncell), + conductance = fill(0.0, connectivity.nconnection), + volume = fill(0.0, ncell), ) + aquifer = Wflow.ConfinedAquifer(; parameters, variables) cell_index = reshape(collect(range(1, ncell; step = 1)), shape) indices = vcat(cell_index[1, :], cell_index[end, :])# , cell_index[:, 1], cell_index[:, end],) - constanthead = Wflow.ConstantHead(fill(10.0, size(indices)), indices) + variables = Wflow.ConstantHeadVariables(; head = fill(10.0, size(indices))) + constanthead = Wflow.ConstantHead(; variables, index = indices) # Place a well in the middle of the domain - well = Wflow.Well([discharge], [0.0], [reverse_indices[wellrow, wellrow]]) + variables = Wflow.WellVariables(; volumetric_rate = [discharge], flux = [0.0]) + well = Wflow.Well(; variables, index = [reverse_indices[wellrow, wellrow]]) gwf = Wflow.GroundwaterFlow{Wflow.Float}(; aquifer = aquifer, connectivity = connectivity, @@ -576,7 +616,7 @@ end end # test for symmetry on x and y axes - head = reshape(gwf.aquifer.head, shape) + head = reshape(gwf.aquifer.variables.head, shape) @test head[1:halfnrow, :] ≈ head[end:-1:(halfnrow + 2), :] @test head[:, 1:halfnrow] ≈ head[:, end:-1:(halfnrow + 2)] diff --git a/test/io.jl b/test/io.jl index aa9b0e37d..744d952f2 100644 --- a/test/io.jl +++ b/test/io.jl @@ -487,7 +487,7 @@ end @test (:vertical, :soil, :variables, :satwaterdepth) in required_states @test (:vertical, :soil, :variables, :ustorelayerdepth) in required_states @test (:vertical, :interception, :variables, :canopy_storage) in required_states - @test (:lateral, :subsurface, :flow, :aquifer, :head) in required_states + @test (:lateral, :subsurface, :flow, :aquifer, :variables, :head) in required_states @test (:lateral, :river, :variables, :q) in required_states @test (:lateral, :river, :variables, :h_av) in required_states @test (:lateral, :land, :variables, :h_av) in required_states diff --git a/test/run_sbm_gwf.jl b/test/run_sbm_gwf.jl index 4428562eb..9dd58e73d 100644 --- a/test/run_sbm_gwf.jl +++ b/test/run_sbm_gwf.jl @@ -56,17 +56,17 @@ end @testset "groundwater" begin gw = model.lateral.subsurface - @test gw.river.stage[1] ≈ 1.2123636929067039f0 - @test gw.flow.aquifer.head[17:21] ≈ [ + @test gw.river.variables.stage[1] ≈ 1.2123636929067039f0 + @test gw.flow.aquifer.variables.head[17:21] ≈ [ 1.2866380350225155f0, 1.3477853512604643f0, 1.7999999523162842f0, 1.6225103807809076f0, 1.4053590307668113f0, ] - @test gw.river.flux[1] ≈ -51.34674583702381f0 - @test gw.drain.flux[1] ≈ 0.0 - @test gw.recharge.rate[19] ≈ -0.0014241196552847502f0 + @test gw.river.variables.flux[1] ≈ -51.34674583702381f0 + @test gw.drain.variables.flux[1] ≈ 0.0 + @test gw.recharge.variables.rate[19] ≈ -0.0014241196552847502f0 end @testset "no drains" begin @@ -174,33 +174,17 @@ end @testset "groundwater warm start" begin gw = model.lateral.subsurface - @test gw.river.stage[1] ≈ 1.2031171676781156f0 - @test gw.flow.aquifer.head[17:21] ≈ [ + @test gw.river.variables.stage[1] ≈ 1.2031171676781156f0 + @test gw.flow.aquifer.variables.head[17:21] ≈ [ 1.2277456867225283f0, 1.286902494792006f0, 1.7999999523162842f0, 1.5901747932190804f0, 1.2094238817776854f0, ] - @test gw.river.flux[1] ≈ -6.692884222603261f0 - @test gw.drain.flux[1] ≈ 0.0 - @test gw.recharge.rate[19] ≈ -0.0014241196552847502f0 -end - -@testset "Exchange and grid location aquifer, recharge and constant head" begin - aquifer = model.lateral.subsurface.flow.aquifer - @test Wflow.exchange(aquifer.head) == true - @test Wflow.exchange(aquifer.k) == true - @test Wflow.grid_loc(aquifer, :head) == "node" - @test Wflow.grid_loc(aquifer, :k) == "node" - recharge = model.lateral.subsurface.recharge - @test Wflow.exchange(recharge.rate) == true - @test Wflow.exchange(recharge.flux) == true - @test Wflow.grid_loc(recharge, :rate) == "node" - @test Wflow.grid_loc(recharge, :flux) == "node" - constanthead = model.lateral.subsurface.flow.constanthead - @test Wflow.exchange(constanthead) == false - @test Wflow.grid_loc(constanthead, :head) == "node" + @test gw.river.variables.flux[1] ≈ -6.692884222603261f0 + @test gw.drain.variables.flux[1] ≈ 0.0 + @test gw.recharge.variables.rate[19] ≈ -0.0014241196552847502f0 end Wflow.close_files(model; delete_output = false) diff --git a/test/sbm_gwf_config.toml b/test/sbm_gwf_config.toml index 0fdf6df53..c917b409f 100644 --- a/test/sbm_gwf_config.toml +++ b/test/sbm_gwf_config.toml @@ -35,7 +35,7 @@ h = "h_land" h_av = "h_av_land" q = "q_land" -[state.lateral.subsurface.flow.aquifer] +[state.lateral.subsurface.flow.aquifer.variables] head = "head" [input] @@ -122,16 +122,16 @@ soilthickness = "soilthickness" [output.lateral.river.variables] q = "q" -[output.lateral.subsurface.flow.aquifer] +[output.lateral.subsurface.flow.aquifer.variables] head = "head" -[output.lateral.subsurface.recharge] +[output.lateral.subsurface.recharge.variables] rate = "rate" -[output.lateral.subsurface.drain] +[output.lateral.subsurface.drain.variables] flux = "drain_flux" -[output.lateral.subsurface.river] +[output.lateral.subsurface.river.variables] flux = "flux" [csv] @@ -145,4 +145,4 @@ parameter = "lateral.river.variables.q_av" [[csv.column]] header = "head" index = 5 -parameter = "lateral.subsurface.flow.aquifer.head" +parameter = "lateral.subsurface.flow.aquifer.variables.head" diff --git a/test/sbm_gwf_piave_demand_config.toml b/test/sbm_gwf_piave_demand_config.toml index a7338f28c..44d87cdba 100644 --- a/test/sbm_gwf_piave_demand_config.toml +++ b/test/sbm_gwf_piave_demand_config.toml @@ -71,7 +71,7 @@ glacier_store = "glacierstore" [state.vertical.demand.paddy.variables] h = "h_paddy" -[state.lateral.subsurface.flow.aquifer] +[state.lateral.subsurface.flow.aquifer.variables] head = "head" [input.vertical.glacier.parameters] @@ -203,7 +203,7 @@ q_av = "q_river" [output.vertical.soil.variables] zi = "zi" -[output.lateral.subsurface.flow.aquifer] +[output.lateral.subsurface.flow.aquifer.variables] head = "head" [csv]