Model type sbm_gwf: limit river flux to groundwater

docs/changelog.qmd

@@ -24,6 +24,8 @@ project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 - Wflow ZMQ server: allow JSON reading and writing of `NaN` and `Inf` values to avoid a JSON
   spec error. For example, during initialization of a wflow model, some (diagnostic) model
   variables are initialized with `NaN` values.
+- Limit flux from `River` boundary of `GroundwaterFlow` to unconfined aquifer based on
+  available river storage (from river flow routing).
 
 ### Changed
 - Removed vertical concepts `HBV` and `FLEXTopo`.

src/groundwater/aquifer.jl

@@ -420,7 +420,7 @@ function update!(gwf, Q, dt, conductivity_profile)
     Q .= 0.0  # TODO: Probably remove this when linking with other components
     flux!(Q, gwf.aquifer, gwf.connectivity, conductivity_profile)
     for boundary in gwf.boundaries
-        flux!(Q, boundary, gwf.aquifer)
+        flux!(Q, boundary, gwf.aquifer; dt)
     end
     gwf.aquifer.variables.head .+=
         (Q ./ gwf.aquifer.parameters.area .* dt ./ storativity(gwf.aquifer))

src/groundwater/boundary_conditions.jl

@@ -1,15 +1,15 @@
-function check_flux(flux, aquifer::UnconfinedAquifer, index::Int)
-    # Check if cell is dry
-    if aquifer.variables.head[index] <= aquifer.parameters.bottom[index]
-        # If cell is dry, no negative flux is allowed
-        return max(0, flux)
+function check_flux(flux, aquifer::UnconfinedAquifer, index::Int; dt = 1.0)
+    # Limit flux, not smaller than negative flux based on aquifer storage
+    if flux < 0.0
+        max_outflux = -aquifer.variables.storage[index] / dt
+        return max(flux, max_outflux)
     else
         return flux
     end
 end
 
 # Do nothing for a confined aquifer: aquifer can always provide flux
-check_flux(flux, aquifer::ConfinedAquifer, index::Int) = flux
+check_flux(flux, aquifer::ConfinedAquifer, index::Int; dt = 1.0) = flux
 
 @get_units @grid_loc @with_kw struct RiverParameters{T}
     infiltration_conductance::Vector{T} | "m2 d-1"
@@ -19,11 +19,16 @@ end
 
 @get_units @grid_loc @with_kw struct RiverVariables{T}
     stage::Vector{T} | "m"
+    storage::Vector{T} | "m3"
     flux::Vector{T} | "m3 d-1"
 end
 
 function RiverVariables(n)
-    variables = RiverVariables{Float}(; stage = fill(mv, n), flux = fill(mv, n))
+    variables = RiverVariables{Float}(;
+        stage = fill(mv, n),
+        storage = fill(mv, n),
+        flux = fill(mv, n),
+    )
     return variables
 end
 
@@ -51,19 +56,21 @@ function River(dataset, config, indices, index)
     return river
 end
 
-function flux!(Q, river::River, aquifer)
+function flux!(Q, river::River, aquifer; dt = 1.0)
     for (i, index) in enumerate(river.index)
         head = aquifer.variables.head[index]
         stage = river.variables.stage[i]
         if stage > head
             cond = river.parameters.infiltration_conductance[i]
             delta_head = min(stage - river.parameters.bottom[i], stage - head)
+            flux = min(cond * delta_head, river.variables.storage[i] / dt)
         else
             cond = river.parameters.exfiltration_conductance[i]
             delta_head = stage - head
+            flux = check_flux(cond * delta_head, aquifer, index; dt)
         end
-        river.variables.flux[i] = check_flux(cond * delta_head, aquifer, index)
-        Q[index] += river.variables.flux[i]
+        river.variables.flux[i] = flux
+        Q[index] += flux
     end
     return Q
 end
@@ -99,12 +106,12 @@ function Drainage(dataset, config, indices, index)
     return drains
 end
 
-function flux!(Q, drainage::Drainage, aquifer)
+function flux!(Q, drainage::Drainage, aquifer; dt = 1.0)
     for (i, index) in enumerate(drainage.index)
         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)
+        drainage.variables.flux[i] = check_flux(cond * delta_head, aquifer, index; dt)
         Q[index] += drainage.variables.flux[i]
     end
     return Q
@@ -125,11 +132,11 @@ end
     index::Vector{Int} | "-"
 end
 
-function flux!(Q, headboundary::HeadBoundary, aquifer)
+function flux!(Q, headboundary::HeadBoundary, aquifer; dt = 1.0)
     for (i, index) in enumerate(headboundary.index)
         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)
+        headboundary.variables.flux[i] = check_flux(cond * delta_head, aquifer, index; dt)
         Q[index] += headboundary.variables.flux[i]
     end
     return Q
@@ -151,12 +158,13 @@ function Recharge(rate, flux, index)
     return recharge
 end
 
-function flux!(Q, recharge::Recharge, aquifer)
+function flux!(Q, recharge::Recharge, aquifer; dt = 1.0)
     for (i, index) in enumerate(recharge.index)
         recharge.variables.flux[i] = check_flux(
             recharge.variables.rate[i] * aquifer.parameters.area[index],
             aquifer,
-            index,
+            index;
+            dt,
         )
         Q[index] += recharge.variables.flux[i]
     end
@@ -173,10 +181,10 @@ end
     index::Vector{Int} | "-"
 end
 
-function flux!(Q, well::Well, aquifer)
+function flux!(Q, well::Well, aquifer; dt = 1.0)
     for (i, index) in enumerate(well.index)
         well.variables.flux[i] =
-            check_flux(well.variables.volumetric_rate[i], aquifer, index)
+            check_flux(well.variables.volumetric_rate[i], aquifer, index; dt)
         Q[index] += well.variables.flux[i]
     end
     return Q

src/sbm_gwf_model.jl

@@ -434,9 +434,11 @@ function update!(model::AbstractModel{<:SbmGwfModel})
 
     update!(land, routing, network, config, dt)
 
-    # set river stage (groundwater) to average h from kinematic wave
+    # set river stage and storage (groundwater boundary) based on river flow routing
+    # variables
     boundaries.river.variables.stage .=
         routing.river_flow.variables.h_av .+ boundaries.river.parameters.bottom
+    boundaries.river.variables.storage .= routing.river_flow.variables.storage
 
     # determine stable time step for groundwater flow
     conductivity_profile = get(config.model, "conductivity_profile", "uniform")

test/groundwater.jl

@@ -300,13 +300,18 @@ end
                 exfiltration_conductance = [200.0, 200.0],
                 bottom = [1.0, 1.0],
             )
-            variables = Wflow.RiverVariables(; stage = [2.0, 2.0], flux = [0.0, 0.0])
+            variables = Wflow.RiverVariables(;
+                stage = [2.0, 2.0],
+                storage = [20.0, 20.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
-            @test Q[1] == (2.0 - 1.0) * 100.0
-            # drainage, flux is () * exf_cond
+            # infiltration, below bottom, flux is (stage - bottom) * inf_cond, limited by
+            # river storage (20.0)
+            @test Q[1] == 20.0
+            # drainage, flux is (stage - head) * exf_cond
             @test Q[3] == (2.0 - 20.0) * 200.0
         end