diff --git a/autotest/test_gwt_ist02.py b/autotest/test_gwt_ist02.py index 21f8029c299..20fa00d0b4e 100644 --- a/autotest/test_gwt_ist02.py +++ b/autotest/test_gwt_ist02.py @@ -249,7 +249,7 @@ def build_models(idx, test): cim_filerecord = f"{gwtname}.ist.ucn" ist = flopy.mf6.ModflowGwtist( gwt, - sorption=True, + sorption="LINEAR", save_flows=True, cim_filerecord=cim_filerecord, cim=0.0, diff --git a/autotest/test_gwt_ist03.py b/autotest/test_gwt_ist03.py new file mode 100644 index 00000000000..d4389ab4908 --- /dev/null +++ b/autotest/test_gwt_ist03.py @@ -0,0 +1,384 @@ +""" +Test the immobile domain isotherms (Freundlich and Langmuir) in the IST Package +using a one-dimensional flow problem. Compare the simulated sorbate +concentrations output from the model with sorbate concentrations calculated +in this script using python. Problem is adapted from test_gwt_ist02.py. +""" + +import os +import pathlib as pl + +import flopy +import numpy as np +import pytest +from framework import TestFramework + +cases = ["ist03a", "ist03b", "ist03c"] +sorption_idx = ["LINEAR", "FREUNDLICH", "LANGMUIR"] +distcoef_idx = [0.1, 0.1, 0.1] +sp2_idx = [None, 0.7, 0.7] + +nlay, nrow, ncol = 1, 1, 300 + + +def build_models(idx, test): + sorption = sorption_idx[idx] + distcoef = distcoef_idx[idx] + sp2 = sp2_idx[idx] + + perlen = [20.0, 30.0] + nper = len(perlen) + nstp = [100, 100] + tsmult = [1.0, 1.0] + delr = 0.5 + delc = 1.0 + top = 1.0 + botm = [0.0] + strt = 9.5 + hk = 1000.0 + + nouter, ninner = 100, 300 + hclose, rclose, relax = 1e-6, 1e-6, 0.97 + + tdis_rc = [] + for id in range(nper): + tdis_rc.append((perlen[id], nstp[id], tsmult[id])) + + name = cases[idx] + + # build MODFLOW 6 files + ws = test.workspace + sim = flopy.mf6.MFSimulation( + sim_name=name, version="mf6", exe_name="mf6", sim_ws=ws + ) + # create tdis package + tdis = flopy.mf6.ModflowTdis( + sim, time_units="DAYS", nper=nper, perioddata=tdis_rc + ) + + # create gwf model + gwfname = "gwf_" + name + newtonoptions = None + gwf = flopy.mf6.ModflowGwf( + sim, + modelname=gwfname, + newtonoptions=newtonoptions, + ) + + # create iterative model solution and register the gwf model with it + imsgwf = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="DBD", + under_relaxation_theta=0.7, + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + filename=f"{gwfname}.ims", + ) + sim.register_ims_package(imsgwf, [gwf.name]) + + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + ) + + # initial conditions + ic = flopy.mf6.ModflowGwfic(gwf, strt=strt) + + # node property flow + npf = flopy.mf6.ModflowGwfnpf( + gwf, save_flows=False, icelltype=0, k=hk, k33=hk + ) + + # chd files + chdspdict = { + 0: [[(0, 0, 0), 10.0], [(0, 0, ncol - 1), 9.0]], + } + chd = flopy.mf6.ModflowGwfchd( + gwf, + print_input=True, + print_flows=True, + stress_period_data=chdspdict, + save_flows=False, + pname="CHD-1", + ) + + # output control + oc = flopy.mf6.ModflowGwfoc( + gwf, + budget_filerecord=f"{gwfname}.cbc", + head_filerecord=f"{gwfname}.hds", + headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], + saverecord=[("HEAD", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) + + # create gwt model + gwtname = "gwt_" + name + gwt = flopy.mf6.ModflowGwt(sim, modelname=gwtname, save_flows=True) + + # create iterative model solution and register the gwt model with it + imsgwt = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="NONE", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + filename=f"{gwtname}.ims", + ) + sim.register_ims_package(imsgwt, [gwt.name]) + + dis = flopy.mf6.ModflowGwtdis( + gwt, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + filename=f"{gwtname}.dis", + ) + + # initial conditions + ic = flopy.mf6.ModflowGwtic(gwt, strt=0.0) + + # advection + adv = flopy.mf6.ModflowGwtadv(gwt, scheme="UPSTREAM") + + # dispersion + xt3d_off = True + dsp = flopy.mf6.ModflowGwtdsp( + gwt, + xt3d_off=xt3d_off, + diffc=0.0, + alh=0.5, + ath1=0.0, + ) + + # mass storage and transfer + bulk_density = 1.6 + thetam = 0.14 + thetaim = 0.12 + volfracm = thetam / (thetam + thetaim) + volfracim = 1.0 - volfracm + porositym = thetam / volfracm + porosityim = thetaim / volfracim + mst = flopy.mf6.ModflowGwtmst( + gwt, + sorbate_filerecord=f"{gwtname}.mst.csrb", + sorption=sorption, + porosity=porositym, + bulk_density=bulk_density, + distcoef=distcoef, + sp2=sp2, + ) + + # immobile storage and transfer + ist = flopy.mf6.ModflowGwtist( + gwt, + sorption=sorption, + save_flows=True, + cim_filerecord=f"{gwtname}.ist.ucn", + sorbate_filerecord=f"{gwtname}.ist.csrb", + cim=0.0, + porosity=porosityim, + volfrac=volfracim, + bulk_density=bulk_density, + zetaim=0.1, + distcoef=distcoef, + sp2=sp2, + ) + + # cnc + cncspdict = { + 0: [[(0, 0, 0), 1.0]], + 1: [[(0, 0, 0), 0.0]], + } + cnc = flopy.mf6.ModflowGwtcnc( + gwt, + print_input=True, + print_flows=True, + maxbound=1, + stress_period_data=cncspdict, + save_flows=False, + pname="CNC-1", + ) + + # sources + sourcerecarray = [()] + ssm = flopy.mf6.ModflowGwtssm(gwt, sources=sourcerecarray) + + # output control + oc = flopy.mf6.ModflowGwtoc( + gwt, + budget_filerecord=f"{gwtname}.cbc", + concentration_filerecord=f"{gwtname}.ucn", + concentrationprintrecord=[ + ("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL") + ], + saverecord=[("CONCENTRATION", "ALL"), ("BUDGET", "ALL")], + printrecord=[("CONCENTRATION", "ALL"), ("BUDGET", "ALL")], + ) + + obs_data = { + f"{gwtname}.obs.csv": [ + ("(1-1-300)", "CONCENTRATION", (0, 0, ncol - 1)), + ], + } + obs_package = flopy.mf6.ModflowUtlobs( + gwt, + # pname="conc_obs", + filename=f"{gwtname}.obs", + digits=10, + print_input=True, + continuous=obs_data, + ) + + # GWF GWT exchange + gwfgwt = flopy.mf6.ModflowGwfgwt( + sim, + exgtype="GWF6-GWT6", + exgmnamea=gwfname, + exgmnameb=gwtname, + filename=f"{name}.gwfgwt", + ) + + return sim, None + + +def make_plot(sim, plot_title): + print("making plots...") + name = sim.name + ws = sim.workspace + sim = flopy.mf6.MFSimulation.load(sim_ws=ws) + gwfname = "gwf_" + name + gwtname = "gwt_" + name + gwf = sim.get_model(gwfname) + gwt = sim.get_model(gwtname) + + output = gwt.obs.output + obs_names = output.obs_names + data = output.obs(f=obs_names[0]).data + + # mobile aqueous sorbed concentrations + fpth = pl.Path(ws) / f"{gwtname}.mst.csrb" + cobj = flopy.utils.HeadFile(fpth, precision="double", text="SORBATE") + csrb = cobj.get_alldata().reshape(200, 300) + + # immobile aqueous + fpth = pl.Path(ws) / f"{gwtname}.ist.ucn" + cobj = flopy.utils.HeadFile(fpth, precision="double", text="CIM") + cim = cobj.get_alldata().reshape(200, 300) + + # immobile sorbed + fpth = pl.Path(ws) / f"{gwtname}.ist.csrb" + cobj = flopy.utils.HeadFile(fpth, precision="double", text="SORBATE") + cimsrb = cobj.get_alldata().reshape(200, 300) + + import matplotlib.pyplot as plt + + fig = plt.figure(figsize=(6, 3)) + ax = fig.add_subplot(1, 1, 1) + ax.plot(data["totim"], data["(1-1-300)"], "k-", label="mobile aqueous") + ax.plot(data["totim"], csrb[:, 299], "b-", label="mobile sorbate") + ax.plot(data["totim"], cim[:, 299], "k--", label="immobile aqueous") + ax.plot(data["totim"], cimsrb[:, 299], "b--", label="immobile sorbate") + plt.xlabel("time, in days") + plt.ylabel("concentration") + plt.legend() + plt.title(plot_title) + fname = os.path.join(ws, gwtname + ".png") + plt.savefig(fname) + + +def check_output(idx, test): + sorption = sorption_idx[idx] + distcoef = distcoef_idx[idx] + sp2 = sp2_idx[idx] + + makeplot = False + if makeplot: + plot_title = sorption + make_plot(test, plot_title) + + name = test.name + gwtname = "gwt_" + name + gwfname = "gwf_" + name + + # load the observed concentrations in column 300 + fname = os.path.join(test.workspace, gwtname + ".obs.csv") + assert os.path.isfile(fname), f"file not found: {fname}" + simvals = np.genfromtxt(fname, names=True, delimiter=",", deletechars="") + + sim = test.sims[0] + gwt = sim.gwt[0] + mst = gwt.mst + ist = gwt.ist + conc = gwt.output.concentration().get_alldata().reshape(200, 300) + csrb = mst.output.sorbate().get_alldata().reshape(200, 300) + cim = ist.output.cim().get_alldata().reshape(200, 300) + cimsrb = ist.output.sorbate().get_alldata().reshape(200, 300) + + # check conc and csrb + if sorption == "LINEAR": + csrb_answer = np.where(conc > 0, distcoef * conc, 0) + if sorption == "FREUNDLICH": + csrb_answer = np.where(conc > 0, distcoef * conc**sp2, 0) + if sorption == "LANGMUIR": + csrb_answer = np.where( + conc > 0, distcoef * sp2 * conc / (1 + distcoef * conc), 0 + ) + if not np.allclose(csrb[:, 1:], csrb_answer[:, 1:]): + diff = csrb - csrb_answer + print("min and max difference") + print(diff) + print(diff.min(), diff.max()) + assert False, "csrb not consistent with known answer" + + # check cim and cimsrb + if sorption == "LINEAR": + cimsrb_answer = np.where(cim > 0, distcoef * cim, 0) + if sorption == "FREUNDLICH": + cimsrb_answer = np.where(cim > 0, distcoef * cim**sp2, 0) + if sorption == "LANGMUIR": + cimsrb_answer = np.where( + cim > 0, distcoef * sp2 * cim / (1 + distcoef * cim), 0 + ) + if not np.allclose(cimsrb[:, 1:], cimsrb_answer[:, 1:]): + diff = cimsrb - cimsrb_answer + print("min and max difference") + print(diff.min(), diff.max()) + assert False, "cimsrb not consistent with known answer" + + +@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() diff --git a/autotest/test_gwt_mt3dms_p01.py b/autotest/test_gwt_mt3dms_p01.py index 1bdcea8e491..06ab7ffd30c 100644 --- a/autotest/test_gwt_mt3dms_p01.py +++ b/autotest/test_gwt_mt3dms_p01.py @@ -417,7 +417,7 @@ def p01mf6( if zeta is not None: ist = flopy.mf6.ModflowGwtist( gwt, - sorption=True, + sorption="LINEAR", first_order_decay=first_order_decay, zero_order_decay=zero_order_decay, bulk_density=rhob, diff --git a/doc/ReleaseNotes/develop.tex b/doc/ReleaseNotes/develop.tex index 6230295a3ae..5ffdd68834d 100644 --- a/doc/ReleaseNotes/develop.tex +++ b/doc/ReleaseNotes/develop.tex @@ -6,15 +6,16 @@ \underline{NEW FUNCTIONALITY} \begin{itemize} \item A new adaptive time stepping (ATS) capability was added to the Advection (ADV) Package of the Groundwater Transport (GWT) Model. A new input option, called ATS\_PERCEL, specifies the fractional cell distance that a particle of water can travel within one time step. When ATS\_PERCEL is specified by the user, and the ATS utility is activated in the TDIS Package, the ADV Package will calculate the largest time step that will meet this fractional cell distance constraint, and will submit this time step to the ATS utility. This option may improve time stepping for solute transport models and for variable-density flow and transport models by allowing step lengths to be calculated as a function of the flow system rather than being specified as input by the user. - \item Added the capability to write sorbate concentrations to a binary output file. A new SORBATE option is now available in the Mass Storage and Transfer (MST) Package of the GWT Model to provide the name of the binary output file for sorbate concentrations. Sorbate concentrations will be written to the binary output file whenever concentrations for the GWT model are saved, as determined by the GWT Output Control option. + \item Added the capability to write sorbate concentrations to binary output files. A new SORBATE option is now available in the Mobile Storage and Transfer (MST) Package and the Immobile Storage and Transfer (IST) Package of the GWT Model to provide the name of the binary output file for sorbate concentrations. Sorbate concentrations will be written to the binary output file whenever concentrations for the GWT model are saved, as determined by the GWT Output Control option. \item Add kinematic-wave routing option for the Streamflow Routing (SFR) Package. Prior to this change, the SFR Package could only simulate unidirectional, steady, uniform flow. With kinematic-wave routing, unidirectional waves can now propagate through the SFR network by explicitly including a storage term in the reach continuity equation. The kinematic-wave routing option is based on the ``TRANSROUTE'' option available in the SFR Package in MODFLOW-NWT. The kinematic-wave routing option is enabled by specifying ``STORAGE'' in the SFR Package OPTIONS block. - % \item xxx + \item The Freundlich and Langmuir isotherms can now be used to represent sorption in the immobile domain when simulating solute transport using the Groundwater Transport (GWT) Model. Prior to this change, only linear sorption could be used with the Immobile Storage and Transfer (IST) Package for the GWT Model. With this change, the SORPTION keyword specified in the OPTIONS block of the IST Package must be followed by LINEAR, FREUNDLICH, or LANGMUIR. If the FREUNDLICH or LANGMUIR isotherms are specified, then the user must also include the SP2 array in the GRIDDATA block of the IST Package. This change breaks backward compatibility for GWT Models that used the SORPTION option in the IST Package, which previously did not require the type of sorption to be specified after the SORPTION keyword. Details on the implementation of the nonlinear Freundlich and Langmuir sorption isotherms are described in a new chapter in the MODFLOW 6 --- Supplemental Technical Information guide, which is included with the MODFLOW 6 distribution. \end{itemize} \underline{EXAMPLES} \begin{itemize} \item A new Toth example was added to show how the classic groundwater problem consisting of local, intermediate, and regional flow paths can be simulated with MODFLOW. \item A new Streamflow Routing (SFR) Package example based on the Pinder-Sauer problem \citep{pinder1971numerical} modified by \cite{lal2001modification} was added to demonstrate the capabilities of the kinematic-wave approximation option. + \item A new Barends example was added to demonstrate the use of the Groundwater Energy (GWE) Model for simulating heat transport in a groundwater reservoir overlain by a low-permeability overburden. % \item xxx \end{itemize} diff --git a/doc/ReleaseNotes/previous/v6.2.1.tex b/doc/ReleaseNotes/previous/v6.2.1.tex index e9fb9ac0350..ca914259ab6 100644 --- a/doc/ReleaseNotes/previous/v6.2.1.tex +++ b/doc/ReleaseNotes/previous/v6.2.1.tex @@ -27,7 +27,7 @@ \item Budget terms for the Immobile domain Storage and Transfer (IST) Package were not being written to the binary budget file for the GWT Model. The package was modified to write these rate terms to the GWT binary budget file using the settings specified in the GWT Output Control file. \item Bulk density does not need to be specified for the Immobile domain Storage and Transfer (IST) Package if sorption is not active; however the program was trying to access bulk density even though it is not needed, which resulted in an access violation. Program was fixed so that bulk density does not need to be specified by the user unless sorption is active for the IST Package. \item Budget tables printed to the listing file had numeric values that were missing the `E` character if the exponent had three digits (e.g. 1.e-100 or 1.e100). Writing of the budget table was modified to include the `E` character in this case. This change should make it easier for programs written in other languages to parse these tables. - \item In the Mass Transfer and Storage (MST) and the Immobile Storage and Transfer (IST) Packages, the keyword to activate sorption was changed from SORBTION to SORPTION. The program will still accept SORBTION, but this keyword will be deprecated in the future. + \item In the Mobile Transfer and Storage (MST) and the Immobile Storage and Transfer (IST) Packages, the keyword to activate sorption was changed from SORBTION to SORPTION. The program will still accept SORBTION, but this keyword will be deprecated in the future. \item Revised several of the text strings written to the headers within the GWT binary budget file. A table of possible text strings for the GWT binary budget file are now included in the mf6io.pdf document. \end{itemize} diff --git a/doc/SuppTechInfo/body.tex b/doc/SuppTechInfo/body.tex index f49c7261456..19e14c3bb2c 100644 --- a/doc/SuppTechInfo/body.tex +++ b/doc/SuppTechInfo/body.tex @@ -65,6 +65,12 @@ \customlabel{ch:prtmodel}{\thechapno} \input{particle-tracking-model.tex} +\newpage +\incchap +\SECTION{Chapter \thechapno. Nonlinear Sorption for the Immobile Domain} +\customlabel{ch:istsorption}{\thechapno} +\input{ist-sorption.tex} + \newpage \ifx\usgsdirector\undefined \addcontentsline{toc}{section}{\hspace{1.5em}\bibname} diff --git a/doc/SuppTechInfo/ist-sorption.tex b/doc/SuppTechInfo/ist-sorption.tex new file mode 100644 index 00000000000..a5270e338d9 --- /dev/null +++ b/doc/SuppTechInfo/ist-sorption.tex @@ -0,0 +1,97 @@ +The Groundwater Transport (GWT) Model \citep{modflow6gwt} in \mf was designed to simulate a range of solute-transport processes including the sorption and desorption of solute within the mobile and immobile domains. For the mobile domain, sorption can be represented using a linear isotherm as well as the nonlinear Freundlich and Langmuir isotherms. The type of sorption isotherm used for the mobile domain is specified by the user in the input file for the Mobile Storage and Transfer (MST) Package. For the immobile domain, which is activated using the Immobile Storage and Transfer (IST) Package, the GWT model was limited to representing linear isotherms only. The program checked to make sure that if sorption was represented in the immobile domain that sorption in the mobile domain was represented using a linear isotherm. This was the behavior of \mf up to and including version 6.5.0. + +As reported by \mf users, some problems may require use of nonlinear sorption isotherms for mobile-immobile domain simulations. In response to this need, the IST Package was extended to support the nonlinear Freundlich and Langmuir isotherms, in addition to the linear isotherm. The purpose of this chapter is to describe the mathematical approach that was used to implement these nonlinear isotherms in the IST Package. The \mf program still requires use of consistent sorption approaches in both the mobile and immobile domains, but with these changes, any of the three sorption isotherms can now be used as long as the same type of isotherm is used for both the mobile and immobile domains. + +Chapter \ref{ch:sorption} of this document presents a revised parameterization of the mobile and immobile domains. Included in Chapter 9 is the revised form of the partial differential equation governing solute mass within an immobile domain (equation \ref{eqn:gwtistpde}). Equation~\ref{eqn:gwtistpde} is reproduced here as, + +\begin{equation} +\label{eqn:gwtistpde2} +\begin{split} +\theta_{im} \frac{\partial C_{im} }{\partial t} + \hat{f}_{im} \rho_{b,im} \frac{\partial \overline{C}_{im}}{\partial t} = +- \lambda_{1,im} \theta_{im} C_{im} - \lambda_{2,im} \hat{f}_{im} \rho_{b,im} \overline{C}_{im} \\ +- \gamma_{1,im} \theta_{im} - \gamma_{2,im} \hat{f}_{im} \rho_{b,im} ++ \zeta_{im} S_w \left ( C - C_{im} \right ), +\end{split} +\end{equation} + +\noindent where +$\theta_{im}$ is the effective porosity in the immobile domain defined as volume of voids participating in immobile-domain transport per unit volume of immobile domain $im$ ($L^3/L^3$), +$C_{im}$ is an immobile-domain volumetric concentration of solute expressed as mass of dissolved solute per unit volume of fluid in immobile domain $im$ ($M/L^3$), +$t$ is time ($T$), +$\hat{f}_{im}$ is the volume fraction of the immobile domain defined as the volume of mobile domain $im$ per volume of aquifer ($L^3/L^3$), +$\rho_{b,im}$ is the bulk density of aquifer material in the immobile domain defined as mass of solid aquifer material per unit volume of immobile domain $im$ ($M/L^3$), +$\overline{C}_{im}$ is the concentration of sorbate (sorbed solute) expressed as mass of sorbate per unit mass of solid aquifer material in the immobile domain ($M/M$), +$\lambda_{1,im}$ is the first-order reaction rate coefficient for dissolved solute in the immobile domain ($1/T$), +$\lambda_{2,im}$ is the first-order reaction rate coefficient for sorbate in the immobile domain ($1/T$), +$\gamma_{1,im}$ is the zero-order reaction rate coefficient for dissolved solute in the immobile domain ($ML^{-3}T^{-1}$), +$\gamma_{2,im}$ is the zero-order reaction rate coefficient for sorbate in the immobile domain ($M M^{-1}T^{-1}$), +$\zeta_{im}$ is the rate coefficient for the transfer of mass between the mobile domain and immobile domain $im$ ($1/T$), +$S_w$ is the water saturation defined as the volume of water per volume of voids ($L^3/L^3$), +and $C$ is the mobile-domain volumetric concentration of solute expressed as mass of dissolved solute per unit volume of mobile-domain fluid ($M/L^3$). + +The GWT Model documentation report \citep{modflow6gwt} describes the approach for including the effects of an immobile domain on solute transport. The approach is based on the method described in \cite{zheng2002} and involves discretizing the solute balance equation for the immobile domain. \cite{modflow6gwt} present a form of the discretized balance equation based on a linear isotherm and the original parameterization, which has since been revised according to the description in Chapter 9. The following discretized form of equation~\ref{eqn:gwtistpde2} includes a more general representation of the sorption isotherm, in order to support both linear and nonlinear expressions, as well as the updated parameterization described in Chapter \ref{ch:sorption}: + +\begin{equation} +\label{eqn:imdomain2} +\begin{split} +\frac{\theta_{im} V_{cell}}{\Delta t} C_{im}^{t + \Delta t} - +\frac{\theta_{im} V_{cell}}{\Delta t} C_{im}^t ++ \frac{\hat{f}_{im} \rho_{b,im} V_{cell}}{\Delta t} \overline{K}_{im}^{t + \Delta t} C_{im}^{t + \Delta t} +- \frac{\hat{f}_{im} \rho_{b,im} V_{cell}}{\Delta t} \overline{K}_{im}^t C_{im}^{t} \\ += \\ +- \lambda_{1,im} \theta_{im} V_{cell} C_{im}^{t + \Delta t} +- \lambda_{2,im} \hat{f}_{im} V_{cell} \rho_{b,im} \overline{K}_{im}^{t + \Delta t} C_{im}^{t + \Delta t} \\ +- \gamma_{1,im} \theta_{im} V_{cell} +- \gamma_{2,im} \hat{f}_{im} \rho_{b,im} V_{cell} \\ ++ V_{cell} S_w^{t + \Delta t} \zeta_{im} C^{t + \Delta t} +- V_{cell} S_w^{t + \Delta t} \zeta_{im} C_{im}^{t + \Delta t} , +\end{split} +\end{equation} + +\noindent where $V_{cell}$ is the volume of the model cell, and $\overline{K}_{im}$ is a linearized and effective distribution coefficient used to relate the immobile-domain sorbate concentration to the dissolved immobile-domain concentration using $\overline{C}_{im} = \overline{K}_{im} C_{im}$. + +The $\overline{K}_{im}$ term in equation~\ref{eqn:imdomain2} can be determined from the mathematical expression used for the different isotherms. Mathematical expressions for the linear, Freundlich, and Langmuir isotherms, respectively, written in terms of immobile domain concentrations with subscript $im$, are + +\begin{equation} +\label{eqn:linear} +\overline{C}_{im} = K_d C_{im}, +\end{equation} + +\begin{equation} +\label{eqn:freundlich} +\overline{C}_{im} = K_f C_{im}^a, +\end{equation} + +\noindent and + +\begin{equation} +\label{eqn:langmuir} +\overline{C}_{im} = \frac{K_l \overline{S} C_{im}}{1 + K_l C_{im}}, +\end{equation} + +\noindent where $K_d$ is the linear distribution coefficient ($L^3/M$), $K_f$ is the Freundlich constant $(L^3 / M)^a$, $a$ is the dimensionless Freundlich exponent, $K_l$ is the Langmuir constant $(L^3 / M)$ and $\overline{S}$ is the total concentration of sorption sites available $(M/M)$. + +It can be seen from equations \ref{eqn:linear} to \ref{eqn:langmuir} that expressions for linearized distribution coefficients can be written as, + +\begin{equation} +\label{eqn:linearkd} +\overline{K}_{im} = K_d, +\end{equation} + +\begin{equation} +\label{eqn:freundlichkd} +\overline{K}_{im} = K_f C_{im}^{a - 1}, +\end{equation} + +\noindent and + +\begin{equation} +\label{eqn:langmuirkd} +\overline{K}_{im} = \frac{K_l \overline{S}}{1 + K_l C_{im}}, +\end{equation} + +\noindent for the linear, Freundlich, and Langmuir isotherms, respectively. + +Unlike the linear isotherm, the Freundlich and Langmuir isotherms have a nonlinear relation between $\overline{C}_{im}$ and $C_{im}$. In order to calculate a value for $\overline{K}_{im}$ at the $t + \Delta t$ time level in equation~\ref{eqn:imdomain2}, the value for $C_{im}^{t + \Delta t}$ from the previous outer iteration is used in equations~\ref{eqn:freundlichkd} and \ref{eqn:langmuirkd}. Use of a previous iterate in these nonlinear calculations for $\overline{K}_{im}$ may cause convergence problems for some applications, compared to use of the linear expression, which does not depend on $C_{im}$. + +The remaining details of the implementation are identical to those described by \cite{modflow6gwt}. Equation~\ref{eqn:imdomain2} is solved for $C_{im}^{t + \Delta t}$, which is then substituted into the equation for the transfer of mass between the mobile and immobile domains. The resulting terms are then added to the system of equations and solved numerically. \ No newline at end of file diff --git a/doc/SuppTechInfo/mf6suptechinfo.bbl b/doc/SuppTechInfo/mf6suptechinfo.bbl index c94fb8b58c9..cfe85d410cf 100644 --- a/doc/SuppTechInfo/mf6suptechinfo.bbl +++ b/doc/SuppTechInfo/mf6suptechinfo.bbl @@ -1,4 +1,4 @@ -\begin{thebibliography}{30} +\begin{thebibliography}{31} \providecommand{\natexlab}[1]{#1} \expandafter\ifx\csname urlstyle\endcsname\relax \providecommand{\doiagency}[1]{doi:\discretionary{}{}{}#1}\else @@ -191,6 +191,10 @@ Zhang, Y., King, M.J., and Datta-Gupta, A., 2012, Robust streamline tracing Zheng, C., 2010, MT3DMS v5.3, Supplemental User's Guide: {Technical Report Prepared for the U.S. Army Corps of Engineers, 51 p.} +\bibitem[{Zheng and Bennett(2002)}]{zheng2002} +Zheng, C., and Bennett, G.D., 2002, Applied contaminant transport modeling, 2nd + edition: New York, New York, Wiley-Interscience, 560~p. + \bibitem[{Zheng and Wang(1999)}]{zheng1999mt3dms} Zheng, C., and Wang, P.P., 1999, MT3DMS---A modular three-dimensional multi-species transport model for simulation of advection, dispersion and diff --git a/doc/SuppTechInfo/sorption.tex b/doc/SuppTechInfo/sorption.tex index 43fb517d0f1..a50e7e5371c 100644 --- a/doc/SuppTechInfo/sorption.tex +++ b/doc/SuppTechInfo/sorption.tex @@ -25,7 +25,7 @@ \subsection{Review of the Original Parameterization} \label{sec:origparamreview} \end{aligned} \end{equation} -\noindent where $S_w$ is the water saturation defined as the volume of water per volume of voids ($L^3/L^3$), $C$ is the mobile-domain volumetric concentration of solute expressed as mass of dissolved solute per unit volume of mobile-domain fluid ($M/L^3$), $t$ is time ($T$), $\matr{q}$ is the vector of specific discharge ($L/T$), $\matr{D}$ is the second-order tensor of hydrodynamic dispersion coefficients ($L^2/T$), $q'_s$ is the volumetric flow rate per unit volume of aquifer (defined as positive for flow into the aquifer) for mass sources and sinks ($1/T$), $C_s$ is the volumetric solute concentration of the source or sink fluid ($M/L^3$), $M_s$ is rate of solute mass loading per unit volume of aquifer ($M/L^3T$), $\lambda_1$ is the first-order decay rate coefficient for dissolved solute in the mobile domain ($1/T$), $\gamma_1$ is the zero-order decay rate coefficient for dissolved solute in the mobile domain ($M/L^3T$), $\overline{C}$ is the mass-fraction concentration of sorbate (sorbed solute) expressed as mass of sorbate per unit mass of solid aquifer material in the mobile domain ($M/M$), $\lambda_2$ is the first-order decay rate coefficient ($1/T$) for sorbate in the mobile domain, $\gamma_2$ is the zero-order decay rate coefficient ($M/MT$) for sorbate in the mobile domain, $nim$ is the number of immobile domains, $\zeta_{im}$ is the rate coefficient for the transfer of mass between the mobile domain and immobile domain $im$ ($1/T$), $C_{im}$ is an immobile-domain volumetric concentration of solute expressed as mass of dissolved solute per unit volume of fluid in immobile domain $im$ ($M/L^3$), $\rho_{b}$ is the bulk density of the aquifer material defined as the mass of solid aquifer material per unit volume of aquifer ($M/L^3$), $\theta_m$ is the mobile-domain effective porosity defined as defined as volume of voids participating in mobile-domain transport per unit volume of aquifer ($L^3/L^3$), and $f_m$ is the fraction of the mass of aquifer solid material that is in the mobile domain ($M/M$). To avoid potential confusion with the total porosity, in equation~\ref{eqn:gwtpdeorig} the mobile porosity is intentionally given the symbol $\theta_m$, which is different than the symbol $\theta$ used by \cite{modflow6gwt}. Also, note that \cite{modflow6gwt} define $f_m$ as ``the fraction of aquifer solid material available for sorptive exchange with the mobile phase under fully saturated conditions," which is correct only if all the aquifer solid material in the mobile domain is available for sorptive exchange and the fraction is understood to be a mass fraction. +\noindent where $S_w$ is the water saturation defined as the volume of water per volume of voids ($L^3/L^3$), $C$ is the mobile-domain volumetric concentration of solute expressed as mass of dissolved solute per unit volume of mobile-domain fluid ($M/L^3$), $t$ is time ($T$), $\matr{q}$ is the vector of specific discharge ($L/T$), $\matr{D}$ is the second-order tensor of hydrodynamic dispersion coefficients ($L^2/T$), $q'_s$ is the volumetric flow rate per unit volume of aquifer (defined as positive for flow into the aquifer) for mass sources and sinks ($1/T$), $C_s$ is the volumetric solute concentration of the source or sink fluid ($M/L^3$), $M_s$ is rate of solute mass loading per unit volume of aquifer ($M/L^3T$), $\lambda_1$ is the first-order decay rate coefficient for dissolved solute in the mobile domain ($1/T$), $\gamma_1$ is the zero-order decay rate coefficient for dissolved solute in the mobile domain ($M/L^3T$), $\overline{C}$ is the concentration of sorbate (sorbed solute) expressed as mass of sorbate per unit mass of solid aquifer material in the mobile domain ($M/M$), $\lambda_2$ is the first-order decay rate coefficient ($1/T$) for sorbate in the mobile domain, $\gamma_2$ is the zero-order decay rate coefficient ($M/MT$) for sorbate in the mobile domain, $nim$ is the number of immobile domains, $\zeta_{im}$ is the rate coefficient for the transfer of mass between the mobile domain and immobile domain $im$ ($1/T$), $C_{im}$ is an immobile-domain volumetric concentration of solute expressed as mass of dissolved solute per unit volume of fluid in immobile domain $im$ ($M/L^3$), $\rho_{b}$ is the bulk density of the aquifer material defined as the mass of solid aquifer material per unit volume of aquifer ($M/L^3$), $\theta_m$ is the mobile-domain effective porosity defined as defined as volume of voids participating in mobile-domain transport per unit volume of aquifer ($L^3/L^3$), and $f_m$ is the fraction of the mass of aquifer solid material that is in the mobile domain ($M/M$). To avoid potential confusion with the total porosity, in equation~\ref{eqn:gwtpdeorig} the mobile porosity is intentionally given the symbol $\theta_m$, which is different than the symbol $\theta$ used by \cite{modflow6gwt}. Also, note that \cite{modflow6gwt} define $f_m$ as ``the fraction of aquifer solid material available for sorptive exchange with the mobile phase under fully saturated conditions," which is correct only if all the aquifer solid material in the mobile domain is available for sorptive exchange and the fraction is understood to be a mass fraction. The Immobile Storage and Transfer (IST) Package for the GWT Model allows users to designate a fraction of a model cell as immobile. With the original parameterization used in versions of \mf up to and including version 6.4.1, the partial differential equation that represents the conservation of solute mass at points in an immobile domain is equation 7--2 of \cite{modflow6gwt}, which is reproduced here as @@ -39,7 +39,7 @@ \subsection{Review of the Original Parameterization} \label{sec:origparamreview} \end{split} \end{equation} -\noindent where $\overline{C}_{im}$ is the mass-fraction concentration of sorbate (sorbed solute) expressed as mass of sorbate per unit mass of solid aquifer material in the immobile domain ($M/M$), $\lambda_{1,im}$ is the first-order reaction rate coefficient for dissolved solute in the immobile domain ($1/T$), $\lambda_{2,im}$ is the first-order reaction rate coefficient for sorbate in the immobile domain ($1/T$), $\gamma_{1,im}$ is the zero-order reaction rate coefficient for dissolved solute in the immobile domain ($ML^{-3}T^{-1}$), $\gamma_{2,im}$ is the zero-order reaction rate coefficient for sorbate in the immobile domain ($M M^{-1}T^{-1}$), $\theta_{im}$ is the effective porosity in the immobile domain defined as defined as volume of voids participating in immobile-domain transport per unit volume of immobile domain $im$ ($L^3/L^3$), and $f_{im}$ is the fraction of the mass of aquifer solid material that is in immobile domain $im$ ($M/M$). Note that \cite{modflow6gwt} define $f_{im}$ as ``the fraction of aquifer solid material available for sorptive exchange with the immobile domain under fully saturated conditions," which is correct only if all the aquifer solid material in the immobile domain is available for sorptive exchange and the fraction is understood to be a mass fraction. +\noindent where $\overline{C}_{im}$ is the concentration of sorbate (sorbed solute) expressed as mass of sorbate per unit mass of solid aquifer material in the immobile domain ($M/M$), $\lambda_{1,im}$ is the first-order reaction rate coefficient for dissolved solute in the immobile domain ($1/T$), $\lambda_{2,im}$ is the first-order reaction rate coefficient for sorbate in the immobile domain ($1/T$), $\gamma_{1,im}$ is the zero-order reaction rate coefficient for dissolved solute in the immobile domain ($ML^{-3}T^{-1}$), $\gamma_{2,im}$ is the zero-order reaction rate coefficient for sorbate in the immobile domain ($M M^{-1}T^{-1}$), $\theta_{im}$ is the effective porosity in the immobile domain defined as volume of voids participating in immobile-domain transport per unit volume of immobile domain $im$ ($L^3/L^3$), and $f_{im}$ is the fraction of the mass of aquifer solid material that is in immobile domain $im$ ($M/M$). Note that \cite{modflow6gwt} define $f_{im}$ as ``the fraction of aquifer solid material available for sorptive exchange with the immobile domain under fully saturated conditions," which is correct only if all the aquifer solid material in the immobile domain is available for sorptive exchange and the fraction is understood to be a mass fraction. The original model parameters that are affected by the revised parameterization are listed in table~\ref{table:origparam}. Note that in the original parameterization, division of the aquifer into mobile and immobile domains is conceptualized in terms solid mass fractions, and porosities and bulk densities are defined on a per-aquifer-volume basis. Note also that the user is required (in versions of \mf up to and including version 6.4.1) to provide the value of the overall bulk density, $\rho_b$, in the input for each package that represents a domain in which sorption is active (the MST Package for the mobile domain, and the IST Package for immobile domains). diff --git a/doc/mf6io/mf6ivar/dfn/gwt-ist.dfn b/doc/mf6io/mf6ivar/dfn/gwt-ist.dfn index c39d95cc26e..77b10f01c06 100644 --- a/doc/mf6io/mf6ivar/dfn/gwt-ist.dfn +++ b/doc/mf6io/mf6ivar/dfn/gwt-ist.dfn @@ -87,11 +87,12 @@ description name of the comma-separated value (CSV) output file to write budget block options name sorption -type keyword +type string +valid linear freundlich langmuir reader urword optional true longname activate sorption -description is a text keyword to indicate that sorption will be activated. Use of this keyword requires that BULK\_DENSITY and DISTCOEF are specified in the GRIDDATA block. The linear sorption isotherm is the only isotherm presently supported in the IST Package. +description is a text keyword to indicate that sorption will be activated. Valid sorption options include LINEAR, FREUNDLICH, and LANGMUIR. Use of this keyword requires that BULK\_DENSITY and DISTCOEF are specified in the GRIDDATA block. If sorption is specified as FREUNDLICH or LANGMUIR then SP2 is also required in the GRIDDATA block. block options name first_order_decay @@ -217,6 +218,39 @@ optional false longname write format description write format can be EXPONENTIAL, FIXED, GENERAL, or SCIENTIFIC. +block options +name sorbate_filerecord +type record sorbate fileout sorbatefile +shape +reader urword +tagged true +optional true +longname +description + +block options +name sorbate +type keyword +shape +in_record true +reader urword +tagged true +optional false +longname sorbate keyword +description keyword to specify that record corresponds to immobile sorbate concentration. + +block options +name sorbatefile +type string +preserve_case true +shape +in_record true +reader urword +tagged false +optional false +longname file keyword +description name of the output file to write immobile sorbate concentration information. Immobile sorbate concentrations will be written whenever aqueous immobile concentrations are saved, as determined by settings in the Output Control option. + # --------------------- gwt ist griddata --------------------- block griddata @@ -296,3 +330,12 @@ layered true longname distribution coefficient description is the distribution coefficient for the equilibrium-controlled linear sorption isotherm in dimensions of length cubed per mass. distcoef is not required unless the SORPTION keyword is specified in the options block. If the SORPTION keyword is not specified in the options block, distcoef will have no effect on simulation results. +block griddata +name sp2 +type double precision +shape (nodes) +reader readarray +layered true +optional true +longname second sorption parameter +description is the exponent for the Freundlich isotherm and the sorption capacity for the Langmuir isotherm. sp2 is not required unless the SORPTION keyword is specified in the options block and sorption is specified as FREUNDLICH or LANGMUIR. If the SORPTION keyword is not specified in the options block, or if sorption is specified as LINEAR, sp2 will have no effect on simulation results. diff --git a/doc/mf6io/mf6ivar/dfn/gwt-mst.dfn b/doc/mf6io/mf6ivar/dfn/gwt-mst.dfn index dbaaf5c6fe9..bc0bd1709bb 100644 --- a/doc/mf6io/mf6ivar/dfn/gwt-mst.dfn +++ b/doc/mf6io/mf6ivar/dfn/gwt-mst.dfn @@ -126,7 +126,7 @@ reader readarray layered true optional true longname distribution coefficient -description is the distribution coefficient for the equilibrium-controlled linear sorption isotherm in dimensions of length cubed per mass. distcoef is not required unless the SORPTION keyword is specified. +description is the distribution coefficient for the equilibrium-controlled linear sorption isotherm in dimensions of length cubed per mass. If the Freunchlich isotherm is specified, then discoef is the Freundlich constant. If the Langmuir isotherm is specified, then distcoef is the Langmuir constant. distcoef is not required unless the SORPTION keyword is specified. block griddata name sp2 @@ -136,5 +136,5 @@ reader readarray layered true optional true longname second sorption parameter -description is the exponent for the Freundlich isotherm and the sorption capacity for the Langmuir isotherm. +description is the exponent for the Freundlich isotherm and the sorption capacity for the Langmuir isotherm. sp2 is not required unless the SORPTION keyword is specified in the options block. If the SORPTION keyword is not specified in the options block, sp2 will have no effect on simulation results. diff --git a/src/Model/GroundWaterTransport/gwt-ist.f90 b/src/Model/GroundWaterTransport/gwt-ist.f90 index a1ed8154800..0164f5c4a7b 100644 --- a/src/Model/GroundWaterTransport/gwt-ist.f90 +++ b/src/Model/GroundWaterTransport/gwt-ist.f90 @@ -14,13 +14,16 @@ module GwtIstModule use KindModule, only: DP, I4B - use ConstantsModule, only: DONE, DZERO, LENFTYPE, & + use ConstantsModule, only: DONE, DZERO, DHALF, LENFTYPE, & LENPACKAGENAME, & LENBUDTXT, DHNOFLO use BndModule, only: BndType use BudgetModule, only: BudgetType use TspFmiModule, only: TspFmiType - use GwtMstModule, only: GwtMstType, get_zero_order_decay + use GwtMstModule, only: GwtMstType, get_zero_order_decay, & + get_freundlich_conc, get_freundlich_derivative, & + get_langmuir_conc, get_langmuir_derivative, & + get_freundlich_kd, get_langmuir_kd use OutputControlDataModule, only: OutputControlDataType use MatrixBaseModule ! @@ -57,17 +60,20 @@ module GwtIstModule integer(I4B), pointer :: icimout => null() !< unit number for binary cim output integer(I4B), pointer :: ibudgetout => null() !< binary budget output file integer(I4B), pointer :: ibudcsv => null() !< unit number for csv budget output file + integer(I4B), pointer :: ioutsorbate => null() !< unit number for sorbate concentration output integer(I4B), pointer :: idcy => null() !< order of decay rate (0:none, 1:first, 2:zero) integer(I4B), pointer :: isrb => null() !< sorption active flag (0:off, 1:on); only linear is supported in ist integer(I4B), pointer :: kiter => null() !< picard iteration counter real(DP), pointer, contiguous :: cim(:) => null() !< concentration for immobile domain real(DP), pointer, contiguous :: cimnew(:) => null() !< immobile concentration at end of current time step real(DP), pointer, contiguous :: cimold(:) => null() !< immobile concentration at end of last time step + real(DP), pointer, contiguous :: cimsrb(:) => null() !< sorbate concentration in immobile domain real(DP), pointer, contiguous :: zetaim(:) => null() !< mass transfer rate to immobile domain real(DP), pointer, contiguous :: porosity(:) => null() !< immobile domain porosity defined as volume of immobile voids per volume of immobile domain real(DP), pointer, contiguous :: volfrac(:) => null() !< volume fraction of the immobile domain defined as volume of immobile domain per aquifer volume real(DP), pointer, contiguous :: bulk_density(:) => null() !< bulk density of immobile domain defined as mass of solids in immobile domain per volume of immobile domain real(DP), pointer, contiguous :: distcoef(:) => null() !< distribution coefficient + real(DP), pointer, contiguous :: sp2(:) => null() !< second sorption parameter real(DP), pointer, contiguous :: decay(:) => null() !< first or zero order rate constant for liquid real(DP), pointer, contiguous :: decaylast(:) => null() !< decay rate used for last iteration (needed for zero order decay) real(DP), pointer, contiguous :: decayslast(:) => null() !< sorbed decay rate used for last iteration (needed for zero order decay) @@ -85,12 +91,15 @@ module GwtIstModule procedure :: bnd_bd => ist_bd procedure :: bnd_ot_model_flows => ist_ot_model_flows procedure :: bnd_ot_dv => ist_ot_dv + procedure :: output_immobile_concentration + procedure :: output_immobile_sorbate_concentration procedure :: bnd_ot_bdsummary => ist_ot_bdsummary procedure :: bnd_da => ist_da procedure :: allocate_scalars procedure :: read_dimensions => ist_read_dimensions procedure :: read_options procedure :: get_thetaim + procedure :: ist_calc_csrb procedure, private :: ist_allocate_arrays procedure, private :: read_data @@ -197,7 +206,8 @@ subroutine ist_ar(this) if (this%isrb /= this%mst%isrb) then call store_error('Sorption is active for the IST Package but it is not & &compatible with the sorption option selected for the MST Package. & - &If sorption is active for the IST Package, then SORPTION LINEAR must & + &If sorption is active for the IST Package, then the same type of & + &sorption (LINEAR, FREUNDLICH, or LANGMUIR) must & &be specified in the options block of the MST Package.') end if if (count_errors() > 0) then @@ -259,20 +269,17 @@ subroutine ist_ad(this) end subroutine ist_ad !> @ brief Fill coefficient method for package - !! - !! Method to calculate and fill coefficients for the package. - !! !< subroutine ist_fc(this, rhs, ia, idxglo, matrix_sln) - ! -- modules + ! modules use TdisModule, only: delt - ! -- dummy + ! dummy class(GwtIstType) :: this !< GwtIstType object real(DP), dimension(:), intent(inout) :: rhs !< right-hand side vector for model integer(I4B), dimension(:), intent(in) :: ia !< solution CRS row pointers integer(I4B), dimension(:), intent(in) :: idxglo !< mapping vector for model (local) to solution (global) class(MatrixBaseType), pointer :: matrix_sln !< solution coefficient matrix - ! -- local + ! local integer(I4B) :: n, idiag real(DP) :: tled real(DP) :: hhcof, rrhs @@ -280,7 +287,8 @@ subroutine ist_fc(this, rhs, ia, idxglo, matrix_sln) real(DP) :: vcell real(DP) :: thetaim real(DP) :: zetaim - real(DP) :: kd + real(DP) :: kdnew + real(DP) :: kdold real(DP) :: volfracim real(DP) :: rhobim real(DP) :: lambda1im @@ -292,37 +300,40 @@ subroutine ist_fc(this, rhs, ia, idxglo, matrix_sln) real(DP) :: cimsrbold real(DP) :: cimsrbnew real(DP), dimension(9) :: ddterm - ! - ! -- set variables + + ! set variables tled = DONE / delt this%kiter = this%kiter + 1 - ! - ! -- loop through and calculate immobile domain contribution to hcof and rhs + + ! loop through each node and calculate immobile domain contribution + ! to hcof and rhs do n = 1, this%dis%nodes - ! - ! -- skip if transport inactive + + ! skip if transport inactive if (this%ibound(n) <= 0) cycle - ! - ! -- calculate new and old water volumes + + ! calculate new and old water volumes vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n)) swtpdt = this%fmi%gwfsat(n) swt = this%fmi%gwfsatold(n, delt) thetaim = this%get_thetaim(n) idiag = ia(n) - ! - ! -- set exchange coefficient + + ! set exchange coefficient zetaim = this%zetaim(n) - ! - ! -- Add dual domain mass transfer contributions to rhs and hcof - kd = DZERO + + ! Add dual domain mass transfer contributions to rhs and hcof + ! dcimsrbdc = DZERO + kdnew = DZERO + kdold = DZERO volfracim = DZERO rhobim = DZERO lambda1im = DZERO lambda2im = DZERO gamma1im = DZERO gamma2im = DZERO - ! - ! -- setup decay variables + + ! set variables for decay of aqueous solute if (this%idcy == 1) lambda1im = this%decay(n) if (this%idcy == 2) then gamma1im = get_zero_order_decay(this%decay(n), this%decaylast(n), & @@ -330,16 +341,48 @@ subroutine ist_fc(this, rhs, ia, idxglo, matrix_sln) this%cimnew(n), delt) this%decaylast(n) = gamma1im end if - ! - ! -- setup sorption variables + + ! setup sorption variables if (this%isrb > 0) then - kd = this%distcoef(n) + + ! initialize sorption variables volfracim = this%volfrac(n) rhobim = this%bulk_density(n) - if (this%idcy == 1) lambda2im = this%decay_sorbed(n) - if (this%idcy == 2) then - cimsrbold = this%cimold(n) * kd - cimsrbnew = this%cimnew(n) * kd + + ! set isotherm dependent sorption variables + select case (this%isrb) + case (1) + ! linear + kdnew = this%distcoef(n) + kdold = this%distcoef(n) + cimsrbnew = this%cimnew(n) * kdnew + cimsrbold = this%cimold(n) * kdold + case (2) + ! freundlich + kdnew = get_freundlich_kd(this%cimnew(n), this%distcoef(n), & + this%sp2(n)) + kdold = get_freundlich_kd(this%cimold(n), this%distcoef(n), & + this%sp2(n)) + cimsrbnew = get_freundlich_conc(this%cimnew(n), this%distcoef(n), & + this%sp2(n)) + cimsrbold = get_freundlich_conc(this%cimold(n), this%distcoef(n), & + this%sp2(n)) + case (3) + ! langmuir + kdnew = get_langmuir_kd(this%cimnew(n), this%distcoef(n), & + this%sp2(n)) + kdold = get_langmuir_kd(this%cimold(n), this%distcoef(n), & + this%sp2(n)) + cimsrbnew = get_langmuir_conc(this%cimnew(n), this%distcoef(n), & + this%sp2(n)) + cimsrbold = get_langmuir_conc(this%cimold(n), this%distcoef(n), & + this%sp2(n)) + end select + + ! set decay of sorbed solute parameters + if (this%idcy == 1) then + lambda2im = this%decay_sorbed(n) + else if (this%idcy == 2) then gamma2im = get_zero_order_decay(this%decay_sorbed(n), & this%decayslast(n), & this%kiter, cimsrbold, & @@ -347,39 +390,34 @@ subroutine ist_fc(this, rhs, ia, idxglo, matrix_sln) this%decayslast(n) = gamma2im end if end if - ! - ! -- calculate the terms and then get the hcof and rhs contributions + + ! calculate dual domain terms and get the hcof and rhs contributions call get_ddterm(thetaim, vcell, delt, swtpdt, & - volfracim, rhobim, kd, lambda1im, lambda2im, & + volfracim, rhobim, kdnew, kdold, lambda1im, lambda2im, & gamma1im, gamma2im, zetaim, ddterm, f) cimold = this%cimold(n) call get_hcofrhs(ddterm, f, cimold, hhcof, rrhs) - ! - ! -- update solution accumulators + + ! update solution accumulators call matrix_sln%add_value_pos(idxglo(idiag), hhcof) rhs(n) = rhs(n) + rrhs - ! + end do - ! - ! -- Return - return + end subroutine ist_fc !> @ brief Calculate package flows. - !! - !! Calculate the flow between connected package control volumes. - !! !< subroutine ist_cq(this, x, flowja, iadv) - ! -- modules + ! modules use TdisModule, only: delt use ConstantsModule, only: DZERO - ! -- dummy + ! dummy class(GwtIstType), intent(inout) :: this !< GwtIstType object real(DP), dimension(:), intent(in) :: x !< current dependent-variable value real(DP), dimension(:), contiguous, intent(inout) :: flowja !< flow between two connected control volumes integer(I4B), optional, intent(in) :: iadv !< flag that indicates if this is an advance package - ! -- local + ! local integer(I4B) :: idiag integer(I4B) :: n real(DP) :: rate @@ -388,7 +426,8 @@ subroutine ist_cq(this, x, flowja, iadv) real(DP) :: vcell real(DP) :: thetaim real(DP) :: zetaim - real(DP) :: kd + real(DP) :: kdnew + real(DP) :: kdold real(DP) :: volfracim real(DP) :: rhobim real(DP) :: lambda1im @@ -401,33 +440,34 @@ subroutine ist_cq(this, x, flowja, iadv) real(DP) :: cimsrbold real(DP) :: cimsrbnew real(DP), dimension(9) :: ddterm - ! -- formats - ! - ! -- initialize + ! formats + + ! initialize this%budterm(:, :) = DZERO - ! - ! -- Calculate immobile domain transfer rate + + ! Calculate immobile domain transfer rate do n = 1, this%dis%nodes - ! - ! -- skip if transport inactive + + ! skip if transport inactive rate = DZERO cimnew = DZERO if (this%ibound(n) > 0) then - ! - ! -- calculate new and old water volumes + + ! calculate new and old water volumes vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n)) swtpdt = this%fmi%gwfsat(n) swt = this%fmi%gwfsatold(n, delt) thetaim = this%get_thetaim(n) - ! - ! -- set exchange coefficient + + ! set exchange coefficient zetaim = this%zetaim(n) - ! - ! -- Calculate exchange with immobile domain + + ! Calculate exchange with immobile domain rate = DZERO hhcof = DZERO rrhs = DZERO - kd = DZERO + kdnew = DZERO + kdold = DZERO volfracim = DZERO rhobim = DZERO lambda1im = DZERO @@ -439,51 +479,119 @@ subroutine ist_cq(this, x, flowja, iadv) gamma1im = get_zero_order_decay(this%decay(n), this%decaylast(n), 0, & this%cimold(n), this%cimnew(n), delt) end if + + ! setup sorption variables if (this%isrb > 0) then - kd = this%distcoef(n) + + ! initialize sorption variables volfracim = this%volfrac(n) rhobim = this%bulk_density(n) - if (this%idcy == 1) lambda2im = this%decay_sorbed(n) - if (this%idcy == 2) then - cimsrbold = this%cimold(n) * kd - cimsrbnew = this%cimnew(n) * kd + + ! set isotherm dependent sorption variables + select case (this%isrb) + case (1) + ! linear + kdnew = this%distcoef(n) + kdold = this%distcoef(n) + cimsrbnew = this%cimnew(n) * kdnew + cimsrbold = this%cimold(n) * kdold + case (2) + ! freundlich + kdnew = get_freundlich_kd(this%cimnew(n), this%distcoef(n), & + this%sp2(n)) + kdold = get_freundlich_kd(this%cimold(n), this%distcoef(n), & + this%sp2(n)) + cimsrbnew = get_freundlich_conc(this%cimnew(n), this%distcoef(n), & + this%sp2(n)) + cimsrbold = get_freundlich_conc(this%cimold(n), this%distcoef(n), & + this%sp2(n)) + case (3) + ! langmuir + kdnew = get_langmuir_kd(this%cimnew(n), this%distcoef(n), & + this%sp2(n)) + kdold = get_langmuir_kd(this%cimold(n), this%distcoef(n), & + this%sp2(n)) + cimsrbnew = get_langmuir_conc(this%cimnew(n), this%distcoef(n), & + this%sp2(n)) + cimsrbold = get_langmuir_conc(this%cimold(n), this%distcoef(n), & + this%sp2(n)) + end select + + ! set decay of sorbed solute parameters + if (this%idcy == 1) then + lambda2im = this%decay_sorbed(n) + else if (this%idcy == 2) then gamma2im = get_zero_order_decay(this%decay_sorbed(n), & this%decayslast(n), & 0, cimsrbold, & cimsrbnew, delt) end if end if - ! - ! -- calculate the terms and then get the hcof and rhs contributions + + ! calculate the terms and then get the hcof and rhs contributions call get_ddterm(thetaim, vcell, delt, swtpdt, & - volfracim, rhobim, kd, lambda1im, lambda2im, & + volfracim, rhobim, kdnew, kdold, lambda1im, lambda2im, & gamma1im, gamma2im, zetaim, ddterm, f) cimold = this%cimold(n) call get_hcofrhs(ddterm, f, cimold, hhcof, rrhs) - ! - ! -- calculate rate from hcof and rhs + + ! calculate rate from hcof and rhs rate = hhcof * x(n) - rrhs - ! - ! -- calculate immobile domain concentration + + ! calculate immobile domain concentration cimnew = get_ddconc(ddterm, f, cimold, x(n)) - ! - ! -- accumulate the budget terms + + ! accumulate the budget terms call accumulate_budterm(this%budterm, ddterm, cimnew, cimold, x(n), & this%idcy) end if - ! - ! -- store rate and add to flowja + + ! store rate and add to flowja this%strg(n) = rate idiag = this%dis%con%ia(n) flowja(idiag) = flowja(idiag) + rate - ! - ! -- store immobile domain concentration + + ! store immobile domain concentration this%cimnew(n) = cimnew - ! + end do - return + + ! calculate csrb + if (this%isrb /= 0) then + call this%ist_calc_csrb(this%cimnew) + end if + end subroutine ist_cq + !> @ brief Calculate immobile sorbed concentration + !< + subroutine ist_calc_csrb(this, cim) + ! -- dummy + class(GwtIstType) :: this !< GwtMstType object + real(DP), intent(in), dimension(:) :: cim !< immobile domain aqueous concentration at end of this time step + ! -- local + integer(I4B) :: n + real(DP) :: distcoef + real(DP) :: csrb + + ! Calculate sorbed concentration + do n = 1, size(cim) + csrb = DZERO + if (this%ibound(n) > 0) then + distcoef = this%distcoef(n) + if (this%isrb == 1) then + csrb = cim(n) * distcoef + else if (this%isrb == 2) then + csrb = get_freundlich_conc(cim(n), distcoef, this%sp2(n)) + else if (this%isrb == 3) then + csrb = get_langmuir_conc(cim(n), distcoef, this%sp2(n)) + end if + end if + this%cimsrb(n) = csrb + end do + + end subroutine ist_calc_csrb + !> @ brief Add package flows to model budget. !! !! Add the flow between IST package and the model (ratin and ratout) to the @@ -578,23 +686,34 @@ subroutine ist_ot_model_flows(this, icbcfl, ibudfl, icbcun, imap) return end subroutine ist_ot_model_flows - !> @ brief Output immobile domain concentration. - !! + !> @ brief Output dependent variables. !< subroutine ist_ot_dv(this, idvsave, idvprint) - ! -- modules + ! dummy variables + class(GwtIstType) :: this !< BndType object + integer(I4B), intent(in) :: idvsave !< flag and unit number for dependent-variable output + integer(I4B), intent(in) :: idvprint !< flag indicating if dependent-variable should be written to the model listing file + + call this%output_immobile_concentration(idvsave, idvprint) + call this%output_immobile_sorbate_concentration(idvsave, idvprint) + + end subroutine ist_ot_dv + + !> @ brief Output immobile domain aqueous concentration. + !< + subroutine output_immobile_concentration(this, idvsave, idvprint) + ! modules use TdisModule, only: kstp, endofperiod - ! -- dummy variables + ! dummy variables class(GwtIstType) :: this !< BndType object integer(I4B), intent(in) :: idvsave !< flag and unit number for dependent-variable output integer(I4B), intent(in) :: idvprint !< flag indicating if dependent-variable should be written to the model listing file - ! -- local + ! local integer(I4B) :: ipflg integer(I4B) :: ibinun ! - ! -- Save cim to a binary file. ibinun is a flag where 1 indicates that - ! cim should be written to a binary file if a binary file is open - ! for it. + ! Save cim to a binary file. ibinun is a flag where 1 indicates that + ! cim should be written to a binary file if a binary file is open for it. ipflg = 0 ibinun = 1 if (idvsave == 0) ibinun = 0 @@ -603,12 +722,51 @@ subroutine ist_ot_dv(this, idvsave, idvprint) iprint_opt=0, isav_opt=ibinun) end if ! - ! -- Print immobile domain concentrations to listing file + ! Print immobile domain concentrations to listing file if (idvprint /= 0) then call this%ocd%ocd_ot(ipflg, kstp, endofperiod, this%iout, & iprint_opt=idvprint, isav_opt=0) end if - end subroutine ist_ot_dv + + end subroutine output_immobile_concentration + + !> @ brief Output immobile domain sorbate concentration. + !< + subroutine output_immobile_sorbate_concentration(this, idvsave, idvprint) + ! modules + ! dummy + class(GwtIstType) :: this !< BndType object + integer(I4B), intent(in) :: idvsave !< flag and unit number for dependent-variable output + integer(I4B), intent(in) :: idvprint !< flag indicating if dependent-variable should be written to the model listing file + ! local + character(len=1) :: cdatafmp = ' ', editdesc = ' ' + integer(I4B) :: ibinun + integer(I4B) :: iprint, nvaluesp, nwidthp + real(DP) :: dinact + + ! Save cimsrb to a binary file. ibinun is a flag where 1 indicates that + ! cim should be written to a binary file if a binary file is open for it. + ! Set unit number for sorbate output + if (this%ioutsorbate /= 0) then + ibinun = 1 + else + ibinun = 0 + end if + if (idvsave == 0) ibinun = 0 + + ! save sorbate concentration array + if (ibinun /= 0) then + iprint = 0 + dinact = DHNOFLO + if (this%ioutsorbate /= 0) then + ibinun = this%ioutsorbate + call this%dis%record_array(this%cimsrb, this%iout, iprint, ibinun, & + ' SORBATE', cdatafmp, nvaluesp, & + nwidthp, editdesc, dinact) + end if + end if + + end subroutine output_immobile_sorbate_concentration !> @ brief Output IST package budget summary. !! @@ -660,17 +818,20 @@ subroutine ist_da(this) call mem_deallocate(this%icimout) call mem_deallocate(this%ibudgetout) call mem_deallocate(this%ibudcsv) + call mem_deallocate(this%ioutsorbate) call mem_deallocate(this%idcy) call mem_deallocate(this%isrb) call mem_deallocate(this%kiter) call mem_deallocate(this%cim) call mem_deallocate(this%cimnew) call mem_deallocate(this%cimold) + call mem_deallocate(this%cimsrb) call mem_deallocate(this%zetaim) call mem_deallocate(this%porosity) call mem_deallocate(this%volfrac) call mem_deallocate(this%bulk_density) call mem_deallocate(this%distcoef) + call mem_deallocate(this%sp2) call mem_deallocate(this%decay) call mem_deallocate(this%decaylast) call mem_deallocate(this%decayslast) @@ -715,6 +876,7 @@ subroutine allocate_scalars(this) call mem_allocate(this%icimout, 'ICIMOUT', this%memoryPath) call mem_allocate(this%ibudgetout, 'IBUDGETOUT', this%memoryPath) call mem_allocate(this%ibudcsv, 'IBUDCSV', this%memoryPath) + call mem_allocate(this%ioutsorbate, 'IOUTSORBATE', this%memoryPath) call mem_allocate(this%isrb, 'ISRB', this%memoryPath) call mem_allocate(this%idcy, 'IDCY', this%memoryPath) call mem_allocate(this%kiter, 'KITER', this%memoryPath) @@ -723,6 +885,7 @@ subroutine allocate_scalars(this) this%icimout = 0 this%ibudgetout = 0 this%ibudcsv = 0 + this%ioutsorbate = 0 this%isrb = 0 this%idcy = 0 this%kiter = 0 @@ -763,11 +926,20 @@ subroutine ist_allocate_arrays(this) if (this%isrb == 0) then call mem_allocate(this%bulk_density, 1, 'BULK_DENSITY', this%memoryPath) call mem_allocate(this%distcoef, 1, 'DISTCOEF', this%memoryPath) + call mem_allocate(this%sp2, 1, 'SP2', this%memoryPath) + call mem_allocate(this%cimsrb, 1, 'SORBATE', this%memoryPath) else call mem_allocate(this%bulk_density, this%dis%nodes, 'BULK_DENSITY', & this%memoryPath) call mem_allocate(this%distcoef, this%dis%nodes, 'DISTCOEF', & this%memoryPath) + call mem_allocate(this%cimsrb, this%dis%nodes, 'SORBATE', & + this%memoryPath) + if (this%isrb == 1) then + call mem_allocate(this%sp2, 1, 'SP2', this%memoryPath) + else + call mem_allocate(this%sp2, this%dis%nodes, 'SP2', this%memoryPath) + end if end if if (this%idcy == 0) then call mem_allocate(this%decay, 1, 'DECAY', this%memoryPath) @@ -795,6 +967,14 @@ subroutine ist_allocate_arrays(this) this%zetaim(n) = DZERO this%volfrac(n) = DZERO end do + do n = 1, size(this%bulk_density) + this%bulk_density(n) = DZERO + this%distcoef(n) = DZERO + this%cimsrb(n) = DZERO + end do + do n = 1, size(this%sp2) + this%sp2(n) = DZERO + end do do n = 1, size(this%decay) this%decay(n) = DZERO this%decaylast(n) = DZERO @@ -825,6 +1005,7 @@ subroutine read_options(this) class(GwtIstType), intent(inout) :: this !< GwtIstType object ! -- local character(len=LINELENGTH) :: errmsg, keyword + character(len=LINELENGTH) :: sorption character(len=LINELENGTH) :: fname character(len=:), allocatable :: keyword2 integer(I4B) :: ierr @@ -834,8 +1015,12 @@ subroutine read_options(this) character(len=*), parameter :: fmtisvflow = & "(4x,'CELL-BY-CELL FLOW INFORMATION WILL BE SAVED TO BINARY FILE & &WHENEVER ICBCFL IS NOT ZERO.')" - character(len=*), parameter :: fmtisrb = & + character(len=*), parameter :: fmtlinear = & &"(4x,'LINEAR SORPTION IS SELECTED. ')" + character(len=*), parameter :: fmtfreundlich = & + &"(4x,'FREUNDLICH SORPTION IS ACTIVE. ')" + character(len=*), parameter :: fmtlangmuir = & + &"(4x,'LANGMUIR SORPTION IS ACTIVE. ')" character(len=*), parameter :: fmtidcy1 = & &"(4x,'FIRST-ORDER DECAY IS ACTIVE. ')" character(len=*), parameter :: fmtidcy2 = & @@ -890,15 +1075,49 @@ subroutine read_options(this) call store_error('Optional BUDGETCSV keyword must be followed by & &FILEOUT') end if - case ('SORBTION', 'SORPTION') - this%isrb = 1 - write (this%iout, fmtisrb) + case ('SORBTION') + call store_error('SORBTION is not a valid option. Use & + &SORPTION instead.') + call this%parser%StoreErrorUnit() + case ('SORPTION') + call this%parser%GetStringCaps(sorption) + select case (sorption) + case ('LINEAR', '') + this%isrb = 1 + write (this%iout, fmtlinear) + case ('FREUNDLICH') + this%isrb = 2 + write (this%iout, fmtfreundlich) + case ('LANGMUIR') + this%isrb = 3 + write (this%iout, fmtlangmuir) + case default + call store_error('Unknown sorption type was specified ' & + & //'('//trim(adjustl(sorption))//').'// & + &' Sorption must be specified as LINEAR, & + &FREUNDLICH, or LANGMUIR.') + call this%parser%StoreErrorUnit() + end select case ('FIRST_ORDER_DECAY') this%idcy = 1 write (this%iout, fmtidcy1) case ('ZERO_ORDER_DECAY') this%idcy = 2 write (this%iout, fmtidcy2) + case ('SORBATE') + call this%parser%GetStringCaps(keyword) + if (keyword == 'FILEOUT') then + call this%parser%GetString(fname) + this%ioutsorbate = getunit() + call openfile(this%ioutsorbate, this%iout, fname, 'DATA(BINARY)', & + form, access, 'REPLACE', mode_opt=MNORMAL) + write (this%iout, fmtistbin) & + 'SORBATE', fname, this%ioutsorbate + else + errmsg = 'Optional SORBATE keyword must be '// & + 'followed by FILEOUT.' + call store_error(errmsg) + end if case default write (errmsg, '(a,a)') 'Unknown IST option: ', & trim(keyword) @@ -946,8 +1165,8 @@ subroutine read_data(this) character(len=:), allocatable :: line integer(I4B) :: istart, istop, lloc, ierr logical :: isfound, endOfBlock - logical, dimension(8) :: lname - character(len=24), dimension(8) :: aname + logical, dimension(9) :: lname + character(len=24), dimension(9) :: aname ! -- formats ! -- data data aname(1)/' BULK DENSITY'/ @@ -958,6 +1177,7 @@ subroutine read_data(this) data aname(6)/' FIRST ORDER TRANS RATE'/ data aname(7)/'IMMOBILE DOMAIN POROSITY'/ data aname(8)/'IMMOBILE VOLUME FRACTION'/ + data aname(9)/' SECOND SORPTION PARAM'/ ! ! -- initialize isfound = .false. @@ -1025,6 +1245,14 @@ subroutine read_data(this) this%parser%iuactive, this%volfrac, & aname(8)) lname(8) = .true. + case ('SP2') + if (this%isrb < 2) & + call mem_reallocate(this%sp2, this%dis%nodes, 'SP2', & + trim(this%memoryPath)) + call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, & + this%parser%iuactive, this%sp2, & + aname(9)) + lname(9) = .true. case ('THETAIM') write (errmsg, '(a)') & 'THETAIM is no longer supported. See Chapter 9 in & @@ -1058,6 +1286,12 @@ subroutine read_data(this) &GRIDDATA block.' call store_error(errmsg) end if + if (this%isrb > 1 .and. .not. lname(9)) then + write (errmsg, '(a)') 'Nonlinear sorption is active but SP2 & + ¬ specified. SP2 must be specified in GRIDDATA block when & + &FREUNDLICH or LANGMUIR sorption is active.' + call store_error(errmsg) + end if else if (lname(1)) then write (this%iout, '(1x,a)') 'Warning. Sorption is not active but & @@ -1069,6 +1303,10 @@ subroutine read_data(this) &distribution coefficient was specified. DISTCOEF will have & &no affect on simulation results.' end if + if (lname(9)) then + write (this%iout, '(1x,a)') 'Warning. Sorption is not active but & + &SP2 was specified. SP2 will have no affect on simulation results.' + end if end if ! ! -- Check for required decay/production rate coefficients @@ -1169,7 +1407,8 @@ end function get_thetaim !! !< subroutine get_ddterm(thetaim, vcell, delt, swtpdt, & - volfracim, rhobim, kd, lambda1im, lambda2im, & + volfracim, rhobim, kdnew, kdold, & + lambda1im, lambda2im, & gamma1im, gamma2im, zetaim, ddterm, f) ! -- dummy real(DP), intent(in) :: thetaim !< immobile domain porosity @@ -1178,7 +1417,8 @@ subroutine get_ddterm(thetaim, vcell, delt, swtpdt, & real(DP), intent(in) :: swtpdt !< cell saturation at end of time step real(DP), intent(in) :: volfracim !< volume fraction of immobile domain real(DP), intent(in) :: rhobim !< bulk density for the immobile domain (fim * rhob) - real(DP), intent(in) :: kd !< distribution coefficient for linear isotherm + real(DP), intent(in) :: kdnew !< effective distribution coefficient for new time + real(DP), intent(in) :: kdold !< effective distribution coefficient for old time real(DP), intent(in) :: lambda1im !< first-order decay rate in aqueous phase real(DP), intent(in) :: lambda2im !< first-order decay rate in sorbed phase real(DP), intent(in) :: gamma1im !< zero-order decay rate in aqueous phase @@ -1198,10 +1438,10 @@ subroutine get_ddterm(thetaim, vcell, delt, swtpdt, & ! information guide (mf6suptechinfo.pdf) ddterm(1) = thetaim * vcell * tled ddterm(2) = thetaim * vcell * tled - ddterm(3) = volfracim * rhobim * vcell * kd * tled - ddterm(4) = volfracim * rhobim * vcell * kd * tled + ddterm(3) = volfracim * rhobim * vcell * kdnew * tled + ddterm(4) = volfracim * rhobim * vcell * kdold * tled ddterm(5) = thetaim * lambda1im * vcell - ddterm(6) = lambda2im * volfracim * rhobim * kd * vcell + ddterm(6) = lambda2im * volfracim * rhobim * kdnew * vcell ddterm(7) = thetaim * gamma1im * vcell ddterm(8) = gamma2im * volfracim * rhobim * vcell ddterm(9) = vcell * swtpdt * zetaim diff --git a/src/Model/GroundWaterTransport/gwt-mst.f90 b/src/Model/GroundWaterTransport/gwt-mst.f90 index d20fb7d15ed..effd4f9baa6 100644 --- a/src/Model/GroundWaterTransport/gwt-mst.f90 +++ b/src/Model/GroundWaterTransport/gwt-mst.f90 @@ -962,7 +962,6 @@ subroutine mst_ot_flow(this, icbcfl, icbcun) integer(I4B), intent(in) :: icbcun !< flag indication if cell-by-cell data should be saved ! -- local integer(I4B) :: ibinun - !character(len=16), dimension(2) :: aname integer(I4B) :: iprint, nvaluesp, nwidthp character(len=1) :: cdatafmp = ' ', editdesc = ' ' real(DP) :: dinact @@ -1224,7 +1223,8 @@ subroutine read_options(this) ! -- dummy class(GwtMstType) :: this !< GwtMstType object ! -- local - character(len=LINELENGTH) :: keyword, keyword2 + character(len=LINELENGTH) :: keyword + character(len=LINELENGTH) :: sorption character(len=MAXCHARLEN) :: fname integer(I4B) :: ierr logical :: isfound, endOfBlock @@ -1232,7 +1232,7 @@ subroutine read_options(this) character(len=*), parameter :: fmtisvflow = & "(4x,'CELL-BY-CELL FLOW INFORMATION WILL BE SAVED TO BINARY FILE & &WHENEVER ICBCFL IS NOT ZERO.')" - character(len=*), parameter :: fmtisrb = & + character(len=*), parameter :: fmtlinear = & &"(4x,'LINEAR SORPTION IS ACTIVE. ')" character(len=*), parameter :: fmtfreundlich = & &"(4x,'FREUNDLICH SORPTION IS ACTIVE. ')" @@ -1261,19 +1261,28 @@ subroutine read_options(this) case ('SAVE_FLOWS') this%ipakcb = -1 write (this%iout, fmtisvflow) - case ('SORBTION', 'SORPTION') - this%isrb = 1 - call this%parser%GetStringCaps(keyword2) - if (trim(adjustl(keyword2)) == 'LINEAR') this%isrb = 1 - if (trim(adjustl(keyword2)) == 'FREUNDLICH') this%isrb = 2 - if (trim(adjustl(keyword2)) == 'LANGMUIR') this%isrb = 3 - select case (this%isrb) - case (1) - write (this%iout, fmtisrb) - case (2) + case ('SORBTION') + call store_error('SORBTION is not a valid option. Use & + &SORPTION instead.') + call this%parser%StoreErrorUnit() + case ('SORPTION') + call this%parser%GetStringCaps(sorption) + select case (sorption) + case ('LINEAR', '') + this%isrb = 1 + write (this%iout, fmtlinear) + case ('FREUNDLICH') + this%isrb = 2 write (this%iout, fmtfreundlich) - case (3) + case ('LANGMUIR') + this%isrb = 3 write (this%iout, fmtlangmuir) + case default + call store_error('Unknown sorption type was specified ' & + & //'('//trim(adjustl(sorption))//').'// & + &' Sorption must be specified as LINEAR, & + &FREUNDLICH, or LANGMUIR.') + call this%parser%StoreErrorUnit() end select case ('FIRST_ORDER_DECAY') this%idcy = 1 @@ -1640,6 +1649,42 @@ function get_langmuir_derivative(conc, kl, sbar) result(derv) return end function + !> @ brief Get effective Freundlich distribution coefficient + !< + function get_freundlich_kd(conc, kf, a) result(kd) + ! -- dummy + real(DP), intent(in) :: conc !< solute concentration + real(DP), intent(in) :: kf !< freundlich constant + real(DP), intent(in) :: a !< freundlich exponent + ! -- return + real(DP) :: kd !< effective distribution coefficient + ! + if (conc > DZERO) then + kd = kf * conc**(a - DONE) + else + kd = DZERO + end if + return + end function get_freundlich_kd + + !> @ brief Get effective Langmuir distribution coefficient + !< + function get_langmuir_kd(conc, kl, sbar) result(kd) + ! -- dummy + real(DP), intent(in) :: conc !< solute concentration + real(DP), intent(in) :: kl !< langmuir constant + real(DP), intent(in) :: sbar !< langmuir sorption sites + ! -- return + real(DP) :: kd !< effective distribution coefficient + ! + if (conc > DZERO) then + kd = (kl * sbar) / (DONE + kl * conc) + else + kd = DZERO + end if + return + end function get_langmuir_kd + !> @ brief Calculate zero-order decay rate and constrain if necessary !! !! Function to calculate the zero-order decay rate from the user specified