fix(gwe-est): add support for energy decay in the solid phase

autotest/test_gwe_decay01.py

@@ -0,0 +1,184 @@
+"""
+Test problem for decay of energy in EST package of GWE.  Uses a single-cell
+model.  Test contributed by Cas Neyens.
+"""
+
+# Imports
+
+import flopy
+import numpy as np
+import pytest
+from framework import TestFramework
+
+# Base simulation and model name and workspace
+
+cases = ["decay-aqe", "decay-sld", "decay-both"]
+
+# Model units
+length_units = "meters"
+time_units = "seconds"
+
+rho_w = 1000  # kg/m^3
+rho_s = 2500  # kg/m^3
+n = 0.2  # -
+gamma_w = -1000  # J/s/m^3, arbitrary value for zero-order aqueous heat production
+gamma_s = -0.1  # J/s/kg, arbitrary value for zero-order solid heat production
+c_w = 4000  # J/kg/degC
+c_s = 1000  # J/kg/decC
+T0 = 0  # degC
+
+nrow = 1
+ncol = 1
+nlay = 1
+delr = 1  # m
+delc = 1  # m
+top = 1  # m
+botm = 0  # m
+
+perlen = 86400  # s
+nstp = 20
+
+parameters = {
+    # aqueous
+    "decay-aqe": {
+        "zero_order_decay_water": True,
+        "zero_order_decay_solid": False,
+        "decay_water": gamma_w,
+        "decay_solid": gamma_s,
+    },
+    # solid
+    "decay-sld": {
+        "zero_order_decay_water": False,
+        "zero_order_decay_solid": True,
+        "decay_water": gamma_w,
+        "decay_solid": gamma_s,
+    },
+    # combined
+    "decay-both": {
+        "zero_order_decay_water": True,
+        "zero_order_decay_solid": True,
+        "decay_water": gamma_w,
+        "decay_solid": gamma_s,
+    },
+}
+
+
+def temp_z_decay(t, rho_w, rho_s, c_w, c_s, gamma_w, gamma_s, n, T0):
+    t = np.atleast_1d(t)
+    coeff = (-gamma_w * n - gamma_s * (1 - n) * rho_s) / (
+        rho_w * c_w * n + rho_s * c_s * (1 - n)
+    )
+    return coeff * t + T0
+
+
+def build_models(idx, test):
+    # Base MF6 GWF model type
+    ws = test.workspace
+    name = cases[idx]
+    gwename = "gwe-" + name
+
+    decay_kwargs = parameters[name]
+
+    print(f"Building MF6 model...{name}")
+
+    sim = flopy.mf6.MFSimulation(
+        sim_name="heat",
+        sim_ws=ws,
+        exe_name="mf6",
+        version="mf6",
+    )
+    tdis = flopy.mf6.ModflowTdis(sim, nper=1, perioddata=[(perlen, nstp, 1.0)])
+    ims = flopy.mf6.ModflowIms(
+        sim, complexity="SIMPLE", inner_dvclose=0.001
+    )  # T can not become negative in this model
+
+    gwe = flopy.mf6.ModflowGwe(
+        sim,
+        modelname=gwename,
+        save_flows=True,
+        model_nam_file=f"{gwename}.nam",
+    )
+    dis = flopy.mf6.ModflowGwedis(
+        gwe, nrow=nrow, ncol=ncol, nlay=nlay, delr=delr, delc=delc, top=top, botm=botm
+    )
+    ic = flopy.mf6.ModflowGweic(gwe, strt=T0)
+    est = flopy.mf6.ModflowGweest(
+        gwe,
+        zero_order_decay_water=decay_kwargs["zero_order_decay_water"],
+        zero_order_decay_solid=decay_kwargs["zero_order_decay_solid"],
+        density_water=rho_w,
+        density_solid=rho_s,
+        heat_capacity_water=c_w,
+        heat_capacity_solid=c_s,
+        porosity=n,
+        decay_water=decay_kwargs["decay_water"],
+        decay_solid=decay_kwargs["decay_solid"],
+    )
+
+    oc = flopy.mf6.ModflowGweoc(
+        gwe,
+        budget_filerecord=f"{gwe.name}.bud",
+        temperature_filerecord=f"{gwe.name}.ucn",
+        printrecord=[("BUDGET", "ALL"), ("TEMPERATURE", "ALL")],
+        saverecord=[("BUDGET", "ALL"), ("TEMPERATURE", "ALL")],
+    )
+
+    return sim, None
+
+
+def check_output(idx, test):
+    print("evaluating results...")
+    msg = (
+        "Differences detected between the simulated results for zeroth-order "
+        "energy decay and the expected solution for decay specified in "
+    )
+    msg0 = msg + "the aqueous phase."
+    msg1 = msg + "the solid phase."
+    msg2 = msg + "both the aqueous and solid phases."
+
+    # read transport results from GWE model
+    name = cases[idx]
+    gwename = "gwe-" + name
+
+    # Get the MF6 temperature output
+    sim = test.sims[0]
+    gwe = sim.get_model(gwename)
+    temp_ts = gwe.output.temperature().get_ts((0, 0, 0))
+    t = temp_ts[:, 0]
+
+    temp_analy_w = temp_z_decay(
+        t, rho_w, rho_s, c_w, c_s, gamma_w, 0, n, T0
+    )  # aqueous decay only
+    temp_analy_s = temp_z_decay(
+        t, rho_w, rho_s, c_w, c_s, 0, gamma_s, n, T0
+    )  # aqueous decay only
+    temp_analy_ws = temp_z_decay(
+        t, rho_w, rho_s, c_w, c_s, gamma_w, gamma_s, n, T0
+    )  # aqueous + solid decay
+
+    print("temperature evaluation: " + str(temp_analy_w))
+
+    if "aqe" in name:
+        assert np.isclose(temp_ts[-1, 1], temp_analy_w[-1], atol=1e-10), msg0
+
+    if "sld" in name:
+        assert np.isclose(temp_ts[-1, 1], temp_analy_s[-1], atol=1e-10), msg1
+
+    if "both" in name:
+        assert np.isclose(temp_ts[-1, 1], temp_analy_ws[-1], atol=1e-10), msg2
+
+
+# - No need to change any code below
+@pytest.mark.parametrize(
+    "idx, name",
+    list(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()

doc/ReleaseNotes/develop.tex

@@ -23,6 +23,7 @@
 	\item Previously, GWF Model exchanges with other model types (GWT, GWE, PRT) verified that the flow model and the coupled model had the same number of active nodes, but did not check that the models' IDOMAIN arrays selected precisely the same set. Exchanges will now check that model IDOMAIN arrays are identical.
 	\item A regression was recently introduced into the PRT model's generalized tracking method, in which a coordinate transformation was carried out prematurely. This could result in incorrect particle positions in and near quad-refined cells. This bug has been fixed.
 	\item The PRT model previously allowed particles to be released at any time. Release times falling outside the bounds of the simulation's time discretization could produce undefined behavior. Any release times occurring before the simulation begins (i.e. negative times) will now be skipped with a warning message. If EXTEND\_TRACKING is not enabled, release times occurring after the end of the simulation will now be skipped with a warning message as well. If EXTEND\_TRACKING is enabled, release times after the end of the simulation are allowed.
+	\item Energy decay in the solid phase was added to the EST Package.  This capability is described in the Supplemental Technical Information document, but no actual support was provided in the source code.  Users can now specify distinct zeroth-order decay rates in either the aqueous or solid phases.
 \end{itemize}
 
 \underline{INTERNAL FLOW PACKAGES}

doc/mf6io/mf6ivar/dfn/gwe-est.dfn

@@ -9,12 +9,20 @@ longname save calculated flows to budget file
 description REPLACE save_flows {'{#1}': 'EST'}
 
 block options
-name zero_order_decay
+name zero_order_decay_water
 type keyword
 reader urword
 optional true
-longname activate zero-order decay
-description is a text keyword to indicate that zero-order decay will occur.  Use of this keyword requires that DECAY is specified in the GRIDDATA block.
+longname activate zero-order decay in aqueous phase
+description is a text keyword to indicate that zero-order decay will occur in the aqueous phase.  Use of this keyword requires that DECAY\_WATER is specified in the GRIDDATA block.
+
+block options
+name zero_order_decay_solid
+type keyword
+reader urword
+optional true
+longname activate zero-order decay in solid phase
+description is a text keyword to indicate that zero-order decay will occur in the solid phase.  Use of this keyword requires that DECAY\_SOLID is specified in the GRIDDATA block.
 
 block options
 name density_water
@@ -58,14 +66,24 @@ longname porosity
 description is the mobile domain porosity, defined as the mobile domain pore volume per mobile domain volume.  The GWE model does not support the concept of an immobile domain in the context of heat transport. 
 
 block griddata
-name decay
+name decay_water
 type double precision
 shape (nodes)
 reader readarray
 layered true
 optional true
 longname aqueous phase decay rate coefficient
-description is the rate coefficient for zero-order decay for the aqueous phase of the mobile domain.  A negative value indicates heat (energy) production.  The dimensions of decay for zero-order decay is energy per length cubed per time.  Zero-order decay will have no effect on simulation results unless zero-order decay is specified in the options block.
+description is the rate coefficient for zero-order decay for the aqueous phase of the mobile domain.  A negative value indicates heat (energy) production.  The dimensions of decay for zero-order decay in the aqueous phase is energy per length cubed per time.  Zero-order decay in the aqueous phase will have no effect on simulation results unless ZERO\_ORDER\_DECAY\_WATER is specified in the options block.
+
+block griddata
+name decay_solid
+type double precision
+shape (nodes)
+reader readarray
+layered true
+optional true
+longname solid phase decay rate coefficient
+description is the rate coefficient for zero-order decay for the solid phase.  A negative value indicates heat (energy) production.  The dimensions of decay for zero-order decay in the solid phase is energy per length cubed per time.  Zero-order decay in the solid phase will have no effect on simulation results unless ZERO\_ORDER\_DECAY\_SOLID is specified in the options block.
 
 block griddata
 name heat_capacity_solid