Quick boundary fix (#430)

* Quick boundary fix * Unpack `aquifer` and `constanthead` for readability * Update changelog --------- Co-authored-by: Willem van Verseveld <[email protected]>
Deltares · Jun 11, 2024 · 4becf7b · 4becf7b
1 parent daa96df
commit 4becf7b
Show file tree

Hide file tree

Showing 2 changed files with 13 additions and 9 deletions.
diff --git a/docs/src/changelog.md b/docs/src/changelog.md
@@ -38,6 +38,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   `Polyester.jl`.
 - Fixed required states of the model type `sbm_gwf`: added `h_av` for the river and land
   domain.
+- For the `constanthead` boundary of `GroundwaterFlow` the `head` should not be changed (was
+  set to `top` elevation of the `aquifer` if `head` > `top`), and `exfiltwater` should be 0
+  for these boundary cells.
 
 ### Changed
 - Stop exposing scalar variables through BMI. The `BMI.get_value_ptr` function was not

diff --git a/src/sbm_gwf_model.jl b/src/sbm_gwf_model.jl
@@ -483,6 +483,8 @@ function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfModel}
     @unpack lateral, vertical, network, clock, config = model
 
     inds_riv = network.index_river
+    aquifer = lateral.subsurface.flow.aquifer
+    constanthead = lateral.subsurface.flow.constanthead
 
     # extract water levels h_av [m] from the land and river domains
     # this is used to limit open water evaporation
@@ -512,7 +514,7 @@ function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfModel}
     # determine stable time step for groundwater flow
     conductivity_profile =
         get(config.input.lateral.subsurface, "conductivity_profile", "uniform")
-    dt_gw = stable_timestep(lateral.subsurface.flow.aquifer, conductivity_profile) # time step in day (Float64)
+    dt_gw = stable_timestep(aquifer, conductivity_profile) # time step in day (Float64)
     dt_sbm = (vertical.dt / tosecond(basetimestep)) # vertical.dt is in seconds (Float64)
     if dt_gw < dt_sbm
         @warn(
@@ -529,18 +531,17 @@ function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfModel}
 
     # determine excess water depth [m] (exfiltwater) in groundwater domain (head > surface)
     # and reset head
-    exfiltwater =
-        (
-            lateral.subsurface.flow.aquifer.head .-
-            min.(lateral.subsurface.flow.aquifer.head, lateral.subsurface.flow.aquifer.top)
-        ) .* storativity(lateral.subsurface.flow.aquifer)
-    lateral.subsurface.flow.aquifer.head .=
-        min.(lateral.subsurface.flow.aquifer.head, lateral.subsurface.flow.aquifer.top)
+    exfiltwater = (aquifer.head .- min.(aquifer.head, aquifer.top)) .* storativity(aquifer)
+    aquifer.head .= min.(aquifer.head, aquifer.top)
+
+    # Adjust for constant head boundary of groundwater domain
+    exfiltwater[constanthead.index] .= 0
+    aquifer.head[constanthead.index] .= constanthead.head
 
     # update vertical sbm concept (runoff, ustorelayerdepth and satwaterdepth)
     update_after_subsurfaceflow(
         vertical,
-        (network.land.altitude .- lateral.subsurface.flow.aquifer.head) .* 1000.0, # zi [mm] in vertical concept SBM
+        (network.land.altitude .- aquifer.head) .* 1000.0, # zi [mm] in vertical concept SBM
         exfiltwater .* 1000.0,
     )