New sorption formulation

autotest/test_gwt_mst05.py

@@ -276,9 +276,9 @@ def check_output(idx, test):
         assert False, f'could not load data from "{fpth}"'
 
     cnorm = obs["X008"] / 0.05
-    cnorm_max = [0.32842034, 0.875391418]
+    cnorm_max = [0.324119806, 0.873678873]
     msg = f"{cnorm_max[idx]} /= {cnorm.max()}"
-    assert np.allclose(cnorm_max[idx], cnorm.max(), atol=0.001), msg
+    assert np.allclose(cnorm_max[idx], cnorm.max(), atol=1e-6), msg
 
     savefig = False
     if savefig:

autotest/test_gwt_mst08.py

@@ -0,0 +1,222 @@
+import os
+
+import flopy
+import numpy as np
+import pytest
+from framework import TestFramework
+from scipy import optimize
+
+cases = ["mst08a", "mst08b", "mst08c"]
+sorption_type = ["linear", "freundlich", "langmuir"]
+
+# Test description
+#
+# This test builds and runs a single-cell with Sorption and a Source.
+# The MST package is tested with several equilibrium isotherm.
+#
+# We inject a known total mass over a finite time and validate the
+# end-of-period concentration against the **analytical mass-balance**
+# solution: $$M_{total} = \theta V C - f_m \rho_b V \overline{C}$$.
+#
+# This solution doesn't have an exact closed form for all isotherms,
+# so we solve it algebraically using a root-finding method.
+
+# Parameters
+## Time parameters
+nper = 1
+perlen = [10.0]
+nstp = [10]
+tsmult = [1.0]
+
+## Solver parameters
+nouter, ninner = 100, 300
+hclose, rclose, relax = 1e-6, 1e-6, 1.0
+
+## Model domain and grid definition
+nlay, nrow, ncol = 1, 1, 1
+delr, delc = 1.0, 1.0
+top = 1.0
+botm = 0.0
+
+## Mass source parameters
+Msrc = 0.05  # mass source term (M/T)
+
+## MST parameters
+porosity = 0.25  # porosity (dimensionless)
+fm = 1  # fraction of mass in the solid phase (dimensionless)
+rho_bulk = 1700.0  # bulk density of the solid phase (M/L^3)
+
+## Sorption parameters
+### Linear sorption parameters
+Kd = 5.0e-4  # linear sorption coefficient (L^3/M)
+# Freundlich sorption parameters
+Kf = 1.0e-4  # Freundlich sorption coefficient (L^3M^(-1) )**n
+n_exp = 0.8  # Freundlich exponent (dimensionless)
+### Langmuir sorption parameters
+Kl = 5.0  # Langmuir sorption coefficient (L^3/M)
+Sbar = 1e-4  # Langmuir maximum solid phase concentration
+
+
+# Algebraic solution for equilibrium concentration
+# solved using a root-finding method
+def algebraic_solution(sorption_type):
+    Mtot = Msrc * perlen[0]  # total mass added over the period (M)
+    V = delr * delc * (top - botm)  # volume of the cell (L^3)
+
+    C0 = Mtot / (porosity * V)
+
+    # Find exact concentration at equilibrium for the listed isotherm
+    def C_equilibrium_linear():
+        def f(C):
+            return porosity * V * C + fm * rho_bulk * Kd * C * V - Mtot
+
+        def fprime(C):
+            return porosity * V + fm * rho_bulk * Kd * V
+
+        c_exact = optimize.newton(f, C0, fprime=fprime)
+        return c_exact
+
+    # Find exact concentration at equilibrium for the Freundlich  isotherm
+    def C_equilibrium_freundlich():
+        def f(C):
+            return porosity * V * C + fm * rho_bulk * Kf * C**n_exp * V - Mtot
+
+        def fprime(C):
+            return porosity * V + fm * rho_bulk * Kf * n_exp * C ** (n_exp - 1) * V
+
+        c_exact = optimize.newton(f, C0, fprime=fprime)
+        return c_exact
+
+    # Find exact concentration at equilibrium for the Langmuir isotherm
+    def C_equilibrium_langmuir():
+        def f(C):
+            return (
+                porosity * V * C
+                + fm * rho_bulk * (Kl * Sbar * C) / (1 + Kl * C) * V
+                - Mtot
+            )
+
+        def fprime(C):
+            return porosity * V + fm * rho_bulk * (Kl * Sbar) / ((1 + Kl * C) ** 2) * V
+
+        c_exact = optimize.newton(f, C0, fprime=fprime)
+        return c_exact
+
+    if sorption_type == "linear":
+        c_exact = C_equilibrium_linear()
+    elif sorption_type == "freundlich":
+        c_exact = C_equilibrium_freundlich()
+    elif sorption_type == "langmuir":
+        c_exact = C_equilibrium_langmuir()
+    else:
+        raise ValueError(f"Unknown sorption type: {sorption_type}")
+
+    return c_exact
+
+
+def build_models(idx, test):
+    name = cases[idx]
+    sorption = sorption_type[idx]
+
+    # build MODFLOW 6 files
+    ws = test.workspace
+    sim = flopy.mf6.MFSimulation(sim_name=name, sim_ws=ws, exe_name="mf6")
+
+    # create tdis package
+    tdis_rc = [perlen[0], nstp[0], tsmult[0]]
+    tdis = flopy.mf6.ModflowTdis(sim, time_units="DAYS", nper=nper, perioddata=tdis_rc)
+
+    # create gwt model
+    gwt = flopy.mf6.ModflowGwt(sim, modelname=name)
+
+    # create iterative model solution
+    ims = flopy.mf6.ModflowIms(
+        sim,
+        print_option="SUMMARY",
+        outer_dvclose=hclose,
+        outer_maximum=nouter,
+        inner_maximum=ninner,
+        inner_dvclose=hclose,
+        rcloserecord=rclose,
+        linear_acceleration="BICGSTAB",
+    )
+
+    # create discretization package
+    dis = flopy.mf6.ModflowGwtdis(
+        gwt, nlay=nlay, nrow=nrow, ncol=ncol, delr=delr, delc=delc, top=top, botm=botm
+    )
+
+    # initial conditions
+    ic = flopy.mf6.ModflowGwtic(gwt, strt=0.0)
+
+    # source term
+    spd = {0: [((0, 0, 0), Msrc)]}
+    src = flopy.mf6.ModflowGwtsrc(gwt, stress_period_data=spd)
+
+    # create mst package
+    if sorption == "linear":
+        mst = flopy.mf6.ModflowGwtmst(
+            gwt,
+            sorption="linear",
+            porosity=porosity,
+            bulk_density=rho_bulk,
+            distcoef=Kd,
+        )
+    elif sorption == "freundlich":
+        mst = flopy.mf6.ModflowGwtmst(
+            gwt,
+            sorption="freundlich",
+            porosity=porosity,
+            bulk_density=rho_bulk,
+            distcoef=Kf,
+            sp2=n_exp,
+        )
+    elif sorption == "langmuir":
+        mst = flopy.mf6.ModflowGwtmst(
+            gwt,
+            sorption="langmuir",
+            porosity=porosity,
+            bulk_density=rho_bulk,
+            distcoef=Kl,
+            sp2=Sbar,
+        )
+    else:
+        raise ValueError(f"Unknown sorption type: {sorption}")
+
+    # output control
+    oc = flopy.mf6.ModflowGwtoc(
+        gwt,
+        budget_filerecord=f"{name}.cbc",
+        concentration_filerecord=f"{name}.ucn",
+        saverecord=[("CONCENTRATION", "ALL")],
+        printrecord=[("CONCENTRATION", "LAST"), ("BUDGET", "LAST")],
+    )
+
+    return sim, None
+
+
+# Algebraic solution for equilibrium concentration
+
+
+def check_output(idx, test):
+    sorption = sorption_type[idx]
+    c_exact = algebraic_solution(sorption)
+
+    name = cases[idx]
+    fpth = os.path.join(test.workspace, f"{name}.ucn")
+    cobj = flopy.utils.HeadFile(fpth, precision="double", text="CONCENTRATION")
+    conc = cobj.get_data(idx=-1)
+
+    assert np.isclose(conc[0, 0, 0], c_exact, rtol=1e-6), "Concentrations do not match!"
+
+
+@pytest.mark.parametrize("idx, name", enumerate(cases))
+def test_mf6model(idx, name, function_tmpdir, targets):
+    test = TestFramework(
+        name=name,
+        workspace=function_tmpdir,
+        targets=targets,
+        build=lambda t: build_models(idx, t),
+        check=lambda t: check_output(idx, t),
+    )
+    test.run()

make/makefile

@@ -41,19 +41,20 @@ SOURCEDIR34=../src/Utilities/Export
 SOURCEDIR35=../src/Utilities/Idm
 SOURCEDIR36=../src/Utilities/Idm/mf6blockfile
 SOURCEDIR37=../src/Utilities/Idm/netcdf
-SOURCEDIR38=../src/Utilities/Libraries
-SOURCEDIR39=../src/Utilities/Libraries/blas
-SOURCEDIR40=../src/Utilities/Libraries/daglib
-SOURCEDIR41=../src/Utilities/Libraries/rcm
-SOURCEDIR42=../src/Utilities/Libraries/sparsekit
-SOURCEDIR43=../src/Utilities/Libraries/sparskit2
-SOURCEDIR44=../src/Utilities/Matrix
-SOURCEDIR45=../src/Utilities/Memory
-SOURCEDIR46=../src/Utilities/Observation
-SOURCEDIR47=../src/Utilities/OutputControl
-SOURCEDIR48=../src/Utilities/Performance
-SOURCEDIR49=../src/Utilities/TimeSeries
-SOURCEDIR50=../src/Utilities/Vector
+SOURCEDIR38=../src/Utilities/InterpolationScheme
+SOURCEDIR39=../src/Utilities/Libraries
+SOURCEDIR40=../src/Utilities/Libraries/blas
+SOURCEDIR41=../src/Utilities/Libraries/daglib
+SOURCEDIR42=../src/Utilities/Libraries/rcm
+SOURCEDIR43=../src/Utilities/Libraries/sparsekit
+SOURCEDIR44=../src/Utilities/Libraries/sparskit2
+SOURCEDIR45=../src/Utilities/Matrix
+SOURCEDIR46=../src/Utilities/Memory
+SOURCEDIR47=../src/Utilities/Observation
+SOURCEDIR48=../src/Utilities/OutputControl
+SOURCEDIR49=../src/Utilities/Performance
+SOURCEDIR50=../src/Utilities/TimeSeries
+SOURCEDIR51=../src/Utilities/Vector
 
 VPATH = \
 ${SOURCEDIR1} \
@@ -105,7 +106,8 @@ ${SOURCEDIR46} \
 ${SOURCEDIR47} \
 ${SOURCEDIR48} \
 ${SOURCEDIR49} \
-${SOURCEDIR50} 
+${SOURCEDIR50} \
+${SOURCEDIR51} 
 
 .SUFFIXES: .f90 .F90 .o
 

src/Model/GroundWaterTransport/gwt-mst.f90

@@ -340,14 +340,10 @@ subroutine mst_fc_srb(this, nodes, cold, nja, matrix_sln, idxglo, rhs, &
     real(DP) :: vcell
     real(DP) :: volfracm
     real(DP) :: rhobm
-    real(DP), dimension(nodes) :: c_half
-    real(DP) :: cbar_derv_half
-    real(DP) :: cbar_new, cbar_half, cbar_old
-    real(DP) :: sat_new, sat_half, sat_old
+    real(DP) :: sat_new, sat_old
     !
     ! -- set variables
     tled = DONE / delt
-    c_half = 0.5_DP * (cold + cnew)
     !
     ! -- loop through and calculate sorption contribution to hcof and rhs
     do n = 1, this%dis%nodes
@@ -362,23 +358,24 @@ subroutine mst_fc_srb(this, nodes, cold, nja, matrix_sln, idxglo, rhs, &
       sat_new = this%fmi%gwfsat(n)
       sat_old = this%fmi%gwfsatold(n, delt)
 
-      ! -- Midpoint formulation using average values
-      cbar_new = this%isotherm%value(cnew, n)
-      cbar_old = this%isotherm%value(cold, n)
-      cbar_half = 0.5 * (cbar_new + cbar_old)
-
-      sat_half = 0.5 * (sat_new + sat_old)
-      cbar_derv_half = this%isotherm%derivative(c_half, n)
-
-      hhcof = -volfracm * rhobm * cbar_derv_half * sat_half * Vcell * tled
+      ! -- Matrix contribution for sorption term
+      hhcof = -volfracm * rhobm * sat_new * this%isotherm%derivative(cnew, n) &
+              * Vcell * tled
       idiag = this%dis%con%ia(n)
       call matrix_sln%add_value_pos(idxglo(idiag), hhcof)
 
-      rrhs = -volfracm * rhobm * cbar_derv_half * sat_half * cold(n) &
-             * Vcell * tled
+      ! -- Right-hand side contribution due to linearization
+      rrhs = -volfracm * rhobm * sat_new * this%isotherm%derivative(cnew, n) &
+             * cnew(n) * Vcell * tled
       rhs(n) = rhs(n) + rrhs
 
-      rrhs = volfracm * rhobm * cbar_half * (sat_new - sat_old) * Vcell * tled
+      rrhs = volfracm * rhobm * sat_new * this%isotherm%value(cnew, n) * &
+             Vcell * tled
+      rhs(n) = rhs(n) + rrhs
+
+      ! -- Right-hand side contribution from previous time step
+      rrhs = -volfracm * rhobm * sat_old * this%isotherm%value(cold, n) * &
+             Vcell * tled
       rhs(n) = rhs(n) + rrhs
 
     end do
@@ -639,14 +636,11 @@ subroutine mst_cq_srb(this, nodes, cnew, cold, flowja)
     real(DP) :: vcell
     real(DP) :: volfracm
     real(DP) :: rhobm
-    real(DP), dimension(nodes) :: c_half
-    real(DP) :: cbar_derv_half
-    real(DP) :: cbar_new, cbar_half, cbar_old
-    real(DP) :: sat_new, sat_half, sat_old
+    real(DP) :: sat_new, sat_old
+    real(DP) :: contribution
     !
     ! -- initialize
     tled = DONE / delt
-    c_half = 0.5_DP * (cold + cnew)
     !
     ! -- Calculate sorption change
     do n = 1, nodes
@@ -666,20 +660,25 @@ subroutine mst_cq_srb(this, nodes, cnew, cold, flowja)
       sat_new = this%fmi%gwfsat(n)
       sat_old = this%fmi%gwfsatold(n, delt)
 
-      ! -- Midpoint formulation using average values
-      cbar_new = this%isotherm%value(cnew, n)
-      cbar_old = this%isotherm%value(cold, n)
-      cbar_half = 0.5 * (cbar_new + cbar_old)
+      ! -- Matrix contribution for sorption term
+      contribution = -volfracm * rhobm * sat_new &
+                     * this%isotherm%derivative(cnew, n) * cnew(n) * Vcell * tled
+      rate = rate + contribution
+
+      ! -- Right-hand side contribution due to linearization
+      ! -- Note: this contribution should cancel with the matrix contribution when the outer loop is converged
+      contribution = -volfracm * rhobm * sat_new * &
+                     this%isotherm%derivative(cnew, n) * cnew(n) * Vcell * tled
+      rate = rate - contribution
 
-      sat_half = 0.5 * (sat_new + sat_old)
-      cbar_derv_half = this%isotherm%derivative(c_half, n)
+      contribution = volfracm * rhobm * sat_new * this%isotherm%value(cnew, n) &
+                     * Vcell * tled
+      rate = rate - contribution
 
-      rate = -volfracm * rhobm * cbar_derv_half * sat_half * cnew(n) &
-             * Vcell * tled
-      rate = rate + volfracm * rhobm * cbar_derv_half * sat_half * cold(n) &
-             * Vcell * tled
-      rate = rate - volfracm * rhobm * cbar_half * (sat_new - sat_old) &
-             * Vcell * tled
+      ! -- Right-hand side contribution from previous time step
+      contribution = -volfracm * rhobm * sat_old * this%isotherm%value(cold, n) &
+                     * Vcell * tled
+      rate = rate - contribution
 
       this%ratesrb(n) = rate
       idiag = this%dis%con%ia(n)