diff --git a/autotest/data/vsc04-laktab/stg-vol-surfarea.dat b/autotest/data/vsc04-laktab/stg-vol-surfarea.dat new file mode 100644 index 00000000000..edc44af27e1 --- /dev/null +++ b/autotest/data/vsc04-laktab/stg-vol-surfarea.dat @@ -0,0 +1,152 @@ +stage,Volume,Surf Area +29.5656,0,0 +29.602176,7645.54858,3870.96 +29.638752,15291.09716,7741.92 +29.675328,22936.64574,11612.88 +29.711904,30582.19432,15483.84 +29.74848,38227.7429,19354.8 +29.785056,45873.29148,23225.76 +29.821632,53518.84006,27096.72 +29.858208,61164.38864,30967.68 +29.894784,68809.93722,34838.64 +29.93136,76455.4858,38709.6 +29.967936,84101.03438,42580.56 +30.004512,91746.58296,46451.52 +30.041088,99392.13154,50322.48 +30.077664,107037.6801,54193.44 +30.11424,114683.2287,58064.4 +30.150816,122328.7773,61935.36 +30.187392,129974.3259,65806.32 +30.223968,137619.8744,69677.28 +30.260544,145265.423,73548.24 +30.29712,152910.9716,77419.2 +30.333696,160556.5202,81290.16 +30.370272,168202.0688,85161.12 +30.406848,175847.6173,89032.08 +30.443424,183493.1659,92903.04 +30.48,191138.7145,96774 +30.516576,198784.2631,100644.96 +30.553152,206429.8117,104515.92 +30.589728,214075.3602,108386.88 +30.626304,221720.9088,112257.84 +30.66288,229366.4574,116128.8 +30.699456,237012.006,119999.76 +30.736032,244657.5546,123870.72 +30.772608,252303.1031,127741.68 +30.809184,259948.6517,131612.64 +30.84576,267594.2003,135483.6 +30.882336,275239.7489,139354.56 +30.918912,282885.2975,143225.52 +30.955488,290530.846,147096.48 +30.992064,298176.3946,150967.44 +31.02864,305821.9432,154838.4 +31.065216,313467.4918,158709.36 +31.101792,321113.0404,162580.32 +31.138368,328758.5889,166451.28 +31.174944,336404.1375,170322.24 +31.21152,344049.6861,174193.2 +31.248096,351695.2347,178064.16 +31.284672,359340.7833,181935.12 +31.321248,366986.3318,185806.08 +31.357824,374631.8804,189677.04 +31.3944,382277.429,193548 +31.430976,389922.9776,197418.96 +31.467552,397568.5262,201289.92 +31.504128,405214.0747,205160.88 +31.540704,412859.6233,209031.84 +31.57728,420505.1719,212902.8 +31.613856,428150.7205,216773.76 +31.650432,435796.2691,220644.72 +31.687008,443441.8176,224515.68 +31.723584,451087.3662,228386.64 +31.76016,458732.9148,232257.6 +31.796736,466378.4634,236128.56 +31.833312,474024.012,239999.52 +31.869888,481669.5605,243870.48 +31.906464,489315.1091,247741.44 +31.94304,496960.6577,251612.4 +31.979616,504606.2063,255483.36 +32.016192,512251.7548,259354.32 +32.052768,519897.3034,263225.28 +32.089344,527542.852,267096.24 +32.12592,535188.4006,270967.2 +32.162496,542833.9492,274838.16 +32.199072,550479.4977,278709.12 +32.235648,558125.0463,282580.08 +32.272224,565770.5949,286451.04 +32.3088,573416.1435,290322 +32.345376,581061.6921,294192.96 +32.381952,588707.2406,298063.92 +32.418528,596352.7892,301934.88 +32.455104,603998.3378,305805.84 +32.49168,611643.8864,309676.8 +32.528256,619289.435,313547.76 +32.564832,626934.9835,317418.72 +32.601408,634580.5321,321289.68 +32.637984,651287.4716,325160.64 +32.67456,672525.1066,329031.6 +32.711136,693762.7415,332902.56 +32.747712,715000.3764,336773.52 +32.784288,736238.0114,340644.48 +32.820864,757475.6463,344515.44 +32.85744,778713.2813,348386.4 +32.894016,799950.9162,352257.36 +32.930592,821188.5512,356128.32 +32.967168,842426.1861,359999.28 +33.003744,863663.8211,363870.24 +33.04032,884901.456,367741.2 +33.076896,906139.0909,371612.16 +33.113472,927376.7259,375483.12 +33.150048,948614.3608,379354.08 +33.186624,969851.9958,383225.04 +33.2232,991089.6307,387096 +33.259776,1012327.266,390966.96 +33.296352,1033564.901,394837.92 +33.332928,1054802.536,398708.88 +33.369504,1076040.17,402579.84 +33.40608,1097277.805,406450.8 +33.442656,1118515.44,410321.76 +33.479232,1139753.075,414192.72 +33.515808,1160990.71,418063.68 +33.552384,1182228.345,421934.64 +33.58896,1203465.98,425805.6 +33.625536,1224703.615,429676.56 +33.662112,1245941.25,433547.52 +33.698688,1267178.885,437418.48 +33.735264,1288416.52,441289.44 +33.77184,1309654.155,445160.4 +33.808416,1330891.79,449031.36 +33.844992,1352129.425,452902.32 +33.881568,1373367.06,456773.28 +33.918144,1394604.695,460644.24 +33.95472,1415842.33,464515.2 +33.991296,1437079.965,468386.16 +34.027872,1458317.599,472257.12 +34.064448,1479555.234,476128.08 +34.101024,1500792.869,479999.04 +34.1376,1522030.504,483870 +34.174176,1543268.139,487740.96 +34.210752,1564505.774,491611.92 +34.247328,1585743.409,495482.88 +34.283904,1606981.044,499353.84 +34.32048,1628218.679,503224.8 +34.357056,1649456.314,507095.76 +34.393632,1670693.949,510966.72 +34.430208,1691931.584,514837.68 +34.466784,1713169.219,518708.64 +34.50336,1734406.854,522579.6 +34.539936,1755644.489,526450.56 +34.576512,1776882.124,530321.52 +34.613088,1798119.759,534192.48 +34.649664,1819357.394,538063.44 +34.68624,1840595.028,541934.4 +34.722816,1861832.663,545805.36 +34.759392,1883070.298,549676.32 +34.795968,1904307.933,553547.28 +34.832544,1925545.568,557418.24 +34.86912,1946783.203,561289.2 +34.905696,1968020.838,565160.16 +34.942272,1989258.473,569031.12 +34.978848,2010496.108,572902.08 +35.015424,2031733.743,576773.04 +35.052,2052971.378,580644 diff --git a/autotest/test_gwf_returncodes.py b/autotest/test_gwf_returncodes.py index daea5cec788..9886798e0b8 100644 --- a/autotest/test_gwf_returncodes.py +++ b/autotest/test_gwf_returncodes.py @@ -19,7 +19,7 @@ name = "gwf_ret_codes01" base_ws = os.path.join("temp", name) if not os.path.isdir(base_ws): - os.makedirs(base_ws) + os.makedirs(base_ws, exist_ok=True) app = "mf6" if sys.platform.lower() == "win32": app += ".exe" diff --git a/autotest/test_gwf_vsc01.py b/autotest/test_gwf_vsc01.py new file mode 100644 index 00000000000..0c547db2e77 --- /dev/null +++ b/autotest/test_gwf_vsc01.py @@ -0,0 +1,359 @@ +# ## Test problem for VSC +# +# Uses constant head and general-head boundaries on the left and right +# sides of the model domain, respectively, to drive flow from left to +# right. Tests that head-dependent boundary conditions are properly +# accounting for viscosity when VSC is active. +# + +# Imports + +import os +import sys + +import numpy as np +import pytest + +try: + import flopy +except: + msg = "Error. FloPy package is not available.\n" + msg += "Try installing using the following command:\n" + msg += " pip install flopy" + raise Exception(msg) + +from framework import testing_framework +from simulation import Simulation + +hyd_cond = [1205.49396942506, 864.0] # Hydraulic conductivity (m/d) +ex = ["no-vsc01-bnd", "vsc01-bnd", "no-vsc01-k"] +viscosity_on = [False, True, False] +hydraulic_conductivity = [hyd_cond[0], hyd_cond[1], hyd_cond[1]] +exdirs = [] +for s in ex: + exdirs.append(os.path.join("temp", s)) + +# Model units + +length_units = "cm" +time_units = "seconds" + +# Table of model parameters + +nper = 1 # Number of periods +nstp = 500 # Number of time steps +perlen = 0.5 # Simulation time length ($d$) +nlay = 1 # Number of layers +nrow = 10 # Number of rows +ncol = 80 # Number of columns +system_length = 2.0 # Length of system ($m$) +delr = 1.0 # Column width ($m$) +delc = 1.0 # Row width ($m$) +delv = 1.0 # Layer thickness +top = 1.0 # Top of the model ($m$) +initial_temperature = 35.0 # Initial temperature (unitless) +porosity = 0.26 # porosity (unitless) +K_therm = 2.0 # Thermal conductivity # ($W/m/C$) +rho_water = 1000 # Density of water ($kg/m^3$) +rho_solids = 2650 # Density of the aquifer material ($kg/m^3$) +C_p_w = 4180 # Heat Capacity of water ($J/kg/C$) +C_s = 880 # Heat capacity of the solids ($J/kg/C$) +D_m = K_therm / (porosity * rho_water * C_p_w) +rhob = (1 - porosity) * rho_solids # Bulk density ($kg/m^3$) +K_d = C_s / (rho_water * C_p_w) # Partitioning coefficient ($m^3/kg$) +inflow = 5.7024 # ($m^3/d$) + +botm = [top - k * delv for k in range(1, nlay + 1)] + +nouter, ninner = 100, 300 +hclose, rclose, relax = 1e-10, 1e-6, 0.97 + +# +# MODFLOW 6 flopy GWF simulation object (sim) is returned +# + + +def build_model(idx, dir): + # Base simulation and model name and workspace + ws = dir + name = ex[idx] + + print("Building model...{}".format(name)) + + # generate names for each model + gwfname = "gwf-" + name + gwtname = "gwt-" + name + + sim = flopy.mf6.MFSimulation( + sim_name=name, version="mf6", exe_name="mf6", sim_ws=ws + ) + + # Instantiating time discretization + tdis_ds = ((perlen, nstp, 1.0),) + flopy.mf6.ModflowTdis( + sim, nper=nper, perioddata=tdis_ds, time_units=time_units + ) + gwf = flopy.mf6.ModflowGwf(sim, modelname=gwfname, save_flows=True) + + # Instantiating solver + ims = flopy.mf6.ModflowIms( + sim, + print_option="ALL", + 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="{}.ims".format(gwfname), + ) + sim.register_ims_package(ims, [gwfname]) + + # Instantiating DIS + flopy.mf6.ModflowGwfdis( + gwf, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + ) + + # Instantiating NPF + flopy.mf6.ModflowGwfnpf( + gwf, + save_specific_discharge=True, + icelltype=0, + k=hydraulic_conductivity[idx], + ) + flopy.mf6.ModflowGwfic(gwf, strt=0.0) + + # Instantiating VSC + if viscosity_on[idx]: + # Instantiate viscosity (VSC) package + vsc_filerecord = "{}.vsc.bin".format(gwfname) + vsc_pd = [(0, 0.0, 20.0, gwtname, "temperature")] + flopy.mf6.ModflowGwfvsc( + gwf, + viscref=8.904e-4, + viscosity_filerecord=vsc_filerecord, + thermal_formulation="nonlinear", + thermal_a2=10.0, + thermal_a3=248.37, + thermal_a4=133.16, + nviscspecies=len(vsc_pd), + packagedata=vsc_pd, + pname="vsc", + filename="{}.vsc".format(gwfname), + ) + + # Instantiating GHB + ghbcond = hydraulic_conductivity[idx] * delv * delc / (0.5 * delr) + ghbspd = [ + [(0, i, ncol - 1), top, ghbcond, initial_temperature] + for i in range(nrow) + ] + flopy.mf6.ModflowGwfghb( + gwf, + stress_period_data=ghbspd, + pname="GHB-1", + auxiliary="temperature", + ) + + # Instantiating CHD + chdspd = [[(0, i, 0), 2.0, initial_temperature] for i in range(nrow)] + flopy.mf6.ModflowGwfchd( + gwf, + stress_period_data=chdspd, + pname="CHD-1", + auxiliary="temperature", + ) + + # Instatiating OC + head_filerecord = "{}.hds".format(gwfname) + budget_filerecord = "{}.bud".format(gwfname) + flopy.mf6.ModflowGwfoc( + gwf, + head_filerecord=head_filerecord, + budget_filerecord=budget_filerecord, + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) + + # Setup the GWT model for simulating heat transport + # ------------------------------------------------- + gwt = flopy.mf6.ModflowGwt(sim, modelname=gwtname) + + # Instantiating solver for GWT + imsgwt = flopy.mf6.ModflowIms( + sim, + print_option="ALL", + 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="{}.ims".format(gwtname), + ) + sim.register_ims_package(imsgwt, [gwtname]) + + # Instantiating DIS for GWT + flopy.mf6.ModflowGwtdis( + gwt, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + ) + + # Instantiating MST for GWT + flopy.mf6.ModflowGwtmst( + gwt, + porosity=porosity, + sorption="linear", + bulk_density=rhob, + distcoef=K_d, + pname="MST-1", + filename="{}.mst".format(gwtname), + ) + + # Instantiating IC for GWT + flopy.mf6.ModflowGwtic(gwt, strt=initial_temperature) + + # Instantiating ADV for GWT + flopy.mf6.ModflowGwtadv(gwt, scheme="UPSTREAM") + + # Instantiating DSP for GWT + flopy.mf6.ModflowGwtdsp(gwt, xt3d_off=True, diffc=D_m) + + # Instantiating SSM for GWT + sourcerecarray = [ + ("CHD-1", "AUX", "TEMPERATURE"), + ("GHB-1", "AUX", "TEMPERATURE"), + ] + flopy.mf6.ModflowGwtssm(gwt, sources=sourcerecarray) + + # Instantiating OC for GWT + flopy.mf6.ModflowGwtoc( + gwt, + concentration_filerecord="{}.ucn".format(gwtname), + saverecord=[("CONCENTRATION", "ALL")], + printrecord=[("CONCENTRATION", "LAST"), ("BUDGET", "LAST")], + ) + + # Instantiating GWF/GWT Exchange + flopy.mf6.ModflowGwfgwt( + sim, exgtype="GWF6-GWT6", exgmnamea=gwfname, exgmnameb=gwtname + ) + + return sim, None + + +def eval_results(sim): + print("evaluating results...") + + # read flow results from model + name = ex[sim.idxsim] + gwfname = "gwf-" + name + + fname = gwfname + ".bud" + fname = os.path.join(sim.simpath, fname) + assert os.path.isfile(fname) + budobj = flopy.utils.CellBudgetFile(fname, precision="double") + outbud = budobj.get_data(text=" GHB") + + # Establish known answer: + stored_ans = -151.63446156218242 + + if sim.idxsim == 0: + no_vsc_bud_last = np.array(outbud[-1].tolist()) + sim_val_1 = no_vsc_bud_last[:, 2].sum() + + # Ensure latest simulated value hasn't changed from stored answer + assert np.allclose( + sim_val_1, stored_ans, atol=1e-4 + ), "Flow in the " + exdirs[ + 0 + ] + " test problem (doesn't simulate " "viscosity) has changed,\n should be " + str( + stored_ans + ) + " but instead is " + str( + sim_val_1 + ) + + elif sim.idxsim == 1: + with_vsc_bud_last = np.array(outbud[-1].tolist()) + sim_val_2 = with_vsc_bud_last[:, 2].sum() + + # Ensure latest simulated value hasn't changed from stored answer + assert np.allclose( + sim_val_2, stored_ans, atol=1e-4 + ), "Flow in the " + exdirs[ + 1 + ] + " test problem (simulates " "viscosity) has changed,\n should be " + str( + stored_ans + ) + " but instead is " + str( + sim_val_2 + ) + + elif sim.idxsim == 2: + no_vsc_low_k_bud_last = np.array(outbud[-1].tolist()) + sim_val_3 = no_vsc_low_k_bud_last[:, 2].sum() + + # Ensure the flow leaving model 3 is less than what leaves model 2 + assert abs(stored_ans) > abs(sim_val_3), ( + "Exit flow from model " + + exdirs[1] + + " should be greater than flow exiting " + + exdirs[2] + + ", but it is not." + ) + + +# - No need to change any code below +@pytest.mark.parametrize( + "idx, dir", + list(enumerate(exdirs)), +) +def test_mf6model(idx, dir): + # initialize testing framework + test = testing_framework() + + # build the model + test.build_mf6_models(build_model, idx, dir) + + # run the test model + test.run_mf6(Simulation(dir, exfunc=eval_results, idxsim=idx)) + + +def main(): + # initialize testing framework + test = testing_framework() + + # run the test model + for idx, dir in enumerate(exdirs): + test.build_mf6_models(build_model, idx, dir) + sim = Simulation(dir, exfunc=eval_results, idxsim=idx) + test.run_mf6(sim) + + +if __name__ == "__main__": + # print message + print(f"standalone run of {os.path.basename(__file__)}") + + # run main routine + main() diff --git a/autotest/test_gwf_vsc02.py b/autotest/test_gwf_vsc02.py new file mode 100644 index 00000000000..2e6039d5116 --- /dev/null +++ b/autotest/test_gwf_vsc02.py @@ -0,0 +1,363 @@ +# ## Test problem for VSC +# +# Uses general-head and drain boundaries on the left and right +# sides of the model domain, respectively, to drive flow from left to +# right. Tests that head-dependent boundary conditions are properly +# accounting for viscosity when VSC is active. Similar to gwf-vsc01-bnd +# but employs head-dependent boundary on the left and right side of the +# model + +# Imports + +import os +import sys + +import numpy as np +import pytest + +try: + import flopy +except: + msg = "Error. FloPy package is not available.\n" + msg += "Try installing using the following command:\n" + msg += " pip install flopy" + raise Exception(msg) + +# Import common functionality +from framework import testing_framework +from simulation import Simulation + +# Setup scenario input +hyd_cond = [1205.49396942506, 864.0] # Hydraulic conductivity (m/d) +ex = ["no-vsc02-bnd", "vsc02-bnd", "no-vsc02-k"] +viscosity_on = [False, True, False] +hydraulic_conductivity = [hyd_cond[0], hyd_cond[1], hyd_cond[1]] +exdirs = [] +for s in ex: + exdirs.append(os.path.join("temp", s)) + +# Model units + +length_units = "cm" +time_units = "seconds" + +# Table of model parameters + +nper = 1 # Number of periods +nstp = 500 # Number of time steps +perlen = 0.5 # Simulation time length ($d$) +nlay = 1 # Number of layers +nrow = 10 # Number of rows +ncol = 80 # Number of columns +system_length = 2.0 # Length of system ($m$) +delr = 1.0 # Column width ($m$) +delc = 1.0 # Row width ($m$) +delv = 1.0 # Layer thickness +top = 1.0 # Top of the model ($m$) +initial_temperature = 35.0 # Initial temperature (unitless) +porosity = 0.26 # porosity (unitless) +K_therm = 2.0 # Thermal conductivity # ($W/m/C$) +rho_water = 1000 # Density of water ($kg/m^3$) +rho_solids = 2650 # Density of the aquifer material ($kg/m^3$) +C_p_w = 4180 # Heat Capacity of water ($J/kg/C$) +C_s = 880 # Heat capacity of the solids ($J/kg/C$) +D_m = K_therm / (porosity * rho_water * C_p_w) +rhob = (1 - porosity) * rho_solids # Bulk density ($kg/m^3$) +K_d = C_s / (rho_water * C_p_w) # Partitioning coefficient ($m^3/kg$) +inflow = 5.7024 # ($m^3/d$) + +botm = [top - k * delv for k in range(1, nlay + 1)] + +nouter, ninner = 100, 300 +hclose, rclose, relax = 1e-10, 1e-6, 0.97 + +# +# MODFLOW 6 flopy GWF simulation object (sim) is returned +# + + +def build_model(idx, dir): + # Base simulation and model name and workspace + ws = dir + name = ex[idx] + + print("Building model...{}".format(name)) + + # generate names for each model + gwfname = "gwf-" + name + gwtname = "gwt-" + name + + sim = flopy.mf6.MFSimulation( + sim_name=name, sim_ws=ws, exe_name="mf6", version="mf6" + ) + + # Instantiating time discretization + tdis_ds = ((perlen, nstp, 1.0),) + flopy.mf6.ModflowTdis( + sim, nper=nper, perioddata=tdis_ds, time_units=time_units + ) + gwf = flopy.mf6.ModflowGwf(sim, modelname=gwfname, save_flows=True) + + # Instantiating solver + ims = flopy.mf6.ModflowIms( + sim, + print_option="ALL", + 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="{}.ims".format(gwfname), + ) + sim.register_ims_package(ims, [gwfname]) + + # Instantiating DIS + flopy.mf6.ModflowGwfdis( + gwf, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + ) + + # Instantiating NPF + flopy.mf6.ModflowGwfnpf( + gwf, + save_specific_discharge=True, + icelltype=0, + k=hydraulic_conductivity[idx], + ) + flopy.mf6.ModflowGwfic(gwf, strt=0.0) + + # Instantiating VSC + if viscosity_on[idx]: + # Instantiate viscosity (VSC) package + vsc_filerecord = "{}.vsc.bin".format(gwfname) + vsc_pd = [(0, 0.0, 20.0, gwtname, "temperature")] + flopy.mf6.ModflowGwfvsc( + gwf, + viscref=8.904e-4, + viscosity_filerecord=vsc_filerecord, + thermal_formulation="nonlinear", + thermal_a2=10.0, + thermal_a3=248.37, + thermal_a4=133.16, + nviscspecies=len(vsc_pd), + packagedata=vsc_pd, + pname="vsc", + filename="{}.vsc".format(gwfname), + ) + + # Instantiating GHB + ghbcond = hydraulic_conductivity[idx] * delv * delc / (0.5 * delr) + ghbspd = [ + [(0, i, 0), top + 3, ghbcond, initial_temperature] for i in range(nrow) + ] + flopy.mf6.ModflowGwfghb( + gwf, + stress_period_data=ghbspd, + pname="GHB-1", + auxiliary="temperature", + ) + + # Instantiating DRN + drnspd = [ + [(0, i, ncol - 1), top, 1.2 * ghbcond, initial_temperature] + for i in range(nrow) + ] + flopy.mf6.ModflowGwfdrn( + gwf, + stress_period_data=drnspd, + pname="DRN-1", + auxiliary="temperature", + ) + + # Instatiatingi OC + head_filerecord = "{}.hds".format(gwfname) + budget_filerecord = "{}.bud".format(gwfname) + flopy.mf6.ModflowGwfoc( + gwf, + head_filerecord=head_filerecord, + budget_filerecord=budget_filerecord, + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) + + # Setup the GWT model for simulating heat transport + # ------------------------------------------------- + gwt = flopy.mf6.ModflowGwt(sim, modelname=gwtname) + + # Instantiating solver for GWT + imsgwt = flopy.mf6.ModflowIms( + sim, + print_option="ALL", + 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="{}.ims".format(gwtname), + ) + sim.register_ims_package(imsgwt, [gwtname]) + + # Instantiating DIS for GWT + flopy.mf6.ModflowGwtdis( + gwt, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + ) + + # Instantiating MST for GWT + flopy.mf6.ModflowGwtmst( + gwt, + porosity=porosity, + sorption="linear", + bulk_density=rhob, + distcoef=K_d, + pname="MST-1", + filename="{}.mst".format(gwtname), + ) + + # Instantiating IC for GWT + flopy.mf6.ModflowGwtic(gwt, strt=initial_temperature) + + # Instantiating ADV for GWT + flopy.mf6.ModflowGwtadv(gwt, scheme="UPSTREAM") + + # Instantiating DSP for GWT + flopy.mf6.ModflowGwtdsp(gwt, xt3d_off=True, diffc=D_m) + + # Instantiating SSM for GWT + sourcerecarray = [ + ("GHB-1", "AUX", "TEMPERATURE"), + ("DRN-1", "AUX", "TEMPERATURE"), + ] + flopy.mf6.ModflowGwtssm(gwt, sources=sourcerecarray) + + # Instantiating OC for GWT + flopy.mf6.ModflowGwtoc( + gwt, + concentration_filerecord="{}.ucn".format(gwtname), + saverecord=[("CONCENTRATION", "ALL")], + printrecord=[("CONCENTRATION", "LAST"), ("BUDGET", "LAST")], + ) + + # Instantiating GWF/GWT Exchange + flopy.mf6.ModflowGwfgwt( + sim, exgtype="GWF6-GWT6", exgmnamea=gwfname, exgmnameb=gwtname + ) + return sim, None + + +def eval_results(sim): + print("evaluating results...") + + # read flow results from model + name = ex[sim.idxsim] + gwfname = "gwf-" + name + + fname = gwfname + ".bud" + fname = os.path.join(sim.simpath, fname) + assert os.path.isfile(fname) + budobj = flopy.utils.CellBudgetFile(fname, precision="double") + outbud = budobj.get_data(text=" GHB") + + # Establish known answer: + stored_ans = 452.5316256451224 + + if sim.idxsim == 0: + no_vsc_bud_last = np.array(outbud[-1].tolist()) + sim_val_1 = no_vsc_bud_last[:, 2].sum() + + # Ensure latest simulated value hasn't changed from stored answer + assert np.allclose( + sim_val_1, stored_ans, atol=1e-4 + ), "Flow in the " + exdirs[ + 0 + ] + " test problem (doesn't simulate " "viscosity) has changed,\n should be " + str( + stored_ans + ) + " but instead is " + str( + sim_val_1 + ) + + elif sim.idxsim == 1: + with_vsc_bud_last = np.array(outbud[-1].tolist()) + sim_val_2 = with_vsc_bud_last[:, 2].sum() + + # Ensure latest simulated value hasn't changed from stored answer + assert np.allclose( + sim_val_2, stored_ans, atol=1e-4 + ), "Flow in the " + exdirs[ + 1 + ] + " test problem (simulates " "viscosity) has changed,\n should be " + str( + stored_ans + ) + " but instead is " + str( + sim_val_2 + ) + + elif sim.idxsim == 2: + no_vsc_low_k_bud_last = np.array(outbud[-1].tolist()) + sim_val_3 = no_vsc_low_k_bud_last[:, 2].sum() + + # Ensure the flow leaving model 3 is less than what leaves model 2 + assert abs(stored_ans) > abs(sim_val_3), ( + "Exit flow from model " + + exdirs[1] + + " should be greater than flow existing " + + exdirs[2] + + ", but it is not." + ) + + +# - No need to change any code below +@pytest.mark.parametrize( + "idx, dir", + list(enumerate(exdirs)), +) +def test_mf6model(idx, dir): + # initialize testing framework + test = testing_framework() + + # build the model + test.build_mf6_models(build_model, idx, dir) + + # run the test model + test.run_mf6(Simulation(dir, exfunc=eval_results, idxsim=idx)) + + +def main(): + # initialize testing framework + test = testing_framework() + + # run the test model + for idx, dir in enumerate(exdirs): + test.build_mf6_models(build_model, idx, dir) + sim = Simulation(dir, exfunc=eval_results, idxsim=idx) + test.run_mf6(sim) + + +if __name__ == "__main__": + # print message + print(f"standalone run of {os.path.basename(__file__)}") + + # run main routine + main() diff --git a/autotest/test_gwf_vsc03_sfr.py b/autotest/test_gwf_vsc03_sfr.py new file mode 100644 index 00000000000..a3c79f0620b --- /dev/null +++ b/autotest/test_gwf_vsc03_sfr.py @@ -0,0 +1,576 @@ +# Scenario envisioned by this test is a river running through a V-shaped +# valley that loses water to the aquifer at the upper end until it goes +# dry, then begins to gain flow again in the lower reaches. River water +# enters the simulation at 8 deg C. Aquifer water starts out at 35 deg C. +# Reference viscosity temperature is 20 deg C. With the VSC package active, +# the simulation should predict less loss of river water to the aquifer +# and more discharge of gw to the stream, compared to the same simulation +# with the VSC package inactive. + +# Imports + +import os +import sys + +import numpy as np +import pytest + +try: + import flopy +except: + msg = "Error. FloPy package is not available.\n" + msg += "Try installing using the following command:\n" + msg += " pip install flopy" + raise Exception(msg) + +from framework import testing_framework +from simulation import Simulation + +ex = ["no-vsc-sfr01", "vsc-sfr01"] +viscosity_on = [False, True] +exdirs = [] +for s in ex: + exdirs.append(os.path.join("temp", s)) + +# Equation for determining land surface elevation with a stream running down the middle +def topElev_sfrCentered(x, y): + return ((-0.003 * x) + 260.0) + ( + ((-2e-9 * (x - 5000.0)) + 1e-5) * (y + 1500.0) ** 2 + ) + + +# Model units +length_units = "m" +time_units = "days" + +# model domain and grid definition +Lx = 10000.0 +Ly = 3000.0 +nrow = 60 +ncol = 200 +nlay = 1 +delr = Lx / ncol +delc = Ly / nrow +xmax = ncol * delr +ymax = nrow * delc +X, Y = np.meshgrid( + np.linspace(delr / 2, xmax - delr / 2, ncol), + np.linspace(ymax - delc / 2, 0 + delc / 2, nrow), +) +ibound = np.ones((nlay, nrow, ncol)) +# Because eqn uses negative values in the Y direction, need to do a little manipulation +Y_m = -1 * np.flipud(Y) +top = topElev_sfrCentered(X, Y_m) +botm = np.zeros(top.shape) +strthd = top - 10.0 + +# NPF parameters +k11 = 1 +ss = 0.00001 +sy = 0.20 +hani = 1 +laytyp = 1 + +# Package boundary conditions +viscref = 8.904e-4 + +# time params +steady = {0: True, 1: False} +transient = {0: False, 1: True} +nstp = [1, 20] +tsmult = [1, 1] +perlen = [1, 20] + +nouter, ninner = 1000, 300 +hclose, rclose, relax = 1e-3, 1e-4, 0.97 + +# Transport related parameters +initial_temperature = 35.0 # Initial temperature (unitless) +porosity = 0.20 # porosity (unitless) +K_therm = 2.0 # Thermal conductivity # ($W/m/C$) +rho_water = 1000 # Density of water ($kg/m^3$) +rho_solids = 2650 # Density of the aquifer material ($kg/m^3$) +C_p_w = 4180 # Heat Capacity of water ($J/kg/C$) +C_s = 880 # Heat capacity of the solids ($J/kg/C$) +D_m = K_therm / (porosity * rho_water * C_p_w) +rhob = (1 - porosity) * rho_solids # Bulk density ($kg/m^3$) +K_d = C_s / (rho_water * C_p_w) # Partitioning coefficient ($m^3/kg$) + +# +# MODFLOW 6 flopy GWF & GWT simulation object (sim) is returned +# + + +def build_model(idx, dir): + # Base simulation and model name and workspace + ws = dir + name = ex[idx] + + print("Building model...{}".format(name)) + + # generate names for each model + gwfname = "gwf-" + name + gwtname = "gwt-" + name + + sim = flopy.mf6.MFSimulation( + sim_name=name, sim_ws=ws, exe_name="mf6", version="mf6" + ) + + # Instantiating time discretization + tdis_rc = [] + for i in range(len(nstp)): + tdis_rc.append((perlen[i], nstp[i], tsmult[i])) + + flopy.mf6.ModflowTdis( + sim, nper=len(nstp), perioddata=tdis_rc, time_units=time_units + ) + + gwf = flopy.mf6.ModflowGwf( + sim, modelname=gwfname, save_flows=True, newtonoptions="newton" + ) + + # Instantiating solver + ims = flopy.mf6.ModflowIms( + sim, + print_option="ALL", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="cooley", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + filename="{}.ims".format(gwfname), + ) + sim.register_ims_package(ims, [gwfname]) + + # Instantiate discretization package + flopy.mf6.ModflowGwfdis( + gwf, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + ) + + # Instantiate node property flow package + flopy.mf6.ModflowGwfnpf( + gwf, + save_specific_discharge=True, + icelltype=1, # >0 means saturated thickness varies with computed head + k=k11, + ) + + # Instantiate storage package + flopy.mf6.ModflowGwfsto( + gwf, + save_flows=False, + iconvert=laytyp, + ss=ss, + sy=sy, + steady_state=steady, + transient=transient, + ) + + # Instantiate initial conditions package + flopy.mf6.ModflowGwfic(gwf, strt=strthd) + + # Instantiate viscosity package + if viscosity_on[idx]: + vsc_filerecord = "{}.vsc.bin".format(gwfname) + vsc_pd = [(0, 0.0, 20.0, gwtname, "TEMPERATURE")] + flopy.mf6.ModflowGwfvsc( + gwf, + viscref=viscref, + viscosity_filerecord=vsc_filerecord, + thermal_formulation="nonlinear", + thermal_a2=10.0, + thermal_a3=248.37, + thermal_a4=133.16, + nviscspecies=len(vsc_pd), + packagedata=vsc_pd, + pname="vsc", + filename="{}.vsc".format(gwfname), + ) + + # Instantiate output control package + 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"), ("BUDGET", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "LAST")], + ) + + # Instantiate recharge package + # total inflow 2000.0 on each side (4,000 total) + rech = np.zeros_like(top) + rech_rate_lo = 0.001 + rech_rate_hi = 0.015 + for i in np.arange(ncol): + rech[0, i] = rech_rate_lo + (rech_rate_hi - rech_rate_lo) / ncol * i + + rech[-1, :] = rech[0, :] + irch = np.zeros_like(rech) + irch = irch.astype(int) + temperature_array = np.ones_like(irch) * 15.0 + aux = {0: [temperature_array]} + flopy.mf6.ModflowGwfrcha( + gwf, + print_flows=True, + recharge=rech, + irch=irch, + auxiliary=["TEMPERATURE"], + aux=aux, + pname="RCHA-1", + filename="{}.rcha".format(gwfname), + ) + + # Instantiate evapotranspiration package + # ET rate is 0.003 everywhere in the model + evtr_lo = 0.0001 + evtr_hi = 0.012 + extdp_hi = 30 + extdp_lo = 10 + evtspd = [] + for i in np.arange(nrow): + for j in np.arange(ncol): + evtr = evtr_hi - (evtr_hi - evtr_lo) / ncol * j + extdp = extdp_hi - (extdp_hi - extdp_lo) / ncol * j + # cellid, surface, rate, depth, [pxdp], [petm], [petm0], [aux] + evtspd.append([(0, i, j), top[i, j], evtr, extdp, 1.0, 0.0]) + surf_rate_specified = True + flopy.mf6.ModflowGwfevt( + gwf, + print_flows=False, + surf_rate_specified=surf_rate_specified, + maxbound=nrow * ncol, + nseg=1, + stress_period_data=evtspd, + auxiliary="TEMPERATURE", + pname="EVT-1", + filename="{}.evt".format(gwfname), + ) + + # Instantiate streamflow routing package + + # Determine the middle row and store in rMid (account for 0-base) + rMid = nrow // 2 - 1 + # sfr data + nreaches = ncol + rlen = delr + rwid = 7.0 + roughness = 0.035 + rbth = 1.0 + rhk = 1.0 + strm_up = 254.899750 + strm_dn = 225.150250 + slope = (strm_up - strm_dn) / ((ncol - 1) * delr) + ustrf = 1.0 + ndv = 0 + strm_incision = 10 + viscaux = 1.111111111 + temperatureaux = 8.0 + + packagedata = [] + for irch in range(nreaches): + nconn = 1 + if 0 < irch < nreaches - 1: + nconn += 1 + rp = [ + irch, + (0, rMid, irch), + rlen, + rwid, + slope, + top[rMid, irch] - strm_incision, + rbth, + rhk, + roughness, + nconn, + ustrf, + ndv, + viscaux, + temperatureaux, + ] + packagedata.append(rp) + + connectiondata = [] + for irch in range(nreaches): + rc = [irch] + if irch > 0: + rc.append(irch - 1) + if irch < nreaches - 1: + rc.append(-(irch + 1)) + connectiondata.append(rc) + + inflow_loc = 0 + sfr_perioddata = [ + [inflow_loc, "inflow", 25000.0], + ] + sfr_perioddata = {0: sfr_perioddata} + + budpth = f"{gwfname}.sfr.cbc" + flopy.mf6.ModflowGwfsfr( + gwf, + save_flows=True, + print_stage=True, + print_flows=True, + print_input=False, + auxiliary=["VDUMMY", "TEMPERATURE"], + unit_conversion=1.486 * 86400, + budget_filerecord=budpth, + mover=False, + nreaches=nreaches, + packagedata=packagedata, + connectiondata=connectiondata, + perioddata=sfr_perioddata, + pname="SFR-1", + filename="{}.sfr".format(gwfname), + ) + + # Setup the GWT model for simulating heat transport + # ------------------------------------------------- + gwt = flopy.mf6.ModflowGwt(sim, modelname=gwtname) + + # Instantiating solver for GWT + imsgwt = flopy.mf6.ModflowIms( + sim, + print_option="ALL", + 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="{}.ims".format(gwtname), + ) + sim.register_ims_package(imsgwt, [gwtname]) + + # Instantiating DIS for GWT + flopy.mf6.ModflowGwtdis( + gwt, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + ) + + # Instantiate Mobile Storage and Transfer package + flopy.mf6.ModflowGwtmst( + gwt, + porosity=porosity, + sorption="linear", + bulk_density=rhob, + distcoef=K_d, + pname="MST-1", + filename="{}.mst".format(gwtname), + ) + + # Instantiate Transport Initial Conditions package + flopy.mf6.ModflowGwtic(gwt, strt=initial_temperature) + + # Instantiate Advection package + flopy.mf6.ModflowGwtadv(gwt, scheme="UPSTREAM") + + # Instantiate Dispersion package (also handles conduction) + flopy.mf6.ModflowGwtdsp(gwt, xt3d_off=True, diffc=D_m) + + # Instantiate Source/Sink Mixing package + sourcerecarray = [ + ("RCHA-1", "AUX", "TEMPERATURE"), + ("EVT-1", "AUX", "TEMPERATURE"), + ] + flopy.mf6.ModflowGwtssm(gwt, sources=sourcerecarray) + + # Instantiate Streamflow Transport package + sftpackagedata = [] + for irno in range(ncol): + t = (irno, 0.0) + sftpackagedata.append(t) + + sftperioddata = [(0, "STATUS", "CONSTANT"), (0, "CONCENTRATION", 8.0)] + + flopy.mf6.modflow.ModflowGwtsft( + gwt, + boundnames=False, + save_flows=True, + print_input=False, + print_flows=True, + print_concentration=True, + concentration_filerecord=gwtname + ".sft.bin", + budget_filerecord=gwtname + ".sft.bud", + packagedata=sftpackagedata, + reachperioddata=sftperioddata, + flow_package_auxiliary_name="TEMPERATURE", + flow_package_name="SFR-1", + pname="SFT-1", + filename="{}.sft".format(gwtname), + ) + + # Instantiate Output Control package for transport + flopy.mf6.ModflowGwtoc( + gwt, + concentration_filerecord="{}.ucn".format(gwtname), + saverecord=[("CONCENTRATION", "ALL")], + printrecord=[("CONCENTRATION", "LAST"), ("BUDGET", "LAST")], + filename="{}.oc".format(gwtname), + ) + + # Instantiate Gwf-Gwt Exchange package + flopy.mf6.ModflowGwfgwt( + sim, + exgtype="GWF6-GWT6", + exgmnamea=gwfname, + exgmnameb=gwtname, + filename="{}.gwfgwt".format(gwtname), + ) + + return sim, None + + +def eval_results(sim): + print("evaluating results...") + + # read flow results from model + name = ex[sim.idxsim] + gwfname = "gwf-" + name + + fname = gwfname + ".sfr.cbc" + fname = os.path.join(sim.simpath, fname) + assert os.path.isfile(fname) + budobj = flopy.utils.CellBudgetFile(fname, precision="double") + outbud = budobj.get_data(text=" GWF") + + # Establish known answer: + stored_ans_up = np.array( + [ + -381.07448455, + -380.78732368, + -380.49858508, + -380.20823421, + -379.91623521, + -379.62255085, + -379.32714247, + -379.02996987, + -378.73099123, + -378.43016302, + ] + ) + + stored_ans_dn = np.array( + [ + 87.34154005, + 92.49173569, + 98.11443969, + 104.34176415, + 111.36392723, + 119.46886139, + 129.12349417, + 141.16795665, + 157.38797791, + 182.65575269, + ] + ) + + if sim.idxsim == 0: + # convert np.array to list + no_vsc_bud_last = np.array(outbud[-1].tolist()) + + # sum up total losses and total gains in the first 10 reaches + # and the last 10 reaches + for i in np.arange(10): + # upper reaches + assert np.allclose( + abs(no_vsc_bud_last[i, 2]), abs(stored_ans_up[i]), atol=1e-4 + ), ( + "GW/SW not as expected in upper reaches of viscosity test " + "problem that uses SFR. This test keeps viscosity inactive " + "prior to activating VSC Package in next variant of this test " + "problem." + ) + + # lower reaches + assert np.allclose( + abs(no_vsc_bud_last[-(i + 1), 2]), + abs(stored_ans_dn[-(i + 1)]), + atol=1e-4, + ), ( + "GW/SW not as expected in lower reaches of viscosity test " + "problem that uses SFR. This test keeps viscosity inactive " + "prior to activating VSC Package in next variant of this test " + " problem." + ) + + elif sim.idxsim == 1: + with_vsc_bud_last = np.array(outbud[-1].tolist()) + + # sum up total losses and total gains in the first 10 reaches + # and the last 10 reaches + for i in np.arange(10): + # upper reaches + assert abs(stored_ans_up[i]) > abs(with_vsc_bud_last[i, 2]), ( + "GW/SW not as expected in upper reaches of viscosity test " + "problem that uses SFR. This test activates the VSC package that " + "should elicit a known relative change in the GW/SW exchange" + ) + + # lower reaches + assert abs(stored_ans_dn[-(i + 1)]) < abs( + with_vsc_bud_last[-(i + 1), 2] + ), ( + "GW/SW not as expected in lower reaches of viscosity test " + "problem that uses SFR. This test activates the VSC package that " + "should elicit a known relative change in the GW/SW exchange" + ) + + +# - No need to change any code below +@pytest.mark.parametrize( + "idx, dir", + list(enumerate(exdirs)), +) +def test_mf6model(idx, dir): + # initialize testing framework + test = testing_framework() + + # build the model + test.build_mf6_models(build_model, idx, dir) + + # run the test model + test.run_mf6(Simulation(dir, exfunc=eval_results, idxsim=idx)) + + +def main(): + # initialize testing framework + test = testing_framework() + + # run the test model + for idx, dir in enumerate(exdirs): + test.build_mf6_models(build_model, idx, dir) + sim = Simulation(dir, exfunc=eval_results, idxsim=idx) + test.run_mf6(sim) + + +if __name__ == "__main__": + # print message + print(f"standalone run of {os.path.basename(__file__)}") + + # run main routine + main() diff --git a/autotest/test_gwf_vsc04_lak.py b/autotest/test_gwf_vsc04_lak.py new file mode 100644 index 00000000000..54757eca6e5 --- /dev/null +++ b/autotest/test_gwf_vsc04_lak.py @@ -0,0 +1,822 @@ +# Simple single lake model. Lake cut into top two layers of a 5 layer +# model. Model is loosely based on the first example problem in +# Merritt and Konikow (2000) which also is one of the MT3D-USGS test +# problems. This test developed to isolate lake-aquifer interaction; +# no SFR or other advanced packages. Problem set up to have groundwater +# pass through the lake: gw inflow on the left side, gw outflow on the +# right side of the lake. Uses constant stage boundary in the lake to +# ensure desired flow conditions for testing budget changes with and +# without VSC active. +# +# starting groundwater temperature: 30.0 +# left chd boundary inflow temperature: 30.0 +# starting lake temperature: 4.0 +# + +import os +import sys + +import numpy as np +import pytest + +try: + import flopy +except: + msg = "Error. FloPy package is not available.\n" + msg += "Try installing using the following command:\n" + msg += " pip install flopy" + raise Exception(msg) + +from framework import testing_framework +from simulation import Simulation + +ex = ["no-vsc04-lak", "vsc04-lak"] +viscosity_on = [False, True] +exdirs = [] +for s in ex: + exdirs.append(os.path.join("temp", s)) + +# Model units +length_units = "m" +time_units = "days" + +# model domain and grid definition +delr = [ + 76.2, + 304.8, + 304.8, + 304.8, + 304.8, + 304.8, + 152.4, + 152.4, + 152.4, + 152.4, + 152.4, + 304.8, + 304.8, + 304.8, + 304.8, + 304.8, + 76.2, +] + +delc = [ + 76.2, + 304.8, + 304.8, + 304.8, + 304.8, + 304.8, + 152.4, + 152.4, + 152.4, + 152.4, + 152.4, + 304.8, + 304.8, + 304.8, + 304.8, + 304.8, + 76.2, +] + +fixedstrthds = [ + 35.052, + 34.9267, + 34.7216, + 34.5062, + 34.2755, + 34.0237, + 33.8143, + 33.6657, + 33.5077, + 33.3394, + 33.1599, + 32.8728, + 32.4431, + 31.9632, + 31.4353, + 30.8627, + 30.48, +] + +nrow = len(delc) +ncol = len(delr) +top = np.ones((nrow, ncol)) * 35.6616 +bot1 = np.ones_like(top) * 32.6136 +bot2 = np.ones_like(top) * 29.5656 +bot3 = np.ones_like(top) * 26.5176 +bot4 = np.ones_like(top) * 23.4696 +bot5 = np.ones_like(top) * 20.4216 +botm = np.array([bot1, bot2, bot3, bot4, bot5]) +nlay = botm.shape[0] +ibound = np.ones_like(botm) + +# deactive gw cells where lake cells are active +ibound[0, 6:11, 6:11] = 0 # layer 1 +ibound[1, 7:10, 7:10] = 0 # layer 2 + +strthd = np.zeros_like(ibound) +for j in np.arange(ncol): + strthd[:, :, j] = fixedstrthds[j] + +# setup lake array +lakibnd = np.zeros_like(ibound) +lakibnd[0] = 1 - ibound[0] # layer 1 +lakibnd[1] = 1 - ibound[1] # layer 2 + +# NPF parameters +k11 = 9.144 # = 30 ft/day +k33 = 0.9144 # = 30 ft/day +ss = 3e-4 +sy = 0.20 +hani = 1 +laytyp = 1 + +# Package boundary conditions +chdl = 35.052 +chdr = 30.48 +viscref = 8.904e-4 + +# time params +transient = {0: True} +nstp = [100] +tsmult = [1.02] +perlen = [5000] + +# solver params +nouter, ninner = 1000, 300 +hclose, rclose, relax = 1e-3, 1e-4, 0.97 + +# Transport related parameters +al = 1 # longitudinal dispersivity ($m$) +ath1 = al # horizontal transverse dispersivity +atv = al # vertical transverse dispersivity +mixelm = 0 # Upstream vs TVD (Upstream selected) +initial_temperature = 35.0 # Initial temperature (unitless) +porosity = 0.20 # porosity (unitless) +K_therm = 2.0 # Thermal conductivity # ($W/m/C$) +rho_water = 1000 # Density of water ($kg/m^3$) +rho_solids = 2650 # Density of the aquifer material ($kg/m^3$) +C_p_w = 4180 # Heat Capacity of water ($J/kg/C$) +C_s = 880 # Heat capacity of the solids ($J/kg/C$) +D_m = K_therm / (porosity * rho_water * C_p_w) +rhob = (1 - porosity) * rho_solids # Bulk density ($kg/m^3$) +K_d = C_s / (rho_water * C_p_w) # Partitioning coefficient ($m^3/kg$) +leftTemp = 30.0 # Temperature of inflow from left constant head ($C$) + +# Viscosity related parameters +tviscref = 20.0 + +# +# MODFLOW 6 flopy GWF & GWT simulation object (sim) is returned +# + + +def build_model(idx, dir): + global lak_lkup_dict + + # Base simulation and model name and workspace + ws = dir + name = ex[idx] + + print("Building model...{}".format(name)) + + # generate names for each model + gwfname = "gwf-" + name + gwtname = "gwt-" + name + + sim = flopy.mf6.MFSimulation( + sim_name=name, sim_ws=ws, exe_name="mf6", version="mf6" + ) + + tdis_rc = [] + for i in range(len(nstp)): + tdis_rc.append((perlen[i], nstp[i], tsmult[i])) + + flopy.mf6.ModflowTdis( + sim, nper=len(nstp), perioddata=tdis_rc, time_units=time_units + ) + + gwf = flopy.mf6.ModflowGwf( + sim, modelname=gwfname, save_flows=True, newtonoptions="newton" + ) + + # Instantiating solver + ims = flopy.mf6.ModflowIms( + sim, + print_option="ALL", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="cooley", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + filename="{}.ims".format(gwfname), + ) + sim.register_ims_package(ims, [gwfname]) + + # Instantiate discretization package + flopy.mf6.ModflowGwfdis( + gwf, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + idomain=ibound, + filename="{}.dis".format(gwfname), + ) + + # Instantiate node property flow package + flopy.mf6.ModflowGwfnpf( + gwf, + save_specific_discharge=True, + icelltype=1, # >0 means saturated thickness varies with computed head + k=k11, + k33=k33, + ) + + # Instantiate storage package + flopy.mf6.ModflowGwfsto( + gwf, + save_flows=False, + iconvert=laytyp, + ss=ss, + sy=sy, + transient=transient, + ) + + # Instantiate initial conditions package + flopy.mf6.ModflowGwfic(gwf, strt=strthd) + + # Instantiate viscosity package + if viscosity_on[idx]: + vsc_filerecord = "{}.vsc.bin".format(gwfname) + vsc_pd = [(0, 0.0, tviscref, gwtname, "TEMPERATURE")] + flopy.mf6.ModflowGwfvsc( + gwf, + viscref=viscref, + viscosity_filerecord=vsc_filerecord, + thermal_formulation="nonlinear", + thermal_a2=10.0, + thermal_a3=248.37, + thermal_a4=133.16, + nviscspecies=len(vsc_pd), + packagedata=vsc_pd, + pname="vsc", + filename="{}.vsc".format(gwfname), + ) + + # Instantiate output control package + flopy.mf6.ModflowGwfoc( + gwf, + budget_filerecord=f"{gwfname}.cbc", + head_filerecord=f"{gwfname}.hds", + headprintrecord=[("COLUMNS", 17, "WIDTH", 15, "DIGITS", 6, "GENERAL")], + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "LAST")], + ) + + # Instantiate constant head package + # (for driving gw flow from left to right) + chdlistl = [] + chdlistr = [] + for k in np.arange(nlay): + for i in np.arange(nrow): + # left side + if botm[k, i, 0] <= chdl: + chdlistl.append([(k, i, 0), chdl, leftTemp]) + # right side + if botm[k, i, -1] <= chdr: + chdlistr.append([(k, i, ncol - 1), chdr, 10.0]) + + flopy.mf6.ModflowGwfchd( + gwf, + stress_period_data=chdlistl, + print_input=True, + print_flows=True, + save_flows=False, + pname="CHD-L", + auxiliary="TEMPERATURE", + filename=f"{gwfname}.left.chd", + ) + + flopy.mf6.ModflowGwfchd( + gwf, + stress_period_data=chdlistr, + print_input=True, + print_flows=True, + save_flows=False, + pname="CHD-R", + auxiliary="TEMPERATURE", + filename=f"{gwfname}.right.chd", + ) + + # Instantiate lake package + lakeconnectiondata = [] + nlakecon = [0] # Expand this to [0, 0, ...] for each additional lake + ilakconn = -1 + lak_leakance = 0.1 + lak_lkup_dict = {} + for k in [0, 1]: + for i in range(nrow): + for j in range(ncol): + if lakibnd[k, i, j] == 0: + continue + else: + ilak = int(lakibnd[k, i, j] - 1) + # back + if i > 0: + if ( + lakibnd[k, i - 1, j] == 0 + and ibound[k, i - 1, j] == 1 + ): + ilakconn += 1 + # by setting belev==telev, MF6 will automatically + # re-assign elevations based on cell dimensions + h = [ + ilak, # + ilakconn, # + (k, i - 1, j), # + "horizontal", # + lak_leakance, # + 0.0, # + 0.0, # + delc[i] / 2.0, # + delr[j], # + ] + lakeconnectiondata.append(h) + lak_lkup_dict.update({ilakconn: (k, i, j)}) + + # left + if j > 0: + if ( + lakibnd[k, i, j - 1] == 0 + and ibound[k, i, j - 1] == 1 + ): + ilakconn += 1 + h = [ + ilak, + ilakconn, + (k, i, j - 1), + "horizontal", + lak_leakance, + 0.0, + 0.0, + delr[j] / 2.0, + delc[i], + ] + lakeconnectiondata.append(h) + lak_lkup_dict.update({ilakconn: (k, i, j)}) + + # right + if j < ncol - 1: + if ( + lakibnd[k, i, j + 1] == 0 + and ibound[k, i, j + 1] == 1 + ): + ilakconn += 1 + h = [ + ilak, + ilakconn, + (k, i, j + 1), + "horizontal", + lak_leakance, + 0.0, + 0.0, + delr[j] / 2.0, + delc[i], + ] + lakeconnectiondata.append(h) + lak_lkup_dict.update({ilakconn: (k, i, j)}) + + # front + if i < nrow - 1: + if ( + lakibnd[k, i + 1, j] == 0 + and ibound[k, i + 1, j] == 1 + ): + ilakconn += 1 + h = [ + ilak, + ilakconn, + (k, i + 1, j), + "horizontal", + lak_leakance, + 0.0, + 0.0, + delc[i] / 2.0, + delr[j], + ] + lakeconnectiondata.append(h) + lak_lkup_dict.update({ilakconn: (k, i, j)}) + + # vertical + if lakibnd[k, i, j] == 1 and ibound[k + 1, i, j] == 1: + ilakconn += 1 + v = [ + ilak, + ilakconn, + (k + 1, i, j), + "vertical", + lak_leakance, + 0.0, + 0.0, + 0.0, + 0.0, + ] + lakeconnectiondata.append(v) + lak_lkup_dict.update({ilakconn: (k, i, j)}) + + strtStg = 33.75 + lakpackagedata = [[0, strtStg, len(lakeconnectiondata), 4.0, "lake1"]] + lak_pkdat_dict = {"filename": "lak_pakdata.in", "data": lakpackagedata} + + lakeperioddata = { + 0: [ + (0, "STATUS", "CONSTANT"), # RAINFALL 0.005 & 0.00504739035 + (0, "STAGE", 33.5), + ] + } + + lak_obs = { + "{}.lakeobs".format(gwfname): [ + ("lakestage", "stage", "lake1"), + ("gwexchng", "lak", "lake1"), + ] + } + lak = flopy.mf6.ModflowGwflak( + gwf, + auxiliary="TEMPERATURE", + time_conversion=86400.0, + print_stage=True, + print_flows=True, + budget_filerecord=gwfname + ".lak.bud", + length_conversion=1.0, + mover=False, + pname="LAK-1", + boundnames=True, + nlakes=len(lakpackagedata), + noutlets=0, + packagedata=lak_pkdat_dict, + connectiondata=lakeconnectiondata, + perioddata=lakeperioddata, + observations=lak_obs, + filename="{}.lak".format(gwfname), + ) + + # pull in the tabfile defining the lake stage, vol, & surface area + fname = os.path.join("data", "vsc04-laktab", "stg-vol-surfarea.dat") + tabinput = [] + with open(fname, "r") as f: + # peel off the hdr line + hdr = next(f) + for line in f: + m_arr = line.strip().split(",") + # , , , + tabinput.append([float(m_arr[0]), m_arr[1], m_arr[2]]) + + tab6_filename = "{}.laktab".format(gwfname) + flopy.mf6.ModflowUtllaktab( + gwf, + nrow=len(tabinput), + ncol=3, + table=tabinput, + filename=tab6_filename, + pname="LAK_tab", + parent_file=lak, + ) + + # create gwt model + # ---------------- + gwt = flopy.mf6.ModflowGwt( + sim, modelname=gwtname, model_nam_file="{}.nam".format(gwtname) + ) + gwt.name_file.save_flows = True + + imsgwt = flopy.mf6.ModflowIms( + sim, + print_option="ALL", + 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]) + + # Instantiating MODFLOW 6 transport discretization package + flopy.mf6.ModflowGwtdis( + gwt, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + idomain=ibound, + filename="{}.dis".format(gwtname), + ) + + # Instantiating MODFLOW 6 transport initial concentrations + strtconc = leftTemp + flopy.mf6.ModflowGwtic( + gwt, strt=strtconc, filename="{}.ic".format(gwtname) + ) + + # Instantiate mobile storage and transfer package + sto = flopy.mf6.ModflowGwtmst( + gwt, porosity=porosity, filename=f"{gwtname}.sto" + ) + + # Instantiating MODFLOW 6 transport advection package + if mixelm == 0: + scheme = "UPSTREAM" + elif mixelm == -1: + scheme = "TVD" + else: + raise Exception() + + # Instantiate advection package + flopy.mf6.ModflowGwtadv( + gwt, scheme=scheme, filename="{}.adv".format(gwtname) + ) + + # Instantiate dispersion package + flopy.mf6.ModflowGwtdsp( + gwt, alh=al, ath1=ath1, atv=atv, filename="{}.dsp".format(gwtname) + ) + + # Instantiate source/sink mixing package + sourcerecarray = [ + ("CHD-L", "AUX", "TEMPERATURE"), + ("CHD-R", "AUX", "TEMPERATURE"), + ] + flopy.mf6.ModflowGwtssm( + gwt, sources=sourcerecarray, filename=f"{gwtname}.ssm" + ) + + # Instantiating MODFLOW 6 transport output control package + flopy.mf6.ModflowGwtoc( + gwt, + budget_filerecord="{}.cbc".format(gwtname), + concentration_filerecord="{}.ucn".format(gwtname), + concentrationprintrecord=[ + ("COLUMNS", 17, "WIDTH", 15, "DIGITS", 6, "GENERAL") + ], + saverecord=[("CONCENTRATION", "ALL"), ("BUDGET", "ALL")], + printrecord=[("CONCENTRATION", "ALL"), ("BUDGET", "ALL")], + filename="{}.oc".format(gwtname), + ) + + # Instantiating MODFLOW 6 lake transport (lkt) package + lktpackagedata = [(0, 4.0, "lake1")] + + lktperioddata = {0: [(0, "STATUS", "CONSTANT"), (0, "CONCENTRATION", 4.0)]} + + # note: for specifying lake number, use fortran indexing! + lkt_obs = { + "{}.lakobs".format(gwtname): [ + ("resTemp", "concentration", 1), + ("resGwMassExchng", "lkt", "lake1"), + ] + } + + flopy.mf6.ModflowGwtlkt( + gwt, # Set time_conversion for use with Manning's eqn. + flow_package_name="LAK-1", + flow_package_auxiliary_name="TEMPERATURE", + budget_filerecord=gwtname + ".lkt.bud", + boundnames=True, + save_flows=True, + print_input=True, + print_flows=False, + print_concentration=True, + packagedata=lktpackagedata, + lakeperioddata=lktperioddata, + observations=lkt_obs, + pname="LKT-1", + filename="{}.lkt".format(gwtname), + ) + + # GWF-GWT exchange + flopy.mf6.ModflowGwfgwt( + sim, + exgtype="GWF6-GWT6", + exgmnamea=gwfname, + exgmnameb=gwtname, + filename=f"{name}.gwfgwt", + ) + + return sim, None + + +def eval_results(sim): + print("evaluating results...") + + # read flow results from model + name = ex[sim.idxsim] + gwfname = "gwf-" + name + + fname = gwfname + ".lak.bud" + fname = os.path.join(sim.simpath, fname) + assert os.path.isfile(fname) + budobj = flopy.utils.CellBudgetFile(fname, precision="double") + outbud = budobj.get_data(text=" GWF") + + # Establish known answer: + stored_ans = np.array( + [ + [1.0, 9.20e1, 7.34424130e-01, 1.15181189e06], + [1.0, 1.08e2, 2.44117249e00, 1.15181189e06], + [1.0, 3.98e2, 2.19216490e01, 1.15181189e06], + [1.0, 9.30e1, -6.78268488e-02, 1.15181189e06], + [1.0, 3.99e2, 3.17042406e-01, 1.15181189e06], + [1.0, 9.40e1, -7.47709994e-01, 1.15181189e06], + [1.0, 4.00e2, -6.88938281e00, 1.15181189e06], + [1.0, 9.50e1, -1.51336530e00, 1.15181189e06], + [1.0, 4.01e2, -1.57614834e01, 1.15181189e06], + [1.0, 9.60e1, -2.54095715e00, 1.15181189e06], + [1.0, 1.14e2, -3.95961161e00, 1.15181189e06], + [1.0, 4.02e2, -5.60853100e01, 1.15181189e06], + [1.0, 1.25e2, 2.35138538e00, 1.15181189e06], + [1.0, 4.15e2, 1.69275311e01, 1.15181189e06], + [1.0, 1.31e2, -3.66648779e00, 1.15181189e06], + [1.0, 4.19e2, -3.45225854e01, 1.15181189e06], + [1.0, 1.42e2, 2.32550672e00, 1.15181189e06], + [1.0, 4.32e2, 1.65405908e01, 1.15181189e06], + [1.0, 1.48e2, -3.58087615e00, 1.15181189e06], + [1.0, 4.36e2, -3.27154545e01, 1.15181189e06], + [1.0, 1.59e2, 2.35138505e00, 1.15181189e06], + [1.0, 4.49e2, 1.69275277e01, 1.15181189e06], + [1.0, 1.65e2, -3.66648777e00, 1.15181189e06], + [1.0, 4.53e2, -3.45225840e01, 1.15181189e06], + [1.0, 1.76e2, 2.44117266e00, 1.15181189e06], + [1.0, 1.94e2, 7.34424819e-01, 1.15181189e06], + [1.0, 4.66e2, 2.19216459e01, 1.15181189e06], + [1.0, 1.95e2, -6.78264346e-02, 1.15181189e06], + [1.0, 4.67e2, 3.17038149e-01, 1.15181189e06], + [1.0, 1.96e2, -7.47709250e-01, 1.15181189e06], + [1.0, 4.68e2, -6.88938656e00, 1.15181189e06], + [1.0, 1.97e2, -1.51336458e00, 1.15181189e06], + [1.0, 4.69e2, -1.57614826e01, 1.15181189e06], + [1.0, 1.82e2, -3.95961151e00, 1.15181189e06], + [1.0, 1.98e2, -2.54095654e00, 1.15181189e06], + [1.0, 4.70e2, -5.60853022e01, 1.15181189e06], + [1.0, 3.99e2, 4.03508517e-03, 1.15181189e06], + [1.0, 4.15e2, 2.15441304e-01, 1.15181189e06], + [1.0, 7.05e2, 8.25215117e-01, 1.15181189e06], + [1.0, 4.00e2, -8.76830539e-02, 1.15181189e06], + [1.0, 7.06e2, -3.90793309e-01, 1.15181189e06], + [1.0, 4.01e2, -2.00600698e-01, 1.15181189e06], + [1.0, 4.19e2, -4.39378360e-01, 1.15181189e06], + [1.0, 7.07e2, -2.43955302e00, 1.15181189e06], + [1.0, 4.32e2, 2.10516610e-01, 1.15181189e06], + [1.0, 7.22e2, 8.37390920e-01, 1.15181189e06], + [1.0, 7.23e2, -6.85716153e-02, 1.15181189e06], + [1.0, 4.36e2, -4.16378511e-01, 1.15181189e06], + [1.0, 7.24e2, -1.73090130e00, 1.15181189e06], + [1.0, 4.49e2, 2.15441262e-01, 1.15181189e06], + [1.0, 4.67e2, 4.03503099e-03, 1.15181189e06], + [1.0, 7.39e2, 8.25212943e-01, 1.15181189e06], + [1.0, 4.68e2, -8.76831016e-02, 1.15181189e06], + [1.0, 7.40e2, -3.90795105e-01, 1.15181189e06], + [1.0, 4.53e2, -4.39378341e-01, 1.15181189e06], + [1.0, 4.69e2, -2.00600688e-01, 1.15181189e06], + [1.0, 7.41e2, -2.43955388e00, 1.15181189e06], + ] + ) + + # talley some flows on the left and right sides of the lake for comparison + # test + left_chk_ans = [] + right_chk_ans = [] + left_chk_no_vsc = [] + right_chk_no_vsc = [] + left_chk_with_vsc = [] + right_chk_with_vsc = [] + + if sim.idxsim == 0: + no_vsc_bud_last = np.array(outbud[-1].tolist()) + no_vsc_bud_np = np.array(no_vsc_bud_last.tolist()) + + for idx in np.arange(stored_ans.shape[0]): + k, i, j = lak_lkup_dict[idx] + + # left side of lake + if j < 7: + if no_vsc_bud_np[idx, 2] > 0 and stored_ans[idx, 2] > 0: + left_chk_no_vsc.append(no_vsc_bud_np[idx, 2]) + left_chk_ans.append(stored_ans[idx, 2]) + + # right side of lake + if j > 9: + if no_vsc_bud_np[idx, 2] < 0 and stored_ans[idx, 2] < 0: + right_chk_no_vsc.append(no_vsc_bud_np[idx, 2]) + right_chk_ans.append(stored_ans[idx, 2]) + + # Check that all the flows entering the lak in the 'with vsc' model are greater + # than their 'no vsc' counterpart + assert np.allclose( + np.array(left_chk_ans), np.array(left_chk_no_vsc), atol=1e-3 + ), ( + "Lake inflow in no-VSC LAK simulation do not match established " + "solution." + ) + + # Check that all the flows leaving the lak in the 'with vsc' model are less + # than their 'no vsc' counterpart (keep in mind values are negative, which + # affects how the comparison is made) + assert np.allclose( + np.array(right_chk_ans), np.array(right_chk_no_vsc), atol=1e-3 + ), ( + "Lake outflow in no-VSC LAK simulation do not match established " + "solution." + ) + + elif sim.idxsim == 1: + with_vsc_bud_last = np.array(outbud[-1].tolist()) + with_vsc_bud_np = np.array(with_vsc_bud_last.tolist()) + + for idx in np.arange(stored_ans.shape[0]): + k, i, j = lak_lkup_dict[idx] + + # left side of lake + if j < 7: + if stored_ans[idx, 2] > 0 and with_vsc_bud_np[idx, 2] > 0: + left_chk_no_vsc.append(stored_ans[idx, 2]) + left_chk_with_vsc.append(with_vsc_bud_np[idx, 2]) + + # right side of lake + if j > 9: + if stored_ans[idx, 2] < 0 and with_vsc_bud_np[idx, 2] < 0: + right_chk_no_vsc.append(stored_ans[idx, 2]) + right_chk_with_vsc.append(with_vsc_bud_np[idx, 2]) + + # Check that all the flows entering the lak in the 'with vsc' model are greater + # than their 'no vsc' counterpart + assert np.greater( + np.array(left_chk_with_vsc), np.array(left_chk_no_vsc) + ).all(), ( + "Lake inflow did no increase with VSC turned on and should have." + ) + + # Check that all the flows leaving the lak in the 'with vsc' model are less + # than their 'no vsc' counterpart (keep in mind values are negative, which + # affects how the comparison is made) + assert np.greater( + np.array(right_chk_with_vsc), np.array(right_chk_no_vsc) + ).all(), ( + "Lake outflow did no decrease with VSC turned on and should have." + ) + + +# - No need to change any code below +@pytest.mark.parametrize( + "idx, dir", + list(enumerate(exdirs)), +) +def test_mf6model(idx, dir): + # initialize testing framework + test = testing_framework() + + # build the model + test.build_mf6_models(build_model, idx, dir) + + # run the test model + test.run_mf6(Simulation(dir, exfunc=eval_results, idxsim=idx)) + + +def main(): + # initialize testing framework + test = testing_framework() + + # run the test model + for idx, dir in enumerate(exdirs): + test.build_mf6_models(build_model, idx, dir) + sim = Simulation(dir, exfunc=eval_results, idxsim=idx) + test.run_mf6(sim) + + +if __name__ == "__main__": + # print message + print(f"standalone run of {os.path.basename(__file__)}") + + # run main routine + main() diff --git a/autotest/test_gwf_vsc05_hfb.py b/autotest/test_gwf_vsc05_hfb.py new file mode 100644 index 00000000000..6846226fd78 --- /dev/null +++ b/autotest/test_gwf_vsc05_hfb.py @@ -0,0 +1,415 @@ +# ## Test problem for VSC and HFB +# +# Uses constant head and general-head boundaries on the left and right +# sides of a 10 row by 10 column by 1 layer model to drive flow from left to +# right. Tests that a horizontal flow barrier accounts for changes in +# viscosity when temperature is simulated. Barrier is between middle two +# columns, but only cuts across the bottom 5 rows. +# Model 1: VSC inactive, uses a higher speified K that matches what the VSC +# package will come up with +# Model 2: VSC active, uses a lower K so that when VSC is applied, resulting +# K's match model 1 and should result in the same flows across the +# model domain +# Model 3: VSC inactive, uses the lower K of model 2 and checks that flows +# in model 3 are indeed lower than in model 2 when turning VSC off. +# Model simulates hot groundwater with lower viscosity resulting in +# more gw flow through the model domain.Flows that are checked are +# the row-wise flows between columns 5 and 6 (e.g., cell 5 to 6, 15 +# to 16, etc.) +# + +# Imports + +import os +import sys + +import numpy as np +import pytest + +try: + import flopy +except: + msg = "Error. FloPy package is not available.\n" + msg += "Try installing using the following command:\n" + msg += " pip install flopy" + raise Exception(msg) + +from framework import testing_framework +from simulation import Simulation + +hyd_cond = [1205.49396942506, 864.0] # Hydraulic conductivity (m/d) +ex = ["no-vsc05-hfb", "vsc05-hfb", "no-vsc05-k"] +viscosity_on = [False, True, False] +hydraulic_conductivity = [hyd_cond[0], hyd_cond[1], hyd_cond[1]] +exdirs = [] +for s in ex: + exdirs.append(os.path.join("temp", s)) + +# Model units + +length_units = "cm" +time_units = "seconds" + +# Table of model parameters + +nper = 1 # Number of periods +nstp = 10 # Number of time steps +perlen = 10 # Simulation time length ($d$) +nlay = 1 # Number of layers +nrow = 10 # Number of rows +ncol = 10 # Number of columns +delr = 1.0 # Column width ($m$) +delc = 1.0 # Row width ($m$) +delv = 1.0 # Layer thickness +top = 1.0 # Top of the model ($m$) +initial_temperature = 35.0 # Initial temperature (unitless) +porosity = 0.26 # porosity (unitless) +K_therm = 2.0 # Thermal conductivity # ($W/m/C$) +rho_water = 1000 # Density of water ($kg/m^3$) +rho_solids = 2650 # Density of the aquifer material ($kg/m^3$) +C_p_w = 4180 # Heat Capacity of water ($J/kg/C$) +C_s = 880 # Heat capacity of the solids ($J/kg/C$) +D_m = K_therm / (porosity * rho_water * C_p_w) +rhob = (1 - porosity) * rho_solids # Bulk density ($kg/m^3$) +K_d = C_s / (rho_water * C_p_w) # Partitioning coefficient ($m^3/kg$) +inflow = 5.7024 # ($m^3/d$) + +botm = [top - k * delv for k in range(1, nlay + 1)] + +nouter, ninner = 100, 300 +hclose, rclose, relax = 1e-10, 1e-6, 0.97 + +# +# MODFLOW 6 flopy GWF simulation object (sim) is returned +# + + +def build_model(idx, dir): + # Base simulation and model name and workspace + ws = dir + name = ex[idx] + + print("Building model...{}".format(name)) + + # generate names for each model + gwfname = "gwf-" + name + gwtname = "gwt-" + name + + sim = flopy.mf6.MFSimulation( + sim_name=name, version="mf6", exe_name="mf6", sim_ws=ws + ) + + # Instantiating time discretization + tdis_ds = ((perlen, nstp, 1.0),) + flopy.mf6.ModflowTdis( + sim, nper=nper, perioddata=tdis_ds, time_units=time_units + ) + gwf = flopy.mf6.ModflowGwf(sim, modelname=gwfname, save_flows=True) + + # Instantiating solver + ims = flopy.mf6.ModflowIms( + sim, + print_option="ALL", + 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="{}.ims".format(gwfname), + ) + sim.register_ims_package(ims, [gwfname]) + + # Instantiating DIS + flopy.mf6.ModflowGwfdis( + gwf, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + ) + + # Instantiating NPF + flopy.mf6.ModflowGwfnpf( + gwf, + save_specific_discharge=True, + icelltype=0, + k=hydraulic_conductivity[idx], + ) + flopy.mf6.ModflowGwfic(gwf, strt=0.0) + + # Instantiating VSC + if viscosity_on[idx]: + # Instantiate viscosity (VSC) package + vsc_filerecord = "{}.vsc.bin".format(gwfname) + vsc_pd = [(0, 0.0, 20.0, gwtname, "temperature")] + flopy.mf6.ModflowGwfvsc( + gwf, + viscref=8.904e-4, + viscosity_filerecord=vsc_filerecord, + thermal_formulation="nonlinear", + thermal_a2=10.0, + thermal_a3=248.37, + thermal_a4=133.16, + nviscspecies=len(vsc_pd), + packagedata=vsc_pd, + pname="vsc", + filename="{}.vsc".format(gwfname), + ) + + # Instantiating CHD (leftside, "inflow" boundary) + chdspd = [[(0, i, 0), 2.0, initial_temperature] for i in range(nrow)] + flopy.mf6.ModflowGwfchd( + gwf, + stress_period_data=chdspd, + pname="CHD-1", + auxiliary="temperature", + ) + + # Instantiating GHB (rightside, "outflow" boundary) + ghbcond = hydraulic_conductivity[idx] * delv * delc / (0.5 * delr) + ghbspd = [ + [(0, i, ncol - 1), top, ghbcond, initial_temperature] + for i in range(nrow) + ] + flopy.mf6.ModflowGwfghb( + gwf, + stress_period_data=ghbspd, + pname="GHB-1", + auxiliary="temperature", + ) + + # Instantiate Horizontal Flow-Barrier (HFB) package + # Barrier present between middle two columns of the model domain, but only + # in rows 6-10. Remember that the hydraulic characteristic is the barrier + # hydraulic conductivity divided by the width of the horizontal-flow + # barrier. Assuming a barrier width of 10 cm (0.1 m) and desire to have + # the barrier's K be 1/10th of the aquifer hydraulic conductivity. + hfbspd = [] + K = 0.1 * hydraulic_conductivity[idx] + for i in np.arange(5, 10, 1): + hfbspd.append(((0, i, 4), (0, i, 5), K)) + flopy.mf6.ModflowGwfhfb( + gwf, + print_input=True, + maxhfb=len(hfbspd), + stress_period_data=hfbspd, + pname="HFB-1", + ) + + # Instatiating OC + head_filerecord = "{}.hds".format(gwfname) + budget_filerecord = "{}.bud".format(gwfname) + flopy.mf6.ModflowGwfoc( + gwf, + head_filerecord=head_filerecord, + budget_filerecord=budget_filerecord, + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) + + # Setup the GWT model for simulating heat transport + # ------------------------------------------------- + gwt = flopy.mf6.ModflowGwt(sim, modelname=gwtname) + + # Instantiating solver for GWT + imsgwt = flopy.mf6.ModflowIms( + sim, + print_option="ALL", + 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="{}.ims".format(gwtname), + ) + sim.register_ims_package(imsgwt, [gwtname]) + + # Instantiating DIS for GWT + flopy.mf6.ModflowGwtdis( + gwt, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + ) + + # Instantiating MST for GWT + flopy.mf6.ModflowGwtmst( + gwt, + porosity=porosity, + sorption="linear", + bulk_density=rhob, + distcoef=K_d, + pname="MST-1", + filename="{}.mst".format(gwtname), + ) + + # Instantiating IC for GWT + flopy.mf6.ModflowGwtic(gwt, strt=initial_temperature) + + # Instantiating ADV for GWT + flopy.mf6.ModflowGwtadv(gwt, scheme="UPSTREAM") + + # Instantiating DSP for GWT + flopy.mf6.ModflowGwtdsp(gwt, xt3d_off=True, diffc=D_m) + + # Instantiating SSM for GWT + sourcerecarray = [ + ("CHD-1", "AUX", "TEMPERATURE"), + ("GHB-1", "AUX", "TEMPERATURE"), + ] + flopy.mf6.ModflowGwtssm(gwt, sources=sourcerecarray) + + # Instantiating OC for GWT + flopy.mf6.ModflowGwtoc( + gwt, + concentration_filerecord="{}.ucn".format(gwtname), + saverecord=[("CONCENTRATION", "ALL")], + printrecord=[("CONCENTRATION", "LAST"), ("BUDGET", "LAST")], + ) + + # Instantiating GWF/GWT Exchange + flopy.mf6.ModflowGwfgwt( + sim, exgtype="GWF6-GWT6", exgmnamea=gwfname, exgmnameb=gwtname + ) + + return sim, None + + +def eval_results(sim): + print("evaluating results...") + + # read flow results from model + name = ex[sim.idxsim] + gwfname = "gwf-" + name + sim1 = flopy.mf6.MFSimulation.load(sim_ws=sim.simpath, load_only=["dis"]) + gwf = sim1.get_model(gwfname) + + # Get grid data + grdname = gwfname + ".dis.grb" + bgf = flopy.mf6.utils.MfGrdFile(os.path.join(sim.simpath, grdname)) + ia, ja = bgf.ia, bgf.ja + + fname = gwfname + ".bud" + fname = os.path.join(sim.simpath, fname) + assert os.path.isfile(fname) + budobj = flopy.utils.CellBudgetFile(fname, precision="double") + outbud = budobj.get_data(text=" FLOW-JA-FACE")[-1].squeeze() + + # Establish known answer for the "with viscosity" variant: + stored_ans = np.array( + [ + [4, 5, 131.03196884892344], + [14, 15, 133.1834658429856], + [24, 25, 139.31716925610493], + [34, 35, 156.14497435040056], + [44, 45, 209.1055337693415], + [54, 55, 36.91267872240113], + [64, 65, 46.16474722642168], + [74, 75, 51.2708505192076], + [84, 85, 54.04369740428511], + [94, 95, 55.27469944201896], + ] + ) + + # Look at flow entering the left face for the cells in the 6th (1-based) column + cells = [gwf.modelgrid.get_node([(0, i, 5)])[0] for i in np.arange(nrow)] + + vals_to_store = [] # Will always have 10 vals, 1 per row + + # Note that the layer, row, column indices will be zero-based + for celln in cells: + for ipos in range(ia[celln] + 1, ia[celln + 1]): + cellm = ja[ipos] + if cellm == celln - 1: + vals_to_store.append([cellm, celln, outbud[ipos]]) + + if sim.idxsim == 0: + no_vsc_bud_last = np.array(vals_to_store) + + # Ensure with and without VSC simulations give nearly identical flow results + # for each cell-to-cell exchange between columns 5 and 6 + assert np.allclose( + no_vsc_bud_last[:, 2], stored_ans[:, 2], atol=1e-3 + ), ( + "Flow in models " + + exdirs[0] + + " and the established answer should be approximately " + "equal, but are not." + ) + + elif sim.idxsim == 1: + with_vsc_bud_last = np.array(vals_to_store) + + assert np.allclose( + with_vsc_bud_last[:, 2], stored_ans[:, 2], atol=1e-3 + ), ( + "Flow in models " + + exdirs[1] + + " and the established answer should be approximately " + "equal, but are not." + ) + + elif sim.idxsim == 2: + no_vsc_low_k_bud_last = np.array(vals_to_store) + + # Ensure the cell-to-cell flow between columns 5 and 6 in model + # 3 is less than what's in the "with viscosity" model + assert np.less(no_vsc_low_k_bud_last[:, 2], stored_ans[:, 2]).all(), ( + "Exit flow from model the established answer " + "should be greater than flow existing " + + exdirs[2] + + ", but it is not." + ) + + +# - No need to change any code below +@pytest.mark.parametrize( + "idx, dir", + list(enumerate(exdirs)), +) +def test_mf6model(idx, dir): + # initialize testing framework + test = testing_framework() + + # build the model + test.build_mf6_models(build_model, idx, dir) + + # run the test model + test.run_mf6(Simulation(dir, exfunc=eval_results, idxsim=idx)) + + +def main(): + # initialize testing framework + test = testing_framework() + + # run the test model + for idx, dir in enumerate(exdirs): + test.build_mf6_models(build_model, idx, dir) + sim = Simulation(dir, exfunc=eval_results, idxsim=idx) + test.run_mf6(sim) + + +if __name__ == "__main__": + # print message + print(f"standalone run of {os.path.basename(__file__)}") + + # run main routine + main() diff --git a/doc/MODFLOW6References.bib b/doc/MODFLOW6References.bib index db6303c1e6e..ebd845d7ab0 100644 --- a/doc/MODFLOW6References.bib +++ b/doc/MODFLOW6References.bib @@ -494,7 +494,7 @@ @book{Voss1984sutra Date-Added = {2019-08-29 15:00:00 -0600}, Date-Modified = {2019-08-29 15:00:00 -0600}, Series = {{U.S. Geological Survey Water-Resources Investigations Report 84--4369, 409 p.}}, - Title = {SUTRA---a finite-element simulation model for saturated-unsaturated fluid-density-dependent ground-water flow with energy transport or chemically-reactive single-species solute transport}, + Title = {SUTRA---A finite-element simulation model for saturated-unsaturated fluid-density-dependent ground-water flow with energy transport or chemically-reactive single-species solute transport}, Year = {1984}} @article{VossSouza1987, @@ -597,6 +597,34 @@ @book{langevin2002seawat Year = {2002}, Bdsk-Url-1 = {https://pubs.er.usgs.gov/publication/tm6A22}} +@article{diersch2002viscosity, + Author = {Diersch, Hans-Jorg G and Kolditz, Olaf}, + Date-Added = {2022-09-27 11:58:19 -0500}, + Date-Modified = {2022-09-27 12:00:01 -0500}, + Journal = {Advances in Water Resources}, + Pages = {899--944}, + Title = {Variable-density flow and transport in porous media: Approaches and challenges}, + Url = {https://www.sciencedirect.com/science/article/pii/S0309170802000635}, + doi = {https://doi.org/10.1016/S0309-1708(02)00063-5}, + Urldate = {September 27, 2022}, + Volume = {25}, + number = {8}, + Year = {2002}, + issn = {0309-1708}, + Bdsk-Url-1 = {https://pubs.er.usgs.gov/publication/tm6A22}} + +@book{kipp1987, + Author = {Kipp, Kenneth L}, + Date-Added = {2022-09-27 11:58:19 -0500}, + Date-Modified = {2022-09-27 12:00:01 -0500}, + Institution = {U.S. Geological Survey}, + Series = {{U.S. Geological Survey Water-Resources Investigation Report 86-4095, 517 p.}}, + Title = {HST3D: A Computer Code for Simulation of Heat and Solute Transport in Three-Dimensional Ground-Water Flow Systems}, + Url = {https://pubs.usgs.gov/wri/1986/4095/report.pdf}, + Urldate = {September 27, 2022}, + Year = {1987}, + Bdsk-Url-1 = {https://pubs.er.usgs.gov/publication/wri864095}} + @book{langevin2008seawat, Author = {Langevin, Christian D and Thorne Jr, Daniel T and Dausman, Alyssa M and Sukop, Michael C and Guo, Weixing}, Date-Added = {2019-07-25 11:53:50 -0500}, @@ -2855,3 +2883,38 @@ @book{vs2d Urldate = {June 27, 2017}, Year = {1987}, Bdsk-Url-1 = {https://pubs.er.usgs.gov/publication/wri834099}} + +@book{healy1996, + Author = {Healy, Richard W and Ronan, Ann D}, + Date-Added = {2022-09-27 11:58:19 -0500}, + Date-Modified = {2022-09-27 12:00:01 -0500}, + Institution = {U.S. Geological Survey}, + Series = {{U.S. Geological Survey Water-Resources Investigation Report 96-4230, 36 p.}}, + Title = {Documentation of Computer Program VS2DH for Simulation of Energy Transport in Variably Saturated Porous Media: Modification of the U.S. Geological Survey's Computer Program VS2DT}, + Doi = {10.3133/wri964230}, + Url = {https://pubs.usgs.gov/wri/1996/4230/report.pdf}, + Urldate = {September 27, 2022}, + Year = {1996}, + Bdsk-Url-1 = {https://pubs.er.usgs.gov/publication/wri964230}} + +@book{hughes2004, + Author = {Hughes, Joseph D and Sanford, Ward E}, + Date-Added = {2019-10-11 15:42:15 -0400}, + Date-Modified = {2019-10-11 15:46:29 -0400}, + Doi = {10.3133/ofr20041207}, + Url = {https://water.usgs.gov/nrp/gwsoftware/SutraMS/OFR2004-1207.pdf}, + Series = {{U.S. Geological Survey Open File Report 2004--1207, 141 p.}}, + Title = {SUTRA-MS: A version of SUTRA modified to simulate heat and multiple-solute transport}, + Urldate = {September 27, 2022}, + Year = {2004}, + Bdsk-Url-1 = {https://pubs.er.usgs.gov/publication/ofr20041207}} + +@book{maidment1993, + Author = {Maidment, David R}, + Date-Added = {2022-10-11 15:42:15 -0400}, + Date-Modified = {2022-10-11 15:46:29 -0400}, + Isbn = {0-07-039732-5}, + Publisher = {McGraw-Hill}, + Address = {New York, USA}, + Title = {Handbook of Hydrology}, + Year = {1993}} diff --git a/doc/ReleaseNotes/ReleaseNotes.tex b/doc/ReleaseNotes/ReleaseNotes.tex index 7394311db4a..9d8aa6dafdb 100644 --- a/doc/ReleaseNotes/ReleaseNotes.tex +++ b/doc/ReleaseNotes/ReleaseNotes.tex @@ -192,17 +192,17 @@ \section{Changes Introduced in this Release} \underline{NEW FUNCTIONALITY} \begin{itemize} - \item xxx + \item A new Viscosity (VSC) package for the Groundwater Flow (GWF) Model is introduced in this release. The effects of viscosity are accounted for by updates to intercell conductance, as well as the conductance between the aquifer and head-dependent boundaries, based on simulated concentrations and\/or temperatures. The viscosity package is activated by specifying ``VSC6'' as the file type in a GWF name file. Changes to the code and input may be required in response to user needs and testing. % \item xxx % \item xxx \end{itemize} - %\underline{EXAMPLES} - %\begin{itemize} - % \item xxx + \underline{EXAMPLES} + \begin{itemize} + \item A new example called ex-gwt-stallman was added. This new problem uses the GWT Model as a surrogate for simulating heat flow. The example represents heat conduction in subsurface with a periodic temperature boundary condition at the surface. % \item xxx % \item xxx - %\end{itemize} + \end{itemize} \textbf{\underline{BUG FIXES AND OTHER CHANGES TO EXISTING FUNCTIONALITY}} \\ \underline{BASIC FUNCTIONALITY} @@ -213,6 +213,7 @@ \section{Changes Introduced in this Release} \item When a GWF Model and a corresponding GWT model are solved in the same simulation, the GWF Model must be solved before the corresponding GWT model. The GWF Model must also be solved by a different IMS than the GWT Model. There was not a check for this in previous versions and if these conditions were not met, the solution would often not converge or it would give erroneous results. \item The DISV Package would not raise an error if a model cell was defined as a line. The program was modified to check for the case where the calculated cell area is equal to zero. If the calculated cell area is equal to zero, the program terminates with an error. \item When searching for a required block in an input file, the program would not terminate with a sensible error message if the end of file was found instead of the required block. Program now indicates that the required block was not found. + \item This release contains a first step toward implementation of generic input routines to read input files. The new input routines were implemented for the DIS, DISV, and DISU Packages of the GWF and GWT Models, for the NPF Package of the GWF Model, and the DSP Package of the GWT Model. Output summaries written to the GWF and GWT Model listing files are different from summaries written using previous versions. For packages that use the new input data model, the IPRN capability of the READARRAY utility (described in mf6io.pdf) is no longer supported as a way to write input arrays to the model listing file. \end{itemize} \underline{INTERNAL FLOW PACKAGES} diff --git a/doc/SuppTechInfo/Figures/VSCnonlinear.pdf b/doc/SuppTechInfo/Figures/VSCnonlinear.pdf new file mode 100644 index 00000000000..96346ab714c Binary files /dev/null and b/doc/SuppTechInfo/Figures/VSCnonlinear.pdf differ diff --git a/doc/SuppTechInfo/body.tex b/doc/SuppTechInfo/body.tex index 35e96c48ebf..9b4c38dc6ba 100644 --- a/doc/SuppTechInfo/body.tex +++ b/doc/SuppTechInfo/body.tex @@ -41,6 +41,12 @@ \customlabel{ch:gencouple}{\thechapno} \input{gen-coupling.tex} +\newpage +\incchap +\SECTION{Chapter \thechapno. Accounting for Fluid Viscosity} +\customlabel{ch:vscpackage}{\thechapno} +\input{viscosity.tex} + \newpage \ifx\usgsdirector\undefined diff --git a/doc/SuppTechInfo/mf6suptechinfo.bbl b/doc/SuppTechInfo/mf6suptechinfo.bbl index e7bd7db541b..7a42f1fe96b 100644 --- a/doc/SuppTechInfo/mf6suptechinfo.bbl +++ b/doc/SuppTechInfo/mf6suptechinfo.bbl @@ -1,16 +1,38 @@ -\begin{thebibliography}{7} +\begin{thebibliography}{15} \providecommand{\natexlab}[1]{#1} \expandafter\ifx\csname urlstyle\endcsname\relax \providecommand{\doi}[1]{doi:\discretionary{}{}{}#1}\else \providecommand{\doi}{doi:\discretionary{}{}{}\begingroup \urlstyle{rm}\Url}\fi +\bibitem[{Guo and Langevin(2002)}]{langevin2002seawat} +Guo, Weixing, and Langevin, C.D., 2002, User's Guide to SEAWAT: A Computer + Program for Simulation of Three-Dimensional Variable-Density Ground-Water + Flow: {U.S. Geological Survey Techniques of Water-Resources Investigations + book 6, Chapter A7, 77 p.}, accessed July 25, 2019, at + \url{https://pubs.er.usgs.gov/publication/ofr03426}. + \bibitem[{Harbaugh(2005)}]{modflow2005} Harbaugh, A.W., 2005, MODFLOW-2005, the U.S. Geological Survey modular ground-water model---the Ground-Water Flow Process: {U.S. Geological Survey Techniques and Methods, book 6, chap. A16, variously paged}, accessed June 27, 2017, at \url{https://pubs.usgs.gov/tm/2005/tm6A16/}. +\bibitem[{Healy and Ronan(1996)}]{healy1996} +Healy, R.W., and Ronan, A.D., 1996, Documentation of Computer Program VS2DH for + Simulation of Energy Transport in Variably Saturated Porous Media: + Modification of the U.S. Geological Survey's Computer Program VS2DT: {U.S. + Geological Survey Water-Resources Investigation Report 96-4230, 36 p.}, + accessed September 27, 2022, at \url{https://doi.org/10.3133/wri964230}, at + \url{https://pubs.usgs.gov/wri/1996/4230/report.pdf}. + +\bibitem[{Hughes and Sanford(2004)}]{hughes2004} +Hughes, J.D., and Sanford, W.E., 2004, SUTRA-MS: A version of SUTRA modified to + simulate heat and multiple-solute transport: {U.S. Geological Survey Open + File Report 2004--1207, 141 p.}, accessed September 27, 2022, at + \url{https://doi.org/10.3133/ofr20041207}, at + \url{https://water.usgs.gov/nrp/gwsoftware/SutraMS/OFR2004-1207.pdf}. + \bibitem[{Hughes and others(2017)Hughes, Langevin, and Banta}]{modflow6framework} Hughes, J.D., Langevin, C.D., and Banta, E.R., 2017, Documentation for the @@ -24,6 +46,20 @@ Kavetski, Dmitri, and Kuczera, George, 2007, Model smoothing strategies to no.~3, \url{https://doi.org/10.1029/2006WR005195}, \url{https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2006WR005195}. +\bibitem[{Kipp(1987)}]{kipp1987} +Kipp, K.L., 1987, HST3D: A Computer Code for Simulation of Heat and Solute + Transport in Three-Dimensional Ground-Water Flow Systems: {U.S. Geological + Survey Water-Resources Investigation Report 86-4095, 517 p.}, accessed + September 27, 2022, at \url{https://pubs.usgs.gov/wri/1986/4095/report.pdf}. + +\bibitem[{Langevin and others(2008)Langevin, Thorne~Jr, Dausman, Sukop, and + Guo}]{langevin2008seawat} +Langevin, C.D., Thorne~Jr, D.T., Dausman, A.M., Sukop, M.C., and Guo, Weixing, + 2008, {SEAWAT} Version 4---A computer program for simulation of multi-species + solute and heat transport: {U.S. Geological Survey Techniques and Methods, + book 6, chap. A22, 39 p.}, accessed June 27, 2017, at + \url{https://pubs.er.usgs.gov/publication/tm6A22}. + \bibitem[{Langevin and others(2017)Langevin, Hughes, Provost, Banta, Niswonger, and Panday}]{modflow6gwf} Langevin, C.D., Hughes, J.D., Provost, A.M., Banta, E.R., Niswonger, R.G., and @@ -31,6 +67,9 @@ Langevin, C.D., Hughes, J.D., Provost, A.M., Banta, E.R., Niswonger, R.G., and Model: {U.S. Geological Survey Techniques and Methods, book 6, chap. A55, 197 p.}, \url{https://doi.org/10.3133/tm6A55}. +\bibitem[{Maidment(1993)}]{maidment1993} +Maidment, D.R., 1993, Handbook of Hydrology: New York, USA, McGraw-Hill. + \bibitem[{Panday and others(2013)Panday, Langevin, Niswonger, Ibaraki, and Hughes}]{modflowusg} Panday, Sorab, Langevin, C.D., Niswonger, R.G., Ibaraki, Motomu, and Hughes, @@ -46,6 +85,16 @@ Provost, A.M., Langevin, C.D., and Hughes, J.D., 2017, Documentation for the Geological Survey Techniques and Methods, book 6, chap. A56, 46 p.}, \url{https://doi.org/10.3133/tm6A56}. +\bibitem[{Voss(1984)}]{voss1984sutra} +Voss, C.I., 1984, SUTRA---A finite-element simulation model for + saturated-unsaturated fluid-density-dependent ground-water flow with energy + transport or chemically-reactive single-species solute transport: {U.S. + Geological Survey Water-Resources Investigations Report 84--4369, 409 p.} + +\bibitem[{Zheng(2010)}]{zheng2010supplemental} +Zheng, Chunmiao, 2010, MT3DMS v5.3, Supplemental User's Guide: {Technical + Report Prepared for the U.S. Army Corps of Engineers, 51 p.} + \bibitem[{Zheng and Wang(1999)}]{zheng1999mt3dms} Zheng, Chunmiao, 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/python/VSC-nonlinear.ipynb b/doc/SuppTechInfo/python/VSC-nonlinear.ipynb new file mode 100644 index 00000000000..c907b5687f0 --- /dev/null +++ b/doc/SuppTechInfo/python/VSC-nonlinear.ipynb @@ -0,0 +1,235 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.gridspec as gridspec\n", + "import math" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "try:\n", + " import spnspecs\n", + " spnspecs.set_graph_specifications()\n", + "except:\n", + " spnspecs = None" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "figpth = '../Figures'\n", + "width = 6.8\n", + "dpi = 300" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### A common nonlinear function relating viscosity and temperature\n", + "\n", + "The equation that follows is for temperature expressed in $^\\circ$C" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dtype = np.float\n", + "step = 0.01\n", + "T0, T1 = 0, 40\n", + "T = np.arange(T0, T1, step=step, dtype=dtype)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def nonlin_vsc(T):\n", + " a1 = 2.394e-5\n", + " a2 = 10\n", + " a3 = 248.37\n", + " a4 = 133.15\n", + " visc = np.zeros(T.shape, dtype=np.float)\n", + " for idx, val in enumerate(T):\n", + " mu_T = a1 * a2 ** (a3 / (T[idx] + a4))\n", + " visc[idx] = mu_T\n", + " \n", + " return visc" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### A linear approximation for the equation given above\n", + "\n", + "That is, contrive a linear approximation to the nonlinear case above" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "x_lin_pts = np.array([7.5, 30])\n", + "y_lin_pts = nonlin_vsc(x_lin_pts)\n", + "\n", + "# calc rise/run\n", + "dviscdT = (y_lin_pts[1] - y_lin_pts[0]) / (x_lin_pts[1] - x_lin_pts[0])\n", + "\n", + "# tempearture of reference viscosity\n", + "Tviscref = 20\n", + "\n", + "# reference viscosity\n", + "mu_0 = nonlin_vsc(np.array([Tviscref]))\n", + "\n", + "def lin_vsc(T, dviscdT, Tviscref):\n", + " # use the linear formula\n", + " visc_lin = np.zeros(T.shape, dtype=np.float)\n", + " for idx, val in enumerate(T):\n", + " mu_T = dviscdT * (val - Tviscref) + mu_0\n", + " visc_lin[idx] = mu_T\n", + " \n", + " return visc_lin" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "visc = nonlin_vsc(T)\n", + "visc_lin = lin_vsc(T, dviscdT, Tviscref)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(T,visc)\n", + "plt.plot(T,visc_lin)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plot data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "x_increment = 5\n", + "x_ticks = np.arange(T0, T1 + x_increment, x_increment)\n", + "y_increment = 2e-4\n", + "y_ticks = np.arange(math.floor(visc.min() * 10000) / 10000, math.ceil(visc.max() * 10000) / 10000 + y_increment, y_increment)\n", + "print(x_ticks)\n", + "print(y_ticks)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#fig, axes = plt.subplots(nrows=1, ncols=1, tight_layout=True, \n", + "# figsize=(width/1.5, width/2))\n", + "\n", + "#if not isinstance(axes, list):\n", + "# axes = [axes]\n", + "\n", + "#for idx, ax in enumerate(axes):\n", + "# ax.set_xlim(T0, T1)\n", + "# ax.set_ylim(y_ticks[0], y_ticks[-2])\n", + "# ax.set_xticks(x_ticks)\n", + "# ax.set_xticklabels(x_ticks)\n", + "# #ax.set_yticks(y_ticks)\n", + "# #ax.set_yticklabels(y_ticks, rotation=90, va='center')\n", + "# if spnspecs is not None:\n", + "# spnspecs.remove_edge_ticks(ax)\n", + "# spnspecs.heading(ax, letter=letters[idx])\n", + " \n", + "\n", + "#ax = axes[0]\n", + "##ax.set_aspect('square', adjustable='box')\n", + "##ax.axhline(0, lw=0.5, ls='-.', color='black')\n", + "##ax.axvline(0, lw=0.5, ls='-.', color='black')\n", + "\n", + "#ax.plot(T, visc, lw=1.25, color='#00BFFF', label=r'$\\mu')\n", + "#ax.plot(T, visc_lin, lw=1., color='black', ls='-', label=r'$\\frac {\\partial \\mu} {\\partial T}')\n", + "#ax.plot([Tviscref], [mu_0], marker=\"o\", markersize=4, markeredgecolor=\"red\", markerfacecolor=\"blue\")\n", + "#text = r'$\\left ( T_0, \\mu_0 \\right )$'\n", + "##ax.text(21, 0.00102, text, style='italic') #fontsize=6)\n", + "#spnspecs.add_text(ax, text=text, x=0.5, y=0.5, transform=False, bold=False, ha='left', va='top')\n", + "\n", + "#if spnspecs is not None:\n", + "# handles, labels = ax.get_legend_handles_labels()\n", + "# spnspecs.graph_legend(ax, handles=handles[::-1], labels=labels[::-1],\n", + "# loc='upper right', bbox_to_anchor=(1.00,0.95),\n", + "# handlelength=1.5)\n", + "#ax.set_xlabel(r'Temperature, in $^\\circ$C')\n", + "#ax.set_ylabel(r'$\\mu$, in $\\frac{kg}{m \\cdot s}$')\n", + "\n", + "#fpth = os.path.join(figpth, 'VSCnonlinear.pdf')\n", + "#fig.savefig(fpth, dpi=dpi)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.1" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/doc/SuppTechInfo/viscosity.tex b/doc/SuppTechInfo/viscosity.tex new file mode 100644 index 00000000000..49fdebe18c8 --- /dev/null +++ b/doc/SuppTechInfo/viscosity.tex @@ -0,0 +1,74 @@ + +The dynamic viscosity of a fluid is a function of fluid composition and temperature and typically has a weak dependence on pressure. In some applications, particularly those that involve large variations in temperature, viscosity can vary significantly in space and over time. Variations in viscosity, in turn, can affect groundwater flow by changing the hydraulic conductivity, which is inversely proportional to viscosity. In the original release of \mf, variations in viscosity and their effects on hydraulic conductivity were not taken into account. The \mf Viscosity (VSC) Package allows viscosity to vary as a function of solute concentrations and fluid temperature and accounts for the effect of variable viscosity on hydraulic conductivity in the Node-Property Flow (NPF) Package and in stress packages that involve head-dependent groundwater flow. Temperature-dependent viscosity has been implemented in other groundwater modeling codes, including HST3D \citep{kipp1987}, VS2DH \citep{healy1996}, SUTRA-MS \citep{hughes2004}, and SEAWAT Version 4 \citep{langevin2008seawat}. Like these other codes, the VSC Package does not account for the dependence of viscosity on pressure, which is negligible in most applications. + +For cases in which the rate of groundwater flow is proportional to the head gradient in the direction of flow, as in an isotropic porous medium, the flow between two adjacent locations can be approximated using a one-dimensional form of Darcy's Law \citep[equation 4-15]{modflow6gwf}: + +\begin{equation} +\label{eqn:modflow6gwf-4-15} +Q = C \left( h_{1} - h_{2} \right), +\end{equation} + +\noindent where $h_{1}$ and $h_{2}$ are the hydraulic heads at locations 1 and 2, respectively, $Q$ is the flow from location 1 to location 2, and $C$ is a hydraulic conductance defined by~\citep[equation 4-14]{modflow6gwf}: + +\begin{equation} +\label{eqn:modflowgwf-4-14} +C = \frac{K A} {L_{1,2}}, +\end{equation} + +\noindent where $K$ is the hydraulic conductivity of the porous medium that connects the two locations, $A$ is the cross-sectional area perpendicular to flow, and $L_{1,2}$ is the distance between two locations. A form of equation~\ref{eqn:modflow6gwf-4-15} is used in the \mf ``conductance-based'' formulation for flow between two adjacent model cells \citep{modflow6gwf}. In that case, the conductance is based on an effective hydraulic conductivity between the two cells, the area of the cell-cell interface, and the distance between the cell centers. A form of equation~\ref{eqn:modflow6gwf-4-15} is also used in \mf stress packages in which the flow into or out of the model at a model boundary is head-dependent. For example, in the General-Head Boundary (GHB) Package, flow between an external source and a groundwater model cell is proportional to the difference between the head assigned to the external source and the simulated head in the cell, and the conductance is representative of the material between the external source and the cell center. + +Because hydraulic conductivity is inversely proportional to viscosity, the conductivity $K$ for a simulated fluid composition and temperature is related to the conductivity $K_{0}$ for some reference composition and temperature by + +\begin{equation} +\label{eqn:conductivity-ratio} +K = \frac{\mu_{0}}{\mu} K_{0}. +\end{equation} + +\noindent It follows from equation~\ref{eqn:modflowgwf-4-14} that the conductance $C$ for a simulated fluid composition and temperature is related to the conductance $C_{0}$ for the reference composition and temperature by + +\begin{equation} +\label{eqn:conductance-ratio} +C = \frac{\mu_{0}}{\mu} C_{0}, +\end{equation} + +\noindent where $\mu$ and $\mu_{0}$ are the viscosities at the simulated and reference conditions, respectively. The VSC Package uses equations~\ref{eqn:conductivity-ratio} and~\ref{eqn:conductance-ratio} to update hydraulic conductivities and conductances to account for changes in solute concentrations and temperature during a flow and transport simulation based on the Groundwater Flow (GWF) and Groundwater Transport (GWT) Models. The GWT Model is designed to simulate solute transport but can be used to simulate heat transport in some applications by appropriately scaling the input parameters to render the solute-transport equation mathematically equivalent to the heat-transport equation \citep{zheng2010supplemental}. In such cases, the ``solute concentration'' calculated by a GWT-based heat-transport model can by identified in the VSC Package input file as representing temperature. Although the VSC Package allows viscosity to vary with solute concentration, variations in viscosity are typically most significant in heat-transport applications. + +\subsection{Dependence of Viscosity on Concentration and Temperature} \label{sec:fluidvsc} + +The VSC Package calculates viscosity as a function of concentration and temperature using the following equation \citep[equation 17]{langevin2008seawat}: + +\begin{equation} +\label{eqn:viscfunc} +\mu = \mu_T \left ( T \right ) + \sum_{k=1}^{NS} \frac{\partial \mu} {\partial C^k} \left ( C^k - C_{0}^k \right ) +\end{equation} + +\noindent where $NS$ is the number of chemical species (solutes) whose concentrations affect viscosity, $C^k$ and $C_{0}^k$ are the simulated and user-specified reference concentrations of species $k$, respectively, and $\partial \mu / \partial C^k$ is the user-specified rate at which viscosity changes linearly with the concentration of species $k$. (Symbols $C^k$ and $C_{0}^k$ are not to be confused with the symbols for conducance, $C$ and $C_0$, introduced earlier.) When all concentrations are equal to their reference values, the viscosity is equal to $\mu_T \left ( T \right )$, which embodies the dependence of viscosity on temperature according to + +\begin{align} + \label{eqn:viscT} + \mu_T \left ( T \right ) = \begin{dcases} + \mu_{0} + \frac{\partial \mu} {\partial T} \left ( T - T_{0}\right ) & linear \\ + \mu_{0} A_2^{\frac{A_3 \left ( T_0 - T \right )} {\left ( T + A_4 \right ) \left ( T_0 + A_4 \right ) }} & nonlinear + \end{dcases} . +\end{align} + +\noindent where $T$ and $T_{0}$ are the simulated and user-specified reference temperature, respectively, and $\mu_{0}$ is the user-specified reference viscosity. The reference viscosity is commonly set to 0.001 kg~m$^{-1}$~s$^{-1}$, which is representative of freshwater at a reference temperature of 20$^{\circ}$C \citep[Table 11.1]{maidment1993}. When $T$ equals $T_{0}$, or if temperature is not designated by the user as affecting viscosity, $\mu_T \left ( T \right )$ equals $\mu_{0}$. Equation~\ref{eqn:viscT} includes both linear and nonlinear variation of viscosity with temperature. Linear variation is the default, in which case $\partial \mu / \partial T$ is the user-specified rate at which viscosity changes linearly with temperature. Nonlinear variation is a user-selectable option, in which case the variation of viscosity from the reference value is determined by user-specified parameters $A_2$, $A_3$, and $A_4$. The nonlinear option in equation~\ref{eqn:viscT} is mathematically equivalent to one of the options available in SEAWAT Version 4 \citep[equation 18]{langevin2002seawat} but is written in a form that explicitly sets $\mu_T \left ( T \right )$ equal to $\mu_{0}$ when $T$ equals $T_{0}$. Setting $\mu_0$, $A_2$, $A_3$, and $A_4$ to 0.001002 kg~m$^{-1}$~s$^{-1}$, 10.0, 248.37, and 133.15, respectively, and rounding the lead coefficient to one decimal place gives the temperature dependence of viscosity used in SUTRA \citep{voss1984sutra} and SUTRA-MS \citep{hughes2004}. VS2DH \citep{healy1996} also calculate the temperature dependence of viscosity using functions that are, respectively, mathematically equivalent to and a special case of equation~\ref{eqn:viscT}. + +Over the temperature range $0-40^{\circ}$C, the viscosity of freshwater varies between approximately 0.0007 and 0.00175 kg~m$^{-1}$~s$^{-1}$, as indicated by the nonlinear (blue) curve in figure~\ref{fig:viscosityrelation}. The black line in figure~\ref{fig:viscosityrelation} depicts a linear approximation that coincides with the nonlinear curve at the reference temperature and viscosity, $\left ( T_0, \mu_0 \right ) = \left(\right.$20$^{\circ}$C, 0.001 kg~m$^{-1}$~s$^{-1}\left.\right)$. + +\begin{figure} + \begin{center} + \includegraphics{./Figures/VSCnonlinear.pdf} + \caption[Graph showing the nonlinear response in the viscosity as temperature changes]{Nonlinear dependence of viscosity on temperature (blue curve) that is representative of freshwater over the temperature range $0-40^{\circ}$C, and a linear approximation (black curve) that has the same value $\left ( \mu_0 \right )$ and slope at the reference temperature $\left ( T_0 \right )$. The blue curve is generated using 0.001002 kg~m$^{-1}$~s$^{-1}$,10.0, 248.37, and 133.15 for $\mu_0$, $A_2$, $A_3$, and $A_4$, respectively} + \label{fig:viscosityrelation} + \end{center} +\end{figure} + +\subsection{Accounting for Variable Viscosity in Flows Between Cells} \label{sec:gwfvsc} + +The VSC Package uses concentrations and temperatures calculated for each cell by a GWT model on the previous time step or outer solution iteration to calculate viscosities for the current time step. These viscosities are used to adjust cell hydraulic conductivity values in the NPF Package using equation~\ref{eqn:conductivity-ratio}. The reference values of conductivity are the values specified by the user in the NPF Package and, optionally, the Time-Varying Conductivity (TVK) Package (Chapter 6 of this document). The conductivity adjustment is performed after the user-specified conductivities are read in but before conductivity values are passed to program units that use cell conductivities to formulate expressions for flow between cells, which include the ``conductance-based'' formulation for flow and the XT3D capability \citep{modflow6xt3d}. If the VSC Package is not active, cell conductivities are not adjusted for variable viscosity, and the user-specified values are used. + +\subsection{Accounting for Variable Viscosity in Boundary Flows} \label{sec:gwfvsc} + +The VSC Package uses concentrations and temperatures calculated for each cell by a GWT model on the previous time step or outer solution iteration in calculating viscosities for the current time step. These viscosities are used to adjust hydraulic conductances, using equation~\ref{eqn:conductance-ratio}, in stress packages that involve head-dependent boundary flows, which include the River (RIV), General-Head Boundary (GHB), Drain (DRN), Streamflow Routing (SFR), Lake (LAK), and Multi-Aquifer Well (MAW) Packages \citep{modflow6gwf}. For a boundary flow out of the model, the viscosity is based on the simulated concentration or temperature in the cell from which the boundary flow originates. For a boundary flow into the model, the viscosity is based on the concentration or temperature of the water entering the cell from the boundary. The reference values of conductance are the values normally set or calculated by the stress package based on user input. The conductance adjustment is performed after the stress package completes its normal conductance calculation but before the conductance value is used to formulate the expression for the boundary flow. If the VSC Package is not active, stress-package conductances are not adjusted for variable viscosity. + diff --git a/doc/mf6io/gwf/gwf.tex b/doc/mf6io/gwf/gwf.tex index 98cd9b37554..452122446ef 100644 --- a/doc/mf6io/gwf/gwf.tex +++ b/doc/mf6io/gwf/gwf.tex @@ -109,6 +109,10 @@ \subsection{Skeletal Storage, Compaction, and Subsidence (CSUB) Package} \subsection{Buoyancy (BUY) Package} \input{gwf/buy} +\newpage +\subsection{Viscosity (VSC) Package} +\input{gwf/vsc} + \newpage \subsection{Constant-Head (CHD) Package} \input{gwf/chd} diff --git a/doc/mf6io/gwf/namefile.tex b/doc/mf6io/gwf/namefile.tex index 4d8724801ca..f54f1ee1e5f 100644 --- a/doc/mf6io/gwf/namefile.tex +++ b/doc/mf6io/gwf/namefile.tex @@ -31,6 +31,7 @@ \subsubsection{Explanation of Variables} STO6 & Storage Package \\ CSUB6 & Compaction and Subsidence Package \\ BUY6 & Buoyancy Package \\ +VSC6 & Viscosity Package \\ HFB6 & Horizontal Flow Barrier Package\\ CHD6 & Time-Variant Specified Head Option & * \\ WEL6 & Well Package & * \\ diff --git a/doc/mf6io/gwf/vsc.tex b/doc/mf6io/gwf/vsc.tex new file mode 100644 index 00000000000..9a1f3857481 --- /dev/null +++ b/doc/mf6io/gwf/vsc.tex @@ -0,0 +1,73 @@ +Input to the Viscosity (VSC) Package is read from the file that has type ``VSC6'' in the Name File. If the VSC Package is active within a groundwater flow model, then the model will account for the dependence of fluid viscosity on solute concentration and the resulting changes in hydraulic conductivity and stress-package conductances, which vary inversely with viscosity. Viscosity can be calculated as a function of one or more groundwater solute transport (GWT) species using an approach described in the Supplemental Technical Information document distributed with MODFLOW 6 (Chapter 8). Only one VSC Package can be specified for a GWF model. The VSC Package can be coupled with one or more GWT Models so that the fluid viscosity is updated dynamically with one or more simulated concentration fields. + +The VSC Package calculates fluid viscosity using the following equation from \cite{langevin2008seawat}: + +\begin{equation} +\label{eqn:visclinear} +\mu = VISCREF + \sum_{i=1}^{NVISCSPECIES} DVISCDC_i \left ( CONCENTRATION_i - CVISCREF_i \right ) +\end{equation} + +\noindent where $\mu$ is the calculated viscosity, $VISCREF$ is the viscosity of a reference fluid, typically taken to be freshwater at a temperature of 20 degrees Celsius, $NVISCSPECIES$ is the number of chemical species that contribute to the viscosity calculation, $DVISCDC_i$ is the parameter that describes how viscosity changes linearly as a function of concentration for chemical species $i$ (i.e. the slope of a line that relates viscosity to concentration), $CONCENTRATION_i$ is the concentration of species $i$, and $CVISCREF_i$ is the reference concentration for species $i$ corresponding to when the viscosity of the reference fluid is equal to $VISCREF$, which is normally set to a concentration of zero. + +In many applications, variations in temperature have a greater effect on fluid viscosity than variations in solute concentration. When a GWT model is formulated such that one of the transported ``species'' is heat (thermal energy), with ``concentration'' used to represent temperature \citep{zheng2010supplemental}, the viscosity can vary linearly with temperature, as it can with any other ``concentration.'' In that case, $CONCENTRATION_i$ and $CVISCREF_i$ represent the simulated and reference temperatures, respectively, and $DVISCDC_i$ represents the rate at which viscosity changes with temperature. In addition, the viscosity formula can optionally include a nonlinear dependence on temperature. In that case, equation 3 becomes + +\begin{equation} +\label{eqn:viscnonlinear} +\mu = \mu_T(T) + \sum_{i=1}^{NVISCSPECIES} DVISCDC_i \left ( CONCENTRATION_i - CVISCREF_i \right ) +\end{equation} + +\noindent where the first term on the right-hand side, $\mu_T(T)$, is a nonlinear function of temperature, and the summation corresponds to the summation in equation \ref{eqn:visclinear}, in which one of the ``species'' is heat. The nonlinear term in equation \ref{eqn:viscnonlinear} is of the form + +\begin{equation} +\label{eqn:munonlinear} +\mu_T(T) = CVISCREF_i \cdot A_2^{\left [ \frac {-A_3 \left ( CONCENTRATION_i - CVISCREF_i \right ) } {\left ( CONCENTRATION_i + A_4 \right ) \left ( CVISCREF_i + A_4 \right )} \right ]} +\end{equation} + +\noindent where the coefficients $A_2$, $A_3$, and $A_4$ are specified by the user. Values for $A_2$, $A_3$, and $A_4$ are commonly 10, 248.7, and 133.15, respectively \citep{langevin2008seawat, Voss1984sutra}. + +\subsubsection{Stress Packages} + +For head-dependent stress packages, the VSC Package can adjust the conductance used to calculate flow between the boundary and the aquifer to account for variations in viscosity. Conductance is assumed to vary inversely with viscosity. + +By default, the boundary viscosity is set equal to VISCREF, which, for freshwater, is typically set equal to 1.0. However, there are two additional options for setting the viscosity of a boundary package. The first is to assign an auxiliary variable with the name ``VISCOSITY''. If an auxiliary variable named ``VISCOSITY'' is detected, then it will be assigned as the viscosity of the fluid entering from the boundary. Alternatively, a viscosity value can be calculated for each boundary using the viscosity equation described above and one or more concentrations provided as auxiliary variables. In this case, the user must assign one auxiliary variable for each AUXSPECIESNAME listed in the PACKAGEDATA block below. Thus, there must be NVISCSPECIES auxiliary variables, each with the identical name as those specified in PACKAGEDATA. The VSC Package will calculate the viscosity for each boundary using these concentrations and the values specified for VISCREF, DVISCDC, and CVISCREF. If the boundary package contains an auxiliary variable named VISCOSITY and also contains AUXSPECIESNAME auxiliary variables, then the boundary viscosity value will be assigned to the one in the VISCOSITY auxiliary variable. + +A GWT Model can be used to calculate concentrations for the advanced stress packages (LAK, SFR, MAW, and UZF) if corresponding advanced transport packages are specified (LKT, SFT, MWT, and UZT). The advanced stress packages have an input option called FLOW\_PACKAGE\_AUXILIARY\_NAME. When activated, this option will result in the simulated concentration for a lake or other feature being copied from the advanced transport package into the auxiliary variable for the corresponding GWF stress package. This means that the viscosity for a lake or stream, for example, can be dynamically updated during the simulation using concentrations from advanced transport packages that are fed into auxiliary variables in the advanced stress packages, and ultimately used by the VSC Package to calculate a fluid viscosity. This concept also applies when multiple GWT Models are used simultaneously to simulate multiple species. In this case, multiple auxiliary variables are required for an advanced stress package, with each one representing a concentration from a different GWT Model. + + +\begin{longtable}{p{3cm} p{12cm}} +\caption{Description of viscosity terms for stress packages} +\tabularnewline +\hline +\hline +\textbf{Stress Package} & \textbf{Note} \\ +\hline +\endhead +\hline +\endfoot +GHB & A VISCOSITY auxiliary variable or one or more auxiliary variables for calculating viscosity in the equation of state can be specified \\ +RIV & A VISCOSITY auxiliary variable or one or more auxiliary variables for calculating viscosity in the equation of state can be specified \\ +DRN & The drain formulation assumes that the drain boundary contains water of the same viscosity as the discharging water; auxiliary variables have no effect on the drain calculation \\ +LAK & A VISCOSITY auxiliary variable or one or more auxiliary variables for calculating viscosity in the equation of state can be specified \\ +SFR & A VISCOSITY auxiliary variable or one or more auxiliary variables for calculating viscosity in the equation of state can be specified \\ +MAW & A VISCOSITY auxiliary variable or one or more auxiliary variables for calculating viscosity in the equation of state can be specified \\ +UZF & Viscosity variations not implemented \\ +\end{longtable} + +\vspace{5mm} +\subsubsection{Structure of Blocks} + +\vspace{5mm} +\noindent \textit{FOR EACH SIMULATION} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/gwf-vsc-options.dat} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/gwf-vsc-dimensions.dat} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/gwf-vsc-packagedata.dat} + +\vspace{5mm} +\subsubsection{Explanation of Variables} +\begin{description} +\input{./mf6ivar/tex/gwf-vsc-desc.tex} +\end{description} + +\vspace{5mm} +\subsubsection{Example Input File} +\lstinputlisting[style=inputfile]{./mf6ivar/examples/gwf-vsc-example.dat} diff --git a/doc/mf6io/mf6io.bbl b/doc/mf6io/mf6io.bbl index 05fffcda052..c4d7f4a09b6 100644 --- a/doc/mf6io/mf6io.bbl +++ b/doc/mf6io/mf6io.bbl @@ -1,4 +1,4 @@ -\begin{thebibliography}{34} +\begin{thebibliography}{35} \providecommand{\natexlab}[1]{#1} \expandafter\ifx\csname urlstyle\endcsname\relax \providecommand{\doi}[1]{doi:\discretionary{}{}{}#1}\else @@ -228,6 +228,12 @@ Prudic, D.E., Konikow, L.F., and Banta, E.R., 2004, A New Streamflow-Routing {U.S. Geological Survey Open File Report 2004--1042, 104 p.}, accessed June 27, 2017, at \url{https://pubs.er.usgs.gov/publication/ofr20041042}. +\bibitem[{Voss(1984)}]{Voss1984sutra} +Voss, C.I., 1984, SUTRA---A finite-element simulation model for + saturated-unsaturated fluid-density-dependent ground-water flow with energy + transport or chemically-reactive single-species solute transport: {U.S. + Geological Survey Water-Resources Investigations Report 84--4369, 409 p.} + \bibitem[{Zheng(2010)}]{zheng2010supplemental} Zheng, Chunmiao, 2010, MT3DMS v5.3, Supplemental User's Guide: {Technical Report Prepared for the U.S. Army Corps of Engineers, 51 p.} diff --git a/doc/mf6io/mf6ivar/dfn/gwf-vsc.dfn b/doc/mf6io/mf6ivar/dfn/gwf-vsc.dfn new file mode 100644 index 00000000000..9e3ab6d0083 --- /dev/null +++ b/doc/mf6io/mf6ivar/dfn/gwf-vsc.dfn @@ -0,0 +1,174 @@ +# --------------------- gwf vsc options --------------------- + +block options +name viscref +type double +reader urword +optional true +longname reference viscosity +description fluid reference viscosity used in the equation of state. This value is set to 1.0 if not specified as an option. +default_value 1.0 + +block options +name temperature_species_name +type string +shape +reader urword +optional true +longname auxspeciesname that corresponds to temperature +description string used to identify the auxspeciesname in PACKAGEDATA that corresponds to the temperature species. There can be only one occurrence of this temperature species name in the PACKAGEDATA block or the program will terminate with an error. This value has no effect if viscosity does not depend on temperature. + +block options +name thermal_formulation +type string +shape +reader urword +optional true +valid linear nonlinear +longname keyword to specify viscosity formulation for the temperature species +description may be used for specifying which viscosity formulation to use for the temperature species. Can be either LINEAR or NONLINEAR. The LINEAR viscosity formulation is the default. + +block options +name thermal_a2 +type double +reader urword +optional true +longname coefficient used in nonlinear viscosity function +description is an empirical parameter specified by the user for calculating viscosity using a nonlinear formulation. If A2 is not specified, a default value of 10.0 is assigned (Voss, 1984). +default_value 10. + +block options +name thermal_a3 +type double +reader urword +optional true +longname coefficient used in nonlinear viscosity function +description is an empirical parameter specified by the user for calculating viscosity using a nonlinear formulation. If A3 is not specified, a default value of 248.37 is assigned (Voss, 1984). +default_value 248.37 + +block options +name thermal_a4 +type double precision +reader urword +optional true +longname coefficient used in nonlinear viscosity function +description is an empirical parameter specified by the user for calculating viscosity using a nonlinear formulation. If A4 is not specified, a default value of 133.15 is assigned (Voss, 1984). +default_value 133.15 + +block options +name viscosity_filerecord +type record viscosity fileout viscosityfile +shape +reader urword +tagged true +optional true +longname +description + +block options +name viscosity +type keyword +shape +in_record true +reader urword +tagged true +optional false +longname viscosity keyword +description keyword to specify that record corresponds to viscosity. + +block options +name fileout +type keyword +shape +in_record true +reader urword +tagged true +optional false +longname file keyword +description keyword to specify that an output filename is expected next. + +block options +name viscosityfile +type string +preserve_case true +shape +in_record true +reader urword +tagged false +optional false +longname file keyword +description name of the binary output file to write viscosity information. The viscosity file has the same format as the head file. Viscosity values will be written to the viscosity file whenever heads are written to the binary head file. The settings for controlling head output are contained in the Output Control option. + + +# --------------------- gwf vsc dimensions --------------------- + +block dimensions +name nviscspecies +type integer +reader urword +optional false +longname number of species used in viscosity equation of state +description number of species used in the viscosity equation of state. If either concentrations or temperature (or both) are used to update viscosity then then nrhospecies needs to be at least one. + + +# --------------------- gwf vsc packagedata --------------------- + +block packagedata +name packagedata +type recarray iviscspec dviscdc cviscref modelname auxspeciesname +shape (nrhospecies) +reader urword +longname +description + +block packagedata +name iviscspec +type integer +shape +tagged false +in_record true +reader urword +longname species number for this entry +description integer value that defines the species number associated with the specified PACKAGEDATA data entered on each line. IVISCSPECIES must be greater than zero and less than or equal to NVISCSPECIES. Information must be specified for each of the NVISCSPECIES species or the program will terminate with an error. The program will also terminate with an error if information for a species is specified more than once. +numeric_index true + +block packagedata +name dviscdc +type double precision +shape +tagged false +in_record true +reader urword +longname slope of the line that defines the linear relationship between viscosity and temperature or between viscosity and concentration, depending on the type of species entered on each line. +description real value that defines the slope of the line defining the linear relationship between viscosity and temperature or between viscosity and concentration, depending on the type of species entered on each line. If the value of AUXSPECIESNAME entered on a line corresponds to TEMPERATURE\_SPECIES\_NAME (in the OPTIONS block), this value will be used when VISCOSITY\_FUNC is equal to LINEAR (the default) in the OPTIONS block. When VISCOSITY\_FUNC is set to NONLINEAR, a value for DVISCDC must be specified though it is not used. + +block packagedata +name cviscref +type double precision +shape +tagged false +in_record true +reader urword +longname reference temperature value or reference concentration value +description real value that defines the reference temperature or reference concentration value used for this species in the viscosity equation of state. If AUXSPECIESNAME entered on a line corresponds to TEMPERATURE\_SPECIES\_NAME (in the OPTIONS block), then CVISCREF refers to a reference temperature, otherwise it refers to a reference concentration. + +block packagedata +name modelname +type string +in_record true +tagged false +shape +reader urword +longname modelname +description name of a GWT model used to simulate a species that will be used in the viscosity equation of state. This name will have no effect if the simulation does not include a GWT model that corresponds to this GWF model. + +block packagedata +name auxspeciesname +type string +in_record true +tagged false +shape +reader urword +longname auxspeciesname +description name of an auxiliary variable in a GWF stress package that will be used for this species to calculate the viscosity values. If a viscosity value is needed by the Viscosity Package then it will use the temperature or concentration values associated with this AUXSPECIESNAME in the viscosity equation of state. For advanced stress packages (LAK, SFR, MAW, and UZF) that have an associated advanced transport package (LKT, SFT, MWT, and UZT), the FLOW\_PACKAGE\_AUXILIARY\_NAME option in the advanced transport package can be used to transfer simulated temperature or concentration(s) into the flow package auxiliary variable. In this manner, the Viscosity Package can calculate viscosity values for lakes, streams, multi-aquifer wells, and unsaturated zone flow cells using simulated concentrations. + diff --git a/doc/mf6io/mf6ivar/examples/gwf-vsc-example.dat b/doc/mf6io/mf6ivar/examples/gwf-vsc-example.dat new file mode 100644 index 00000000000..48d21bd801b --- /dev/null +++ b/doc/mf6io/mf6ivar/examples/gwf-vsc-example.dat @@ -0,0 +1,18 @@ +BEGIN OPTIONS + VISCREF 8.904E-04 + THERMAL_FORMULATION NONLINEAR + THERMAL_A2 10.0 + THERMAL_A3 248.37 + THERMAL_A4 133.15 + VISCOSITY FILEOUT GWF-VSC.vsc.bin +END OPTIONS + +BEGIN DIMENSIONS + NVISCSPECIES 2 +END DIMENSIONS + +BEGIN PACKAGEDATA +# ISPEC DVISCDC CVISCREF MODELNAME AUXSPECIESNAME + 1 1.92e-6 0.0 GWT-SALT SALINITY + 2 0.00 25.0 GWT-TEMP TEMPERATURE +END PACKAGEDATA diff --git a/doc/mf6io/mf6ivar/md/mf6ivar.md b/doc/mf6io/mf6ivar/md/mf6ivar.md index 172ce2b91ed..4423337ddb1 100644 --- a/doc/mf6io/mf6ivar/md/mf6ivar.md +++ b/doc/mf6io/mf6ivar/md/mf6ivar.md @@ -148,8 +148,8 @@ | GWF | DISV | DIMENSIONS | NCPL | INTEGER | is the number of cells per layer. This is a constant value for the grid and it applies to all layers. | | GWF | DISV | DIMENSIONS | NVERT | INTEGER | is the total number of (x, y) vertex pairs used to characterize the horizontal configuration of the model grid. | | GWF | DISV | GRIDDATA | TOP | DOUBLE PRECISION (NCPL) | is the top elevation for each cell in the top model layer. | -| GWF | DISV | GRIDDATA | BOTM | DOUBLE PRECISION (NLAY, NCPL) | is the bottom elevation for each cell. | -| GWF | DISV | GRIDDATA | IDOMAIN | INTEGER (NLAY, NCPL) | is an optional array that characterizes the existence status of a cell. If the IDOMAIN array is not specified, then all model cells exist within the solution. If the IDOMAIN value for a cell is 0, the cell does not exist in the simulation. Input and output values will be read and written for the cell, but internal to the program, the cell is excluded from the solution. If the IDOMAIN value for a cell is 1 or greater, the cell exists in the simulation. If the IDOMAIN value for a cell is -1, the cell does not exist in the simulation. Furthermore, the first existing cell above will be connected to the first existing cell below. This type of cell is referred to as a ``vertical pass through'' cell. | +| GWF | DISV | GRIDDATA | BOTM | DOUBLE PRECISION (NCPL, NLAY) | is the bottom elevation for each cell. | +| GWF | DISV | GRIDDATA | IDOMAIN | INTEGER (NCPL, NLAY) | is an optional array that characterizes the existence status of a cell. If the IDOMAIN array is not specified, then all model cells exist within the solution. If the IDOMAIN value for a cell is 0, the cell does not exist in the simulation. Input and output values will be read and written for the cell, but internal to the program, the cell is excluded from the solution. If the IDOMAIN value for a cell is 1 or greater, the cell exists in the simulation. If the IDOMAIN value for a cell is -1, the cell does not exist in the simulation. Furthermore, the first existing cell above will be connected to the first existing cell below. This type of cell is referred to as a ``vertical pass through'' cell. | | GWF | DISV | VERTICES | IV | INTEGER | is the vertex number. Records in the VERTICES block must be listed in consecutive order from 1 to NVERT. | | GWF | DISV | VERTICES | XV | DOUBLE PRECISION | is the x-coordinate for the vertex. | | GWF | DISV | VERTICES | YV | DOUBLE PRECISION | is the y-coordinate for the vertex. | @@ -187,6 +187,7 @@ | GWF | DISU | CELL2D | ICVERT | INTEGER (NCVERT) | is an array of integer values containing vertex numbers (in the VERTICES block) used to define the cell. Vertices must be listed in clockwise order. | | GWF | IC | GRIDDATA | STRT | DOUBLE PRECISION (NODES) | is the initial (starting) head---that is, head at the beginning of the GWF Model simulation. STRT must be specified for all simulations, including steady-state simulations. One value is read for every model cell. For simulations in which the first stress period is steady state, the values used for STRT generally do not affect the simulation (exceptions may occur if cells go dry and (or) rewet). The execution time, however, will be less if STRT includes hydraulic heads that are close to the steady-state solution. A head value lower than the cell bottom can be provided if a cell should start as dry. | | GWF | NPF | OPTIONS | SAVE_FLOWS | KEYWORD | keyword to indicate that budget flow terms will be written to the file specified with ``BUDGET SAVE FILE'' in Output Control. | +| GWF | NPF | OPTIONS | PRINT_FLOWS | KEYWORD | keyword to indicate that calculated flows between cells will be printed to the listing file for every stress period time step in which ``BUDGET PRINT'' is specified in Output Control. If there is no Output Control option and ``PRINT\_FLOWS'' is specified, then flow rates are printed for the last time step of each stress period. This option can produce extremely large list files because all cell-by-cell flows are printed. It should only be used with the NPF Package for models that have a small number of cells. | | GWF | NPF | OPTIONS | ALTERNATIVE_CELL_AVERAGING | STRING | is a text keyword to indicate that an alternative method will be used for calculating the conductance for horizontal cell connections. The text value for ALTERNATIVE\_CELL\_AVERAGING can be ``LOGARITHMIC'', ``AMT-LMK'', or ``AMT-HMK''. ``AMT-LMK'' signifies that the conductance will be calculated using arithmetic-mean thickness and logarithmic-mean hydraulic conductivity. ``AMT-HMK'' signifies that the conductance will be calculated using arithmetic-mean thickness and harmonic-mean hydraulic conductivity. If the user does not specify a value for ALTERNATIVE\_CELL\_AVERAGING, then the harmonic-mean method will be used. This option cannot be used if the XT3D option is invoked. | | GWF | NPF | OPTIONS | THICKSTRT | KEYWORD | indicates that cells having a negative ICELLTYPE are confined, and their cell thickness for conductance calculations will be computed as STRT-BOT rather than TOP-BOT. | | GWF | NPF | OPTIONS | VARIABLECV | KEYWORD | keyword to indicate that the vertical conductance will be calculated using the saturated thickness and properties of the overlying cell and the thickness and properties of the underlying cell. If the DEWATERED keyword is also specified, then the vertical conductance is calculated using only the saturated thickness and properties of the overlying cell if the head in the underlying cell is below its top. If these keywords are not specified, then the default condition is to calculate the vertical conductance at the start of the simulation using the initial head and the cell properties. The vertical conductance remains constant for the entire simulation. | @@ -204,7 +205,12 @@ | GWF | NPF | OPTIONS | K33OVERK | KEYWORD | keyword to indicate that specified K33 is a ratio of K33 divided by K. If this option is specified, then the K33 array entered in the NPF Package will be multiplied by K after being read. | | GWF | NPF | OPTIONS | TVK6 | KEYWORD | keyword to specify that record corresponds to a time-varying hydraulic conductivity (TVK) file. The behavior of TVK and a description of the input file is provided separately. | | GWF | NPF | OPTIONS | FILEIN | KEYWORD | keyword to specify that an input filename is expected next. | -| GWF | NPF | OPTIONS | TVK_FILENAME | STRING | defines a time-varying hydraulic conductivity (TVK) input file. Records in the TVK file can be used to change hydraulic conductivity properties at specified times or stress periods. | +| GWF | NPF | OPTIONS | TVK6_FILENAME | STRING | defines a time-varying hydraulic conductivity (TVK) input file. Records in the TVK file can be used to change hydraulic conductivity properties at specified times or stress periods. | +| GWF | NPF | OPTIONS | DEV_NO_NEWTON | KEYWORD | turn off Newton for unconfined cells | +| GWF | NPF | OPTIONS | DEV_MODFLOWUSG_UPSTREAM_WEIGHTED_SATURATION | KEYWORD | use MODFLOW-USG upstream-weighted saturation approach | +| GWF | NPF | OPTIONS | DEV_MODFLOWNWT_UPSTREAM_WEIGHTING | KEYWORD | use MODFLOW-NWT approach for upstream weighting | +| GWF | NPF | OPTIONS | DEV_MINIMUM_SATURATED_THICKNESS | DOUBLE PRECISION | set minimum allowed saturated thickness | +| GWF | NPF | OPTIONS | DEV_OMEGA | DOUBLE PRECISION | set saturation omega value | | GWF | NPF | GRIDDATA | ICELLTYPE | INTEGER (NODES) | flag for each cell that specifies how saturated thickness is treated. 0 means saturated thickness is held constant; $>$0 means saturated thickness varies with computed head when head is below the cell top; $<$0 means saturated thickness varies with computed head unless the THICKSTRT option is in effect. When THICKSTRT is in effect, a negative value of icelltype indicates that saturated thickness will be computed as STRT-BOT and held constant. | | GWF | NPF | GRIDDATA | K | DOUBLE PRECISION (NODES) | is the hydraulic conductivity. For the common case in which the user would like to specify the horizontal hydraulic conductivity and the vertical hydraulic conductivity, then K should be assigned as the horizontal hydraulic conductivity, K33 should be assigned as the vertical hydraulic conductivity, and K22 and the three rotation angles should not be specified. When more sophisticated anisotropy is required, then K corresponds to the K11 hydraulic conductivity axis. All included cells (IDOMAIN $>$ 0) must have a K value greater than zero. | | GWF | NPF | GRIDDATA | K22 | DOUBLE PRECISION (NODES) | is the hydraulic conductivity of the second ellipsoid axis (or the ratio of K22/K if the K22OVERK option is specified); for an unrotated case this is the hydraulic conductivity in the y direction. If K22 is not included in the GRIDDATA block, then K22 is set equal to K. For a regular MODFLOW grid (DIS Package is used) in which no rotation angles are specified, K22 is the hydraulic conductivity along columns in the y direction. For an unstructured DISU grid, the user must assign principal x and y axes and provide the angle for each cell face relative to the assigned x direction. All included cells (IDOMAIN $>$ 0) must have a K22 value greater than zero. | @@ -785,6 +791,21 @@ | GWF | OC | PERIOD | LAST | KEYWORD | keyword to indicate save for last step in period. This keyword may be used in conjunction with other keywords to print or save results for multiple time steps. | | GWF | OC | PERIOD | FREQUENCY | INTEGER | save at the specified time step frequency. This keyword may be used in conjunction with other keywords to print or save results for multiple time steps. | | GWF | OC | PERIOD | STEPS | INTEGER ( +END DIMENSIONS diff --git a/doc/mf6io/mf6ivar/tex/gwf-vsc-options.dat b/doc/mf6io/mf6ivar/tex/gwf-vsc-options.dat new file mode 100644 index 00000000000..0c62b85da3a --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/gwf-vsc-options.dat @@ -0,0 +1,9 @@ +BEGIN OPTIONS + [VISCREF ] + [TEMPERATURE_SPECIES_NAME ] + [THERMAL_FORMULATION ] + [THERMAL_A2 ] + [THERMAL_A3 ] + [THERMAL_A4 ] + [VISCOSITY FILEOUT ] +END OPTIONS diff --git a/doc/mf6io/mf6ivar/tex/gwf-vsc-packagedata.dat b/doc/mf6io/mf6ivar/tex/gwf-vsc-packagedata.dat new file mode 100644 index 00000000000..ff3d65886f9 --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/gwf-vsc-packagedata.dat @@ -0,0 +1,5 @@ +BEGIN PACKAGEDATA + + + ... +END PACKAGEDATA diff --git a/make/makefile b/make/makefile index e17d8abd2bd..2cf7532c530 100644 --- a/make/makefile +++ b/make/makefile @@ -6,28 +6,28 @@ include ./makedefaults # Define the source file directories SOURCEDIR1=../src SOURCEDIR2=../src/Exchange -SOURCEDIR3=../src/Solution -SOURCEDIR4=../src/Solution/LinearMethods -SOURCEDIR5=../src/Timing -SOURCEDIR6=../src/Utilities -SOURCEDIR7=../src/Utilities/Idm -SOURCEDIR8=../src/Utilities/TimeSeries -SOURCEDIR9=../src/Utilities/Memory -SOURCEDIR10=../src/Utilities/OutputControl -SOURCEDIR11=../src/Utilities/ArrayRead -SOURCEDIR12=../src/Utilities/Libraries -SOURCEDIR13=../src/Utilities/Libraries/rcm -SOURCEDIR14=../src/Utilities/Libraries/blas -SOURCEDIR15=../src/Utilities/Libraries/sparskit2 -SOURCEDIR16=../src/Utilities/Libraries/daglib -SOURCEDIR17=../src/Utilities/Libraries/sparsekit -SOURCEDIR18=../src/Utilities/Observation -SOURCEDIR19=../src/Model -SOURCEDIR20=../src/Model/Connection -SOURCEDIR21=../src/Model/GroundWaterTransport -SOURCEDIR22=../src/Model/ModelUtilities -SOURCEDIR23=../src/Model/GroundWaterFlow -SOURCEDIR24=../src/Model/Geometry +SOURCEDIR3=../src/Model +SOURCEDIR4=../src/Model/Connection +SOURCEDIR5=../src/Model/Geometry +SOURCEDIR6=../src/Model/GroundWaterFlow +SOURCEDIR7=../src/Model/GroundWaterTransport +SOURCEDIR8=../src/Model/ModelUtilities +SOURCEDIR9=../src/Solution +SOURCEDIR10=../src/Solution/LinearMethods +SOURCEDIR11=../src/Timing +SOURCEDIR12=../src/Utilities +SOURCEDIR13=../src/Utilities/ArrayRead +SOURCEDIR14=../src/Utilities/Idm +SOURCEDIR15=../src/Utilities/Libraries +SOURCEDIR16=../src/Utilities/Libraries/blas +SOURCEDIR17=../src/Utilities/Libraries/daglib +SOURCEDIR18=../src/Utilities/Libraries/rcm +SOURCEDIR19=../src/Utilities/Libraries/sparsekit +SOURCEDIR20=../src/Utilities/Libraries/sparskit2 +SOURCEDIR21=../src/Utilities/Memory +SOURCEDIR22=../src/Utilities/Observation +SOURCEDIR23=../src/Utilities/OutputControl +SOURCEDIR24=../src/Utilities/TimeSeries VPATH = \ ${SOURCEDIR1} \ @@ -111,47 +111,55 @@ $(OBJDIR)/gwf3disv8idm.o \ $(OBJDIR)/gwf3disu8idm.o \ $(OBJDIR)/gwf3dis8idm.o \ $(OBJDIR)/Integer2dReader.o \ +$(OBJDIR)/BudgetFileReader.o \ $(OBJDIR)/TimeArraySeriesManager.o \ $(OBJDIR)/PackageMover.o \ $(OBJDIR)/Obs3.o \ $(OBJDIR)/NumericalPackage.o \ $(OBJDIR)/Budget.o \ -$(OBJDIR)/BudgetFileReader.o \ $(OBJDIR)/StructVector.o \ $(OBJDIR)/IdmLogger.o \ $(OBJDIR)/InputDefinitionSelector.o \ $(OBJDIR)/Integer1dReader.o \ $(OBJDIR)/Double2dReader.o \ $(OBJDIR)/Double1dReader.o \ +$(OBJDIR)/sort.o \ +$(OBJDIR)/SfrCrossSectionUtils.o \ +$(OBJDIR)/BudgetTerm.o \ $(OBJDIR)/BoundaryPackage.o \ $(OBJDIR)/BaseModel.o \ -$(OBJDIR)/BudgetTerm.o \ $(OBJDIR)/StructArray.o \ $(OBJDIR)/ModflowInput.o \ $(OBJDIR)/LayeredArrayReader.o \ +$(OBJDIR)/SfrCrossSectionManager.o \ +$(OBJDIR)/dag_module.o \ +$(OBJDIR)/BudgetObject.o \ $(OBJDIR)/NumericalModel.o \ $(OBJDIR)/mf6lists.o \ $(OBJDIR)/PackageBudget.o \ $(OBJDIR)/HeadFileReader.o \ -$(OBJDIR)/BudgetObject.o \ -$(OBJDIR)/sort.o \ -$(OBJDIR)/SfrCrossSectionUtils.o \ $(OBJDIR)/PrintSaveManager.o \ $(OBJDIR)/Xt3dAlgorithm.o \ $(OBJDIR)/gwf3tvbase8.o \ $(OBJDIR)/LoadMf6FileType.o \ +$(OBJDIR)/gwf3sfr8.o \ +$(OBJDIR)/gwf3riv8.o \ +$(OBJDIR)/gwf3maw8.o \ +$(OBJDIR)/gwf3lak8.o \ +$(OBJDIR)/GwfVscInputData.o \ +$(OBJDIR)/gwf3ghb8.o \ +$(OBJDIR)/gwf3drn8.o \ $(OBJDIR)/DistributedModel.o \ $(OBJDIR)/BaseExchange.o \ $(OBJDIR)/UzfCellGroup.o \ $(OBJDIR)/gwt1fmi1.o \ -$(OBJDIR)/SfrCrossSectionManager.o \ -$(OBJDIR)/dag_module.o \ $(OBJDIR)/OutputControlData.o \ $(OBJDIR)/gwf3ic8.o \ $(OBJDIR)/Xt3dInterface.o \ $(OBJDIR)/gwf3tvk8.o \ $(OBJDIR)/MemoryManagerExt.o \ $(OBJDIR)/IdmMf6FileLoader.o \ +$(OBJDIR)/gwf3vsc8.o \ $(OBJDIR)/GwfNpfOptions.o \ $(OBJDIR)/CellWithNbrs.o \ $(OBJDIR)/NumericalExchange.o \ @@ -159,11 +167,8 @@ $(OBJDIR)/Iunit.o \ $(OBJDIR)/gwf3uzf8.o \ $(OBJDIR)/gwt1apt1.o \ $(OBJDIR)/GwtSpc.o \ -$(OBJDIR)/gwf3sfr8.o \ $(OBJDIR)/OutputControl.o \ $(OBJDIR)/gwt1ic1.o \ -$(OBJDIR)/gwf3maw8.o \ -$(OBJDIR)/gwf3lak8.o \ $(OBJDIR)/gwt1mst1.o \ $(OBJDIR)/GwtDspOptions.o \ $(OBJDIR)/gwf3npf8.o \ @@ -199,7 +204,6 @@ $(OBJDIR)/gwf3disv8.o \ $(OBJDIR)/gwf3dis8.o \ $(OBJDIR)/gwf3api8.o \ $(OBJDIR)/gwf3wel8.o \ -$(OBJDIR)/gwf3riv8.o \ $(OBJDIR)/gwf3rch8.o \ $(OBJDIR)/gwf3sto8.o \ $(OBJDIR)/gwf3oc8.o \ @@ -209,9 +213,7 @@ $(OBJDIR)/gwf3hfb8.o \ $(OBJDIR)/gwf3csub8.o \ $(OBJDIR)/gwf3buy8.o \ $(OBJDIR)/GhostNode.o \ -$(OBJDIR)/gwf3ghb8.o \ $(OBJDIR)/gwf3evt8.o \ -$(OBJDIR)/gwf3drn8.o \ $(OBJDIR)/gwf3chd8.o \ $(OBJDIR)/ims8reordering.o \ $(OBJDIR)/GridConnection.o \ diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index 547d51d3ee3..15eb52dea77 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -105,6 +105,7 @@ + @@ -135,6 +136,7 @@ + diff --git a/src/Exchange/GwfGwtExchange.f90 b/src/Exchange/GwfGwtExchange.f90 index b8a121c474f..44e040b5986 100644 --- a/src/Exchange/GwfGwtExchange.f90 +++ b/src/Exchange/GwfGwtExchange.f90 @@ -275,12 +275,18 @@ subroutine exg_ar(this) end if end if ! - ! -- Set a pointer to conc + ! -- Set a pointer to conc in buy if (gwfmodel%inbuy > 0) then call gwfmodel%buy%set_concentration_pointer(gwtmodel%name, gwtmodel%x, & gwtmodel%ibound) end if ! + ! -- Set a pointer to conc (which could be a temperature) in vsc + if (gwfmodel%invsc > 0) then + call gwfmodel%vsc%set_concentration_pointer(gwtmodel%name, gwtmodel%x, & + gwtmodel%ibound) + end if + ! ! -- transfer the boundary package information from gwf to gwt call this%gwfbnd2gwtfmi() ! diff --git a/src/Model/Connection/GwfInterfaceModel.f90 b/src/Model/Connection/GwfInterfaceModel.f90 index a4dfe220e1a..b8f71e6cdb3 100644 --- a/src/Model/Connection/GwfInterfaceModel.f90 +++ b/src/Model/Connection/GwfInterfaceModel.f90 @@ -94,7 +94,7 @@ subroutine gwfifm_df(this) ! define NPF package call npfOptions%construct() call this%setNpfOptions(npfOptions) - call this%npf%npf_df(this%dis, this%xt3d, 0, npfOptions) + call this%npf%npf_df(this%dis, this%xt3d, 0, 0, npfOptions) call npfOptions%destroy() ! define BUY package @@ -119,7 +119,7 @@ end subroutine gwfifm_df subroutine gwfifm_ar(this) class(GwfInterfaceModelType) :: this !< the GWF interface model - call this%npf%npf_ar(this%ic, this%ibound, this%x) + call this%npf%npf_ar(this%ic, this%vsc, this%ibound, this%x) if (this%inbuy > 0) call this%buy%buy_ar(this%npf, this%ibound) end subroutine gwfifm_ar @@ -147,6 +147,7 @@ subroutine gwfifm_da(this) call mem_deallocate(this%inobs) call mem_deallocate(this%innpf) call mem_deallocate(this%inbuy) + call mem_deallocate(this%invsc) call mem_deallocate(this%insto) call mem_deallocate(this%incsub) call mem_deallocate(this%inmvr) diff --git a/src/Model/GroundWaterFlow/gwf3.f90 b/src/Model/GroundWaterFlow/gwf3.f90 index 9b519ca16d7..dd4964b62af 100644 --- a/src/Model/GroundWaterFlow/gwf3.f90 +++ b/src/Model/GroundWaterFlow/gwf3.f90 @@ -11,6 +11,7 @@ module GwfModule use GwfNpfModule, only: GwfNpfType use Xt3dModule, only: Xt3dType use GwfBuyModule, only: GwfBuyType + use GwfVscModule, only: GwfVscType use GwfHfbModule, only: GwfHfbType use GwfStoModule, only: GwfStoType use GwfCsubModule, only: GwfCsubType @@ -35,6 +36,7 @@ module GwfModule type(GwfNpfType), pointer :: npf => null() ! node property flow package type(Xt3dType), pointer :: xt3d => null() ! xt3d option for npf type(GwfBuyType), pointer :: buy => null() ! buoyancy package + type(GwfVscType), pointer :: vsc => null() ! viscosity package type(GwfStoType), pointer :: sto => null() ! storage package type(GwfCsubType), pointer :: csub => null() ! subsidence package type(GwfOcType), pointer :: oc => null() ! output control package @@ -47,6 +49,7 @@ module GwfModule integer(I4B), pointer :: inoc => null() ! unit number OC integer(I4B), pointer :: innpf => null() ! unit number NPF integer(I4B), pointer :: inbuy => null() ! unit number BUY + integer(I4B), pointer :: invsc => null() ! unit number VSC integer(I4B), pointer :: insto => null() ! unit number STO integer(I4B), pointer :: incsub => null() ! unit number CSUB integer(I4B), pointer :: inmvr => null() ! unit number MVR @@ -96,7 +99,7 @@ module GwfModule &'GHB6 ', 'RCH6 ', 'EVT6 ', 'OBS6 ', 'GNC6 ', & ! 15 &'API6 ', 'CHD6 ', ' ', ' ', ' ', & ! 20 &' ', 'MAW6 ', 'SFR6 ', 'LAK6 ', 'UZF6 ', & ! 25 - &'DISV6', 'MVR6 ', 'CSUB6', 'BUY6 ', ' ', & ! 30 + &'DISV6', 'MVR6 ', 'CSUB6', 'BUY6 ', 'VSC6 ', & ! 30 &70*' '/ contains @@ -122,6 +125,7 @@ subroutine gwf_cr(filename, id, modelname) use GwfNpfModule, only: npf_cr use Xt3dModule, only: xt3d_cr use GwfBuyModule, only: buy_cr + use GwfVscModule, only: vsc_cr use GwfStoModule, only: sto_cr use GwfCsubModule, only: csub_cr use GwfMvrModule, only: mvr_cr @@ -235,6 +239,7 @@ subroutine gwf_cr(filename, id, modelname) call namefile_obj%get_unitnumber('OC6', this%inoc, 1) call namefile_obj%get_unitnumber('NPF6', this%innpf, 1) call namefile_obj%get_unitnumber('BUY6', this%inbuy, 1) + call namefile_obj%get_unitnumber('VSC6', this%invsc, 1) call namefile_obj%get_unitnumber('STO6', this%insto, 1) call namefile_obj%get_unitnumber('CSUB6', this%incsub, 1) call namefile_obj%get_unitnumber('MVR6', this%inmvr, 1) @@ -261,6 +266,7 @@ subroutine gwf_cr(filename, id, modelname) call npf_cr(this%npf, this%name, this%innpf, this%iout) call xt3d_cr(this%xt3d, this%name, this%innpf, this%iout) call buy_cr(this%buy, this%name, this%inbuy, this%iout) + call vsc_cr(this%vsc, this%name, this%invsc, this%iout) call gnc_cr(this%gnc, this%name, this%ingnc, this%iout) call hfb_cr(this%hfb, this%name, this%inhfb, this%iout) call sto_cr(this%sto, this%name, this%insto, this%iout) @@ -306,10 +312,11 @@ subroutine gwf_df(this) ! ! -- Define packages and utility objects call this%dis%dis_df() - call this%npf%npf_df(this%dis, this%xt3d, this%ingnc) + call this%npf%npf_df(this%dis, this%xt3d, this%ingnc, this%invsc) call this%oc%oc_df() call this%budget%budget_df(niunit, 'VOLUME', 'L**3') if (this%inbuy > 0) call this%buy%buy_df(this%dis) + if (this%invsc > 0) call this%vsc%vsc_df(this%dis) if (this%ingnc > 0) call this%gnc%gnc_df(this) ! ! -- Assign or point model members to dis members @@ -416,9 +423,12 @@ subroutine gwf_ar(this) ! ! -- Allocate and read modules attached to model if (this%inic > 0) call this%ic%ic_ar(this%x) - if (this%innpf > 0) call this%npf%npf_ar(this%ic, this%ibound, this%x) + if (this%innpf > 0) call this%npf%npf_ar(this%ic, this%vsc, this%ibound, & + this%x) + if (this%invsc > 0) call this%vsc%vsc_ar(this%ibound) if (this%inbuy > 0) call this%buy%buy_ar(this%npf, this%ibound) - if (this%inhfb > 0) call this%hfb%hfb_ar(this%ibound, this%xt3d, this%dis) + if (this%inhfb > 0) call this%hfb%hfb_ar(this%ibound, this%xt3d, this%dis, & + this%invsc, this%vsc) if (this%insto > 0) call this%sto%sto_ar(this%dis, this%ibound) if (this%incsub > 0) call this%csub%csub_ar(this%dis, this%ibound) if (this%inmvr > 0) call this%mvr%mvr_ar() @@ -439,6 +449,7 @@ subroutine gwf_ar(this) ! -- Read and allocate package call packobj%bnd_ar() if (this%inbuy > 0) call this%buy%buy_ar_bnd(packobj, this%x) + if (this%invsc > 0) call this%vsc%vsc_ar_bnd(packobj) end do ! ! -- return @@ -466,6 +477,7 @@ subroutine gwf_rp(this) ! -- Read and prepare if (this%innpf > 0) call this%npf%npf_rp() if (this%inbuy > 0) call this%buy%buy_rp() + if (this%invsc > 0) call this%vsc%vsc_rp() if (this%inhfb > 0) call this%hfb%hfb_rp() if (this%inoc > 0) call this%oc%oc_rp() if (this%insto > 0) call this%sto%sto_rp() @@ -515,6 +527,7 @@ subroutine gwf_ad(this) end if ! ! -- Advance + if (this%invsc > 0) call this%vsc%vsc_ad() if (this%innpf > 0) call this%npf%npf_ad(this%dis%nodes, this%xold, & this%x, irestore) if (this%insto > 0) call this%sto%sto_ad() @@ -524,6 +537,7 @@ subroutine gwf_ad(this) do ip = 1, this%bndlist%Count() packobj => GetBndFromList(this%bndlist, ip) call packobj%bnd_ad() + if (this%invsc > 0) call this%vsc%vsc_ad_bnd(packobj, this%x) if (isimcheck > 0) then call packobj%bnd_ck() end if @@ -1123,24 +1137,31 @@ subroutine gwf_ot_dv(this, idvsave, idvprint, ipflag) integer(I4B), intent(inout) :: ipflag class(BndType), pointer :: packobj integer(I4B) :: ip - + ! ! -- Save compaction to binary file if (this%incsub > 0) call this%csub%csub_ot_dv(idvsave, idvprint) - + ! ! -- save density to binary file if (this%inbuy > 0) then call this%buy%buy_ot_dv(idvsave) end if - + ! + ! -- save viscosity to binary file + if (this%invsc > 0) then + call this%vsc%vsc_ot_dv(idvsave) + end if + ! ! -- Print advanced package dependent variables do ip = 1, this%bndlist%Count() packobj => GetBndFromList(this%bndlist, ip) call packobj%bnd_ot_dv(idvsave, idvprint) end do - + ! ! -- save head and print head call this%oc%oc_ot(ipflag) - + ! + ! -- Return + return end subroutine gwf_ot_dv subroutine gwf_ot_bdsummary(this, ibudfl, ipflag) @@ -1207,6 +1228,7 @@ subroutine gwf_da(this) call this%npf%npf_da() call this%xt3d%xt3d_da() call this%buy%buy_da() + call this%vsc%vsc_da() call this%gnc%gnc_da() call this%sto%sto_da() call this%csub%csub_da() @@ -1222,6 +1244,7 @@ subroutine gwf_da(this) deallocate (this%npf) deallocate (this%xt3d) deallocate (this%buy) + deallocate (this%vsc) deallocate (this%gnc) deallocate (this%sto) deallocate (this%csub) @@ -1244,6 +1267,7 @@ subroutine gwf_da(this) call mem_deallocate(this%inobs) call mem_deallocate(this%innpf) call mem_deallocate(this%inbuy) + call mem_deallocate(this%invsc) call mem_deallocate(this%insto) call mem_deallocate(this%incsub) call mem_deallocate(this%inmvr) @@ -1336,6 +1360,7 @@ subroutine allocate_scalars(this, modelname) call mem_allocate(this%inoc, 'INOC', this%memoryPath) call mem_allocate(this%innpf, 'INNPF', this%memoryPath) call mem_allocate(this%inbuy, 'INBUY', this%memoryPath) + call mem_allocate(this%invsc, 'INVSC', this%memoryPath) call mem_allocate(this%insto, 'INSTO', this%memoryPath) call mem_allocate(this%incsub, 'INCSUB', this%memoryPath) call mem_allocate(this%inmvr, 'INMVR', this%memoryPath) @@ -1349,6 +1374,7 @@ subroutine allocate_scalars(this, modelname) this%inoc = 0 this%innpf = 0 this%inbuy = 0 + this%invsc = 0 this%insto = 0 this%incsub = 0 this%inmvr = 0 diff --git a/src/Model/GroundWaterFlow/gwf3hfb8.f90 b/src/Model/GroundWaterFlow/gwf3hfb8.f90 index 2687b00ad3f..ffec99e84bd 100644 --- a/src/Model/GroundWaterFlow/gwf3hfb8.f90 +++ b/src/Model/GroundWaterFlow/gwf3hfb8.f90 @@ -3,6 +3,7 @@ module GwfHfbModule use KindModule, only: DP, I4B use Xt3dModule, only: Xt3dType + use GwfVscModule, only: GwfVscType use NumericalPackageModule, only: NumericalPackageType use BlockParserModule, only: BlockParserType use BaseDisModule, only: DisBaseType @@ -14,27 +15,34 @@ module GwfHfbModule public :: hfb_cr type, extends(NumericalPackageType) :: GwfHfbType - integer(I4B), pointer :: maxhfb => null() !max number of hfb's - integer(I4B), pointer :: nhfb => null() !number of hfb's - integer(I4B), dimension(:), pointer, contiguous :: noden => null() !first cell - integer(I4B), dimension(:), pointer, contiguous :: nodem => null() !second cell - integer(I4B), dimension(:), pointer, contiguous :: idxloc => null() !position in model ja - real(DP), dimension(:), pointer, contiguous :: hydchr => null() !hydraulic characteristic of the barrier - real(DP), dimension(:), pointer, contiguous :: csatsav => null() !value of condsat prior to hfb modification - real(DP), dimension(:), pointer, contiguous :: condsav => null() !saved conductance of combined npf and hfb - type(Xt3dType), pointer :: xt3d => null() !pointer to xt3d object - ! - integer(I4B), dimension(:), pointer, contiguous :: ibound => null() !pointer to model ibound - integer(I4B), dimension(:), pointer, contiguous :: icelltype => null() !pointer to model icelltype - integer(I4B), dimension(:), pointer, contiguous :: ihc => null() !pointer to model ihc - integer(I4B), dimension(:), pointer, contiguous :: ia => null() !pointer to model ia - integer(I4B), dimension(:), pointer, contiguous :: ja => null() !pointer to model ja - integer(I4B), dimension(:), pointer, contiguous :: jas => null() !pointer to model jas - integer(I4B), dimension(:), pointer, contiguous :: isym => null() !pointer to model isym - real(DP), dimension(:), pointer, contiguous :: condsat => null() !pointer to model condsat - real(DP), dimension(:), pointer, contiguous :: top => null() !pointer to model top - real(DP), dimension(:), pointer, contiguous :: bot => null() !pointer to model bot - real(DP), dimension(:), pointer, contiguous :: hwva => null() !pointer to model hwva + + type(GwfVscType), pointer :: vsc => null() !< viscosity object + integer(I4B), pointer :: maxhfb => null() !< max number of hfb's + integer(I4B), pointer :: nhfb => null() !< number of hfb's + integer(I4B), dimension(:), pointer, contiguous :: noden => null() !< first cell + integer(I4B), dimension(:), pointer, contiguous :: nodem => null() !< second cell + integer(I4B), dimension(:), pointer, contiguous :: idxloc => null() !< position in model ja + real(DP), dimension(:), pointer, contiguous :: hydchr => null() !< hydraulic characteristic of the barrier + real(DP), dimension(:), pointer, contiguous :: csatsav => null() !< value of condsat prior to hfb modification + real(DP), dimension(:), pointer, contiguous :: condsav => null() !< saved conductance of combined npf and hfb + type(Xt3dType), pointer :: xt3d => null() !< pointer to xt3d object + ! + integer(I4B), dimension(:), pointer, contiguous :: ibound => null() !< pointer to model ibound + integer(I4B), dimension(:), pointer, contiguous :: icelltype => null() !< pointer to model icelltype + integer(I4B), dimension(:), pointer, contiguous :: ihc => null() !< pointer to model ihc + integer(I4B), dimension(:), pointer, contiguous :: ia => null() !< pointer to model ia + integer(I4B), dimension(:), pointer, contiguous :: ja => null() !< pointer to model ja + integer(I4B), dimension(:), pointer, contiguous :: jas => null() !< pointer to model jas + integer(I4B), dimension(:), pointer, contiguous :: isym => null() !< pointer to model isym + real(DP), dimension(:), pointer, contiguous :: condsat => null() !< pointer to model condsat + real(DP), dimension(:), pointer, contiguous :: top => null() !< pointer to model top + real(DP), dimension(:), pointer, contiguous :: bot => null() !< pointer to model bot + real(DP), dimension(:), pointer, contiguous :: hwva => null() !< pointer to model hwva + real(DP), dimension(:), pointer, contiguous :: hnew => null() !< pointer to model xnew + ! + ! -- viscosity flag + integer(I4B), pointer :: ivsc => null() !< flag indicating if viscosity is active in the model + contains procedure :: hfb_ar procedure :: hfb_rp @@ -49,6 +57,7 @@ module GwfHfbModule procedure, private :: check_data procedure, private :: condsat_reset procedure, private :: condsat_modify + end type GwfHfbType contains @@ -87,7 +96,7 @@ subroutine hfb_cr(hfbobj, name_model, inunit, iout) return end subroutine hfb_cr - subroutine hfb_ar(this, ibound, xt3d, dis) + subroutine hfb_ar(this, ibound, xt3d, dis, invsc, vsc) ! ****************************************************************************** ! hfb_ar -- Allocate and read ! ****************************************************************************** @@ -101,7 +110,10 @@ subroutine hfb_ar(this, ibound, xt3d, dis) class(GwfHfbType) :: this integer(I4B), dimension(:), pointer, contiguous :: ibound type(Xt3dType), pointer :: xt3d - class(DisBaseType), pointer, intent(inout) :: dis + class(DisBaseType), pointer, intent(inout) :: dis !< discretization package + integer(I4B), pointer :: invsc !< indicates if viscosity package is active + type(GwfVscType), pointer, intent(in) :: vsc !< viscosity package + ! -- local ! -- formats character(len=*), parameter :: fmtheader = & "(1x, /1x, 'HFB -- HORIZONTAL FLOW BARRIER PACKAGE, VERSION 8, ', & @@ -133,6 +145,16 @@ subroutine hfb_ar(this, ibound, xt3d, dis) call this%read_dimensions() call this%allocate_arrays() ! + ! -- If vsc package active, set ivsc + if (invsc /= 0) then + this%ivsc = 1 + this%vsc => vsc + ! + ! -- Notify user via listing file viscosity accounted for by HFB package + write (this%iout, '(/1x,a,a)') 'Viscosity active in ', & + trim(this%filtyp)//' Package calculations: '//trim(adjustl(this%packName)) + end if + ! ! -- return return end subroutine hfb_ar @@ -213,7 +235,7 @@ subroutine hfb_fc(this, kiter, njasln, amat, idxglo, rhs, hnew) ! SPECIFICATIONS: ! ------------------------------------------------------------------------------ ! -- modules - use ConstantsModule, only: DHALF, DZERO + use ConstantsModule, only: DHALF, DZERO, DONE ! -- dummy class(GwfHfbType) :: this integer(I4B) :: kiter @@ -231,8 +253,11 @@ subroutine hfb_fc(this, kiter, njasln, amat, idxglo, rhs, hnew) real(DP) :: cond, condhfb, aterm real(DP) :: fawidth, faheight real(DP) :: topn, topm, botn, botm + real(DP) :: viscratio ! ------------------------------------------------------------------------------ ! + ! -- initialize variables + viscratio = DONE nodes = this%dis%nodes nja = this%dis%con%nja if (associated(this%xt3d%ixt3d)) then @@ -248,7 +273,10 @@ subroutine hfb_fc(this, kiter, njasln, amat, idxglo, rhs, hnew) m = max(this%noden(ihfb), this%nodem(ihfb)) ! -- Skip if either cell is inactive. if (this%ibound(n) == 0 .or. this%ibound(m) == 0) cycle - !!! if(this%icelltype(n) == 1 .or. this%icelltype(m) == 1) then + !!! if(this%icelltype(n) == 1 .or. this%icelltype(m) == 1) then + if (this%ivsc /= 0) then + call this%vsc%get_visc_ratio(n, m, hnew(n), hnew(m), viscratio) + end if ! -- Compute scale factor for hfb correction if (this%hydchr(ihfb) > DZERO) then if (this%inewton == 0) then @@ -269,9 +297,10 @@ subroutine hfb_fc(this, kiter, njasln, amat, idxglo, rhs, hnew) faheight = DHALF * ((topn - botn) + (topm - botm)) end if fawidth = this%hwva(this%jas(ipos)) - condhfb = this%hydchr(ihfb) * fawidth * faheight + condhfb = this%hydchr(ihfb) * viscratio * & + fawidth * faheight else - condhfb = this%hydchr(ihfb) + condhfb = this%hydchr(ihfb) * viscratio end if else condhfb = this%hydchr(ihfb) @@ -291,7 +320,13 @@ subroutine hfb_fc(this, kiter, njasln, amat, idxglo, rhs, hnew) n = this%noden(ihfb) m = this%nodem(ihfb) if (this%ibound(n) == 0 .or. this%ibound(m) == 0) cycle - if (this%icelltype(n) == 1 .or. this%icelltype(m) == 1) then + ! + if (this%ivsc /= 0) then + call this%vsc%get_visc_ratio(n, m, hnew(n), hnew(m), viscratio) + end if + ! + if (this%icelltype(n) == 1 .or. this%icelltype(m) == 1 .or. & + this%ivsc /= 0) then ! ! -- Calculate hfb conductance topn = this%top(n) @@ -311,7 +346,8 @@ subroutine hfb_fc(this, kiter, njasln, amat, idxglo, rhs, hnew) end if if (this%hydchr(ihfb) > DZERO) then fawidth = this%hwva(this%jas(ipos)) - condhfb = this%hydchr(ihfb) * fawidth * faheight + condhfb = this%hydchr(ihfb) * viscratio * & + fawidth * faheight cond = aterm * condhfb / (aterm + condhfb) else cond = -aterm * this%hydchr(ihfb) @@ -351,7 +387,7 @@ subroutine hfb_cq(this, hnew, flowja) ! SPECIFICATIONS: ! ------------------------------------------------------------------------------ ! -- modules - use ConstantsModule, only: DHALF, DZERO + use ConstantsModule, only: DHALF, DZERO, DONE ! -- dummy class(GwfHfbType) :: this real(DP), intent(inout), dimension(:) :: hnew @@ -365,8 +401,12 @@ subroutine hfb_cq(this, hnew, flowja) real(DP) :: condhfb real(DP) :: fawidth, faheight real(DP) :: topn, topm, botn, botm + real(DP) :: viscratio ! ------------------------------------------------------------------------------ ! + ! -- initialize viscratio + viscratio = DONE + ! if (associated(this%xt3d%ixt3d)) then ixt3d = this%xt3d%ixt3d else @@ -380,7 +420,11 @@ subroutine hfb_cq(this, hnew, flowja) m = max(this%noden(ihfb), this%nodem(ihfb)) ! -- Skip if either cell is inactive. if (this%ibound(n) == 0 .or. this%ibound(m) == 0) cycle - !!! if(this%icelltype(n) == 1 .or. this%icelltype(m) == 1) then + !!! if(this%icelltype(n) == 1 .or. this%icelltype(m) == 1) then + if (this%ivsc /= 0) then + call this%vsc%get_visc_ratio(n, m, hnew(n), hnew(m), viscratio) + end if + ! ! -- Compute scale factor for hfb correction if (this%hydchr(ihfb) > DZERO) then if (this%inewton == 0) then @@ -401,7 +445,8 @@ subroutine hfb_cq(this, hnew, flowja) faheight = DHALF * ((topn - botn) + (topm - botm)) end if fawidth = this%hwva(this%jas(ipos)) - condhfb = this%hydchr(ihfb) * fawidth * faheight + condhfb = this%hydchr(ihfb) * viscratio * & + fawidth * faheight else condhfb = this%hydchr(ihfb) end if @@ -420,8 +465,11 @@ subroutine hfb_cq(this, hnew, flowja) n = this%noden(ihfb) m = this%nodem(ihfb) if (this%ibound(n) == 0 .or. this%ibound(m) == 0) cycle - if (this%icelltype(n) == 1 .or. this%icelltype(m) == 1) then + if (this%icelltype(n) == 1 .or. this%icelltype(m) == 1 .or. & + this%ivsc /= 0) then ipos = this%dis%con%getjaindex(n, m) + ! + ! -- condsav already accnts for visc adjustment cond = this%condsav(ihfb) qnm = cond * (hnew(m) - hnew(n)) flowja(ipos) = qnm @@ -456,6 +504,7 @@ subroutine hfb_da(this) ! -- Scalars call mem_deallocate(this%maxhfb) call mem_deallocate(this%nhfb) + call mem_deallocate(this%ivsc) ! ! -- Arrays if (this%inunit > 0) then @@ -484,6 +533,7 @@ subroutine hfb_da(this) this%top => null() this%bot => null() this%hwva => null() + this%vsc => null() ! ! -- return return @@ -509,9 +559,13 @@ subroutine allocate_scalars(this) call mem_allocate(this%maxhfb, 'MAXHFB', this%memoryPath) call mem_allocate(this%nhfb, 'NHFB', this%memoryPath) ! + ! -- allocate flag for determining if vsc active + call mem_allocate(this%ivsc, 'IVSC', this%memoryPath) + ! ! -- initialize this%maxhfb = 0 this%nhfb = 0 + this%ivsc = 0 ! ! -- return return @@ -855,6 +909,7 @@ subroutine condsat_modify(this) this%csatsav(ihfb) = cond n = this%noden(ihfb) m = this%nodem(ihfb) + ! if (this%inewton == 1 .or. & (this%icelltype(n) == 0 .and. this%icelltype(m) == 0)) then ! @@ -870,7 +925,8 @@ subroutine condsat_modify(this) end if if (this%hydchr(ihfb) > DZERO) then fawidth = this%hwva(this%jas(ipos)) - condhfb = this%hydchr(ihfb) * fawidth * faheight + condhfb = this%hydchr(ihfb) * & + fawidth * faheight cond = cond * condhfb / (cond + condhfb) else cond = -cond * this%hydchr(ihfb) diff --git a/src/Model/GroundWaterFlow/gwf3lak8.f90 b/src/Model/GroundWaterFlow/gwf3lak8.f90 index d11275dadbf..cfb921248fa 100644 --- a/src/Model/GroundWaterFlow/gwf3lak8.f90 +++ b/src/Model/GroundWaterFlow/gwf3lak8.f90 @@ -190,6 +190,9 @@ module LakModule integer(I4B), pointer :: idense real(DP), dimension(:, :), pointer, contiguous :: denseterms => null() ! + ! -- viscosity variables + real(DP), dimension(:, :), pointer, contiguous :: viscratios => null() !< viscosity ratios (1: lak vsc ratio; 2: gwf vsc ratio) + ! ! -- type bound procedures contains procedure :: lak_allocate_scalars @@ -267,6 +270,8 @@ module LakModule ! -- density procedure :: lak_activate_density procedure, private :: lak_calculate_density_exchange + ! -- viscosity + procedure :: lak_activate_viscosity end type LakType contains @@ -375,6 +380,7 @@ subroutine lak_allocate_scalars(this) this%bditems = 11 this%cbcauxitems = 1 this%idense = 0 + this%ivsc = 0 ! ! -- return return @@ -445,6 +451,9 @@ subroutine lak_allocate_arrays(this) ! -- allocate denseterms to size 0 call mem_allocate(this%denseterms, 3, 0, 'DENSETERMS', this%memoryPath) ! + ! -- allocate viscratios to size 0 + call mem_allocate(this%viscratios, 2, 0, 'VISCRATIOS', this%memoryPath) + ! ! -- return return end subroutine lak_allocate_arrays @@ -2338,9 +2347,11 @@ subroutine lak_calculate_conn_conductance(this, ilak, iconn, stage, head, cond) real(DP) :: botl real(DP) :: sat real(DP) :: wa + real(DP) :: vscratio ! -- formats ! ------------------------------------------------------------------------------ cond = DZERO + vscratio = DONE topl = this%telev(iconn) botl = this%belev(iconn) call this%lak_calculate_cond_head(iconn, stage, head, vv) @@ -2370,7 +2381,18 @@ subroutine lak_calculate_conn_conductance(this, ilak, iconn, stage, head, cond) end if sat = wa end if - cond = sat * this%satcond(iconn) + ! + ! -- account for viscosity effects (if vsc active) + if (this%ivsc == 1) then + ! flow from lake to aquifer + if (stage > head) then + vscratio = this%viscratios(1, iconn) + ! flow from aquifer to lake + else + vscratio = this%viscratios(2, iconn) + end if + end if + cond = sat * this%satcond(iconn) * vscratio ! ! -- return return @@ -4451,6 +4473,7 @@ subroutine lak_da(this) call mem_deallocate(this%qleak) call mem_deallocate(this%qsto) call mem_deallocate(this%denseterms) + call mem_deallocate(this%viscratios) ! ! -- tables if (this%ntables > 0) then @@ -6361,6 +6384,36 @@ subroutine lak_activate_density(this) return end subroutine lak_activate_density + !> @brief Activate viscosity terms + !! + !! Method to activate addition of viscosity terms for a LAK package reach. + !! + !< + subroutine lak_activate_viscosity(this) + ! -- modules + use MemoryManagerModule, only: mem_reallocate + ! -- dummy variables + class(LakType), intent(inout) :: this !< LakType object + ! -- local variables + integer(I4B) :: i + integer(I4B) :: j + ! + ! -- Set ivsc and reallocate viscratios to be of size MAXBOUND + this%ivsc = 1 + call mem_reallocate(this%viscratios, 2, this%MAXBOUND, 'VISCRATIOS', & + this%memoryPath) + do i = 1, this%maxbound + do j = 1, 2 + this%viscratios(j, i) = DONE + end do + end do + write (this%iout, '(/1x,a)') 'VISCOSITY HAS BEEN ACTIVATED FOR LAK & + &PACKAGE: '//trim(adjustl(this%packName)) + ! + ! -- return + return + end subroutine lak_activate_viscosity + subroutine lak_calculate_density_exchange(this, iconn, stage, head, cond, & botl, flow, gwfhcof, gwfrhs) ! ****************************************************************************** diff --git a/src/Model/GroundWaterFlow/gwf3maw8.f90 b/src/Model/GroundWaterFlow/gwf3maw8.f90 index 415d79e6d13..295b84d8472 100644 --- a/src/Model/GroundWaterFlow/gwf3maw8.f90 +++ b/src/Model/GroundWaterFlow/gwf3maw8.f90 @@ -158,6 +158,9 @@ module MawModule integer(I4B), pointer :: idense real(DP), dimension(:, :), pointer, contiguous :: denseterms => null() ! + ! -- viscosity variables + real(DP), dimension(:, :), pointer, contiguous :: viscratios => null() !< viscosity ratios (1: maw vsc ratio; 2: gwf vsc ratio) + ! ! -- type bound procedures contains procedure :: maw_allocate_scalars @@ -213,6 +216,8 @@ module MawModule ! -- MAW reduced flow outputs procedure, private :: maw_redflow_csv_init procedure, private :: maw_redflow_csv_write + ! -- viscosity + procedure :: maw_activate_viscosity end type MawType contains @@ -318,6 +323,7 @@ subroutine maw_allocate_scalars(this) this%kappa = DEM4 this%cbcauxitems = 1 this%idense = 0 + this%ivsc = 0 ! ! -- return return @@ -521,6 +527,9 @@ subroutine maw_allocate_well_conn_arrays(this) ! -- allocate denseterms to size 0 call mem_allocate(this%denseterms, 3, 0, 'DENSETERMS', this%memoryPath) ! + ! -- allocate viscratios to size 0 + call mem_allocate(this%viscratios, 2, 0, 'VISCRATIOS', this%memoryPath) + ! ! -- return return end subroutine maw_allocate_well_conn_arrays @@ -3003,6 +3012,7 @@ subroutine maw_da(this) call mem_deallocate(this%qsto) call mem_deallocate(this%qconst) call mem_deallocate(this%denseterms) + call mem_deallocate(this%viscratios) call mem_deallocate(this%idxlocnode) call mem_deallocate(this%idxdglo) call mem_deallocate(this%idxoffdglo) @@ -3858,10 +3868,12 @@ subroutine maw_calculate_conn_terms(this, n, j, icflow, cmaw, cterm, term, & real(DP) :: hbar real(DP) :: drterm real(DP) :: dhbarterm + real(DP) :: vscratio ! ------------------------------------------------------------------------------ ! ! -- initialize terms cterm = DZERO + vscratio = DONE icflow = 0 if (present(term2)) then inewton = 1 @@ -3877,9 +3889,19 @@ subroutine maw_calculate_conn_terms(this, n, j, icflow, cmaw, cterm, term, & tmaw = this%topscrn(jpos) bmaw = this%botscrn(jpos) ! + ! -- if vsc active, select appropriate viscosity ratio + if (this%ivsc == 1) then + ! flow out of well (flow is negative) + if (flow < 0) then + vscratio = this%viscratios(1, igwfnode) + else + vscratio = this%viscratios(2, igwfnode) + end if + end if + ! ! -- calculate saturation call this%maw_calculate_saturation(n, j, igwfnode, sat) - cmaw = this%satcond(jpos) * sat + cmaw = this%satcond(jpos) * vscratio * sat ! ! -- set upstream head, term, and term2 if returning newton terms if (inewton == 1) then @@ -3928,14 +3950,14 @@ subroutine maw_calculate_conn_terms(this, n, j, icflow, cmaw, cterm, term, & ! -- maw is upstream if (hmaw > hgwf) then hbar = sQuadratic0sp(hgwf, en, this%satomega) - term = drterm * this%satcond(jpos) * (hbar - hmaw) + term = drterm * this%satcond(jpos) * vscratio * (hbar - hmaw) dhbarterm = sQuadratic0spDerivative(hgwf, en, this%satomega) term2 = cmaw * (dhbarterm - DONE) ! ! -- gwf is upstream else hbar = sQuadratic0sp(hmaw, en, this%satomega) - term = -drterm * this%satcond(jpos) * (hgwf - hbar) + term = -drterm * this%satcond(jpos) * vscratio * (hgwf - hbar) dhbarterm = sQuadratic0spDerivative(hmaw, en, this%satomega) term2 = cmaw * (DONE - dhbarterm) end if @@ -3944,7 +3966,7 @@ subroutine maw_calculate_conn_terms(this, n, j, icflow, cmaw, cterm, term, & ! ! -- flow is not corrected, so calculate term for newton formulation if (inewton /= 0) then - term = drterm * this%satcond(jpos) * (hgwf - hmaw) + term = drterm * this%satcond(jpos) * vscratio * (hgwf - hmaw) end if end if ! @@ -4165,20 +4187,32 @@ subroutine maw_calculate_qpot(this, n, qnet) real(DP) :: bmaw real(DP) :: htmp real(DP) :: hv + real(DP) :: vscratio ! -- format ! ------------------------------------------------------------------------------ ! ! -- initialize qnet and htmp qnet = DZERO + vscratio = DONE htmp = this%shutofflevel(n) ! + ! -- if vsc active, select appropriate viscosity ratio + if (this%ivsc == 1) then + ! flow out of well (flow is negative) + if (qnet < 0) then + vscratio = this%viscratios(1, igwfnode) + else + vscratio = this%viscratios(2, igwfnode) + end if + end if + ! ! -- calculate discharge to flowing wells if (this%iflowingwells > 0) then if (this%fwcond(n) > DZERO) then bt = this%fwelev(n) tp = bt + this%fwrlen(n) scale = sQSaturation(tp, bt, htmp) - cfw = scale * this%fwcond(n) + cfw = scale * this%fwcond(n) * this%viscratios(2, n) this%ifwdischarge(n) = 0 if (cfw > DZERO) then this%ifwdischarge(n) = 1 @@ -4203,7 +4237,7 @@ subroutine maw_calculate_qpot(this, n, qnet) jpos = this%get_jpos(n, j) igwfnode = this%get_gwfnode(n, j) call this%maw_calculate_saturation(n, j, igwfnode, sat) - cmaw = this%satcond(jpos) * sat + cmaw = this%satcond(jpos) * vscratio * sat hgwf = this%xnew(igwfnode) bmaw = this%botscrn(jpos) hv = htmp @@ -4843,6 +4877,36 @@ subroutine maw_activate_density(this) return end subroutine maw_activate_density + !> @brief Activate viscosity terms + !! + !! Method to activate addition of viscosity terms for a MAW package reach. + !! + !< + subroutine maw_activate_viscosity(this) + ! -- modules + use MemoryManagerModule, only: mem_reallocate + ! -- dummy variables + class(MawType), intent(inout) :: this !< MawType object + ! -- local variables + integer(I4B) :: i + integer(I4B) :: j + ! + ! -- Set ivsc and reallocate viscratios to be of size MAXBOUND + this%ivsc = 1 + call mem_reallocate(this%viscratios, 2, this%MAXBOUND, 'VISCRATIOS', & + this%memoryPath) + do i = 1, this%maxbound + do j = 1, 2 + this%viscratios(j, i) = DONE + end do + end do + write (this%iout, '(/1x,a)') 'VISCOSITY HAS BEEN ACTIVATED FOR MAW & + &PACKAGE: '//trim(adjustl(this%packName)) + ! + ! -- return + return + end subroutine maw_activate_viscosity + subroutine maw_calculate_density_exchange(this, iconn, hmaw, hgwf, cond, & bmaw, flow, hcofterm, rhsterm) ! ****************************************************************************** diff --git a/src/Model/GroundWaterFlow/gwf3npf8.f90 b/src/Model/GroundWaterFlow/gwf3npf8.f90 index dcfdf99bc31..889c49e3fc9 100644 --- a/src/Model/GroundWaterFlow/gwf3npf8.f90 +++ b/src/Model/GroundWaterFlow/gwf3npf8.f90 @@ -11,6 +11,7 @@ module GwfNpfModule use GwfNpfOptionsModule, only: GwfNpfOptionsType use BaseDisModule, only: DisBaseType use GwfIcModule, only: GwfIcType + use GwfVscModule, only: GwfVscType use Xt3dModule, only: Xt3dType use BlockParserModule, only: BlockParserType use InputOutputModule, only: GetUnit, openfile @@ -33,6 +34,7 @@ module GwfNpfModule type, extends(NumericalPackageType) :: GwfNpfType type(GwfIcType), pointer :: ic => null() !< initial conditions object + type(GwfVscType), pointer :: vsc => null() !< viscosity object type(Xt3dType), pointer :: xt3d => null() !< xt3d pointer integer(I4B), pointer :: iname => null() !< length of variable names character(len=24), dimension(:), pointer :: aname => null() !< variable names @@ -65,6 +67,9 @@ module GwfNpfModule real(DP), dimension(:), pointer, contiguous :: k11 => null() !< hydraulic conductivity; if anisotropic, then this is Kx prior to rotation real(DP), dimension(:), pointer, contiguous :: k22 => null() !< hydraulic conductivity; if specified then this is Ky prior to rotation real(DP), dimension(:), pointer, contiguous :: k33 => null() !< hydraulic conductivity; if specified then this is Kz prior to rotation + real(DP), dimension(:), pointer, contiguous :: k11input => null() !< hydraulic conductivity originally specified by user prior to TVK or VSC modification + real(DP), dimension(:), pointer, contiguous :: k22input => null() !< hydraulic conductivity originally specified by user prior to TVK or VSC modification + real(DP), dimension(:), pointer, contiguous :: k33input => null() !< hydraulic conductivity originally specified by user prior to TVK or VSC modification integer(I4B), pointer :: iavgkeff => null() !< effective conductivity averaging (0: harmonic, 1: arithmetic) integer(I4B), pointer :: ik22 => null() !< flag that k22 is specified integer(I4B), pointer :: ik33 => null() !< flag that k33 is specified @@ -92,6 +97,7 @@ module GwfNpfModule real(DP), dimension(:, :), pointer, contiguous :: propsedge => null() !< edge properties (Q, area, nx, ny, distance) ! integer(I4B), pointer :: intvk => null() ! TVK (time-varying K) unit number (0 if unused) + integer(I4B), pointer :: invsc => null() ! VSC (viscosity) unit number (0 if unused); viscosity leads to time-varying K's type(TvkType), pointer :: tvk => null() ! TVK object integer(I4B), pointer :: kchangeper => null() ! last stress period in which any node K (or K22, or K33) values were changed (0 if unchanged from start of simulation) integer(I4B), pointer :: kchangestp => null() ! last time step in which any node K (or K22, or K33) values were changed (0 if unchanged from start of simulation) @@ -117,6 +123,7 @@ module GwfNpfModule procedure, private :: wd => sgwf_npf_wetdry procedure, private :: wdmsg => sgwf_npf_wdmsg procedure :: allocate_scalars + procedure, private :: store_original_k_arrays procedure, private :: allocate_arrays procedure, private :: source_options procedure, private :: source_griddata @@ -152,7 +159,7 @@ subroutine npf_cr(npfobj, name_model, inunit, iout) use IdmMf6FileLoaderModule, only: input_load use ConstantsModule, only: LENPACKAGETYPE ! -- dummy - type(GwfNpftype), pointer :: npfobj + type(GwfNpfType), pointer :: npfobj character(len=*), intent(in) :: name_model integer(I4B), intent(in) :: inunit integer(I4B), intent(in) :: iout @@ -201,7 +208,7 @@ end subroutine npf_cr !! should be passed. A consistency check is performed, and finally !! xt3d_df is called, when enabled. !< - subroutine npf_df(this, dis, xt3d, ingnc, npf_options) + subroutine npf_df(this, dis, xt3d, ingnc, invsc, npf_options) ! ****************************************************************************** ! npf_df -- Define ! ****************************************************************************** @@ -216,6 +223,7 @@ subroutine npf_df(this, dis, xt3d, ingnc, npf_options) class(DisBaseType), pointer, intent(inout) :: dis !< the pointer to the discretization type(Xt3dType), pointer :: xt3d !< the pointer to the XT3D 'package' integer(I4B), intent(in) :: ingnc !< ghostnodes enabled? (>0 means yes) + integer(I4B), intent(in) :: invsc !< viscosity enabled? (>0 means yes) type(GwfNpfOptionsType), optional, intent(in) :: npf_options !< the optional options, for when not constructing from file ! -- local ! -- data @@ -224,6 +232,9 @@ subroutine npf_df(this, dis, xt3d, ingnc, npf_options) ! -- Set a pointer to dis this%dis => dis ! + ! -- Set flag signifying whether vsc is active + if (invsc > 0) this%invsc = invsc + ! if (.not. present(npf_options)) then ! ! -- source options @@ -310,7 +321,7 @@ end subroutine npf_mc !! Allocate remaining package arrays, preprocess the input data and !! call *_ar on xt3d, when active. !< - subroutine npf_ar(this, ic, ibound, hnew) + subroutine npf_ar(this, ic, vsc, ibound, hnew) ! ****************************************************************************** ! npf_ar -- Allocate and Read ! ****************************************************************************** @@ -322,6 +333,7 @@ subroutine npf_ar(this, ic, ibound, hnew) ! -- dummy class(GwfNpftype) :: this !< instance of the NPF package type(GwfIcType), pointer, intent(in) :: ic !< initial conditions + type(GwfVscType), pointer, intent(in) :: vsc !< viscosity package integer(I4B), dimension(:), pointer, contiguous, intent(inout) :: ibound !< model ibound array real(DP), dimension(:), pointer, contiguous, intent(inout) :: hnew !< pointer to model head array ! -- local @@ -346,6 +358,28 @@ subroutine npf_ar(this, ic, ibound, hnew) end do end if ! + ! -- Store pointer to VSC if active + if (this%invsc /= 0) then + this%vsc => vsc + end if + + ! + ! -- allocate arrays to store original user input in case TVK/VSC modify them + if (this%invsc > 0) then + ! + ! -- Reallocate arrays that store user-input values. + call mem_reallocate(this%k11input, this%dis%nodes, 'K11INPUT', & + this%memoryPath) + call mem_reallocate(this%k22input, this%dis%nodes, 'K22INPUT', & + this%memoryPath) + call mem_reallocate(this%k33input, this%dis%nodes, 'K33INPUT', & + this%memoryPath) + ! Allocate arrays that will store the original K values. When VSC active, + ! the current Kxx arrays carry the viscosity-adjusted K values. + ! This approach leverages existing functionality that makes use of K. + call this%store_original_k_arrays(this%dis%nodes, this%dis%njas) + end if + ! ! -- preprocess data call this%preprocess_input() ! @@ -425,6 +459,12 @@ subroutine npf_ad(this, nodes, hold, hnew, irestore) call this%tvk%ad() end if ! + ! -- VSC + ! -- Hit the TVK-updated K's with VSC correction before calling/updating condsat + if (this%invsc /= 0) then + call this%vsc%update_k_with_vsc() + end if + ! ! -- If any K values have changed, we need to update CONDSAT or XT3D arrays if (this%kchangeper == kper .and. this%kchangestp == kstp) then if (this%ixt3d == 0) then @@ -1063,6 +1103,11 @@ subroutine npf_da(this) deallocate (this%tvk) end if ! + ! -- VSC + if (this%invsc /= 0) then + nullify (this%vsc) + end if + ! ! -- Strings ! ! -- Scalars @@ -1100,6 +1145,7 @@ subroutine npf_da(this) call mem_deallocate(this%ik22overk) call mem_deallocate(this%ik33overk) call mem_deallocate(this%intvk) + call mem_deallocate(this%invsc) call mem_deallocate(this%kchangeper) call mem_deallocate(this%kchangestp) ! @@ -1110,6 +1156,9 @@ subroutine npf_da(this) call mem_deallocate(this%k11) call mem_deallocate(this%k22) call mem_deallocate(this%k33) + call mem_deallocate(this%k11input) + call mem_deallocate(this%k22input) + call mem_deallocate(this%k33input) call mem_deallocate(this%sat) call mem_deallocate(this%condsat) call mem_deallocate(this%wetdry) @@ -1129,13 +1178,13 @@ subroutine npf_da(this) return end subroutine npf_da + !> @ brief Allocate scalars + !! + !! Allocate and initialize scalars for the VSC package. The base model + !! allocate scalars method is also called. + !! + !< subroutine allocate_scalars(this) -! ****************************************************************************** -! allocate_scalars -- Allocate scalar pointer variables -! ****************************************************************************** -! -! SPECIFICATIONS: -! ------------------------------------------------------------------------------ ! -- modules use MemoryHelperModule, only: create_mem_path ! -- dummy @@ -1179,6 +1228,7 @@ subroutine allocate_scalars(this) call mem_allocate(this%nedges, 'NEDGES', this%memoryPath) call mem_allocate(this%lastedge, 'LASTEDGE', this%memoryPath) call mem_allocate(this%intvk, 'INTVK', this%memoryPath) + call mem_allocate(this%invsc, 'INVSC', this%memoryPath) call mem_allocate(this%kchangeper, 'KCHANGEPER', this%memoryPath) call mem_allocate(this%kchangestp, 'KCHANGESTP', this%memoryPath) ! @@ -1220,6 +1270,7 @@ subroutine allocate_scalars(this) this%nedges = 0 this%lastedge = 0 this%intvk = 0 + this%invsc = 0 this%kchangeper = 0 this%kchangestp = 0 ! @@ -1230,6 +1281,40 @@ subroutine allocate_scalars(this) return end subroutine allocate_scalars + !> @ brief Store backup copy of hydraulic conductivity when the VSC + !! package is activate + !! + !! The K arrays (K11, etc.) get multiplied by the viscosity ratio + !! so that subsequent uses of K already take into account the effect + !! of viscosity. Thus the original user-specified K array values are + !! lost unless they are backed up in k11input, for example. In a new + !! stress period/time step, the values in k11input are multiplied by + !! the viscosity ratio, not k11 since it contains viscosity-adjusted + !! hydraulic conductivity values. + !! + !< + subroutine store_original_k_arrays(this, ncells, njas) + ! -- modules + use MemoryManagerModule, only: mem_allocate + ! -- dummy + class(GwfNpftype) :: this + integer(I4B), intent(in) :: ncells + integer(I4B), intent(in) :: njas + ! -- local + integer(I4B) :: n +! ------------------------------------------------------------------------------ + ! + ! -- Retain copy of user-specified K arrays + do n = 1, ncells + this%k11input(n) = this%k11(n) + this%k22input(n) = this%k22(n) + this%k33input(n) = this%k33(n) + end do + ! + ! -- Return + return + end subroutine store_original_k_arrays + subroutine allocate_arrays(this, ncells, njas) ! ****************************************************************************** ! allocate_arrays -- Allocate npf arrays @@ -1267,6 +1352,11 @@ subroutine allocate_arrays(this, ncells, njas) call mem_allocate(this%ihcedge, 0, 'IHCEDGE', this%memoryPath) call mem_allocate(this%propsedge, 0, 0, 'PROPSEDGE', this%memoryPath) ! + ! -- Optional arrays only needed when vsc package is active + call mem_allocate(this%k11input, 0, 'K11INPUT', this%memoryPath) + call mem_allocate(this%k22input, 0, 'K22INPUT', this%memoryPath) + call mem_allocate(this%k33input, 0, 'K33INPUT', this%memoryPath) + ! ! -- Specific discharge is (re-)allocated when nedges is known call mem_allocate(this%spdis, 3, 0, 'SPDIS', this%memoryPath) ! diff --git a/src/Model/GroundWaterFlow/gwf3sfr8.f90 b/src/Model/GroundWaterFlow/gwf3sfr8.f90 index 8f9a6ae28b9..649445fa55e 100644 --- a/src/Model/GroundWaterFlow/gwf3sfr8.f90 +++ b/src/Model/GroundWaterFlow/gwf3sfr8.f90 @@ -146,6 +146,9 @@ module SfrModule integer(I4B), pointer :: idense !< flag indicating if density corrections are active real(DP), dimension(:, :), pointer, contiguous :: denseterms => null() !< density terms ! + ! -- viscosity variables + real(DP), dimension(:, :), pointer, contiguous :: viscratios => null() !< viscosity ratios (1: sfr vsc ratio; 2: gwf vsc ratio) + ! ! -- type bound procedures contains procedure :: sfr_allocate_scalars @@ -208,6 +211,8 @@ module SfrModule ! -- density procedure :: sfr_activate_density procedure, private :: sfr_calculate_density_exchange + ! -- viscosity + procedure :: sfr_activate_viscosity end type SfrType contains @@ -316,6 +321,7 @@ subroutine sfr_allocate_scalars(this) this%icheck = 1 this%iconvchk = 1 this%idense = 0 + this%ivsc = 0 this%ianynone = 0 this%ncrossptstot = 0 ! @@ -499,12 +505,15 @@ subroutine sfr_allocate_arrays(this) this%qauxcbc(i) = DZERO end do ! - !-- fill cauxcbc + ! -- fill cauxcbc this%cauxcbc(1) = 'FLOW-AREA ' ! ! -- allocate denseterms to size 0 call mem_allocate(this%denseterms, 3, 0, 'DENSETERMS', this%memoryPath) ! + ! -- allocate viscratios to size 0 + call mem_allocate(this%viscratios, 2, 0, 'VISCRATIOS', this%memoryPath) + ! ! -- return return end subroutine sfr_allocate_arrays @@ -2471,7 +2480,7 @@ subroutine sfr_ot_dv(this, idvsave, idvprint) call this%stagetab%add_term(stage) call this%stagetab%add_term(depth) call this%stagetab%add_term(w) - call this%sfr_calc_cond(n, depth, cond) + call this%sfr_calc_cond(n, depth, cond, stage, hgwf) if (node > 0) then sbot = this%strtop(n) - this%bthick(n) if (hgwf < sbot) then @@ -2557,6 +2566,7 @@ subroutine sfr_da(this) call mem_deallocate(this%stage0) call mem_deallocate(this%usflow0) call mem_deallocate(this%denseterms) + call mem_deallocate(this%viscratios) ! ! -- deallocate reach order and connection data call mem_deallocate(this%isfrorder) @@ -3405,7 +3415,7 @@ subroutine sfr_solve(this, n, h, hcof, rhs, update) ! ! -- calculate reach conductance for a unit depth of water ! if equal to zero will skip iterations - call this%sfr_calc_cond(n, d1, cstr) + call this%sfr_calc_cond(n, d1, cstr, hsfr, hgwf) ! ! -- set flag to skip iterations isolve = 1 @@ -3968,10 +3978,7 @@ subroutine sfr_calc_qgwf(this, n, depth, hgwf, qgwf, gwfhcof, gwfrhs) ! -- calculate saturation call sChSmooth(depth, sat, derv) ! - ! -- calculate conductance - call this%sfr_calc_cond(n, depth, cond) - ! - ! -- calculate groundwater leakage + ! -- terms for calculating direction of gradient across streambed tp = this%strtop(n) bt = tp - this%bthick(n) hsfr = tp + depth @@ -3979,6 +3986,11 @@ subroutine sfr_calc_qgwf(this, n, depth, hgwf, qgwf, gwfhcof, gwfrhs) if (htmp < bt) then htmp = bt end if + ! + ! -- calculate conductance + call this%sfr_calc_cond(n, depth, cond, hsfr, htmp) + ! + ! -- calculate groundwater leakage qgwf = sat * cond * (htmp - hsfr) gwfrhs0 = -sat * cond * hsfr gwfhcof0 = -sat * cond @@ -4002,25 +4014,41 @@ end subroutine sfr_calc_qgwf !! Method to calculate the reach-aquifer conductance for a SFR package reach. !! !< - subroutine sfr_calc_cond(this, n, depth, cond) + subroutine sfr_calc_cond(this, n, depth, cond, hsfr, htmp) ! -- dummy variables class(SfrType) :: this !< SfrType object integer(I4B), intent(in) :: n !< reach number real(DP), intent(in) :: depth !< reach depth real(DP), intent(inout) :: cond !< reach-aquifer conductance + real(DP), intent(in), optional :: hsfr !< stream stage + real(DP), intent(in), optional :: htmp !< head in gw cell ! -- local variables integer(I4B) :: node real(DP) :: wp + real(DP) :: vscratio ! ! -- initialize conductance cond = DZERO ! + ! -- initial viscosity ratio to 1 + vscratio = DONE + ! ! -- calculate conductance if GWF cell is active node = this%igwfnode(n) if (node > 0) then if (this%ibound(node) > 0) then + ! + ! -- direction of gradient across streambed determines which vsc ratio + if (this%ivsc == 1) then + if (hsfr > htmp) then + ! strm stg > gw head + vscratio = this%viscratios(1, n) + else + vscratio = this%viscratios(2, n) + end if + end if wp = this%calc_perimeter_wet(n, depth) - cond = this%hk(n) * this%length(n) * wp / this%bthick(n) + cond = this%hk(n) * vscratio * this%length(n) * wp / this%bthick(n) end if end if ! @@ -5563,6 +5591,37 @@ subroutine sfr_activate_density(this) return end subroutine sfr_activate_density + !> @brief Activate viscosity terms + !! + !! Method to activate addition of viscosity terms for exhange + !! with groundwater along a SFR package reach. + !! + !< + subroutine sfr_activate_viscosity(this) + ! -- modules + use MemoryManagerModule, only: mem_reallocate + ! -- dummy variables + class(SfrType), intent(inout) :: this !< SfrType object + ! -- local variables + integer(I4B) :: i + integer(I4B) :: j + ! + ! -- Set ivsc and reallocate viscratios to be of size MAXBOUND + this%ivsc = 1 + call mem_reallocate(this%viscratios, 2, this%MAXBOUND, 'VISCRATIOS', & + this%memoryPath) + do i = 1, this%maxbound + do j = 1, 2 + this%viscratios(j, i) = DONE + end do + end do + write (this%iout, '(/1x,a)') 'VISCOSITY HAS BEEN ACTIVATED FOR SFR & + &PACKAGE: '//trim(adjustl(this%packName)) + ! + ! -- return + return + end subroutine sfr_activate_viscosity + !> @brief Calculate density terms !! !! Method to galculate groundwater-reach density exchange terms for a diff --git a/src/Model/GroundWaterFlow/gwf3vsc8.f90 b/src/Model/GroundWaterFlow/gwf3vsc8.f90 new file mode 100644 index 00000000000..e5f3a6fbf20 --- /dev/null +++ b/src/Model/GroundWaterFlow/gwf3vsc8.f90 @@ -0,0 +1,1546 @@ +! Viscosity Package for representing variable-viscosity groundwater flow + +module GwfVscModule + + use KindModule, only: DP, I4B + use SimModule, only: store_error, store_warning, count_errors + use MemoryManagerModule, only: mem_allocate, mem_reallocate, & + mem_deallocate, mem_setptr + use MemoryHelperModule, only: create_mem_path + use ConstantsModule, only: DHALF, DZERO, DONE, LENMODELNAME, LENAUXNAME, & + DHNOFLO, MAXCHARLEN, LINELENGTH, LENMEMPATH + use TdisModule, only: kper, kstp + use NumericalPackageModule, only: NumericalPackageType + use BaseDisModule, only: DisBaseType + use GwfVscInputDataModule, only: GwfVscInputDataType + use ListsModule, only: basemodellist + + implicit none + + private + public :: GwfVscType + public :: vsc_cr + + type :: ConcentrationPointer + real(DP), dimension(:), pointer :: conc => null() !< pointer to concentration array + integer(I4B), dimension(:), pointer :: icbund => null() !< store pointer to gwt ibound array + end type ConcentrationPointer + + type, extends(NumericalPackageType) :: GwfVscType + integer(I4B), pointer :: thermivisc => null() !< viscosity formulation flag (1:Linear, 2:Nonlinear) + integer(I4B), pointer :: idxtmpr => null() !< if greater than 0 then an index for identifying whether the "species" array is temperature + integer(I4B), pointer :: ioutvisc => null() !< unit number for saving viscosity + integer(I4B), pointer :: iconcset => null() !< if 1 then conc points to a gwt (or gwe) model%x array + integer(I4B), pointer :: ireadelev => null() !< if 1 then elev has been allocated and filled + integer(I4B), dimension(:), pointer, contiguous :: ivisc => null() !< viscosity formulation flag for each species (1:Linear, 2:Nonlinear) + real(DP), pointer :: viscref => null() !< reference fluid viscosity + real(DP), dimension(:), pointer, contiguous :: visc => null() !< viscosity + real(DP), dimension(:), pointer, contiguous :: elev => null() !< cell center elevation (optional; if not specified, then use (top+bot)/2) + integer(I4B), dimension(:), pointer :: ibound => null() !< store pointer to ibound + + integer(I4B), pointer :: nviscspecies => null() !< number of concentration species used in viscosity equation + real(DP), dimension(:), pointer, contiguous :: dviscdc => null() !< linear change in viscosity with change in concentration + real(DP), dimension(:), pointer, contiguous :: cviscref => null() !< reference concentration used in viscosity equation + real(DP), dimension(:), pointer, contiguous :: ctemp => null() !< temporary array of size (nviscspec) to pass to calc_visc_x + character(len=LENMODELNAME), dimension(:), allocatable :: cmodelname !< names of gwt (or gwe) models used in viscosity equation + character(len=LENAUXNAME), dimension(:), allocatable :: cauxspeciesname !< names of aux columns used in viscosity equation + character(len=LENAUXNAME) :: name_temp_spec = 'TEMPERATURE' + ! + ! -- Viscosity constants + real(DP), pointer :: a2 => null() !< an empirical parameter specified by the user for calculating viscosity + real(DP), pointer :: a3 => null() !< an empirical parameter specified by the user for calculating viscosity + real(DP), pointer :: a4 => null() !< an empirical parameter specified by the user for calculating viscosity + + type(ConcentrationPointer), allocatable, dimension(:) :: modelconc !< concentration (or temperature) pointer for each solute (or heat) transport model + + real(DP), dimension(:), pointer, contiguous :: k11 => null() !< NPF hydraulic conductivity; if anisotropic, then this is Kx prior to rotation + real(DP), dimension(:), pointer, contiguous :: k22 => null() !< NPF hydraulic conductivity; if specified then this is Ky prior to rotation + real(DP), dimension(:), pointer, contiguous :: k33 => null() !< NPF hydraulic conductivity; if specified then this is Kz prior to rotation + real(DP), dimension(:), pointer, contiguous :: k11input => null() !< NPF hydraulic conductivity as originally specified by the user + real(DP), dimension(:), pointer, contiguous :: k22input => null() !< NPF hydraulic conductivity as originally specified by the user + real(DP), dimension(:), pointer, contiguous :: k33input => null() !< NPF hydraulic conductivity as originally specified by the user + integer(I4B), pointer :: kchangeper => null() ! last stress period in which any node K (or K22, or K33) values were changed (0 if unchanged from start of simulation) + integer(I4B), pointer :: kchangestp => null() ! last time step in which any node K (or K22, or K33) values were changed (0 if unchanged from start of simulation) + integer(I4B), dimension(:), pointer, contiguous :: nodekchange => null() ! grid array of flags indicating for each node whether its K (or K22, or K33) value changed (1) at (kchangeper, kchangestp) or not (0) + + contains + procedure :: vsc_df + procedure :: vsc_ar + procedure, public :: vsc_ar_bnd + procedure :: vsc_rp + procedure :: vsc_ad + procedure, public :: vsc_ad_bnd + procedure :: vsc_ot_dv + procedure :: vsc_da + procedure, private :: vsc_calcvisc + procedure :: allocate_scalars + procedure, private :: allocate_arrays + procedure, private :: read_options + procedure, private :: set_options + procedure, private :: read_dimensions + procedure, private :: read_packagedata + procedure, private :: set_packagedata + procedure, private :: set_npf_pointers + procedure, public :: update_k_with_vsc + procedure, private :: vsc_set_changed_at + procedure, public :: calc_q_visc + procedure, public :: get_visc_ratio + procedure :: set_concentration_pointer + end type GwfVscType + +contains + + function calc_visc(ivisc, viscref, dviscdc, cviscref, conc, & + a2, a3, a4) result(visc) +! ****************************************************************************** +! calc_visc -- generic function to calculate changes in fluid viscosity +! using a linear formulation +! ****************************************************************************** +! +! SPECIFICATIONS: +! ------------------------------------------------------------------------------ + + ! -- dummy + integer(I4B), dimension(:), intent(in) :: ivisc + real(DP), intent(in) :: viscref + real(DP), dimension(:), intent(in) :: dviscdc + real(DP), dimension(:), intent(in) :: cviscref + real(DP), dimension(:), intent(in) :: conc + real(DP), intent(in) :: a2, a3, a4 + ! -- return + real(DP) :: visc + ! -- local + integer(I4B) :: nviscspec + integer(I4B) :: i + real(DP) :: mu_t + real(DP) :: expon +! ------------------------------------------------------------------------------ + ! + nviscspec = size(dviscdc) + visc = viscref + + do i = 1, nviscspec + if (ivisc(i) == 1) then + visc = visc + dviscdc(i) * (conc(i) - cviscref(i)) + else + expon = -1 * a3 * ((conc(i) - cviscref(i)) / & + ((conc(i) + a4) * (cviscref(i) + a4))) + mu_t = viscref * a2**expon + ! Order matters!! (This assumes we apply the temperature correction after + ! accounting for solute concentrations) + ! If a nonlinear correction is applied, then b/c it takes into account + ! viscref, need to subtract it in this case + ! At most, there will only ever be 1 nonlinear correction + visc = (visc - viscref) + mu_t + end if + ! end if + end do + ! + ! -- return + return + end function calc_visc + + !> @ brief Create a new package object + !! + !! Create a new VSC Package object. + !! + !< + subroutine vsc_cr(vscobj, name_model, inunit, iout) + ! -- dummy + type(GwfVscType), pointer :: vscobj + character(len=*), intent(in) :: name_model + integer(I4B), intent(in) :: inunit + integer(I4B), intent(in) :: iout +! ------------------------------------------------------------------------------ + ! + ! -- Create the object + allocate (vscobj) + ! + ! -- create name and memory path + call vscobj%set_names(1, name_model, 'VSC', 'VSC') + ! + ! -- Allocate scalars + call vscobj%allocate_scalars() + ! + ! -- Set variables + vscobj%inunit = inunit + vscobj%iout = iout + ! + ! -- Initialize block parser + call vscobj%parser%Initialize(vscobj%inunit, vscobj%iout) + ! + ! -- Return + return + end subroutine vsc_cr + + !> @ brief Define viscosity package options and dimensions + !! + !! Define viscosity package options and dimensions + !! + !< + subroutine vsc_df(this, dis, vsc_input) + ! -- modules + ! -- dummy + class(GwfVscType) :: this !< this viscosity package + class(DisBaseType), pointer, intent(in) :: dis !< pointer to discretization + type(GwfVscInputDataType), optional, intent(in) :: vsc_input !< optional vsc input data, otherwise read from file + ! -- local + ! -- formats + character(len=*), parameter :: fmtvsc = & + "(1x,/1x,'VSC -- Viscosity Package, version 1, 11/15/2022', & + &' input read from unit ', i0, //)" +! ------------------------------------------------------------------------------ + ! + ! --print a message identifying the viscosity package + write (this%iout, fmtvsc) this%inunit + ! + ! -- store pointers to arguments that were passed in + this%dis => dis + + if (.not. present(vsc_input)) then + ! + ! -- Read viscosity options + call this%read_options() + ! + ! -- Read viscosity dimensions + call this%read_dimensions() + else + ! set from input data instead + call this%set_options(vsc_input) + this%nviscspecies = vsc_input%nviscspecies + end if + ! + ! -- Allocate arrays + call this%allocate_arrays(dis%nodes) + + if (.not. present(vsc_input)) then + ! + ! -- Read viscosity packagedata + call this%read_packagedata() + else + ! set from input data instead + call this%set_packagedata(vsc_input) + end if + ! + ! -- Return + return + end subroutine vsc_df + + !> @ brief Allocate and read method for viscosity package + !! + !! Generic method to allocate and read static data for the viscosity + !! package available within the GWF model type. + !! + !< + subroutine vsc_ar(this, ibound) + ! -- modules + ! -- dummy + class(GwfVscType) :: this + integer(I4B), dimension(:), pointer :: ibound + ! -- local + ! -- formats +! ------------------------------------------------------------------------------ + ! + ! -- store pointers to arguments that were passed in + this%ibound => ibound + ! + ! -- Set pointers to npf variables + call this%set_npf_pointers() + ! + ! -- Return + return + end subroutine vsc_ar + + !> @brief Activate viscosity in advanced packages + !! + !! Viscosity ar_bnd rountine to activate viscosity in the advanced + !! packages. This routine is called from gwf_ar() as it moves through each + !! package + !! + !< + subroutine vsc_ar_bnd(this, packobj) + ! -- modules + use BndModule, only: BndType + use DrnModule, only: DrnType + use GhbModule, only: GhbType + use RivModule, only: RivType + use LakModule, only: LakType + use SfrModule, only: SfrType + use MawModule, only: MawType + ! -- dummy + class(GwfVscType) :: this + class(BndType), pointer :: packobj + ! -- local + ! ---------------------------------------------------------------------------- + ! + ! -- Add density terms based on boundary package type + select case (packobj%filtyp) + case ('DRN') + ! + ! -- activate viscosity for the drain package + select type (packobj) + type is (DrnType) + call packobj%bnd_activate_viscosity() + end select + case ('GHB') + ! + ! -- activate viscosity for the drain package + select type (packobj) + type is (GhbType) + call packobj%bnd_activate_viscosity() + end select + case ('RIV') + ! + ! -- activate viscosity for the drain package + select type (packobj) + type is (RivType) + call packobj%bnd_activate_viscosity() + end select + case ('LAK') + ! + ! -- activate viscosity for lake package + select type (packobj) + type is (LakType) + call packobj%lak_activate_viscosity() + end select + + case ('SFR') + ! + ! -- activate viscosity for sfr package + select type (packobj) + type is (SfrType) + call packobj%sfr_activate_viscosity() + end select + + case ('MAW') + ! + ! -- activate viscosity for maw package + select type (packobj) + type is (MawType) + call packobj%maw_activate_viscosity() + end select + + case default + ! + ! -- nothing + end select + ! + ! -- Return + return + end subroutine vsc_ar_bnd + + !> @brief Set pointers to NPF variables + !! + !! Set array and variable pointers from the NPF + !! package for access by VSC. + !! + !< + subroutine set_npf_pointers(this) + ! -- dummy variables + class(GwfVscType) :: this + ! -- local variables + character(len=LENMEMPATH) :: npfMemoryPath + ! + ! -- Set pointers to other package variables + ! -- NPF + npfMemoryPath = create_mem_path(this%name_model, 'NPF') + call mem_setptr(this%k11, 'K11', npfMemoryPath) + call mem_setptr(this%k22, 'K22', npfMemoryPath) + call mem_setptr(this%k33, 'K33', npfMemoryPath) + call mem_setptr(this%k11input, 'K11INPUT', npfMemoryPath) + call mem_setptr(this%k22input, 'K22INPUT', npfMemoryPath) + call mem_setptr(this%k33input, 'K33INPUT', npfMemoryPath) + call mem_setptr(this%kchangeper, 'KCHANGEPER', npfMemoryPath) + call mem_setptr(this%kchangestp, 'KCHANGESTP', npfMemoryPath) + call mem_setptr(this%nodekchange, 'NODEKCHANGE', npfMemoryPath) + ! + return + end subroutine set_npf_pointers + + !> @ brief Read new period data in viscosity package + !! + !! Method to read and prepare period data for the VSC package. + !! + !< + subroutine vsc_rp(this) + ! -- modules + use TdisModule, only: kstp, kper + ! -- dummy + class(GwfVscType) :: this + ! -- local + character(len=LINELENGTH) :: errmsg + integer(I4B) :: i + ! -- formats + character(len=*), parameter :: fmtc = & + "('Viscosity Package does not have a concentration set & + &for species ',i0,'. One or more model names may be specified & + &incorrectly in the PACKAGEDATA block or a GWF-GWT exchange may need & + &to be activated.')" +! ------------------------------------------------------------------------------ + ! + ! -- Check to make sure all concentration pointers have been set + if (kstp * kper == 1) then + do i = 1, this%nviscspecies + if (.not. associated(this%modelconc(i)%conc)) then + write (errmsg, fmtc) i + call store_error(errmsg) + end if + end do + if (count_errors() > 0) then + call this%parser%StoreErrorUnit() + end if + end if + ! + ! -- return + return + end subroutine vsc_rp + + !> @ brief Advance the viscosity package + !! + !! Advance data in the VSC package. The method sets or + !! advances time series, time array series, and observation + !! data. + !! + !< + subroutine vsc_ad(this) + ! -- dummy + class(GwfVscType) :: this + ! -- local +! ------------------------------------------------------------------------------ + ! + ! -- update viscosity using the latest concentration/temperature + call this%vsc_calcvisc() + ! + ! -- Return + return + end subroutine vsc_ad + + !> @brief Advance the boundary packages when viscosity is active + !! + !! Update the conductance values associate with inflow from a boundary + !! when VSC package is active. + !< + subroutine vsc_ad_bnd(this, packobj, hnew) + ! -- modules + use BndModule, only: BndType + ! -- dummy + class(GwfVscType) :: this + class(BndType), pointer :: packobj + real(DP), intent(in), dimension(:) :: hnew + ! -- local + integer(I4B) :: i, j + integer(I4B) :: n, locvisc, locelev + integer(I4B), dimension(:), allocatable :: locconc + ! + ! -- initialize + locvisc = 0 + locelev = 0 + allocate (locconc(this%nviscspecies)) + locconc(:) = 0 + ! + ! -- Add viscosity terms for conductance-dependent boundaries + do n = 1, packobj%naux + if (packobj%auxname(n) == 'VISCOSITY') then + locvisc = n + else if (packobj%auxname(n) == 'ELEVATION') then + locelev = n + end if + end do + ! + ! -- find aux columns for conc (or temp.) that affect viscosity + do i = 1, this%nviscspecies + locconc(i) = 0 + do j = 1, packobj%naux + if (this%cauxspeciesname(i) == packobj%auxname(j)) then + locconc(i) = j + exit + end if + end do + if (locconc(i) == 0) then + ! -- one not found, so don't use and mark all as 0 + locconc(:) = 0 + exit + end if + end do + ! + ! -- apply viscosity terms to inflow from boundary based on package type + select case (packobj%filtyp) + case ('GHB', 'DRN', 'RIV') + ! + ! -- general head, drain, and river boundary + call vsc_ad_standard_bnd(packobj, hnew, this%visc, this%viscref, & + locelev, locvisc, locconc, this%dviscdc, & + this%cviscref, this%ivisc, this%a2, this%a3, & + this%a4, this%ctemp) + case ('LAK') + ! + ! -- lake + ! Update 'viscratios' internal to lak such that they are + ! automatically applied in the LAK calc_cond() routine + call vsc_ad_lak(packobj, this%visc, this%viscref, this%elev, locvisc, & + locconc, this%dviscdc, this%cviscref, this%ivisc, & + this%a2, this%a3, this%a4, this%ctemp) + case ('SFR') + ! + ! -- streamflow routing + ! Update 'viscratios' internal to sfr such that they are + ! automatically applied in the SFR calc_cond() routine + call vsc_ad_sfr(packobj, this%visc, this%viscref, this%elev, locvisc, & + locconc, this%dviscdc, this%cviscref, this%ivisc, & + this%a2, this%a3, this%a4, this%ctemp) + case ('MAW') + ! + ! -- multi-aquifer well + call vsc_ad_maw(packobj, this%visc, this%viscref, this%elev, locvisc, & + locconc, this%dviscdc, this%cviscref, this%ivisc, & + this%a2, this%a3, this%a4, this%ctemp) + case ('UZF') + ! + ! -- unsaturated-zone flow + case default + ! + ! -- nothing + end select + ! + ! -- deallocate + deallocate (locconc) + ! + ! -- Return + return + end subroutine vsc_ad_bnd + + !> @brief advance ghb while accounting for viscosity + !! + !! When flow enters from ghb boundary type, take into account the effects + !! of viscosity on the user-specified conductance terms + !< + subroutine vsc_ad_standard_bnd(packobj, hnew, visc, viscref, locelev, & + locvisc, locconc, dviscdc, cviscref, & + ivisc, a2, a3, a4, ctemp) + ! -- modules + use BndModule, only: BndType + class(BndType), pointer :: packobj + ! -- dummy + real(DP), intent(in), dimension(:) :: hnew + real(DP), intent(in), dimension(:) :: visc + real(DP), intent(in) :: a2, a3, a4 + real(DP), intent(in) :: viscref + integer(I4B), intent(in) :: locelev + integer(I4B), intent(in) :: locvisc + integer(I4B), dimension(:), intent(in) :: locconc + integer(I4B), dimension(:), intent(in) :: ivisc + real(DP), dimension(:), intent(in) :: dviscdc + real(DP), dimension(:), intent(in) :: cviscref + real(DP), dimension(:), intent(inout) :: ctemp + ! -- local + integer(I4B) :: n + integer(I4B) :: node + real(DP) :: viscbnd +! ------------------------------------------------------------------------------- + ! + ! -- Process density terms for each GHB + do n = 1, packobj%nbound + node = packobj%nodelist(n) + ! + ! -- Check if boundary cell is active, cycle if not + if (packobj%ibound(node) <= 0) cycle + ! + ! -- calculate the viscosity associcated with the boundary + viscbnd = calc_bnd_viscosity(n, locvisc, locconc, viscref, dviscdc, & + cviscref, ctemp, ivisc, a2, a3, a4, & + packobj%auxvar) + ! + ! -- update boundary conductance based on viscosity effects + packobj%bound(2, n) = update_bnd_cond(viscbnd, viscref, & + packobj%condinput(n)) + ! + end do + ! + ! -- Return + return + end subroutine vsc_ad_standard_bnd + + !> @brief Update sfr-related viscosity ratios + !! + !! When the viscosity package is active, update the viscosity ratio that is + !! applied to the hydraulic conductivity specified in the SFR package + !< + subroutine vsc_ad_sfr(packobj, visc, viscref, elev, locvisc, locconc, & + dviscdc, cviscref, ivisc, a2, a3, a4, ctemp) + ! -- modules + use BndModule, only: BndType + use SfrModule, only: SfrType + class(BndType), pointer :: packobj + ! -- dummy + real(DP), intent(in) :: viscref + real(DP), intent(in) :: a2, a3, a4 + integer(I4B), intent(in) :: locvisc + integer(I4B), dimension(:), intent(in) :: locconc + integer(I4B), dimension(:), intent(in) :: ivisc + real(DP), dimension(:), intent(in) :: visc + real(DP), dimension(:), intent(in) :: elev + real(DP), dimension(:), intent(in) :: dviscdc + real(DP), dimension(:), intent(in) :: cviscref + real(DP), dimension(:), intent(inout) :: ctemp + ! -- local + integer(I4B) :: n + integer(I4B) :: node + real(DP) :: viscsfr +! ------------------------------------------------------------------------------- + ! + ! -- update viscosity ratios for updating hyd. cond (and conductance) + select type (packobj) + type is (SfrType) + do n = 1, packobj%nbound + ! + ! -- get gwf node number + node = packobj%nodelist(n) + ! + ! -- Check if boundary cell is active, cycle if not + if (packobj%ibound(node) <= 0) cycle + ! + ! -- + ! + ! -- calculate the viscosity associcated with the boundary + viscsfr = calc_bnd_viscosity(n, locvisc, locconc, viscref, dviscdc, & + cviscref, ctemp, ivisc, a2, a3, a4, & + packobj%auxvar) + ! + ! -- fill sfr relative viscosity into column 1 of viscratios + packobj%viscratios(1, n) = calc_vsc_ratio(viscref, viscsfr) + ! + ! -- fill gwf relative viscosity into column 2 of viscratios + packobj%viscratios(2, n) = calc_vsc_ratio(viscref, visc(node)) + end do + end select + ! + ! -- Return + return + end subroutine vsc_ad_sfr + + !> @brief Update lak-related viscosity ratios + !! + !! When the viscosity package is active, update the viscosity ratio that is + !! applied to the lakebed conductance calculated in the LAK package + !< + subroutine vsc_ad_lak(packobj, visc, viscref, elev, locvisc, locconc, & + dviscdc, cviscref, ivisc, a2, a3, a4, ctemp) + ! -- modules + use BndModule, only: BndType + use LakModule, only: LakType + class(BndType), pointer :: packobj + ! -- dummy + real(DP), intent(in) :: viscref + real(DP), intent(in) :: a2, a3, a4 + integer(I4B), intent(in) :: locvisc + integer(I4B), dimension(:), intent(in) :: locconc + integer(I4B), dimension(:), intent(in) :: ivisc + real(DP), dimension(:), intent(in) :: visc + real(DP), dimension(:), intent(in) :: elev + real(DP), dimension(:), intent(in) :: dviscdc + real(DP), dimension(:), intent(in) :: cviscref + real(DP), dimension(:), intent(inout) :: ctemp + ! -- local + integer(I4B) :: n + integer(I4B) :: node + real(DP) :: visclak +! ------------------------------------------------------------------------------- + ! + ! -- update viscosity ratios for updating hyd. cond (and conductance) + select type (packobj) + type is (LakType) + do n = 1, packobj%nbound + ! + ! -- get gwf node number + node = packobj%nodelist(n) + ! + ! -- Check if boundary cell is active, cycle if not + if (packobj%ibound(node) <= 0) cycle + ! + ! -- + ! + ! -- calculate the viscosity associcated with the boundary + visclak = calc_bnd_viscosity(n, locvisc, locconc, viscref, dviscdc, & + cviscref, ctemp, ivisc, a2, a3, a4, & + packobj%auxvar) + ! + ! -- fill lak relative viscosity into column 1 of viscratios + packobj%viscratios(1, n) = calc_vsc_ratio(viscref, visclak) + ! + ! -- fill gwf relative viscosity into column 2 of viscratios + packobj%viscratios(2, n) = calc_vsc_ratio(viscref, visc(node)) + end do + end select + ! + ! -- Return + return + end subroutine vsc_ad_lak + + !> @brief Update maw-related viscosity ratios + !! + !! When the viscosity package is active, update the viscosity ratio that is + !! applied to the conductance calculated in the MAW package + !< + subroutine vsc_ad_maw(packobj, visc, viscref, elev, locvisc, locconc, & + dviscdc, cviscref, ivisc, a2, a3, a4, ctemp) + ! -- modules + use BndModule, only: BndType + use MawModule, only: MawType + class(BndType), pointer :: packobj + ! -- dummy + real(DP), intent(in) :: viscref + real(DP), intent(in) :: a2, a3, a4 + integer(I4B), intent(in) :: locvisc + integer(I4B), dimension(:), intent(in) :: locconc + integer(I4B), dimension(:), intent(in) :: ivisc + real(DP), dimension(:), intent(in) :: visc + real(DP), dimension(:), intent(in) :: elev + real(DP), dimension(:), intent(in) :: dviscdc + real(DP), dimension(:), intent(in) :: cviscref + real(DP), dimension(:), intent(inout) :: ctemp + ! -- local + integer(I4B) :: n + integer(I4B) :: node + real(DP) :: viscmaw +! ------------------------------------------------------------------------------- + ! + ! -- update viscosity ratios for updating hyd. cond (and conductance) + select type (packobj) + type is (MawType) + do n = 1, packobj%nbound + ! + ! -- get gwf node number + node = packobj%nodelist(n) + ! + ! -- Check if boundary cell is active, cycle if not + if (packobj%ibound(node) <= 0) cycle + ! + ! -- + ! + ! -- calculate the viscosity associcated with the boundary + viscmaw = calc_bnd_viscosity(n, locvisc, locconc, viscref, dviscdc, & + cviscref, ctemp, ivisc, a2, a3, a4, & + packobj%auxvar) + ! + ! -- fill lak relative viscosity into column 1 of viscratios + packobj%viscratios(1, n) = calc_vsc_ratio(viscref, viscmaw) + ! + ! -- fill gwf relative viscosity into column 2 of viscratios + packobj%viscratios(2, n) = calc_vsc_ratio(viscref, visc(node)) + end do + end select + ! + ! -- Return + return + end subroutine vsc_ad_maw + + !> @brief Apply viscosity to the conductance term + !! + !! When the viscosity package is active apply the viscosity ratio to the + !! active boundary package's conductance term. + !< + function update_bnd_cond(bndvisc, viscref, spcfdcond) result(updatedcond) + ! -- modules + ! -- dummy + real(DP), intent(in) :: viscref + real(DP), intent(in) :: bndvisc + real(DP), intent(in) :: spcfdcond + ! -- local + real(DP) :: vscratio + real(DP) :: updatedcond +! ------------------------------------------------------------------------------- + ! + vscratio = calc_vsc_ratio(viscref, bndvisc) + ! + ! -- calculate new conductance here + updatedcond = vscratio * spcfdcond + ! + ! -- Return + return + end function update_bnd_cond + + !> @brief calculate and return the viscosity ratio + !< + function calc_vsc_ratio(viscref, bndvisc) result(viscratio) + ! -- dummy + real(DP), intent(in) :: viscref + real(DP), intent(in) :: bndvisc + ! -- local + real(DP) :: viscratio +! ------------------------------------------------------------------------------- + ! + viscratio = viscref / bndvisc + ! + ! -- Return + return + end function calc_vsc_ratio + + !> @ brief Calculate the boundary viscosity + !! + !! Return the viscosity of the boundary package using one of + !! the options in the following order of priority: + !! 1. Assign as aux variable in column with name 'VISCOSITY' + !! 2. Calculate using viscosity equation and nviscspecies aux columns + !! 3. If neither of those, then assign as viscref !! + !< + function calc_bnd_viscosity(n, locvisc, locconc, viscref, dviscdc, cviscref, & + ctemp, ivisc, a2, a3, a4, auxvar) result(viscbnd) + ! -- modules + ! -- dummy + integer(I4B), intent(in) :: n + integer(I4B), intent(in) :: locvisc + real(DP), intent(in) :: a2, a3, a4 + integer(I4B), dimension(:), intent(in) :: ivisc + integer(I4B), dimension(:), intent(in) :: locconc + real(DP), intent(in) :: viscref + real(DP), dimension(:), intent(in) :: dviscdc + real(DP), dimension(:), intent(in) :: cviscref + real(DP), dimension(:), intent(inout) :: ctemp + real(DP), dimension(:, :), intent(in) :: auxvar + ! -- return + real(DP) :: viscbnd + ! -- local + integer(I4B) :: i +! ------------------------------------------------------------------------------ + ! + ! -- assign boundary viscosity based on one of three options + if (locvisc > 0) then + ! -- assign viscosity to an aux column named 'VISCOSITY' + viscbnd = auxvar(locvisc, n) + else if (locconc(1) > 0) then + ! -- calculate viscosity using one or more concentration auxcolumns + do i = 1, size(locconc) + ctemp(i) = DZERO + if (locconc(i) > 0) then + ctemp(i) = auxvar(locconc(i), n) + end if + end do + viscbnd = calc_visc(ivisc, viscref, dviscdc, cviscref, ctemp, a2, a3, a4) + else + ! -- neither of the above, so assign as viscref + viscbnd = viscref + end if + ! + ! -- return + return + end function calc_bnd_viscosity + + !> @brief Calculate the viscosity ratio + !! + !! Calculate the viscosity ratio applied to the hydraulic characteristic + !! provided by the user. The viscosity ratio is applicable only + !! when the hydraulic characteristic is specified as positive and will not + !! be applied when the hydchr is negative + !< + subroutine get_visc_ratio(this, n, m, gwhdn, gwhdm, viscratio) + ! -- modules + use ConstantsModule, only: DONE + ! -- dummy + class(GwfVscType) :: this + integer(I4B), intent(in) :: n, m + real(DP), intent(in) :: gwhdn, gwhdm + real(DP), intent(inout) :: viscratio + ! -- loca + integer(I4B) :: cellid +! ------------------------------------------------------------------------------ +! + viscratio = DONE + if (gwhdm > gwhdn) then + cellid = m + else if (gwhdn >= gwhdm) then + cellid = n + end if + call this%calc_q_visc(cellid, viscratio) + ! + ! -- return + return + end subroutine get_visc_ratio + + !> @brief Account for viscosity in the aquiferhorizontal flow barriers + !! + !! Will return the viscosity associated with the upgradient node (cell) + !! to the HFB package for adjusting the hydraulic characteristic (hydchr) + !! of the barrier + !< + subroutine calc_q_visc(this, cellid, viscratio) + ! -- dummy variables + class(GwfVscType) :: this + integer(I4B), intent(in) :: cellid + ! -- return + real(DP), intent(inout) :: viscratio + ! -- local + real(DP) :: visc +! ------------------------------------------------------------------------------ + ! + ! -- Retrieve viscosity for the passed node number + visc = this%visc(cellid) + ! + ! -- Calculate the viscosity ratio for the + viscratio = calc_vsc_ratio(this%viscref, visc) + ! + ! -- return + return + end subroutine calc_q_visc + + !> @brief Appled the viscosity ratio (mu_o/mu) to the hydraulic conductivity + !! + !! This routine called after updating the viscosity values using the latest + !! concentration and/or temperature values. The ratio mu_o/mu, reference + !! viscosity divided by the updated viscosity value, is multiplied by K + !! for each cell. + !< + subroutine update_k_with_vsc(this) + ! -- modules + ! -- dummy + class(GwfVscType) :: this + ! -- local + integer(I4B) :: n + real(DP) :: viscratio +! ------------------------------------------------------------------------------ + ! + ! -- For viscosity-based K's, apply change of K to K11 by starting with + ! user-specified K values and not the K's leftover from the last viscosity + ! update. + do n = 1, this%dis%nodes + call this%calc_q_visc(n, viscratio) + this%k11(n) = this%k11input(n) * viscratio + this%k22(n) = this%k22input(n) * viscratio + this%k33(n) = this%k33input(n) * viscratio + this%nodekchange(n) = 1 + end do + ! + ! -- Flag kchange + call this%vsc_set_changed_at(kper, kstp) + ! + ! -- return + return + end subroutine update_k_with_vsc + + !> @brief Mark K changes as having occurred at (kper, kstp) + !! + !! Procedure called by VSC code when K updated due to viscosity changes. + !! K values changed at (kper, kstp). + !! + !< + subroutine vsc_set_changed_at(this, kper, kstp) + ! -- dummy variables + class(GwfVscType) :: this + integer(I4B), intent(in) :: kper + integer(I4B), intent(in) :: kstp + ! + this%kchangeper = kper + this%kchangestp = kstp + ! + return + end subroutine vsc_set_changed_at + + !> @ brief Output viscosity package dependent-variable terms. + !! + !! Save calculated viscosity array to binary file + !! + !< + subroutine vsc_ot_dv(this, idvfl) + ! -- dummy + class(GwfVscType) :: this + integer(I4B), intent(in) :: idvfl + ! -- local + character(len=1) :: cdatafmp = ' ', editdesc = ' ' + integer(I4B) :: ibinun + integer(I4B) :: iprint + integer(I4B) :: nvaluesp + integer(I4B) :: nwidthp + real(DP) :: dinact +! ------------------------------------------------------------------------------ + ! + ! -- Set unit number for viscosity output + if (this%ioutvisc /= 0) then + ibinun = 1 + else + ibinun = 0 + end if + if (idvfl == 0) ibinun = 0 + ! + ! -- save viscosity array + if (ibinun /= 0) then + iprint = 0 + dinact = DHNOFLO + ! + ! -- write viscosity to binary file + if (this%ioutvisc /= 0) then + ibinun = this%ioutvisc + call this%dis%record_array(this%visc, this%iout, iprint, ibinun, & + ' VISCOSITY', cdatafmp, nvaluesp, & + nwidthp, editdesc, dinact) + end if + end if + ! + ! -- Return + return + end subroutine vsc_ot_dv + + !> @ brief Deallocate viscosity package memory + !! + !! Deallocate viscosity package scalars and arrays. + !! + !< + subroutine vsc_da(this) + ! -- modules + ! -- dummy + class(GwfVscType) :: this +! ------------------------------------------------------------------------------ + ! + ! -- Deallocate arrays if package was active + if (this%inunit > 0) then + call mem_deallocate(this%visc) + call mem_deallocate(this%ivisc) + call mem_deallocate(this%dviscdc) + call mem_deallocate(this%cviscref) + call mem_deallocate(this%ctemp) + deallocate (this%cmodelname) + deallocate (this%cauxspeciesname) + deallocate (this%modelconc) + end if + ! + ! -- Scalars + call mem_deallocate(this%thermivisc) + call mem_deallocate(this%idxtmpr) + call mem_deallocate(this%ioutvisc) + call mem_deallocate(this%ireadelev) + call mem_deallocate(this%iconcset) + call mem_deallocate(this%viscref) + call mem_deallocate(this%nviscspecies) + call mem_deallocate(this%a2) + call mem_deallocate(this%a3) + call mem_deallocate(this%a4) + ! + ! -- Nullify pointers to other package variables + nullify (this%k11) + nullify (this%k22) + nullify (this%k33) + nullify (this%k11input) + nullify (this%k22input) + nullify (this%k33input) + nullify (this%kchangeper) + nullify (this%kchangestp) + nullify (this%nodekchange) + ! + ! -- deallocate parent + call this%NumericalPackageType%da() + ! + ! -- Return + return + end subroutine vsc_da + + !> @ brief Read dimensions + !! + !! Read dimensions for the viscosity package + !! + !< + subroutine read_dimensions(this) + ! -- modules + ! -- dummy + class(GwfVscType), intent(inout) :: this + ! -- local + character(len=LINELENGTH) :: errmsg, keyword + integer(I4B) :: ierr + logical :: isfound, endOfBlock + ! -- format +! ------------------------------------------------------------------------------ + ! + ! -- get dimensions block + call this%parser%GetBlock('DIMENSIONS', isfound, ierr, & + supportOpenClose=.true.) + ! + ! -- parse dimensions block if detected + if (isfound) then + write (this%iout, '(/1x,a)') 'Processing VSC DIMENSIONS block' + do + call this%parser%GetNextLine(endOfBlock) + if (endOfBlock) exit + call this%parser%GetStringCaps(keyword) + select case (keyword) + case ('NVISCSPECIES') + this%nviscspecies = this%parser%GetInteger() + write (this%iout, '(4x,a,i0)') 'NVISCSPECIES = ', this%nviscspecies + case default + write (errmsg, '(4x,a,a)') & + 'unknown VSC dimension: ', trim(keyword) + call store_error(errmsg) + call this%parser%StoreErrorUnit() + end select + end do + write (this%iout, '(1x,a)') 'End of VSC DIMENSIONS block' + else + call store_error('Required VSC DIMENSIONS block not found.') + call this%parser%StoreErrorUnit() + end if + ! + ! -- check dimension + if (this%nviscspecies < 1) then + call store_error('NVISCSPECIES must be greater than zero.') + call this%parser%StoreErrorUnit() + end if + ! + ! -- return + return + end subroutine read_dimensions + + !> @ brief Read data for package + !! + !! Method to read data for the viscosity package. + !! + !< + subroutine read_packagedata(this) + ! -- modules + ! -- dummy + class(GwfVscType) :: this + ! -- local + character(len=LINELENGTH) :: errmsg + character(len=LINELENGTH) :: line + integer(I4B) :: ierr + integer(I4B) :: iviscspec + logical :: isfound, endOfBlock + logical :: blockrequired + integer(I4B), dimension(:), allocatable :: itemp + character(len=10) :: c10 + character(len=16) :: c16 + ! -- format + character(len=*), parameter :: fmterr = & + "('Invalid value for IRHOSPEC (',i0,') detected in VSC Package. & + &IRHOSPEC must be > 0 and <= NVISCSPECIES, and duplicate values & + &are not allowed.')" +! ------------------------------------------------------------------------------ + ! + ! -- initialize + allocate (itemp(this%nviscspecies)) + itemp(:) = 0 + ! + ! -- get packagedata block + blockrequired = .true. + call this%parser%GetBlock('PACKAGEDATA', isfound, ierr, & + blockRequired=blockRequired, & + supportOpenClose=.true.) + ! + ! -- parse packagedata block + if (isfound) then + write (this%iout, '(1x,a)') 'Procesing VSC PACKAGEDATA block' + do + call this%parser%GetNextLine(endOfBlock) + if (endOfBlock) exit + iviscspec = this%parser%GetInteger() + if (iviscspec < 1 .or. iviscspec > this%nviscspecies) then + write (errmsg, fmterr) iviscspec + call store_error(errmsg) + end if + if (itemp(iviscspec) /= 0) then + write (errmsg, fmterr) iviscspec + call store_error(errmsg) + end if + itemp(iviscspec) = 1 + ! + this%dviscdc(iviscspec) = this%parser%GetDouble() + this%cviscref(iviscspec) = this%parser%GetDouble() + call this%parser%GetStringCaps(this%cmodelname(iviscspec)) + call this%parser%GetStringCaps(this%cauxspeciesname(iviscspec)) + ! + if (this%cauxspeciesname(iviscspec) == this%name_temp_spec) then + if (this%idxtmpr > 0) then + write (errmsg, '(a)') 'More than one species in VSC input identified & + &as '//trim(this%name_temp_spec)//'. Only one species may be & + &designated to represent temperature.' + call store_error(errmsg) + else + this%idxtmpr = iviscspec + if (this%thermivisc == 2) then + this%ivisc(iviscspec) = 2 + end if + end if + end if + end do + else + call store_error('Required VSC PACKAGEDATA block not found.') + call this%parser%StoreErrorUnit() + end if + ! + ! -- Check for errors. + if (count_errors() > 0) then + call this%parser%StoreErrorUnit() + end if + ! + ! -- write packagedata information + write (this%iout, '(/,1x,a)') 'Summary of species information in VSC Package' + write (this%iout, '(1a11,5a17)') & + 'Species', 'DVISCDC', 'CVISCREF', 'Model', 'AUXSPECIESNAME' + do iviscspec = 1, this%nviscspecies + write (c10, '(i0)') iviscspec + line = ' '//adjustr(c10) + + write (c16, '(g15.6)') this%dviscdc(iviscspec) + line = trim(line)//' '//adjustr(c16) + write (c16, '(g15.6)') this%cviscref(iviscspec) + line = trim(line)//' '//adjustr(c16) + write (c16, '(a)') this%cmodelname(iviscspec) + line = trim(line)//' '//adjustr(c16) + write (c16, '(a)') this%cauxspeciesname(iviscspec) + line = trim(line)//' '//adjustr(c16) + write (this%iout, '(a)') trim(line) + end do + ! + ! -- deallocate + deallocate (itemp) + ! + write (this%iout, '(/,1x,a)') 'End of VSC PACKAGEDATA block' + ! + ! -- return + return + end subroutine read_packagedata + + !> @brief Sets package data instead of reading from file + !< + subroutine set_packagedata(this, input_data) + class(GwfVscType) :: this !< this vscoancy pkg + type(GwfVscInputDataType), intent(in) :: input_data !< the input data to be set + ! local + integer(I4B) :: ispec + + do ispec = 1, this%nviscspecies + this%dviscdc(ispec) = input_data%dviscdc(ispec) + this%cviscref(ispec) = input_data%cviscref(ispec) + this%cmodelname(ispec) = input_data%cmodelname(ispec) + this%cauxspeciesname(ispec) = input_data%cauxspeciesname(ispec) + end do + + end subroutine set_packagedata + + !> @brief Calculate fluid viscosity + !! + !! Calculates fluid viscosity based on concentration or + !! temperature + !! + !< + subroutine vsc_calcvisc(this) + ! -- dummy + class(GwfVscType) :: this + + ! -- local + integer(I4B) :: n + integer(I4B) :: i +! ------------------------------------------------------------------------------ + ! + ! -- Calculate the viscosity using the specified concentration and/or + ! temperature arrays + do n = 1, this%dis%nodes + do i = 1, this%nviscspecies + if (this%modelconc(i)%icbund(n) == 0) then + this%ctemp = DZERO + else + this%ctemp(i) = this%modelconc(i)%conc(n) + end if + end do + ! + this%visc(n) = calc_visc(this%ivisc, this%viscref, this%dviscdc, & + this%cviscref, this%ctemp, this%a2, & + this%a3, this%a4) + end do + ! + ! -- Return + return + end subroutine vsc_calcvisc + + !> @ brief Allocate scalars + !! + !! Allocate and initialize scalars for the VSC package. The base model + !! allocate scalars method is also called. + !! + !< + subroutine allocate_scalars(this) +! ****************************************************************************** +! allocate_scalars +! ****************************************************************************** +! +! SPECIFICATIONS: +! ------------------------------------------------------------------------------ + ! -- modules + use ConstantsModule, only: DZERO, DTEN, DEP3 + ! -- dummy + class(GwfVscType) :: this + ! -- local +! ------------------------------------------------------------------------------ + ! + ! -- allocate scalars in NumericalPackageType + call this%NumericalPackageType%allocate_scalars() + ! + ! -- Allocate + call mem_allocate(this%thermivisc, 'THERMIVISC', this%memoryPath) + call mem_allocate(this%idxtmpr, 'IDXTMPR', this%memoryPath) + call mem_allocate(this%ioutvisc, 'IOUTVISC', this%memoryPath) + call mem_allocate(this%ireadelev, 'IREADELEV', this%memoryPath) + call mem_allocate(this%iconcset, 'ICONCSET', this%memoryPath) + call mem_allocate(this%viscref, 'VISCREF', this%memoryPath) + call mem_allocate(this%a2, 'A2', this%memoryPath) + call mem_allocate(this%a3, 'A3', this%memoryPath) + call mem_allocate(this%a4, 'A4', this%memoryPath) + ! + call mem_allocate(this%nviscspecies, 'NVISCSPECIES', this%memoryPath) + ! + ! -- Initialize + this%thermivisc = 0 + this%idxtmpr = 0 + this%ioutvisc = 0 + this%ireadelev = 0 + this%iconcset = 0 + this%viscref = DEP3 + this%A2 = DTEN + this%A3 = 248.37_DP + this%A4 = 133.15_DP + ! + this%nviscspecies = 0 + ! + ! -- Return + return + end subroutine allocate_scalars + + !> @ brief Allocate arrays + !! + !! Allocate and initialize arrays for the VSC package. + !! + !< + subroutine allocate_arrays(this, nodes) + ! -- modules + ! -- dummy + class(GwfVscType) :: this + integer(I4B), intent(in) :: nodes + ! -- local + integer(I4B) :: i +! ------------------------------------------------------------------------------ + ! + ! -- Allocate + call mem_allocate(this%visc, nodes, 'VISC', this%memoryPath) + call mem_allocate(this%ivisc, this%nviscspecies, 'IVISC', this%memoryPath) + call mem_allocate(this%dviscdc, this%nviscspecies, 'DRHODC', & + this%memoryPath) + call mem_allocate(this%cviscref, this%nviscspecies, 'CRHOREF', & + this%memoryPath) + call mem_allocate(this%ctemp, this%nviscspecies, 'CTEMP', this%memoryPath) + allocate (this%cmodelname(this%nviscspecies)) + allocate (this%cauxspeciesname(this%nviscspecies)) + allocate (this%modelconc(this%nviscspecies)) + ! + ! -- Initialize + do i = 1, nodes + this%visc(i) = this%viscref + end do + ! + ! -- Initialize nviscspecies arrays + do i = 1, this%nviscspecies + this%ivisc(i) = 1 + this%dviscdc(i) = DZERO + this%cviscref(i) = DZERO + this%ctemp(i) = DZERO + this%cmodelname(i) = '' + this%cauxspeciesname(i) = '' + end do + ! + ! -- Return + return + end subroutine allocate_arrays + + !> @ brief Read Options block + !! + !! Reads the options block inside the VSC package. + !! + !< + subroutine read_options(this) + ! -- modules + use OpenSpecModule, only: access, form + use InputOutputModule, only: urword, getunit, urdaux, openfile + ! -- dummy + class(GwfVscType) :: this + ! -- local + character(len=LINELENGTH) :: warnmsg, errmsg, keyword, keyword2 + character(len=MAXCHARLEN) :: fname + integer(I4B) :: ierr + logical :: isfound, endOfBlock + ! -- formats + character(len=*), parameter :: fmtfileout = & + "(1x, 'VSC', 1x, a, 1x, 'Will be saved to file: ', & + &a, /4x, 'opened on unit: ', I7)" + character(len=*), parameter :: fmtlinear = & + "(/,1x,'Viscosity will vary linearly with temperature & + &change ')" + character(len=*), parameter :: fmtnonlinear = & + "(/,1x,'Viscosity will vary non-linearly with temperature & + &change ')" +! ------------------------------------------------------------------------------ + ! + ! -- get options block + call this%parser%GetBlock('OPTIONS', isfound, ierr, & + supportOpenClose=.true., blockRequired=.false.) + ! + ! -- parse options block if detected + if (isfound) then + write (this%iout, '(1x,a)') 'Processing VSC OPTIONS block' + do + call this%parser%GetNextLine(endOfBlock) + if (endOfBlock) exit + call this%parser%GetStringCaps(keyword) + select case (keyword) + case ('VISCREF') + this%viscref = this%parser%GetDouble() + write (this%iout, '(4x,a,1pg15.6)') & + 'Reference viscosity has been set to: ', & + this%viscref + case ('VISCOSITY') + call this%parser%GetStringCaps(keyword) + if (keyword == 'FILEOUT') then + call this%parser%GetString(fname) + this%ioutvisc = getunit() + call openfile(this%ioutvisc, this%iout, fname, 'DATA(BINARY)', & + form, access, 'REPLACE') + write (this%iout, fmtfileout) & + 'VISCOSITY', fname, this%ioutvisc + else + errmsg = 'Optional VISCOSITY keyword must be '// & + 'followed by FILEOUT' + call store_error(errmsg) + end if + case ('TEMPERATURE_SPECIES_NAME') + call this%parser%GetStringCaps(this%name_temp_spec) + write (this%iout, '(4x, a)') 'Temperature species name set to: '// & + trim(this%name_temp_spec) + case ('THERMAL_FORMULATION') + call this%parser%GetStringCaps(keyword2) + if (trim(adjustl(keyword2)) == 'LINEAR') this%thermivisc = 1 + if (trim(adjustl(keyword2)) == 'NONLINEAR') this%thermivisc = 2 + select case (this%thermivisc) + case (1) + write (this%iout, fmtlinear) + case (2) + write (this%iout, fmtnonlinear) + end select + case ('THERMAL_A2') + this%a2 = this%parser%GetDouble() + if (this%thermivisc == 2) then + write (this%iout, '(4x,a,1pg15.6)') & + 'A2 in nonlinear viscosity formulation has been set to: ', & + this%a2 + else + write (this%iout, '(4x,a,/,4x,a,/,4x,a)') 'THERMAL_A2 specified by user & + &in VSC Package input file. LINEAR viscosity ', 'formulation also & + &specified. THERMAL_A2 will not affect ', 'viscosity calculations.' + end if + case ('THERMAL_A3') + this%a3 = this%parser%GetDouble() + if (this%thermivisc == 2) then + write (this%iout, '(4x,a,1pg15.6)') & + 'A3 in nonlinear viscosity formulation has been set to: ', & + this%a3 + else + write (this%iout, '(4x,a,/,4x,a,/,4x,a)') 'THERMAL_A3 specified by user & + &in VSC Package input file. LINEAR viscosity ', 'formulation also & + &specified. THERMAL_A3 will not affect ', 'viscosity calculations.' + end if + case ('THERMAL_A4') + this%a4 = this%parser%GetDouble() + if (this%thermivisc == 2) then + write (this%iout, '(4x,a,1pg15.6)') & + 'A4 in nonlinear viscosity formulation has been set to: ', & + this%a4 + else + write (this%iout, '(4x,a,/,4x,a,/,4x,a)') 'THERMAL_A4 specified by user & + &in VSC Package input file. LINEAR viscosity ', 'formulation also & + &specified. THERMAL_A4 will not affect ', 'viscosity calculations.' + end if + case default + write (errmsg, '(4x,a,a)') '**Error. Unknown VSC option: ', & + trim(keyword) + call store_error(errmsg) + call this%parser%StoreErrorUnit() + end select + end do + ! + if (this%thermivisc == 1) then + if (this%a2 == 0.0) then + write (errmsg, '(a)') 'LINEAR option selected for varying & + &viscosity with temperature, but A1, a surrogate for & + &dVISC/dT, set equal to 0.0' + call store_error(errmsg) + end if + end if + if (this%thermivisc > 1) then + if (this%a2 == 0) then + write (warnmsg, '(a)') 'NONLINEAR option selected for & + &varying viscosity with temperature, but A2 set equal to & + &zero which may lead to unintended values for viscosity' + call store_warning(errmsg) + end if + if (this%a3 == 0) then + write (warnmsg, '(a)') 'NONLINEAR option selected for & + &varying viscosity with temperature,, but A3 set equal to & + &zero which may lead to unintended values for viscosity' + call store_warning(warnmsg) + end if + if (this%a4 == 0) then + write (warnmsg, '(a)') 'NONLINEAR option selected for & + &varying viscosity with temperature, BUT A4 SET EQUAL TO & + &zero which may lead to unintended values for viscosity' + call store_warning(warnmsg) + end if + end if + end if + ! + write (this%iout, '(/,1x,a)') 'end of VSC options block' + ! + ! -- Return + return + end subroutine read_options + + !> @brief Sets options as opposed to reading them from a file + !< + subroutine set_options(this, input_data) + class(GwfVscType) :: this + type(GwfVscInputDataType), intent(in) :: input_data !< the input data to be set + + this%viscref = input_data%viscref + ! + ! -- Return + return + end subroutine set_options + + !> @ brief Set pointers to concentration(s) + !! + !! Pass in a gwt model name, concentration array, and ibound, + !! and store a pointer to these in the VSC package so that + !! viscosity can be calculated from them. This routine is called + !! from the gwfgwt exchange in the exg_ar() method. + !! + !< + subroutine set_concentration_pointer(this, modelname, conc, icbund, istmpr) + ! -- modules + ! -- dummy + class(GwfVscType) :: this + character(len=LENMODELNAME), intent(in) :: modelname + real(DP), dimension(:), pointer :: conc + integer(I4B), dimension(:), pointer :: icbund + integer(I4B), optional, intent(in) :: istmpr + ! -- local + integer(I4B) :: i + logical :: found +! ------------------------------------------------------------------------------ + ! + this%iconcset = 1 + found = .false. + do i = 1, this%nviscspecies + if (this%cmodelname(i) == modelname) then + this%modelconc(i)%conc => conc + this%modelconc(i)%icbund => icbund + found = .true. + exit + end if + end do + ! + ! -- Return + return + end subroutine set_concentration_pointer + +end module GwfVscModule diff --git a/src/Model/ModelUtilities/BoundaryPackage.f90 b/src/Model/ModelUtilities/BoundaryPackage.f90 index c978db1bbf5..e602b2436ed 100644 --- a/src/Model/ModelUtilities/BoundaryPackage.f90 +++ b/src/Model/ModelUtilities/BoundaryPackage.f90 @@ -80,9 +80,13 @@ module BndModule real(DP), dimension(:), pointer, contiguous :: simtomvr => null() !< simulated to mover values ! ! -- water mover flag and object - integer(I4B), pointer :: imover => null() !< flag indicating of the mover is active in the package + integer(I4B), pointer :: imover => null() !< flag indicating if the mover is active in the package type(PackageMoverType), pointer :: pakmvrobj => null() !< mover object for package ! + ! -- viscosity flag and safe-copy of conductance array + integer(I4B), pointer :: ivsc => null() !< flag indicating if viscosity is active in the model + real(DP), dimension(:), pointer, contiguous :: condinput => null() !< stores user-specified conductance values + ! ! -- timeseries type(TimeSeriesManagerType), pointer :: TsManager => null() !< time series manager type(TimeArraySeriesManagerType), pointer :: TasManager => null() !< time array series manager @@ -152,6 +156,12 @@ module BndModule ! -- procedure to support time series procedure, public :: bnd_rp_ts ! + ! -- procedure to inform package that viscosity active + procedure, public :: bnd_activate_viscosity + ! + ! -- procedure to backup user-specified conductance + procedure, private :: bnd_store_user_cond + ! end type BndType contains @@ -362,6 +372,12 @@ subroutine bnd_rp(this) this%packName, this%tsManager, this%iscloc) this%nbound = nlist ! + ! -- save user-specified conductance if vsc package is active + if (this%ivsc == 1) then + call this%bnd_store_user_cond(nlist, this%nodelist, this%bound, & + this%condinput) + end if + ! ! Define the tsLink%Text value(s) appropriately. ! E.g. for WEL package, entry 1, assign tsLink%Text = 'Q' ! For RIV package, entry 1 text = 'STAGE', entry 2 text = 'COND', @@ -902,6 +918,7 @@ subroutine bnd_da(this) call mem_deallocate(this%nodelist, 'NODELIST', this%memoryPath) call mem_deallocate(this%noupdateauxvar, 'NOUPDATEAUXVAR', this%memoryPath) call mem_deallocate(this%bound, 'BOUND', this%memoryPath) + call mem_deallocate(this%condinput, 'CONDINPUT', this%memoryPath) call mem_deallocate(this%hcof, 'HCOF', this%memoryPath) call mem_deallocate(this%rhs, 'RHS', this%memoryPath) call mem_deallocate(this%simvals, 'SIMVALS', this%memoryPath) @@ -958,6 +975,7 @@ subroutine bnd_da(this) call mem_deallocate(this%imover) call mem_deallocate(this%npakeq) call mem_deallocate(this%ioffset) + call mem_deallocate(this%ivsc) ! ! -- deallocate methods on objects call this%obs%obs_da() @@ -1016,6 +1034,9 @@ subroutine allocate_scalars(this) ! -- allocate the object and assign values to object variables call mem_allocate(this%imover, 'IMOVER', this%memoryPath) ! + ! -- allocate flag for determining if vsc active + call mem_allocate(this%ivsc, 'IVSC', this%memoryPath) + ! ! -- allocate scalars for packages that add rows to the matrix (e.g. MAW) call mem_allocate(this%npakeq, 'NPAKEQ', this%memoryPath) call mem_allocate(this%ioffset, 'IOFFSET', this%memoryPath) @@ -1043,6 +1064,7 @@ subroutine allocate_scalars(this) this%imover = 0 this%npakeq = 0 this%ioffset = 0 + this%ivsc = 0 ! ! -- Set pointer to model inewton variable call mem_setptr(imodelnewton, 'INEWTON', create_mem_path(this%name_model)) @@ -1092,6 +1114,10 @@ subroutine allocate_arrays(this, nodelist, auxvar) call mem_allocate(this%bound, this%ncolbnd, this%maxbound, 'BOUND', & this%memoryPath) ! + !-- Allocate array for storing user-specified conductances + ! Will be reallocated to size maxbound if vsc active + call mem_allocate(this%condinput, 0, 'CONDINPUT', this%memoryPath) + ! ! -- Allocate hcof and rhs call mem_allocate(this%hcof, this%maxbound, 'HCOF', this%memoryPath) call mem_allocate(this%rhs, this%maxbound, 'RHS', this%memoryPath) @@ -1472,6 +1498,47 @@ subroutine bnd_read_dimensions(this) return end subroutine bnd_read_dimensions + !> @ brief Store user-specified conductances when vsc is active + !! + !! VSC will update boundary package conductance values. Because + !! viscosity can change every stress period, but user-specified + !! conductances may not, the base user-input should be stored in + !! backup array so that viscosity-updated conductances may be + !! recalculated every stress period/time step + !! + !< + subroutine bnd_store_user_cond(this, nlist, nodelist, rlist, condinput) + ! -- modules + use SimModule, only: store_error + ! -- dummy variables + class(BndType), intent(inout) :: this !< BndType object + integer(I4B), intent(in) :: nlist + integer(I4B), dimension(:), pointer, contiguous, intent(inout) :: nodelist + real(DP), dimension(:, :), pointer, contiguous, intent(in) :: rlist + real(DP), dimension(:), pointer, contiguous, intent(inout) :: condinput + ! -- local variables + integer(I4B) :: l + integer(I4B) :: nodeu, noder + character(len=LINELENGTH) :: nodestr + ! + ! -- store backup copy of conductance values + do l = 1, nlist + nodeu = nodelist(l) + noder = this%dis%get_nodenumber(nodeu, 0) + if (noder <= 0) then + call this%dis%nodeu_to_string(nodeu, nodestr) + write (errmsg, *) & + ' Cell is outside active grid domain: '// & + trim(adjustl(nodestr)) + call store_error(errmsg) + end if + condinput(l) = rlist(2, l) + end do + ! + ! -- return + return + end subroutine + !> @ brief Read initial parameters for package !! !! Read initial parameters for a boundary package. This method is not @@ -2003,4 +2070,38 @@ subroutine save_print_model_flows(icbcfl, ibudfl, icbcun, iprflow, & return end subroutine save_print_model_flows + !> @brief Activate viscosity terms + !! + !! Method to activate addition of viscosity terms when package type + !! is DRN, GHB, or RIV (method not needed by other packages at this point) + !! + !< + subroutine bnd_activate_viscosity(this) + ! -- modules + use MemoryManagerModule, only: mem_reallocate + ! -- dummy variables + class(BndType), intent(inout) :: this !< BndType object + ! -- local variables + integer(I4B) :: i + ! + ! -- Set ivsc and reallocate viscratios to be of size MAXBOUND + this%ivsc = 1 + ! + ! -- Allocate array for storing user-specified conductances + ! modified by updated viscosity values + call mem_reallocate(this%condinput, this%maxbound, 'CONDINPUT', & + this%memoryPath) + do i = 1, this%maxbound + this%condinput(i) = DZERO + end do + ! + ! -- Notify user via listing file viscosity accounted for by standard + ! boundary package. + write (this%iout, '(/1x,a,a)') 'VISCOSITY ACTIVE IN ', & + trim(this%filtyp)//' PACKAGE CALCULATIONS: '//trim(adjustl(this%packName)) + ! + ! -- return + return + end subroutine bnd_activate_viscosity + end module BndModule diff --git a/src/Model/ModelUtilities/GwfVscInputData.f90 b/src/Model/ModelUtilities/GwfVscInputData.f90 new file mode 100644 index 00000000000..fb39c0d1994 --- /dev/null +++ b/src/Model/ModelUtilities/GwfVscInputData.f90 @@ -0,0 +1,57 @@ +module GwfVscInputDataModule + use KindModule, only: I4B, DP + use ConstantsModule, only: LENMODELNAME, LENAUXNAME, DZERO + + implicit none + private + + !> Data structure to transfer input configuration to the + !< VSC package, as opposed to reading from file + type, public :: GwfVscInputDataType + + ! options + real(DP) :: viscref !< see VSC for description + ! dim + integer(I4B) :: nviscspecies !< see VSC for description + ! pkg data + integer(I4B), dimension(:), pointer :: ivisc => null() !< indicates if species uses linear or nonlinear relationship + real(DP), dimension(:), pointer, contiguous :: dviscdc => null() !< see VSC for description + real(DP), dimension(:), pointer, contiguous :: cviscref => null() !< see VSC for description + character(len=LENMODELNAME), dimension(:), allocatable :: cmodelname !< see VSC for description + character(len=LENAUXNAME), dimension(:), allocatable :: cauxspeciesname !< see VSC for description + + contains + procedure, pass(this) :: construct + procedure, pass(this) :: destruct + end type GwfVscInputDataType + +contains + +!> @brief Allocate the input data +!< + subroutine construct(this, nviscspecies) + class(GwfVscInputDataType) :: this !< the input data block + integer(I4B) :: nviscspecies !< the number of species + + allocate (this%ivisc(nviscspecies)) + allocate (this%dviscdc(nviscspecies)) + allocate (this%cviscref(nviscspecies)) + allocate (this%cmodelname(nviscspecies)) + allocate (this%cauxspeciesname(nviscspecies)) + + end subroutine construct + + !> @brief clean up + !< + subroutine destruct(this) + class(GwfVscInputDataType) :: this !< the input data block + + deallocate (this%ivisc) + deallocate (this%dviscdc) + deallocate (this%cviscref) + deallocate (this%cmodelname) + deallocate (this%cauxspeciesname) + + end subroutine destruct + +end module GwfVscInputDataModule diff --git a/src/meson.build b/src/meson.build index 0532308f8cd..301b07c431f 100644 --- a/src/meson.build +++ b/src/meson.build @@ -62,6 +62,7 @@ modflow_sources = files( 'Model' / 'GroundWaterFlow' / 'gwf3tvk8.f90', 'Model' / 'GroundWaterFlow' / 'gwf3tvs8.f90', 'Model' / 'GroundWaterFlow' / 'gwf3uzf8.f90', + 'Model' / 'GroundWaterFlow' / 'gwf3vsc8.f90', 'Model' / 'GroundWaterFlow' / 'gwf3wel8.f90', 'Model' / 'GroundWaterTransport' / 'gwt1.f90', 'Model' / 'GroundWaterTransport' / 'gwt1adv1.f90', @@ -90,6 +91,7 @@ modflow_sources = files( 'Model' / 'ModelUtilities' / 'GwfMvrPeriodData.f90', 'Model' / 'ModelUtilities' / 'GwfNpfOptions.f90', 'Model' / 'ModelUtilities' / 'GwfStorageUtils.f90', + 'Model' / 'ModelUtilities' / 'GwfVscInputData.f90', 'Model' / 'ModelUtilities' / 'GwtAdvOptions.f90', 'Model' / 'ModelUtilities' / 'GwtDspOptions.f90', 'Model' / 'ModelUtilities' / 'GwtSpc.f90',