diff --git a/autotest/test_gwe_sfe_abc00.py b/autotest/test_gwe_sfe_abc00.py new file mode 100644 index 00000000000..571d1f53f5b --- /dev/null +++ b/autotest/test_gwe_sfe_abc00.py @@ -0,0 +1,522 @@ +# Test the use of the atmospheric boundary condition utility used in conjunction +# with the SFE advanced package. This test is a single cell with a single +# reach. Channel flow characteristics are unrealistic: Manning's n is +# unrealistically low and slope is extremely high. These conditions result in +# an extremely high streamflow velocity that results in nearly all of the heat +# being added to the channel exiting at the outlet with very near negligle heat +# storage increases in the channel. This test only uses latent heat flux (lhf). +# The result is a -1 deg C change in temperature in the +# streamflow - an easy result to confirm in this test. + +import math +import os + +import flopy +import numpy as np +import pandas as pd +import pytest +from framework import TestFramework + +cases = ["sfe-abc-lhf"] + +DCTOK = 273.16 + +# Model units +length_units = "m" +time_units = "seconds" + +# model domain and grid definition +nrow = 1 +ncol = 1 +nlay = 1 +delr = 1.0 +delc = 1.0 +xmax = ncol * delr +ymax = nrow * delc + +ibound = 1 +top = 1.0 +botm = 0.0 +strthd = 0.0 +chd_on = False + +# model input parameters +k11 = 500.0 +strt_gw_temp = 99.0 + +ss = 0.00001 +sy = 0.20 +hani = 1 +laytyp = 1 + +# Package boundary conditions +sfr_evaprate = 0.0 +rhk = 0.0 +rwid = 1.0 +strm_temp = 11.0 +surf_Q_in = [ + [10.0], +] +# sensible and latent heat flux parameter values +wspd = 126005.30 # unrealistically high to drive a -1C change +tatm = 5.0 + DCTOK +# shortwave radiation parameter values +solr = 47880870.9 # unrealistically high to drive a 1 deg C rise in stream temperature +shd = 1.0 # 100% shade "turns off" solar flux +swrefl = 0.03 +rh = 30.0 # percent + +# Transport related parameters +porosity = sy # porosity (unitless) +K_therm = 2.0 # Thermal conductivity # ($W/m/C$) +rhow = 1000.0 # Density of water ($kg/m^3$) +rhos = 2650.0 # Density of the aquifer material ($kg/m^3$) +rhoa = 1.225 # Density of the atmosphere ($kg/m^3$) +Cpw = 4180.0 # Heat capacity of water ($J/kg/C$) +Cps = 880.0 # Heat capacity of the solids ($J/kg/C$) +Cpa = 717.0 # Heat capacity of the atmosphere ($J/kg/C$) +lhv = 2454000.0 # Latent heat of vaporization ($J/kg$) +c_d = 0.0 # Drag coefficient ($unitless$) !! +wf_slope = 1.383e-08 # wind function slope ($1/mbar$) +wf_int = 3.445e-09 # wind function intercept ($m/s$) +# Thermal conductivity of the streambed material ($W/m/C$) +K_therm_strmbed = 0.0 +rbthcnd = 0.0001 + +# time params +steady = {0: True, 1: False} +transient = {0: False, 1: True} +nstp = [1] +tsmult = [1] +perlen = [1] + +nouter, ninner = 1000, 300 +hclose, rclose, relax = 1e-10, 1e-10, 0.97 + +# +# MODFLOW 6 flopy GWF object +# + + +def build_models(idx, test): + # Base simulation and model name and workspace + ws = test.workspace + name = cases[idx] + + print(f"Building model...{name}") + + # generate names for each model + gwfname = "gwf-" + name + gwename = "gwe-" + 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=f"{gwfname}.ims", + ) + 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 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", "ALL")], + ) + + # Instantiate streamflow routing package + # Determine the middle row and store in rMid (account for 0-base) + rMid = 1 + # sfr data + nreaches = ncol + rlen = delr + roughness = 1e-10 + rbth = 0.1 + strmbd_hk = rhk + strm_up = 0.95 + strm_dn = 0.94 + # divide by 10 to further reduce slop + slope = 0.04 + nconn = 0 + ustrf = 1.0 + ndv = 0 + strm_incision = 0.05 + + packagedata = [] + for irch in range(nreaches): + rp = [ + irch, + (0, 0, 0), + rlen, + rwid, + slope, + top - strm_incision, + rbth, + strmbd_hk, + roughness, + nconn, + ustrf, + ndv, + ] + packagedata.append(rp) + + connectiondata = [0] + + sfr_perioddata = {} + for t in np.arange(len(surf_Q_in[idx])): + sfrbndx = [] + for i in np.arange(nreaches): + if i == 0: + sfrbndx.append([i, "INFLOW", surf_Q_in[idx][t]]) + # sfrbndx.append([i, "EVAPORATION", sfr_evaprate]) + + sfr_perioddata.update({t: sfrbndx}) + + # Instantiate SFR observation points + sfr_obs = { + f"{gwfname}.sfr.obs.csv": [ + ("rch1_depth", "depth", 1), + ("rch1_outf", "ext-outflow", 1), + ("rch1_wetwidth", "wet-width", 1), + ], + "digits": 8, + "print_input": True, + "filename": gwfname + ".sfr.obs", + } + + budpth = f"{gwfname}.sfr.cbc" + flopy.mf6.ModflowGwfsfr( + gwf, + save_flows=True, + print_stage=True, + print_flows=True, + print_input=True, + length_conversion=1.0, + time_conversion=1.0, + budget_filerecord=budpth, + mover=False, + nreaches=nreaches, + packagedata=packagedata, + connectiondata=connectiondata, + perioddata=sfr_perioddata, + observations=sfr_obs, + pname="SFR", + filename=f"{gwfname}.sfr", + ) + + # -------------------------------------------------- + # Setup the GWE model for simulating heat transport + # -------------------------------------------------- + gwe = flopy.mf6.ModflowGwe(sim, modelname=gwename) + + # Instantiating solver for GWT + imsgwe = flopy.mf6.ModflowIms( + sim, + print_option="ALL", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="NONE", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=f"{rclose} strict", + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + filename=f"{gwename}.ims", + ) + sim.register_ims_package(imsgwe, [gwename]) + + # Instantiating DIS for GWE + flopy.mf6.ModflowGwedis( + gwe, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + pname="DIS", + filename=f"{gwename}.dis", + ) + + # Instantiate Mobile Storage and Transfer package + flopy.mf6.ModflowGweest( + gwe, + save_flows=True, + porosity=porosity, + heat_capacity_water=Cpw, + density_water=rhow, + latent_heat_vaporization=lhv, + heat_capacity_solid=Cps, + density_solid=rhos, + pname="EST", + filename=f"{gwename}.est", + ) + + # Instantiate Energy Transport Initial Conditions package + flopy.mf6.ModflowGweic(gwe, strt=strt_gw_temp) + + # Instantiate Advection package + flopy.mf6.ModflowGweadv(gwe, scheme="UPSTREAM") + + # Instantiate Dispersion package (also handles conduction) + flopy.mf6.ModflowGwecnd( + gwe, + xt3d_off=True, + ktw=0.5918, + kts=0.2700, + pname="CND", + filename=f"{gwename}.cnd", + ) + + # Instantiating MODFLOW 6 transport source-sink mixing package + # [b/c at least one boundary back is active (SFR), ssm must be on] + sourcerecarray = [[]] + flopy.mf6.ModflowGwessm(gwe, sources=sourcerecarray, filename=f"{gwename}.ssm") + + # Instantiate Streamflow Energy Transport package + sfepackagedata = [] + for irno in range(ncol): + t = (irno, strm_temp, K_therm_strmbed, rbthcnd) + sfepackagedata.append(t) + + sfeperioddata = [] + for irno in range(ncol): + if irno == 0: + sfeperioddata.append((irno, "INFLOW", strm_temp)) + + # Instantiate SFE observation points + sfe_obs = { + f"{gwename}.sfe.obs.csv": [ + ("rch1_outftemp", "temperature", 1), + ("rch1_outfener", "ext-outflow", 1), + ], + "digits": 8, + "print_input": True, + "filename": gwename + ".sfe.obs", + } + + abc_filename = f"{gwename}.sfe.abc" + sfe = flopy.mf6.modflow.ModflowGwesfe( + gwe, + boundnames=False, + save_flows=True, + print_input=False, + print_flows=True, + print_temperature=True, + temperature_filerecord=gwename + ".sfe.bin", + budget_filerecord=gwename + ".sfe.bud", + packagedata=sfepackagedata, + reachperioddata=sfeperioddata, + flow_package_name="SFR", + observations=sfe_obs, + pname="SFE", + filename=f"{gwename}.sfe", + ) + + # Abc utility + abc_spd = {} + for kper in range(len(nstp)): + spd = [] + for irno in range(ncol): + spd.append([irno, "WSPD", wspd]) + spd.append([irno, "TATM", tatm]) + spd.append([irno, "SOLR", solr]) + spd.append([irno, "SHD", shd]) + spd.append([irno, "SWREFL", swrefl]) + spd.append([irno, "RH", rh]) + abc_spd[kper] = spd + + abc = flopy.mf6.ModflowUtlabc( + sfe, + print_input=True, + density_air=rhoa, + heat_capacity_air=Cpa, + drag_coefficient=c_d, + wind_func_slope=wf_slope, + wind_func_int=wf_int, + reachperioddata=abc_spd, + filename=abc_filename, + ) + + # Instantiate Output Control package for transport + flopy.mf6.ModflowGweoc( + gwe, + temperature_filerecord=f"{gwename}.ucn", + saverecord=[("TEMPERATURE", "ALL")], + temperatureprintrecord=[("COLUMNS", 3, "WIDTH", 20, "DIGITS", 8, "GENERAL")], + printrecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")], + filename=f"{gwename}.oc", + ) + + # Instantiate Gwf-Gwe Exchange package + flopy.mf6.ModflowGwfgwe( + sim, + exgtype="GWF6-GWE6", + exgmnamea=gwfname, + exgmnameb=gwename, + filename=f"{gwename}.gwfgwe", + ) + + return sim, None + + +# sim, dum = build_models(0, r"c:\temp\_shf00") +# sim.write_simulation() + + +def calc_ener_transfer(updated_strm_temp, mf_strm_wid): + L = 2499.64 - (2.51 * (updated_strm_temp - DCTOK)) + e_w = 6.1275 * math.exp( + 17.2693882 * ((updated_strm_temp - DCTOK) / (updated_strm_temp - 35.86)) + ) + e_s = 6.1275 * math.exp(17.2693882 * ((tatm - DCTOK) / (tatm - 35.86))) + e_a = (rh / 100) * e_s + vap_press_deficit = e_w - e_a + wind_function = wf_int + wf_slope * wspd + Ev = wind_function * vap_press_deficit + lhf_ener_per_sqm = Ev * L * rhow + + ener_transfer = lhf_ener_per_sqm * delr * mf_strm_wid + + return -ener_transfer + + +def check_output(idx, test): + print("evaluating results...") + msg0 = "Stream channel width less than 1.0, should be 1.0 m" + + # read flow results from model + name = cases[idx] + gwfname = "gwf-" + name + gwename = "gwe-" + name + + # calc expected rise in temperature independent of mf6 + + fpth = os.path.join(test.workspace, gwfname + ".sfr.obs.csv") + assert os.path.isfile(fpth) + df = pd.read_csv(fpth) + mf_strm_wid = df.loc[0, "RCH1_WETWIDTH"].copy() + # confirm stream width is 1.0 m + assert np.isclose(mf_strm_wid, 1.0, atol=1e-9), msg0 + + # confirm that the energy added to the stream results in a -1C change in temp + # temperature gradient + + tgrad = (tatm - DCTOK) - strm_temp + shf_ener_per_sqm = c_d * rhoa * Cpa * wspd * tgrad + swr_ener_per_sqm = solr * (1 - shd) * (1 - swrefl) + + # latent calcs + chng = 1 + strt_strm_temp = strm_temp + updated_strm_temp = strm_temp + while chng > hclose: + ener_transfer = calc_ener_transfer(updated_strm_temp + DCTOK, mf_strm_wid) + temp_change = ener_transfer / (surf_Q_in[idx][0] * Cpw * rhow) + updated_temp = strt_strm_temp + temp_change + chng = abs(updated_strm_temp - updated_temp) + updated_strm_temp = updated_temp + + fpth2 = os.path.join(test.workspace, gwename + ".sfe.obs.csv") + assert os.path.isfile(fpth2) + df2 = pd.read_csv(fpth2) + + # confirm 1 deg C decrease in temp + + msg1 = "Python temperature change is = " + str(temp_change) + msg2 = "MODFLOW temperature = " + str(df2.loc[0, "RCH1_OUTFTEMP"]) + msg3 = "MODFLOW temperature change is " + str( + strm_temp - df2.loc[0, "RCH1_OUTFTEMP"] + ) + + assert np.isclose( + df2.loc[0, "RCH1_OUTFTEMP"], strt_strm_temp + temp_change, atol=1e-6 + ), msg2 + ". " + msg3 + ". " + msg1 + + +# - No need to change any code below +@pytest.mark.parametrize( + "idx, name", + list(enumerate(cases)), +) +def test_mf6model(idx, name, function_tmpdir, targets): + test = TestFramework( + name=name, + workspace=function_tmpdir, + targets=targets, + build=lambda t: build_models(idx, t), + check=lambda t: check_output(idx, t), + ) + test.run() diff --git a/autotest/test_gwe_sfe_abc01.py b/autotest/test_gwe_sfe_abc01.py new file mode 100644 index 00000000000..0e163bc6e24 --- /dev/null +++ b/autotest/test_gwe_sfe_abc01.py @@ -0,0 +1,494 @@ +""" + Test the use of the atmospheric boundary condition utility used in conjunction + with the SFE advanced package. This test uses five reaches all hosted within + a single cell. The test is meant to error-out because it specifies a user- + defined evaporation rate but also tries to invoke the ABC package to calculate + latent heat flux. MF6 should prevent the user from invoking both approaches + and instead force the user to choose. If this is successful, it means MF6 + successfully errored out. + + + Reach configuration: + + \ + \ Reach 1 + \ + \ + Reach 2 v Reach 4 Reach 5 + ------>+------------+-----------> + ^ + / + / + / Reach 3 (shf) + / +""" + +import flopy +import numpy as np +import pytest +from framework import TestFramework + +cases = ["sfe-abc-fail"] + +DCTOK = 273.16 + +# Model units +length_units = "m" +time_units = "seconds" + +# model domain and grid definition + +nrow = 1 +ncol = 1 +nlay = 1 +delr = 100.0 +delc = 100.0 +xmax = ncol * delr +ymax = nrow * delc + +ibound = 1 +top = 1.0 +botm = 0.0 +strthd = 0.0 +chd_on = False + +# model input parameters +k11 = 500.0 +strt_gw_temp = 99.0 + +ss = 0.00001 +sy = 0.20 +hani = 1 +laytyp = 1 + +# Package boundary conditions +sfr_evaprate = 0.0 +rhk = 0.0 +rwid = 1.0 +rlen = 1.0 +strm_initial_temp = 11.0 +surf_Q_in = 10.0 + +# sensible and latent heat flux parameter values +wspd = 1.0 +patm = 954.0 +tatm = 20.0 +tatmK = tatm + DCTOK +# shortwave radiation parameter values + +# unrealistically high to drive a 0.1 deg C rise in stream temperature +solr = 4180000.0 +shd = 0.20 +swrefl = 0.10 +rh = 20.0 +lwrefl = 0.03 # Fogg et al 2023 + +# atmosphere composition adjustment factor (using dummy value to drive +# half a degree change) +atmc = 0.80 + +# latent heat flux parameter values (these values are specified in the options +# block and are therefore constant across reaches +c_d = 0.0 # Drag coefficient ($unitless$) +wf_slope = 1.383e-08 # wind function slope ($1/mbar$) +wf_int = 3.445e-09 # wind function intercept ($m/s$) +emiss_water = 0.95 # Fogg et al 2023 +emiss_riparian = 0.97 # Fogg et al 2023 +stephan_boltzmann = 5.670374419e-08 + +# Transport related parameters +porosity = sy # porosity (unitless) +K_therm = 2.0 # Thermal conductivity # ($W/m/C$) +rhow = 1000.0 # Density of water ($kg/m^3$) +rhos = 2650.0 # Density of the aquifer material ($kg/m^3$) +rhoa = 1.225 # Density of the atmosphere ($kg/m^3$) +Cpw = 4180.0 # Heat capacity of water ($J/kg/C$) +Cps = 880.0 # Heat capacity of the solids ($J/kg/C$) +Cpa = 717.0 # Heat capacity of the atmosphere ($J/kg/C$) +lhv = 2454000.0 # Latent heat of vaporization ($J/kg$) + +# Thermal conductivity of the streambed material ($W/m/C$) +K_therm_strmbed = 0.0 +rbthcnd = 0.0001 + +# time params +steady = {0: True, 1: False} +transient = {0: False, 1: True} +nstp = [1] +tsmult = [1] +perlen = [1] + +nouter, ninner = 1000, 300 +hclose, rclose, relax = 1e-10, 1e-10, 0.97 + +# +# MODFLOW 6 flopy GWF object +# + + +def build_models(idx, test): + # Base simulation and model name and workspace + ws = test.workspace + name = cases[idx] + + print(f"Building model...{name}") + + # generate names for each model + gwfname = "gwf-" + name + gwename = "gwe-" + 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=f"{gwfname}.ims", + ) + 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 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", "ALL")], + ) + + # Instantiate streamflow routing package + # Determine the middle row and store in rMid (account for 0-base) + rMid = 1 + # sfr data + roughness = 1e-10 + rbth = 0.1 + strmbd_hk = rhk + strm_up = 0.95 + strm_dn = 0.94 + # divide by 10 to further reduce slop + slope = 0.04 + nconn = 0 + ustrf = 1.0 + ndv = 0 + strm_incision = 0.05 + + # explicitly set connections + conns = [(0, -3), (1, -3), (2, -3), (3, 0, 1, 2, -4), (4, 3)] + nreaches = len(conns) + + packagedata = [] + for irch in range(nreaches): + ncon = len(conns[irch]) - 1 + rp = [ + irch, + (0, 0, 0), + rlen, + rwid, + slope, + top - strm_incision, + rbth, + strmbd_hk, + roughness, + ncon, + ustrf, + ndv, + ] + packagedata.append(rp) + + sfr_perioddata = {} + for t in np.arange(len(perlen)): + sfrbndx = [] + for i in np.arange(nreaches): + # only specify inflow for the first three reaches + if i < 3: + sfrbndx.append([i, "INFLOW", surf_Q_in]) + + sfrbndx.append([4, "EVAPORATION", "0.005"]) + + sfr_perioddata.update({t: sfrbndx}) + + # Instantiate SFR observation points + sfr_obs = { + f"{gwfname}.sfr.obs.csv": [ + ("rch1_depth", "depth", 1), + ("rch1_outf", "ext-outflow", 1), + ("rch1_wetwidth", "wet-width", 1), + ("rch2_wetwidth", "wet-width", 2), + ("rch3_wetwidth", "wet-width", 3), + ("rch4_wetwidth", "wet-width", 4), + ("rch5_wetwidth", "wet-width", 5), + ("rch1_stg", "stage", 1), + ("rch2_stg", "stage", 2), + ("rch3_stg", "stage", 3), + ("rch4_stg", "stage", 4), + ("rch5_stg", "stage", 5), + ], + "digits": 15, + "print_input": True, + "filename": gwfname + ".sfr.obs", + } + + budpth = f"{gwfname}.sfr.cbc" + flopy.mf6.ModflowGwfsfr( + gwf, + save_flows=True, + print_stage=True, + print_flows=True, + print_input=True, + length_conversion=1.0, + time_conversion=1.0, + budget_filerecord=budpth, + mover=False, + nreaches=nreaches, + packagedata=packagedata, + connectiondata=conns, + perioddata=sfr_perioddata, + observations=sfr_obs, + pname="SFR", + filename=f"{gwfname}.sfr", + ) + + # -------------------------------------------------- + # Setup the GWE model for simulating heat transport + # -------------------------------------------------- + gwe = flopy.mf6.ModflowGwe(sim, modelname=gwename) + + # Instantiating solver for GWT + imsgwe = flopy.mf6.ModflowIms( + sim, + print_option="ALL", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="NONE", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=f"{rclose} strict", + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + filename=f"{gwename}.ims", + ) + sim.register_ims_package(imsgwe, [gwename]) + + # Instantiating DIS for GWE + flopy.mf6.ModflowGwedis( + gwe, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + pname="DIS", + filename=f"{gwename}.dis", + ) + + # Instantiate Mobile Storage and Transfer package + flopy.mf6.ModflowGweest( + gwe, + save_flows=True, + porosity=porosity, + heat_capacity_water=Cpw, + density_water=rhow, + latent_heat_vaporization=lhv, + heat_capacity_solid=Cps, + density_solid=rhos, + pname="EST", + filename=f"{gwename}.est", + ) + + # Instantiate Energy Transport Initial Conditions package + flopy.mf6.ModflowGweic(gwe, strt=strt_gw_temp) + + # Instantiate Advection package + flopy.mf6.ModflowGweadv(gwe, scheme="UPSTREAM") + + # Instantiate Dispersion package (also handles conduction) + flopy.mf6.ModflowGwecnd( + gwe, + xt3d_off=True, + ktw=0.5918, + kts=0.2700, + pname="CND", + filename=f"{gwename}.cnd", + ) + + # Instantiating MODFLOW 6 transport source-sink mixing package + # [b/c at least one boundary back is active (SFR), ssm must be on] + sourcerecarray = [[]] + flopy.mf6.ModflowGwessm(gwe, sources=sourcerecarray, filename=f"{gwename}.ssm") + + # Instantiate Streamflow Energy Transport package + sfepackagedata = [] + for irno in range(len(conns)): + t = (irno, strm_initial_temp, K_therm_strmbed, rbthcnd) + sfepackagedata.append(t) + + sfeperioddata = [] + sfeperioddata.append((0, "INFLOW", strm_initial_temp)) + sfeperioddata.append((1, "INFLOW", strm_initial_temp)) + sfeperioddata.append((2, "INFLOW", strm_initial_temp)) + + # Instantiate SFE observation points + sfe_obs = { + f"{gwename}.sfe.obs.csv": [ + ("rch5_evap", "surfevap", 5), + ], + "digits": 20, + "print_input": True, + "filename": gwename + ".sfe.obs", + } + + abc_filename = f"{gwename}.sfe.abc" + sfe = flopy.mf6.modflow.ModflowGwesfe( + gwe, + boundnames=False, + save_flows=True, + print_input=False, + print_flows=True, + print_temperature=True, + temperature_filerecord=gwename + ".sfe.bin", + budget_filerecord=gwename + ".sfe.bud", + packagedata=sfepackagedata, + reachperioddata=sfeperioddata, + flow_package_name="SFR", + observations=sfe_obs, + pname="SFE", + filename=f"{gwename}.sfe", + ) + + # Abc utility + abc_spd = {} + for kper in range(len(nstp)): + spd = [] + for irno in range(len(conns)): + spd.append([irno, "WSPD", wspd]) + spd.append([irno, "TATM", tatm]) + spd.append([irno, "SOLR", solr]) + spd.append([irno, "SHD", shd]) + spd.append([irno, "SWREFL", swrefl]) + spd.append([irno, "RH", rh]) + spd.append([irno, "ATMC", atmc]) + spd.append([irno, "PATM", patm]) + + abc_spd[kper] = spd + + abc = flopy.mf6.ModflowUtlabc( + sfe, + print_input=True, + density_air=rhoa, + heat_capacity_air=Cpa, + drag_coefficient=c_d, + emissivity_water=emiss_water, + emissivity_canopy=emiss_riparian, + wind_func_slope=wf_slope, + wind_func_int=wf_int, + longwave_reflectance=lwrefl, + reachperioddata=abc_spd, + filename=abc_filename, + ) + + # Instantiate Output Control package for transport + flopy.mf6.ModflowGweoc( + gwe, + temperature_filerecord=f"{gwename}.ucn", + saverecord=[("TEMPERATURE", "ALL")], + temperatureprintrecord=[("COLUMNS", 3, "WIDTH", 20, "DIGITS", 8, "GENERAL")], + printrecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")], + filename=f"{gwename}.oc", + ) + + # Instantiate Gwf-Gwe Exchange package + flopy.mf6.ModflowGwfgwe( + sim, + exgtype="GWF6-GWE6", + exgmnamea=gwfname, + exgmnameb=gwename, + filename=f"{gwename}.gwfgwe", + ) + + return sim, None + + +# - No need to change any code below +@pytest.mark.parametrize( + "idx, name", + list(enumerate(cases)), +) +def test_mf6model(idx, name, function_tmpdir, targets): + test = TestFramework( + name=name, + workspace=function_tmpdir, + targets=targets, + build=lambda t: build_models(idx, t), + compare=None, + xfail="fail" in name, + ) + test.run() diff --git a/autotest/test_gwe_sfe_abc02.py b/autotest/test_gwe_sfe_abc02.py new file mode 100644 index 00000000000..004bad6e25e --- /dev/null +++ b/autotest/test_gwe_sfe_abc02.py @@ -0,0 +1,690 @@ +# Test the use of the atmospheric boundary condition utility used in conjunction +# with the SFE advanced package. This test uses four cells that each host a +# single reach. Channel flow characteristics are unrealistic: Manning's n is +# unrealistically low and slope is extremely high. These conditions result in +# an extremely high streamflow velocity that results in nearly all of the heat +# being added to, or subtracted from, the channel at the outlet with near- +# negligle heat storage increases (or decreases) in the channel. +# +# This test applies each of the atmospheric boundary conditions heat fluxes to +# only 1 of the reaches. In other words, the shortwave radiation is applied to +# the first reach, longwave radiation to the second reach, sensible heat flux +# to the third, and latent heat flux to the fourth reach. The idea is that the +# temperature change at the outlet of each reach, before it flows into the next +# downstream reach, is equal to externally calculated changes. + +""" + + + \ + \ Reach 1 (swr) + \ + \ + Reach 2 (lwr) v Reach 4 (precip) Reach 5 (lhf) + ------>+---------------------+---------------> + ^ + / + / + / Reach 3 (shf) + / +""" + +import math +import os + +import flopy +import numpy as np +import pandas as pd +import pytest +from framework import TestFramework + +cases = ["sfe-abc"] + +DCTOK = 273.16 + +# Model units +length_units = "m" +time_units = "seconds" + +# model domain and grid definition + +nrow = 1 +ncol = 1 +nlay = 1 +delr = 100.0 +delc = 100.0 +xmax = ncol * delr +ymax = nrow * delc + +ibound = 1 +top = 1.0 +botm = 0.0 +strthd = 0.0 +chd_on = False + +# model input parameters +k11 = 500.0 +strt_gw_temp = 99.0 + +ss = 0.00001 +sy = 0.20 +hani = 1 +laytyp = 1 + +# Package boundary conditions +sfr_evaprate = 0.0 +rhk = 0.0 +rwid = 1.0 +rlen = 1.0 +# strm_temp = 11.0 +strm_in_init = [11.0, 21.867108141, 25.0, 19.3310145906107, 19.27964737200225] +surf_Q_in = 10.0 + +# sensible and latent heat flux parameter values +wspd = [ + 0.0, + 0.0, + 3919295.19947608, + 0.0, + 16118322.9664089, +] # unrealistically high to drive observable temperature change +patm = [0.0, 0.0, 954.680843658077, 0.0, 954.680843658077] +tatm = [11.0, 40.9752, 11.0, 11.0, 1.0] # used by lwr, shf, lhf +tatm = [t + DCTOK for t in tatm] +# shortwave radiation parameter values + +# unrealistically high to drive a 0.1 deg C rise in stream temperature +solr = [ + 4180000.0, + 0.0, + 0.0, + 0.0, + 0.0, +] +shd = [0.0, 0.0, 0.0, 0.0, 0.0] # 100% shade "turns off" solar flux +swrefl = [0.0, 1.0, 0.0, 1.0, 0.0] +rh = [0.0, 20.0, 1.0, 0.0, 0.01] # percent +lwrefl = 0.03 # Fogg et al 2023 + +# atmosphere composition adjustment factor (using dummy value to drive +# half a degree change) +atmc = [0.0, 9667.121567, 0.0, 0.0, 0.0] + +# latent heat flux parameter values (these values are specified in the options +# block and are therefore constant across reaches +c_d = 0.0 # Drag coefficient ($unitless$) +wf_slope = 1.383e-08 # wind function slope ($1/mbar$) +wf_int = 3.445e-09 # wind function intercept ($m/s$) +emiss_water = 0.95 # Fogg et al 2023 +emiss_riparian = 0.97 # Fogg et al 2023 +stephan_boltzmann = 5.670374419e-08 + +# Transport related parameters +porosity = sy # porosity (unitless) +K_therm = 2.0 # Thermal conductivity # ($W/m/C$) +rhow = 1000.0 # Density of water ($kg/m^3$) +rhos = 2650.0 # Density of the aquifer material ($kg/m^3$) +rhoa = 1.225 # Density of the atmosphere ($kg/m^3$) +Cpw = 4180.0 # Heat capacity of water ($J/kg/C$) +Cps = 880.0 # Heat capacity of the solids ($J/kg/C$) +Cpa = 717.0 # Heat capacity of the atmosphere ($J/kg/C$) +lhv = 2454000.0 # Latent heat of vaporization ($J/kg$) + +# Thermal conductivity of the streambed material ($W/m/C$) +K_therm_strmbed = 0.0 +rbthcnd = 0.0001 + +# time params +steady = {0: True, 1: False} +transient = {0: False, 1: True} +nstp = [1] +tsmult = [1] +perlen = [1] + +nouter, ninner = 1000, 300 +hclose, rclose, relax = 1e-10, 1e-10, 0.97 + +# +# MODFLOW 6 flopy GWF object +# + + +def build_models(idx, test): + # Base simulation and model name and workspace + ws = test.workspace + name = cases[idx] + + print(f"Building model...{name}") + + # generate names for each model + gwfname = "gwf-" + name + gwename = "gwe-" + 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=f"{gwfname}.ims", + ) + 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 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", "ALL")], + ) + + # Instantiate streamflow routing package + # Determine the middle row and store in rMid (account for 0-base) + rMid = 1 + # sfr data + roughness = 1e-10 + rbth = 0.1 + strmbd_hk = rhk + strm_up = 0.95 + strm_dn = 0.94 + # divide by 10 to further reduce slop + slope = 0.04 + nconn = 0 + ustrf = 1.0 + ndv = 0 + strm_incision = 0.05 + + # explicitly set connections + conns = [(0, -3), (1, -3), (2, -3), (3, 0, 1, 2, -4), (4, 3)] + nreaches = len(conns) + + packagedata = [] + for irch in range(nreaches): + ncon = len(conns[irch]) - 1 + rp = [ + irch, + (0, 0, 0), + rlen, + rwid, + slope, + top - strm_incision, + rbth, + strmbd_hk, + roughness, + ncon, + ustrf, + ndv, + ] + packagedata.append(rp) + + sfr_perioddata = {} + for t in np.arange(len(perlen)): + sfrbndx = [] + for i in np.arange(nreaches): + # only specify inflow for the first three reaches + if i < 3: + sfrbndx.append([i, "INFLOW", surf_Q_in]) + + sfr_perioddata.update({t: sfrbndx}) + + # Instantiate SFR observation points + sfr_obs = { + f"{gwfname}.sfr.obs.csv": [ + ("rch1_depth", "depth", 1), + ("rch1_outf", "ext-outflow", 1), + ("rch1_wetwidth", "wet-width", 1), + ("rch2_wetwidth", "wet-width", 2), + ("rch3_wetwidth", "wet-width", 3), + ("rch4_wetwidth", "wet-width", 4), + ("rch5_wetwidth", "wet-width", 5), + ("rch1_stg", "stage", 1), + ("rch2_stg", "stage", 2), + ("rch3_stg", "stage", 3), + ("rch4_stg", "stage", 4), + ("rch5_stg", "stage", 5), + ], + "digits": 15, + "print_input": True, + "filename": gwfname + ".sfr.obs", + } + + budpth = f"{gwfname}.sfr.cbc" + flopy.mf6.ModflowGwfsfr( + gwf, + save_flows=True, + print_stage=True, + print_flows=True, + print_input=True, + length_conversion=1.0, + time_conversion=1.0, + budget_filerecord=budpth, + mover=False, + nreaches=nreaches, + packagedata=packagedata, + connectiondata=conns, + perioddata=sfr_perioddata, + observations=sfr_obs, + pname="SFR", + filename=f"{gwfname}.sfr", + ) + + # -------------------------------------------------- + # Setup the GWE model for simulating heat transport + # -------------------------------------------------- + gwe = flopy.mf6.ModflowGwe(sim, modelname=gwename) + + # Instantiating solver for GWT + imsgwe = flopy.mf6.ModflowIms( + sim, + print_option="ALL", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="NONE", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=f"{rclose} strict", + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + filename=f"{gwename}.ims", + ) + sim.register_ims_package(imsgwe, [gwename]) + + # Instantiating DIS for GWE + flopy.mf6.ModflowGwedis( + gwe, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + pname="DIS", + filename=f"{gwename}.dis", + ) + + # Instantiate Mobile Storage and Transfer package + flopy.mf6.ModflowGweest( + gwe, + save_flows=True, + porosity=porosity, + heat_capacity_water=Cpw, + density_water=rhow, + latent_heat_vaporization=lhv, + heat_capacity_solid=Cps, + density_solid=rhos, + pname="EST", + filename=f"{gwename}.est", + ) + + # Instantiate Energy Transport Initial Conditions package + flopy.mf6.ModflowGweic(gwe, strt=strt_gw_temp) + + # Instantiate Advection package + flopy.mf6.ModflowGweadv(gwe, scheme="UPSTREAM") + + # Instantiate Dispersion package (also handles conduction) + flopy.mf6.ModflowGwecnd( + gwe, + xt3d_off=True, + ktw=0.5918, + kts=0.2700, + pname="CND", + filename=f"{gwename}.cnd", + ) + + # Instantiating MODFLOW 6 transport source-sink mixing package + # [b/c at least one boundary back is active (SFR), ssm must be on] + sourcerecarray = [[]] + flopy.mf6.ModflowGwessm(gwe, sources=sourcerecarray, filename=f"{gwename}.ssm") + + # Instantiate Streamflow Energy Transport package + sfepackagedata = [] + for irno in range(len(conns)): + t = (irno, strm_in_init[irno], K_therm_strmbed, rbthcnd) + sfepackagedata.append(t) + + sfeperioddata = [] + sfeperioddata.append((0, "INFLOW", strm_in_init[0])) + sfeperioddata.append((1, "INFLOW", strm_in_init[1])) + sfeperioddata.append((2, "INFLOW", strm_in_init[2])) + + # Instantiate SFE observation points + sfe_obs = { + f"{gwename}.sfe.obs.csv": [ + ("rch5_evap", "surfevap", 5), + ("rch1_outftemp", "temperature", 1), + ("rch2_outftemp", "temperature", 2), + ("rch3_outftemp", "temperature", 3), + ("rch4_outftemp", "temperature", 4), + ("rch5_outftemp", "temperature", 5), + ("rch1_outfener", "ext-outflow", 1), + ("rch2_outfener", "ext-outflow", 2), + ("rch3_outfener", "ext-outflow", 3), + ("rch4_outfener", "ext-outflow", 4), + ("rch5_outfener", "ext-outflow", 5), + ("rch1_swr", "swr", 1), + ("rch2_swr", "swr", 2), + ("rch3_swr", "swr", 3), + ("rch4_swr", "swr", 4), + ("rch5_swr", "swr", 5), + ("rch1_lwr", "lwr", 1), + ("rch2_lwr", "lwr", 2), + ("rch3_lwr", "lwr", 3), + ("rch4_lwr", "lwr", 4), + ("rch5_lwr", "lwr", 5), + ("rch1_lhf", "lhf", 1), + ("rch2_lhf", "lhf", 2), + ("rch3_lhf", "lhf", 3), + ("rch4_lhf", "lhf", 4), + ("rch5_lhf", "lhf", 5), + ("rch1_shf", "shf", 1), + ("rch2_shf", "shf", 2), + ("rch3_shf", "shf", 3), + ("rch4_shf", "shf", 4), + ("rch5_shf", "shf", 5), + ], + "digits": 20, + "print_input": True, + "filename": gwename + ".sfe.obs", + } + + abc_filename = f"{gwename}.sfe.abc" + sfe = flopy.mf6.modflow.ModflowGwesfe( + gwe, + boundnames=False, + save_flows=True, + print_input=False, + print_flows=True, + print_temperature=True, + temperature_filerecord=gwename + ".sfe.bin", + budget_filerecord=gwename + ".sfe.bud", + packagedata=sfepackagedata, + reachperioddata=sfeperioddata, + flow_package_name="SFR", + observations=sfe_obs, + pname="SFE", + filename=f"{gwename}.sfe", + ) + + # Abc utility + abc_spd = {} + for kper in range(len(nstp)): + spd = [] + for irno in range(len(conns)): + spd.append([irno, "WSPD", wspd[irno]]) + spd.append([irno, "TATM", tatm[irno]]) + spd.append([irno, "SOLR", solr[irno]]) + spd.append([irno, "SHD", shd[irno]]) + spd.append([irno, "SWREFL", swrefl[irno]]) + spd.append([irno, "RH", rh[irno]]) + spd.append([irno, "ATMC", atmc[irno]]) + spd.append([irno, "PATM", patm[irno]]) + + abc_spd[kper] = spd + + abc = flopy.mf6.ModflowUtlabc( + sfe, + print_input=True, + density_air=rhoa, + heat_capacity_air=Cpa, + drag_coefficient=c_d, + emissivity_water=emiss_water, + emissivity_canopy=emiss_riparian, + wind_func_slope=wf_slope, + wind_func_int=wf_int, + longwave_reflectance=lwrefl, + reachperioddata=abc_spd, + filename=abc_filename, + ) + + # Instantiate Output Control package for transport + flopy.mf6.ModflowGweoc( + gwe, + temperature_filerecord=f"{gwename}.ucn", + saverecord=[("TEMPERATURE", "ALL")], + temperatureprintrecord=[("COLUMNS", 3, "WIDTH", 20, "DIGITS", 8, "GENERAL")], + printrecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")], + filename=f"{gwename}.oc", + ) + + # Instantiate Gwf-Gwe Exchange package + flopy.mf6.ModflowGwfgwe( + sim, + exgtype="GWF6-GWE6", + exgmnamea=gwfname, + exgmnameb=gwename, + filename=f"{gwename}.gwfgwe", + ) + + return sim, None + + +# sim, dum = build_models(0, r"c:\temp\_shf00") +# sim.write_simulation() + + +def calc_lwr_ener_transfer(ifno, updated_strm_temp): + Ql_up = emiss_water * stephan_boltzmann * ((updated_strm_temp + DCTOK) ** 4) + e_s = 6.1275 * math.exp(17.2693882 * ((tatm[ifno] - DCTOK) / (tatm[ifno] - 35.86))) + e_a = (rh[ifno] / 100.0) * e_s + emiss_air = (1.24 * (e_a / tatm[ifno]) ** (1.0 / 7.0)) * atmc[ifno] + emiss_down = (1.0 - shd[ifno]) * emiss_air + shd[ + ifno + ] * emiss_riparian # calcs for emiss_riparian + + Ql_dn = emiss_down * stephan_boltzmann * (tatm[ifno] ** 4) + # calc net longwave + lwrflx = Ql_dn * (1.0 - lwrefl) - Ql_up + + return lwrflx + + +def calc_lhf_ener_transfer(ifno, updated_strm_temp, mf_strm_wid): + L = 2499.64 - (2.51 * (updated_strm_temp - DCTOK)) + e_w = 6.1275 * math.exp( + 17.2693882 * ((updated_strm_temp - DCTOK) / (updated_strm_temp - 35.86)) + ) + e_s = 6.1275 * math.exp(17.2693882 * ((tatm[ifno] - DCTOK) / (tatm[ifno] - 35.86))) + e_a = (rh[ifno] / 100) * e_s + vap_press_deficit = e_w - e_a + wind_function = wf_int + wf_slope * wspd[ifno] + Ev = wind_function * vap_press_deficit + lhf_ener_per_sqm = Ev * L * rhow + + ener_transfer = lhf_ener_per_sqm * rlen * mf_strm_wid + + return ener_transfer + + +def calc_shf_ener_transfer(ifno, lhflx, updated_strm_temp): + e_w = 6.1275 * math.exp( + 17.2693882 * ((updated_strm_temp - DCTOK) / (updated_strm_temp - 35.86)) + ) + e_s = 6.1275 * math.exp(17.2693882 * ((tatm[ifno] - DCTOK) / (tatm[ifno] - 35.86))) + e_a = (rh[ifno] / 100) * e_s + br = 0.00061 * patm[ifno] * (updated_strm_temp - tatm[ifno]) / (e_w - e_a) + + shflx = br * lhflx + + return shflx + + +def check_output(idx, test): + print("evaluating results...") + msg0 = "Stream channel width less than 1.0, should be 1.0 m" + msg1 = "SFE reach 1 should have roughly a 0.1 deg C rise from short wave" + msg2 = "SFE reach 2 should have roughly a 0.1 deg C rise from long wave" + msg3 = "SFE reach 3 not as expected" + msg4 = "SFE reach 4 not reflective of simple mixing of reaches 1-3" + msg5 = "SFE reach 5 evaporation rate from latent heat flux not 5 as expected" + + # read flow results from model + name = cases[idx] + gwfname = "gwf-" + name + gwename = "gwe-" + name + + # calc expected rise in temperature independent of mf6 + + fpth = os.path.join(test.workspace, gwfname + ".sfr.obs.csv") + assert os.path.isfile(fpth) + df = pd.read_csv(fpth) + mf_strm_wid = df.loc[0, "RCH1_WETWIDTH"].copy() + # confirm stream width is 1.0 m + assert np.isclose(mf_strm_wid, 1.0, atol=1e-9), msg0 + + # confirm that the stream temperature change fits with expectations + # reach 1 - solar radiation + fpth2 = os.path.join(test.workspace, gwename + ".sfe.obs.csv") + assert os.path.isfile(fpth2) + df2 = pd.read_csv(fpth2) + # input such that a roughly 0.1 degC rise should result + ifno = 0 + swr_ener_input = ( + (1.0 - shd[ifno]) * (1 - swrefl[ifno]) * solr[ifno] * mf_strm_wid * rlen + ) + rch1_Tdelt = swr_ener_input / (Cpw * rhow * surf_Q_in) + rch1Temp = df2.loc[0, "RCH1_OUTFTEMP"] + sim_Tdelt = rch1Temp - strm_in_init[ifno] + # atol needed to account for a small difference associated with a small + # amount of energy storage change in the reach + assert np.isclose(rch1_Tdelt, sim_Tdelt, atol=1e-4), msg1 + + # reach 2 - longwave radiation + chng = 1 + ifno = 1 # which is really 2 in 0-based! + strt_strm_temp = strm_in_init[ifno] + updated_strm_temp = strm_in_init[ifno] + mf_strm_wid = df.loc[0, "RCH2_WETWIDTH"].copy() + while chng > hclose: + lwr_ener_transfer = calc_lwr_ener_transfer(ifno, updated_strm_temp) + temp_change = lwr_ener_transfer / (surf_Q_in * Cpw * rhow) + updated_temp = strm_in_init[ifno] + temp_change + chng = abs(updated_strm_temp - updated_temp) + updated_strm_temp = updated_temp + + rch2Temp = df2.loc[0, "RCH2_OUTFTEMP"].copy() + sim_Tdelt = rch2Temp - strm_in_init[ifno] + rch2_Tdelt = updated_strm_temp - strm_in_init[ifno] + assert np.isclose(rch2_Tdelt, sim_Tdelt, atol=1e-4), msg2 + + # reach 3 - a focus on sensible heat flux (though other fluxes present) + chng = 1 + ifno = 2 # which is really 3 in 0-based! + strt_strm_temp = strm_in_init[ifno] + updated_strm_temp = strm_in_init[ifno] + mf_strm_wid = df.loc[0, "RCH3_WETWIDTH"].copy() + # latent calcs come first + while chng > hclose: + lwr_ener_transfer = calc_lwr_ener_transfer(ifno, updated_strm_temp) + lhf_ener_transfer = calc_lhf_ener_transfer( + ifno, updated_strm_temp + DCTOK, mf_strm_wid + ) + shf_ener_transfer = calc_shf_ener_transfer( + ifno, lhf_ener_transfer, updated_strm_temp + DCTOK + ) + temp_change = (lwr_ener_transfer - lhf_ener_transfer + shf_ener_transfer) / ( + surf_Q_in * Cpw * rhow + ) + updated_temp = strt_strm_temp + temp_change + chng = abs(updated_strm_temp - updated_temp) + updated_strm_temp = updated_temp + + rch3Temp = df2.loc[0, "RCH3_OUTFTEMP"].copy() + assert np.isclose(updated_strm_temp, rch3Temp), msg3 + + # reach 4 - simple mixing of reaches 1, 2, & 3, ABC boundary fluxes kept to a min + rch4Temp = df2.loc[0, "RCH4_OUTFTEMP"].copy() + rch4_expected_Temp = ( + rch1Temp * surf_Q_in + rch2Temp * surf_Q_in + rch3Temp * surf_Q_in + ) / (3 * surf_Q_in) + assert np.isclose(rch4Temp, rch4_expected_Temp, atol=1e-5), msg4 + + # reach 5 - a focus on latent heat flux (though other fluxes present) + # an expected evaporation rate of 5 mm/day + rch5_simEvap = df2.loc[0, "RCH5_EVAP"].copy() + assert np.isclose(rch5_simEvap, 5.0, atol=1e-8), msg5 + + +# - No need to change any code below +@pytest.mark.parametrize( + "idx, name", + list(enumerate(cases)), +) +def test_mf6model(idx, name, function_tmpdir, targets): + test = TestFramework( + name=name, + workspace=function_tmpdir, + targets=targets, + build=lambda t: build_models(idx, t), + check=lambda t: check_output(idx, t), + ) + test.run() diff --git a/autotest/test_gwe_sfe_abc03.py b/autotest/test_gwe_sfe_abc03.py new file mode 100644 index 00000000000..83bcac98a01 --- /dev/null +++ b/autotest/test_gwe_sfe_abc03.py @@ -0,0 +1,545 @@ +# Test the use of the atmospheric boundary condition utility used in conjunction +# with the SFE advanced package. This test is focused on the use of time series +# in setting up ABC input. There is no "check_output" routine, the test is +# simply that the time series input is read in and the model runs to completion. + + +import flopy +import numpy as np +import pytest +from framework import TestFramework + +cases = ["sfe-abc-ts"] + +DCTOK = 273.16 + +# Model units +length_units = "m" +time_units = "seconds" + +# time params +transient = {0: True, 1: True, 2: True, 3: True} +nstp = [1, 1, 1, 1] +tsmult = [1, 1, 1, 1] +perlen = [1, 1, 1, 1] + +# model domain and grid definition +nrow = 1 +ncol = 5 +nlay = 1 +delr = 10.0 +delc = 10.0 + +ibound = 1 +top = 1.0 +botm = 0.0 +strthd = 0.0 +chd_on = False + +# model input parameters +k11 = 500.0 +strt_gw_temp = 99.0 + +ss = 0.00001 +sy = 0.20 +hani = 1 +laytyp = 1 + +# Package boundary conditions +sfr_evaprate = 0.0 +rhk = 0.0 +rwid = 1.0 +strm_temp = 11.0 +surf_Q_in = 10.0 +# ABC parameter values +# wind speed +wspd = np.random.uniform(0.5, 5.0, size=[5, 5]) +wspd_cols = ["wnd_rch1", "wnd_rch2", "wnd_rch3", "wnd_rch4", "wnd_rch5"] +wspd_data = [] +for t in range(len(nstp) + 1): + wspd_data.append((t, wspd[t, 0], wspd[t, 1], wspd[t, 2], wspd[t, 3], wspd[t, 4])) + +# atmospheric temperature +tatm = np.random.uniform(10.0, 15.0, size=[5, 5]) +tatmK = tatm + DCTOK +tatm_cols = ["tatm_rch1", "tatm_rch2", "tatm_rch3", "tatm_rch4", "tatm_rch5"] +tatm_data = [] +for t in range(len(nstp) + 1): + tatm_data.append( + (t, tatmK[t, 0], tatmK[t, 1], tatmK[t, 2], tatmK[t, 3], tatmK[t, 4]) + ) + +# shortwave radiation parameter values +solr = np.random.uniform(800.0, 1200.0, size=[5, 5]) +solr_cols = ["solr_rch1", "solr_rch2", "solr_rch3", "solr_rch4", "solr_rch5"] +solr_data = [] +for t in range(len(nstp) + 1): + solr_data.append((t, solr[t, 0], solr[t, 1], solr[t, 2], solr[t, 3], solr[t, 4])) + +# shading from solar radiation +shd = np.random.uniform(0.3, 0.7, size=[5, 5]) +shd_cols = ["shd_rch1", "shd_rch2", "shd_rch3", "shd_rch4", "shd_rch5"] +shd_data = [] +for t in range(len(nstp) + 1): + shd_data.append((t, shd[t, 0], shd[t, 1], shd[t, 2], shd[t, 3], shd[t, 4])) + +# surface water reflectance +swrefl = np.random.uniform(0.03, 0.15, size=[5, 5]) +swrefl_cols = ["swr_rch1", "swr_rch2", "swr_rch3", "swr_rch4", "swr_rch5"] +swrefl_data = [] +for t in range(len(nstp) + 1): + swrefl_data.append( + (t, swrefl[t, 0], swrefl[t, 1], swrefl[t, 2], swrefl[t, 3], swrefl[t, 4]) + ) + +# relative humidity +rh = np.random.uniform(20.0, 60.0, size=[5, 5]) +rh_cols = ["rh_rch1", "rh_rch2", "rh_rch3", "rh_rch4", "rh_rch5"] +rh_data = [] +for t in range(len(nstp) + 1): + rh_data.append((t, rh[t, 0], rh[t, 1], rh[t, 2], rh[t, 3], rh[t, 4])) + +# Transport related parameters +porosity = sy # porosity (unitless) +K_therm = 2.0 # Thermal conductivity # ($W/m/C$) +rhow = 1000.0 # Density of water ($kg/m^3$) +rhos = 2650.0 # Density of the aquifer material ($kg/m^3$) +rhoa = 1.225 # Density of the atmosphere ($kg/m^3$) +Cpw = 4180.0 # Heat capacity of water ($J/kg/C$) +Cps = 880.0 # Heat capacity of the solids ($J/kg/C$) +Cpa = 717.0 # Heat capacity of the atmosphere ($J/kg/C$) +lhv = 2454000.0 # Latent heat of vaporization ($J/kg$) +c_d = 0.0 # Drag coefficient ($unitless$) !! +wf_slope = 1.383e-08 # wind function slope ($1/mbar$) +wf_int = 3.445e-09 # wind function intercept ($m/s$) +# Thermal conductivity of the streambed material ($W/m/C$) +K_therm_strmbed = 0.0 +rbthcnd = 0.0001 + + +nouter, ninner = 1000, 300 +hclose, rclose, relax = 1e-10, 1e-10, 0.97 + +# +# MODFLOW 6 flopy GWF object +# + + +def build_models(idx, test): + # Base simulation and model name and workspace + ws = test.workspace + name = cases[idx] + + print(f"Building model...{name}") + + # generate names for each model + gwfname = "gwf-" + name + gwename = "gwe-" + 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=f"{gwfname}.ims", + ) + 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, + transient=transient, + ) + + # Instantiate initial conditions package + flopy.mf6.ModflowGwfic(gwf, strt=strthd) + + # 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", "ALL")], + ) + + # Instantiate streamflow routing package + # Determine the middle row and store in rMid (account for 0-base) + rMid = 1 + # sfr data + nreaches = ncol + rlen = delr + roughness = 0.03 + rbth = 0.1 + strmbd_hk = rhk + strm_up = [up for up in np.arange(0.95, 0.91, -0.01)] + strm_dn = [dwn for dwn in np.arange(0.94, 0.90, -0.01)] + # divide by 10 to further reduce slop + slope = 0.001 + ustrf = 1.0 + ndv = 0 + strm_incision = 0.05 + + packagedata = [] + for irch in range(nreaches): + nconn = 1 + if 0 < irch < nreaches - 1: + nconn += 1 + rp = [ + irch, + (0, 0, 0), + rlen, + rwid, + slope, + top - strm_incision, + rbth, + strmbd_hk, + roughness, + nconn, + ustrf, + ndv, + ] + 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) + + sfrbndx = [0, "INFLOW", surf_Q_in] + sfr_perioddata = {0: sfrbndx} + + # Instantiate SFR observation points + sfr_obs = { + f"{gwfname}.sfr.obs.csv": [ + ("rch1_depth", "depth", 1), + ("rch1_outf", "ext-outflow", 1), + ("rch1_wetwidth", "wet-width", 1), + ], + "digits": 8, + "print_input": True, + "filename": gwfname + ".sfr.obs", + } + + budpth = f"{gwfname}.sfr.cbc" + flopy.mf6.ModflowGwfsfr( + gwf, + save_flows=True, + print_stage=True, + print_flows=True, + print_input=True, + length_conversion=1.0, + time_conversion=1.0, + budget_filerecord=budpth, + mover=False, + nreaches=nreaches, + packagedata=packagedata, + connectiondata=connectiondata, + perioddata=sfr_perioddata, + observations=sfr_obs, + pname="SFR", + filename=f"{gwfname}.sfr", + ) + + # -------------------------------------------------- + # Setup the GWE model for simulating heat transport + # -------------------------------------------------- + gwe = flopy.mf6.ModflowGwe(sim, modelname=gwename) + + # Instantiating solver for GWE + imsgwe = flopy.mf6.ModflowIms( + sim, + print_option="ALL", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="NONE", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=f"{rclose} strict", + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + filename=f"{gwename}.ims", + ) + sim.register_ims_package(imsgwe, [gwename]) + + # Instantiating DIS for GWE + flopy.mf6.ModflowGwedis( + gwe, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + pname="DIS", + filename=f"{gwename}.dis", + ) + + # Instantiate Energy Storage and Transfer package + flopy.mf6.ModflowGweest( + gwe, + save_flows=True, + porosity=porosity, + heat_capacity_water=Cpw, + density_water=rhow, + latent_heat_vaporization=lhv, + heat_capacity_solid=Cps, + density_solid=rhos, + pname="EST", + filename=f"{gwename}.est", + ) + + # Instantiate Energy Transport Initial Conditions package + flopy.mf6.ModflowGweic(gwe, strt=strt_gw_temp) + + # Instantiate Advection package + flopy.mf6.ModflowGweadv(gwe, scheme="UPSTREAM") + + # Instantiate Conduction package + flopy.mf6.ModflowGwecnd( + gwe, + xt3d_off=True, + ktw=0.5918, + kts=0.2700, + pname="CND", + filename=f"{gwename}.cnd", + ) + + # Instantiating MODFLOW 6 transport source-sink mixing package + # [b/c at least one boundary back is active (SFR), ssm must be on] + sourcerecarray = [[]] + flopy.mf6.ModflowGwessm(gwe, sources=sourcerecarray, filename=f"{gwename}.ssm") + + # Instantiate Streamflow Energy Transport package + sfepackagedata = [] + for irno in range(ncol): + t = (irno, strm_temp, K_therm_strmbed, rbthcnd) + sfepackagedata.append(t) + + sfeperioddata = [] + for irno in range(ncol): + if irno == 0: + sfeperioddata.append((irno, "INFLOW", strm_temp)) + + # Instantiate SFE observation points + sfe_obs = { + f"{gwename}.sfe.obs.csv": [ + ("rch1_outftemp", "temperature", 1), + ("rch1_outfener", "ext-outflow", 1), + ], + "digits": 12, + "print_input": True, + "filename": gwename + ".sfe.obs", + } + + abc_filename = f"{gwename}.sfe.abc" + sfe = flopy.mf6.modflow.ModflowGwesfe( + gwe, + boundnames=False, + save_flows=True, + print_input=False, + print_flows=True, + print_temperature=True, + temperature_filerecord=gwename + ".sfe.bin", + budget_filerecord=gwename + ".sfe.bud", + packagedata=sfepackagedata, + reachperioddata=sfeperioddata, + flow_package_name="SFR", + observations=sfe_obs, + pname="SFE", + filename=f"{gwename}.sfe", + ) + + # Abc utility + # use time series to specify boundary conditions + abc_spd = [] + abc_spd.append([0, "WSPD", "wnd_rch1"]) + abc_spd.append([0, "TATM", "tatm_rch1"]) + abc_spd.append([0, "SOLR", "solr_rch1"]) + abc_spd.append([0, "SHD", "shd_rch1"]) + abc_spd.append([0, "SWREFL", "swr_rch1"]) + abc_spd.append([0, "RH", "rh_rch1"]) + abc_spd.append([2, "WSPD", "wnd_rch3"]) + abc_spd.append([2, "TATM", "tatm_rch3"]) + abc_spd.append([2, "SOLR", "solr_rch3"]) + abc_spd.append([2, "SHD", "shd_rch3"]) + abc_spd.append([2, "SWREFL", "swr_rch3"]) + abc_spd.append([2, "RH", "rh_rch3"]) + abc_spd.append([1, "WSPD", "wnd_rch2"]) + abc_spd.append([1, "TATM", "tatm_rch2"]) + abc_spd.append([1, "SOLR", "solr_rch2"]) + abc_spd.append([1, "SHD", "shd_rch2"]) + abc_spd.append([1, "SWREFL", "swr_rch2"]) + abc_spd.append([1, "RH", "rh_rch2"]) + abc_spd.append([3, "WSPD", "wnd_rch4"]) + abc_spd.append([3, "TATM", "tatm_rch4"]) + abc_spd.append([3, "SOLR", "solr_rch4"]) + abc_spd.append([3, "SHD", "shd_rch4"]) + abc_spd.append([3, "SWREFL", "swr_rch4"]) + abc_spd.append([3, "RH", "rh_rch4"]) + abc_spd.append([4, "WSPD", "wnd_rch5"]) + abc_spd.append([4, "TATM", "tatm_rch5"]) + abc_spd.append([4, "SOLR", "solr_rch5"]) + abc_spd.append([4, "SHD", "shd_rch5"]) + abc_spd.append([4, "SWREFL", "swr_rch5"]) + abc_spd.append([4, "RH", "rh_rch5"]) + + abc = flopy.mf6.ModflowUtlabc( + sfe, + print_input=True, + density_air=rhoa, + heat_capacity_air=Cpa, + drag_coefficient=c_d, + wind_func_slope=wf_slope, + wind_func_int=wf_int, + reachperioddata=abc_spd, + filename=abc_filename, + ) + + fname1 = f"{gwename}.sfe.abc.wspd.ts" + abc.ts.initialize( + filename=fname1, + timeseries=wspd_data, # wspd_df.reset_index().to_numpy(), + time_series_namerecord=wspd_cols, # wspd_df.columns.tolist(), + interpolation_methodrecord=["linearend"] * len(wspd_cols), + ) + + fname2 = f"{gwename}.sfe.abc.tatm.ts" + abc.ts.append_package( + filename=fname2, + timeseries=tatm_data, # tatm_df.reset_index().to_numpy(), + time_series_namerecord=tatm_cols, # tatm_df.columns.tolist(), + interpolation_methodrecord=["linearend"] * len(tatm_cols), + ) + + fname3 = f"{gwename}.sfe.abc.solr.ts" + abc.ts.append_package( + filename=fname3, + timeseries=solr_data, # solr_df.reset_index().to_numpy(), + time_series_namerecord=solr_cols, # solr_df.columns.tolist(), + interpolation_methodrecord=["linearend"] * len(solr_cols), + ) + + fname4 = f"{gwename}.sfe.abc.shd.ts" + abc.ts.append_package( + filename=fname4, + timeseries=shd_data, # shd_df.reset_index().to_numpy(), + time_series_namerecord=shd_cols, # shd_df.columns.tolist(), + interpolation_methodrecord=["linearend"] * len(shd_cols), + ) + + fname5 = f"{gwename}.sfe.abc.swrefl.ts" + abc.ts.append_package( + filename=fname5, + timeseries=swrefl_data, # swrefl_df.reset_index().to_numpy(), + time_series_namerecord=swrefl_cols, # swrefl_df.columns.tolist(), + interpolation_methodrecord=["linearend"] * len(swrefl_cols), + ) + + fname6 = f"{gwename}.sfe.abc.rh.ts" + abc.ts.append_package( + filename=fname6, + timeseries=rh_data, # rh_df.reset_index().to_numpy(), + time_series_namerecord=rh_cols, # rh_df.columns.tolist(), + interpolation_methodrecord=["linearend"] * len(rh_cols), + ) + + # Instantiate Output Control package for transport + flopy.mf6.ModflowGweoc( + gwe, + temperature_filerecord=f"{gwename}.ucn", + saverecord=[("TEMPERATURE", "ALL")], + temperatureprintrecord=[("COLUMNS", 3, "WIDTH", 20, "DIGITS", 8, "GENERAL")], + printrecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")], + filename=f"{gwename}.oc", + ) + + # Instantiate Gwf-Gwe Exchange package + flopy.mf6.ModflowGwfgwe( + sim, + exgtype="GWF6-GWE6", + exgmnamea=gwfname, + exgmnameb=gwename, + filename=f"{gwename}.gwfgwe", + ) + + return sim, None + + +# - No need to change any code below +@pytest.mark.parametrize( + "idx, name", + list(enumerate(cases)), +) +def test_mf6model(idx, name, function_tmpdir, targets): + test = TestFramework( + name=name, + workspace=function_tmpdir, + targets=targets, + build=lambda t: build_models(idx, t), + xfail=[False], + ) + test.run() diff --git a/autotest/test_gwe_sfe_lhf00.py b/autotest/test_gwe_sfe_lhf00.py new file mode 100644 index 00000000000..fd6492fea66 --- /dev/null +++ b/autotest/test_gwe_sfe_lhf00.py @@ -0,0 +1,530 @@ +# Test the use of the atmospheric boundary condition utility used in conjunction with +# the SFE advanced package. This test is a single cell with a single reach. +# Channel flow characteristics are unrealistic: Manning's n is unrealistically +# low and slope is extremely high. These conditions result in an extremely high +# streamflow velocity that results in nearly all of the heat being added to the +# channel exiting at the outlet with very near negligle heat storage increases +# in the channel. This test only uses latent heat flux (lhf). +# The result is a -1 deg C change in temperature in the +# streamflow - an easy result to confirm in this test. + +import math +import os + +import flopy +import numpy as np +import pandas as pd +import pytest +from framework import TestFramework + +cases = ["sfe-abc"] + +DCTOK = 273.16 + +# Model units +length_units = "m" +time_units = "seconds" + +# model domain and grid definition + +nrow = 1 +ncol = 1 +nlay = 1 +delr = 1.0 +delc = 1.0 +xmax = ncol * delr +ymax = nrow * delc + +ibound = 1 +top = 1.0 +botm = 0.0 +strthd = 0.0 +chd_on = False + +# model input parameters +k11 = 500.0 +strt_gw_temp = 99.0 + +ss = 0.00001 +sy = 0.20 +hani = 1 +laytyp = 1 + +# Package boundary conditions +sfr_evaprate = 0.0 +rhk = 0.0 +rwid = 1.0 +strm_temp = 11.0 +surf_Q_in = [ + [10.0], +] +# sensible and latent heat flux parameter values +wspd = 126005.30 # unrealistically high to drive a -1C change +tatm = 5.0 + DCTOK +# shortwave radiation parameter values +solr = 47880870.9 # unrealistically high to drive a 1 deg C rise in stream temperature +shd = 1.0 # 100% shade "turns off" solar flux +swrefl = 0.03 +rh = 30.0 # percent + +# Transport related parameters +porosity = sy # porosity (unitless) +K_therm = 2.0 # Thermal conductivity # ($W/m/C$) +rhow = 1000.0 # Density of water ($kg/m^3$) +rhos = 2650.0 # Density of the aquifer material ($kg/m^3$) +rhoa = 1.225 # Density of the atmosphere ($kg/m^3$) +Cpw = 4180.0 # Heat capacity of water ($J/kg/C$) +Cps = 880.0 # Heat capacity of the solids ($J/kg/C$) +Cpa = 717.0 # Heat capacity of the atmosphere ($J/kg/C$) +lhv = 2454000.0 # Latent heat of vaporization ($J/kg$) +c_d = 0.0 # Drag coefficient ($unitless$) !! +wf_slope = 1.383e-08 # wind function slope ($1/mbar$) +wf_int = 3.445e-09 # wind function intercept ($m/s$) +# Thermal conductivity of the streambed material ($W/m/C$) +K_therm_strmbed = 0.0 +rbthcnd = 0.0001 + +# time params +steady = {0: True, 1: False} +transient = {0: False, 1: True} +nstp = [1] +tsmult = [1] +perlen = [1] + +nouter, ninner = 1000, 300 +hclose, rclose, relax = 1e-10, 1e-10, 0.97 + +# +# MODFLOW 6 flopy GWF object +# + + +def build_models(idx, test): + # Base simulation and model name and workspace + ws = test.workspace + name = cases[idx] + + print(f"Building model...{name}") + + # generate names for each model + gwfname = "gwf-" + name + gwename = "gwe-" + 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=f"{gwfname}.ims", + ) + 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 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", "ALL")], + ) + + # Instantiate streamflow routing package + # Determine the middle row and store in rMid (account for 0-base) + rMid = 1 + # sfr data + nreaches = ncol + rlen = delr + roughness = 1e-10 + rbth = 0.1 + strmbd_hk = rhk + strm_up = 0.95 + strm_dn = 0.94 + # divide by 10 to further reduce slop + slope = 0.04 + nconn = 0 + ustrf = 1.0 + ndv = 0 + strm_incision = 0.05 + + packagedata = [] + for irch in range(nreaches): + rp = [ + irch, + (0, 0, 0), + rlen, + rwid, + slope, + top - strm_incision, + rbth, + strmbd_hk, + roughness, + nconn, + ustrf, + ndv, + ] + packagedata.append(rp) + + connectiondata = [0] + + sfr_perioddata = {} + for t in np.arange(len(surf_Q_in[idx])): + sfrbndx = [] + for i in np.arange(nreaches): + if i == 0: + sfrbndx.append([i, "INFLOW", surf_Q_in[idx][t]]) + # sfrbndx.append([i, "EVAPORATION", sfr_evaprate]) + + sfr_perioddata.update({t: sfrbndx}) + + # Instantiate SFR observation points + sfr_obs = { + f"{gwfname}.sfr.obs.csv": [ + ("rch1_depth", "depth", 1), + ("rch1_outf", "ext-outflow", 1), + ("rch1_wetwidth", "wet-width", 1), + ], + "digits": 8, + "print_input": True, + "filename": gwfname + ".sfr.obs", + } + + budpth = f"{gwfname}.sfr.cbc" + flopy.mf6.ModflowGwfsfr( + gwf, + save_flows=True, + print_stage=True, + print_flows=True, + print_input=True, + length_conversion=1.0, + time_conversion=1.0, + budget_filerecord=budpth, + mover=False, + nreaches=nreaches, + packagedata=packagedata, + connectiondata=connectiondata, + perioddata=sfr_perioddata, + observations=sfr_obs, + pname="SFR", + filename=f"{gwfname}.sfr", + ) + + # -------------------------------------------------- + # Setup the GWE model for simulating heat transport + # -------------------------------------------------- + gwe = flopy.mf6.ModflowGwe(sim, modelname=gwename) + + # Instantiating solver for GWT + imsgwe = flopy.mf6.ModflowIms( + sim, + print_option="ALL", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="NONE", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=f"{rclose} strict", + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + filename=f"{gwename}.ims", + ) + sim.register_ims_package(imsgwe, [gwename]) + + # Instantiating DIS for GWE + flopy.mf6.ModflowGwedis( + gwe, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + pname="DIS", + filename=f"{gwename}.dis", + ) + + # Instantiate Mobile Storage and Transfer package + flopy.mf6.ModflowGweest( + gwe, + save_flows=True, + porosity=porosity, + heat_capacity_water=Cpw, + density_water=rhow, + latent_heat_vaporization=lhv, + heat_capacity_solid=Cps, + density_solid=rhos, + pname="EST", + filename=f"{gwename}.est", + ) + + # Instantiate Energy Transport Initial Conditions package + flopy.mf6.ModflowGweic(gwe, strt=strt_gw_temp) + + # Instantiate Advection package + flopy.mf6.ModflowGweadv(gwe, scheme="UPSTREAM") + + # Instantiate Dispersion package (also handles conduction) + flopy.mf6.ModflowGwecnd( + gwe, + xt3d_off=True, + ktw=0.5918, + kts=0.2700, + pname="CND", + filename=f"{gwename}.cnd", + ) + + # Instantiating MODFLOW 6 transport source-sink mixing package + # [b/c at least one boundary back is active (SFR), ssm must be on] + sourcerecarray = [[]] + flopy.mf6.ModflowGwessm(gwe, sources=sourcerecarray, filename=f"{gwename}.ssm") + + # Instantiate Streamflow Energy Transport package + sfepackagedata = [] + for irno in range(ncol): + t = (irno, strm_temp, K_therm_strmbed, rbthcnd) + sfepackagedata.append(t) + + sfeperioddata = [] + for irno in range(ncol): + if irno == 0: + sfeperioddata.append((irno, "INFLOW", strm_temp)) + + # Instantiate SFE observation points + sfe_obs = { + f"{gwename}.sfe.obs.csv": [ + ("rch1_outftemp", "temperature", 1), + ("rch1_outfener", "ext-outflow", 1), + ], + "digits": 8, + "print_input": True, + "filename": gwename + ".sfe.obs", + } + + abc_filename = f"{gwename}.sfe.abc" + sfe = flopy.mf6.modflow.ModflowGwesfe( + gwe, + boundnames=False, + save_flows=True, + print_input=False, + print_flows=True, + print_temperature=True, + temperature_filerecord=gwename + ".sfe.bin", + budget_filerecord=gwename + ".sfe.bud", + packagedata=sfepackagedata, + reachperioddata=sfeperioddata, + flow_package_name="SFR", + observations=sfe_obs, + pname="SFE", + filename=f"{gwename}.sfe", + ) + + # Abc utility + swr_optional_off = True + lwr_optional_off = True + lhf_optional_off = False + shf_optional_off = True + + abc_spd = {} + for kper in range(len(nstp)): + spd = [] + for irno in range(ncol): + spd.append([irno, "WSPD", wspd]) + spd.append([irno, "TATM", tatm]) + spd.append([irno, "SOLR", solr]) + spd.append([irno, "SHD", shd]) + spd.append([irno, "SWREFL", swrefl]) + spd.append([irno, "RH", rh]) + abc_spd[kper] = spd + + abc = flopy.mf6.ModflowUtlabc( + sfe, + print_input=True, + swr_off=swr_optional_off, + lwr_off=lwr_optional_off, + lhf_off=lhf_optional_off, + shf_off=shf_optional_off, + density_air=rhoa, + heat_capacity_air=Cpa, + drag_coefficient=c_d, + wind_func_slope=wf_slope, + wind_func_int=wf_int, + reachperioddata=abc_spd, + filename=abc_filename, + ) + + # Instantiate Output Control package for transport + flopy.mf6.ModflowGweoc( + gwe, + temperature_filerecord=f"{gwename}.ucn", + saverecord=[("TEMPERATURE", "ALL")], + temperatureprintrecord=[("COLUMNS", 3, "WIDTH", 20, "DIGITS", 8, "GENERAL")], + printrecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")], + filename=f"{gwename}.oc", + ) + + # Instantiate Gwf-Gwe Exchange package + flopy.mf6.ModflowGwfgwe( + sim, + exgtype="GWF6-GWE6", + exgmnamea=gwfname, + exgmnameb=gwename, + filename=f"{gwename}.gwfgwe", + ) + + return sim, None + + +# sim, dum = build_models(0, r"c:\temp\_shf00") +# sim.write_simulation() + + +def calc_ener_transfer(updated_strm_temp, mf_strm_wid): + L = 2499.64 - (2.51 * (updated_strm_temp - DCTOK)) + e_w = 6.1275 * math.exp( + 17.2693882 * ((updated_strm_temp - DCTOK) / (updated_strm_temp - 35.86)) + ) + e_s = 6.1275 * math.exp(17.2693882 * ((tatm - DCTOK) / (tatm - 35.86))) + e_a = (rh / 100) * e_s + vap_press_deficit = e_w - e_a + wind_function = wf_int + wf_slope * wspd + Ev = wind_function * vap_press_deficit + lhf_ener_per_sqm = Ev * L * rhow + + ener_transfer = lhf_ener_per_sqm * delr * mf_strm_wid + + return -ener_transfer + + +def check_output(idx, test): + print("evaluating results...") + msg0 = "Stream channel width less than 1.0, should be 1.0 m" + + # read flow results from model + name = cases[idx] + gwfname = "gwf-" + name + gwename = "gwe-" + name + + # calc expected rise in temperature independent of mf6 + + fpth = os.path.join(test.workspace, gwfname + ".sfr.obs.csv") + assert os.path.isfile(fpth) + df = pd.read_csv(fpth) + mf_strm_wid = df.loc[0, "RCH1_WETWIDTH"].copy() + # confirm stream width is 1.0 m + assert np.isclose(mf_strm_wid, 1.0, atol=1e-9), msg0 + + # confirm that the energy added to the stream results in a -1C change in temp + # temperature gradient + + tgrad = (tatm - DCTOK) - strm_temp + shf_ener_per_sqm = c_d * rhoa * Cpa * wspd * tgrad + swr_ener_per_sqm = solr * (1 - shd) * (1 - swrefl) + + # latent calcs + chng = 1 + strt_strm_temp = strm_temp + updated_strm_temp = strm_temp + while chng > hclose: + ener_transfer = calc_ener_transfer(updated_strm_temp + DCTOK, mf_strm_wid) + temp_change = ener_transfer / (surf_Q_in[idx][0] * Cpw * rhow) + updated_temp = strt_strm_temp + temp_change + chng = abs(updated_strm_temp - updated_temp) + updated_strm_temp = updated_temp + + fpth2 = os.path.join(test.workspace, gwename + ".sfe.obs.csv") + assert os.path.isfile(fpth2) + df2 = pd.read_csv(fpth2) + + # confirm 1 deg C decrease in temp + + msg1 = "Python temp change: " + str(temp_change) + msg2 = "MODFLOW temperature = " + str(df2.loc[0, "RCH1_OUTFTEMP"]) + msg3 = "MF temp change: " + str(strm_temp - df2.loc[0, "RCH1_OUTFTEMP"]) + + assert np.isclose( + df2.loc[0, "RCH1_OUTFTEMP"], strt_strm_temp + temp_change, atol=1e-6 + ), msg3 + ". " + msg1 + + +# - No need to change any code below +@pytest.mark.parametrize( + "idx, name", + list(enumerate(cases)), +) +def test_mf6model(idx, name, function_tmpdir, targets): + test = TestFramework( + name=name, + workspace=function_tmpdir, + targets=targets, + build=lambda t: build_models(idx, t), + check=lambda t: check_output(idx, t), + ) + test.run() diff --git a/autotest/test_gwe_sfe_lwr00.py b/autotest/test_gwe_sfe_lwr00.py new file mode 100644 index 00000000000..bae086d4b2b --- /dev/null +++ b/autotest/test_gwe_sfe_lwr00.py @@ -0,0 +1,544 @@ +# Test the use of the atmospheric boundary condition utility used in conjunction with +# the SFE advanced package. This test is a single cell with a single reach. +# Channel flow characteristics are unrealistic: Manning's n is unrealistically +# low and slope is extremely high. These conditions result in an extremely high +# streamflow velocity that results in nearly all of the heat being added to the +# channel exiting at the outlet with very near negligle heat storage increases +# in the channel. This test only uses longwave radiation heat flux (lwr). +# The result is a 1 deg C change in temperature in the +# streamflow - an easy result to confirm in this test. + +import math +import os + +import flopy +import numpy as np +import pandas as pd +import pytest +from framework import TestFramework + +cases = ["sfe-abc"] + +DCTOK = 273.16 + +# Model units +length_units = "m" +time_units = "seconds" + +# model domain and grid definition + +nrow = 1 +ncol = 1 +nlay = 1 +delr = 1.0 +delc = 1.0 +xmax = ncol * delr +ymax = nrow * delc + +ibound = 1 +top = 1.0 +botm = 0.0 +strthd = 0.0 +chd_on = False + +# model input parameters +k11 = 500.0 +strt_gw_temp = 99.0 + +ss = 0.00001 +sy = 0.20 +hani = 1 +laytyp = 1 + +# Package boundary conditions +sfr_evaprate = 0.0 +rhk = 0.0 +rwid = 1.0 +strm_temp = 11.0 +surf_Q_in = [ + [10.0], +] +# sensible and latent heat flux parameter values +wspd = 126005.30 # unrealistically high to drive a -1C change +# tatm = 5.0 +atmp = 100.0 +# shortwave radiation parameter values +solr = 47880870.9 # unrealistically high to drive a 1 deg C rise in stream temperature +shd = 1.0 # 100% shade "turns off" solar flux +swrefl = 0.1 # Fogg et al 2023 +rh = 30.0 # percent +# longwave radiation parameter values +lwrefl = 0.03 # Fogg et al 2023 +emiss_riparian = 0.97 # Fogg et al 2023 +emiss_water = 0.95 # Fogg et al 2023 +tatm = 686.8339 # unrealistically high atm temp that results in a 1C increase in +# stream temperature if lwr is only flux +atmc = 0.0 + + +# Transport related parameters +porosity = sy # porosity (unitless) +K_therm = 2.0 # Thermal conductivity # ($W/m/C$) +rhow = 1000.0 # Density of water ($kg/m^3$) +rhos = 2650.0 # Density of the aquifer material ($kg/m^3$) +rhoa = 1.225 # Density of the atmosphere ($kg/m^3$) +Cpw = 4180.0 # Heat capacity of water ($J/kg/C$) +Cps = 880.0 # Heat capacity of the solids ($J/kg/C$) +Cpa = 717.0 # Heat capacity of the atmosphere ($J/kg/C$) +lhv = 2454000.0 # Latent heat of vaporization ($J/kg$) +c_d = 0.0 # Drag coefficient ($unitless$) !! +wf_slope = 0.0 # 1.383e-08 # wind function slope ($1/mbar$) +## changing to 0 to "turn off" evap rate +wf_int = 0.0 # 3.445e-09 # wind function intercept ($m/s$) +## changing to 0 to "turn off" evap rate +# Thermal conductivity of the streambed material ($W/m/C$) +K_therm_strmbed = 0.0 +rbthcnd = 0.0001 + +# Constants +stephan_boltzmann = 5.670374419e-08 + +# time params +steady = {0: True, 1: False} +transient = {0: False, 1: True} +nstp = [1] +tsmult = [1] +perlen = [1] + +nouter, ninner = 1000, 300 +hclose, rclose, relax = 1e-10, 1e-10, 0.97 + +# +# MODFLOW 6 flopy GWF object +# + + +def build_models(idx, test): + # Base simulation and model name and workspace + ws = test.workspace + name = cases[idx] + + print(f"Building model...{name}") + + # generate names for each model + gwfname = "gwf-" + name + gwename = "gwe-" + 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=f"{gwfname}.ims", + ) + 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 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", "ALL")], + ) + + # Instantiate streamflow routing package + # Determine the middle row and store in rMid (account for 0-base) + rMid = 1 + # sfr data + nreaches = ncol + rlen = delr + roughness = 1e-10 + rbth = 0.1 + strmbd_hk = rhk + strm_up = 0.95 + strm_dn = 0.94 + # divide by 10 to further reduce slop + slope = 0.04 + nconn = 0 + ustrf = 1.0 + ndv = 0 + strm_incision = 0.05 + + packagedata = [] + for irch in range(nreaches): + rp = [ + irch, + (0, 0, 0), + rlen, + rwid, + slope, + top - strm_incision, + rbth, + strmbd_hk, + roughness, + nconn, + ustrf, + ndv, + ] + packagedata.append(rp) + + connectiondata = [0] + + sfr_perioddata = {} + for t in np.arange(len(surf_Q_in[idx])): + sfrbndx = [] + for i in np.arange(nreaches): + if i == 0: + sfrbndx.append([i, "INFLOW", surf_Q_in[idx][t]]) + # sfrbndx.append([i, "EVAPORATION", sfr_evaprate]) + + sfr_perioddata.update({t: sfrbndx}) + + # Instantiate SFR observation points + sfr_obs = { + f"{gwfname}.sfr.obs.csv": [ + ("rch1_depth", "depth", 1), + ("rch1_outf", "ext-outflow", 1), + ("rch1_wetwidth", "wet-width", 1), + ], + "digits": 8, + "print_input": True, + "filename": gwfname + ".sfr.obs", + } + + budpth = f"{gwfname}.sfr.cbc" + flopy.mf6.ModflowGwfsfr( + gwf, + save_flows=True, + print_stage=True, + print_flows=True, + print_input=True, + length_conversion=1.0, + time_conversion=1.0, + budget_filerecord=budpth, + mover=False, + nreaches=nreaches, + packagedata=packagedata, + connectiondata=connectiondata, + perioddata=sfr_perioddata, + observations=sfr_obs, + pname="SFR", + filename=f"{gwfname}.sfr", + ) + + # -------------------------------------------------- + # Setup the GWE model for simulating heat transport + # -------------------------------------------------- + gwe = flopy.mf6.ModflowGwe(sim, modelname=gwename) + + # Instantiating solver for GWT + imsgwe = flopy.mf6.ModflowIms( + sim, + print_option="ALL", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="NONE", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=f"{rclose} strict", + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + filename=f"{gwename}.ims", + ) + sim.register_ims_package(imsgwe, [gwename]) + + # Instantiating DIS for GWE + flopy.mf6.ModflowGwedis( + gwe, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + pname="DIS", + filename=f"{gwename}.dis", + ) + + # Instantiate Mobile Storage and Transfer package + flopy.mf6.ModflowGweest( + gwe, + save_flows=True, + porosity=porosity, + heat_capacity_water=Cpw, + density_water=rhow, + latent_heat_vaporization=lhv, + heat_capacity_solid=Cps, + density_solid=rhos, + pname="EST", + filename=f"{gwename}.est", + ) + + # Instantiate Energy Transport Initial Conditions package + flopy.mf6.ModflowGweic(gwe, strt=strt_gw_temp) + + # Instantiate Advection package + flopy.mf6.ModflowGweadv(gwe, scheme="UPSTREAM") + + # Instantiate Dispersion package (also handles conduction) + flopy.mf6.ModflowGwecnd( + gwe, + xt3d_off=True, + ktw=0.5918, + kts=0.2700, + pname="CND", + filename=f"{gwename}.cnd", + ) + + # Instantiating MODFLOW 6 transport source-sink mixing package + # [b/c at least one boundary back is active (SFR), ssm must be on] + sourcerecarray = [[]] + flopy.mf6.ModflowGwessm(gwe, sources=sourcerecarray, filename=f"{gwename}.ssm") + + # Instantiate Streamflow Energy Transport package + sfepackagedata = [] + for irno in range(ncol): + t = (irno, strm_temp, K_therm_strmbed, rbthcnd) + sfepackagedata.append(t) + + sfeperioddata = [] + for irno in range(ncol): + if irno == 0: + sfeperioddata.append((irno, "INFLOW", strm_temp)) + + # Instantiate SFE observation points + sfe_obs = { + f"{gwename}.sfe.obs.csv": [ + ("rch1_outftemp", "temperature", 1), + ("rch1_outfener", "ext-outflow", 1), + ], + "digits": 8, + "print_input": True, + "filename": gwename + ".sfe.obs", + } + + abc_filename = f"{gwename}.sfe.abc" + sfe = flopy.mf6.modflow.ModflowGwesfe( + gwe, + boundnames=False, + save_flows=True, + print_input=False, + print_flows=True, + print_temperature=True, + temperature_filerecord=gwename + ".sfe.bin", + budget_filerecord=gwename + ".sfe.bud", + packagedata=sfepackagedata, + reachperioddata=sfeperioddata, + flow_package_name="SFR", + observations=sfe_obs, + pname="SFE", + filename=f"{gwename}.sfe", + ) + + # Abc utility + abc_spd = {} + for kper in range(len(nstp)): + spd = [] + for irno in range(ncol): + spd.append([irno, "WSPD", wspd]) + spd.append([irno, "TATM", tatm]) + spd.append([irno, "SOLR", solr]) + spd.append([irno, "SHD", shd]) + spd.append([irno, "SWREFL", swrefl]) + spd.append([irno, "RH", rh]) + spd.append([irno, "ATMC", atmc]) + abc_spd[kper] = spd + + abc = flopy.mf6.ModflowUtlabc( + sfe, + print_input=True, + density_air=rhoa, + heat_capacity_air=Cpa, + reachperioddata=abc_spd, + filename=abc_filename, + longwave_reflectance=lwrefl, + emissivity_water=emiss_water, + emissivity_canopy=emiss_riparian, + swr_off=True, + lhf_off=True, + shf_off=True, + ) + + # Instantiate Output Control package for transport + flopy.mf6.ModflowGweoc( + gwe, + temperature_filerecord=f"{gwename}.ucn", + saverecord=[("TEMPERATURE", "ALL")], + temperatureprintrecord=[("COLUMNS", 3, "WIDTH", 20, "DIGITS", 8, "GENERAL")], + printrecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")], + filename=f"{gwename}.oc", + ) + + # Instantiate Gwf-Gwe Exchange package + flopy.mf6.ModflowGwfgwe( + sim, + exgtype="GWF6-GWE6", + exgmnamea=gwfname, + exgmnameb=gwename, + filename=f"{gwename}.gwfgwe", + ) + + return sim, None + + +# sim, dum = build_models(0, r"c:\temp\_shf00") +# sim.write_simulation() + +# Other energy transfers should equal 0 +# tgrad = tatm - strm_temp +# shf_ener_per_sqm = c_d * rhoa * Cpa * wspd * tgrad + + +def calc_ener_transfer(updated_strm_temp, mf_strm_wid): + # longwave + Ql_up = emiss_water * stephan_boltzmann * (updated_strm_temp**4) + + e_s = 6.1275 * math.exp(17.2693882 * ((tatm - DCTOK) / (tatm - 35.86))) + e_a = (rh / 100.0) * e_s + emiss_air = (1.24 * (e_a / tatm) ** (1.0 / 7.0)) * atmc # calcs to 0 + emiss_down = ( + 1.0 - shd + ) * emiss_air + shd * emiss_riparian # calcs to emiss_riparian + + Ql_down = emiss_down * stephan_boltzmann * (tatm**4) + + lwr_ener_per_sqm = Ql_down * (1.0 - lwrefl) - Ql_up + + ener_transfer = lwr_ener_per_sqm * delr * mf_strm_wid + + return ener_transfer + + +def check_output(idx, test): + print("evaluating results...") + msg0 = "Stream channel width less than 1.0, should be 1.0 m" + + # read flow results from model + name = cases[idx] + gwfname = "gwf-" + name + gwename = "gwe-" + name + + # calc expected rise in temperature independent of mf6 + + fpth = os.path.join(test.workspace, gwfname + ".sfr.obs.csv") + assert os.path.isfile(fpth) + df = pd.read_csv(fpth) + mf_strm_wid = df.loc[0, "RCH1_WETWIDTH"].copy() + # confirm stream width is 1.0 m + assert np.isclose(mf_strm_wid, 1.0, atol=1e-9), msg0 + + # confirm that the energy added to the stream results in a 1C change in temp + # temperature gradient + + # lw calcs + chng = 1 + strt_strm_temp = strm_temp + updated_strm_temp = strm_temp + while chng > hclose: + ener_transfer = calc_ener_transfer(updated_strm_temp + DCTOK, mf_strm_wid) + temp_change = ener_transfer / (surf_Q_in[idx][0] * Cpw * rhow) + updated_temp = strt_strm_temp + temp_change + chng = abs(updated_strm_temp - updated_temp) + updated_strm_temp = updated_temp + + fpth2 = os.path.join(test.workspace, gwename + ".sfe.obs.csv") + assert os.path.isfile(fpth2) + df2 = pd.read_csv(fpth2) + + # confirm 1 deg C decrease in temp + + msg1 = "Python temperature change is = " + str(temp_change) + msg2 = "MODFLOW temperature = " + str(df2.loc[0, "RCH1_OUTFTEMP"]) + msg3 = "MODFLOW temperature change is " + str( + strm_temp - df2.loc[0, "RCH1_OUTFTEMP"] + ) + + assert np.isclose( + df2.loc[0, "RCH1_OUTFTEMP"], strt_strm_temp + temp_change, atol=1e-6 + ), msg3 + ". " + msg1 # msg2 + ". " + + + +# - No need to change any code below +@pytest.mark.parametrize( + "idx, name", + list(enumerate(cases)), +) +def test_mf6model(idx, name, function_tmpdir, targets): + test = TestFramework( + name=name, + workspace=function_tmpdir, + targets=targets, + build=lambda t: build_models(idx, t), + check=lambda t: check_output(idx, t), + ) + test.run() diff --git a/autotest/test_gwe_sfe_shf00.py b/autotest/test_gwe_sfe_shf00.py new file mode 100644 index 00000000000..92fc1111f9c --- /dev/null +++ b/autotest/test_gwe_sfe_shf00.py @@ -0,0 +1,601 @@ +# Test the use of the sensible heat flux utility used in conjunction with the +# SFE advanced package. This test is a single cell with a single reach. +# After MODFLOW 6 completes the model run, latent and sensible heat fluxes are +# written to the sfe output file where results are pulled and compared to +# similar calculations herein. + +import math +import os + +import flopy +import numpy as np +import pandas as pd +import pytest +from framework import TestFramework + +cases = ["sfe-shf-opt1", "sfe-shf-opt2"] + +DCTOK = 273.16 + + +def calc_ea(rh, es): + # ambient atmospheric vapor pressure + ea = rh / 100.0 * es + return ea + + +def calc_e(temp): + # temperature must enter as degrees Kelvin + e = 6.1275 * math.exp(17.2693882 * ((temp - DCTOK) / (temp - 35.86))) + return e + + +def calc_lhv(tstrm, rhow): + # calculate latent heat of vaporization + l = 2499.64 - (2.51 * tstrm) + return l + + +def calc_evap(wfint, wfslp, wspd, ew, ea): + # calculate evaporation rate + evap = (wfint + (wfslp * wspd)) * (ew - ea) + return evap + + +def calc_bowen(patm, tstrmK, tatm, ew, ea): + # calculate bowen ratio + br = 0.00061 * patm * ((tstrmK - tatm) / (ew - ea)) + return br + + +def process_list_file(lstfile): + bud = None + # get the sfe budget items + with open(lstfile, "r") as f: + for line in f: + if srch_str in line: + # read two more lines + line = next(f) + line = next(f) + # process the next line + line = next(f) + budnames = line.split(" ") + budnames = [itm for itm in budnames if len(itm) > 0] + # skip a line + line = next(f) + # process another line + line = next(f) + vals = line.strip().split() + vals = [float(v) for v in vals] + # make a dataframe + bud = pd.DataFrame({"buditem": budnames, "value": vals}) + break + + # return the sfe budget items as a pandas dataframe + return bud + + +# Model units +length_units = "m" +time_units = "seconds" + +# model domain and grid definition + +nrow = 1 +ncol = 1 +nlay = 1 +delr = 1.0 +delc = 1.0 +xmax = ncol * delr +ymax = nrow * delc + +ibound = 1 +top = 1.0 +botm = 0.0 +strthd = 0.0 +chd_on = False + +# model input parameters +k11 = 500.0 +strt_gw_temp = 99.0 + +ss = 0.00001 +sy = 0.20 +hani = 1 +laytyp = 1 + +# Package boundary conditions +sfr_evaprate = 0.0 +rhk = 0.0 +rwid = 1.0 +strm_temp = 21.8671310408894 # deg C +surf_Q_in = [ + [10.0], + [10.0], # m^3/s +] + +# sensible heat flux parameter values +wfslp = 1.383e-8 +wfint = 3.445e-9 +wspd = 1.0 # m/s +patm = 954.680843658077 # mbar +tatm = 278.16 # deg K +rh = 20.0 # % (expressed as a percentage) + +# Transport related parameters +porosity = sy # porosity (unitless) +K_therm = 2.0 # Thermal conductivity # ($W/m/C$) +rhow = 1000.0 # Density of water ($kg/m^3$) +rhos = 2650.0 # Density of the aquifer material ($kg/m^3$) +rhoa = 1.225 # Density of the atmosphere ($kg/m^3$) +Cpw = 4180.0 # Heat capacity of water ($J/kg/C$) +Cps = 880.0 # Heat capacity of the solids ($J/kg/C$) +Cpa = 717.0 # Heat capacity of the atmosphere ($J/kg/C$) +lhv = 2454000.0 # Latent heat of vaporization ($J/kg$) +c_d = 0.002 # Drag coefficient ($unitless$) +# Thermal conductivity of the streambed material ($W/m/C$) +K_therm_strmbed = 0.0 +rbthcnd = 0.0001 + +# time params +steady = {0: True, 1: False} +transient = {0: False, 1: True} +nstp = [1] +tsmult = [1] +perlen = [1] + +nouter, ninner = 1000, 300 +hclose, rclose, relax = 1e-3, 1e-4, 0.97 + +srch_str = "SFE PACKAGE - SUMMARY OF FLOWS FOR EACH CONTROL VOLUME" + +# +# MODFLOW 6 flopy GWF object +# + + +def build_models(idx, test): + # Base simulation and model name and workspace + ws = test.workspace + name = cases[idx] + + print(f"Building model...{name}") + + # generate names for each model + gwfname = "gwf-" + name + gwename = "gwe-" + 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=f"{gwfname}.ims", + ) + 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 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", "ALL")], + ) + + # Instantiate streamflow routing package + # Determine the middle row and store in rMid (account for 0-base) + rMid = 1 + # sfr data + nreaches = ncol + rlen = delr + roughness = 1e-10 + rbth = 0.1 + strmbd_hk = rhk + strm_up = 0.95 + strm_dn = 0.94 + # divide by 10 to further reduce slop + slope = 0.04 + nconn = 0 + ustrf = 1.0 + ndv = 0 + strm_incision = 0.05 + + packagedata = [] + for irch in range(nreaches): + rp = [ + irch, + (0, 0, 0), + rlen, + rwid, + slope, + top - strm_incision, + rbth, + strmbd_hk, + roughness, + nconn, + ustrf, + ndv, + ] + packagedata.append(rp) + + connectiondata = [0] + + sfr_perioddata = {} + for t in np.arange(len(surf_Q_in[idx])): + sfrbndx = [] + for i in np.arange(nreaches): + if i == 0: + sfrbndx.append([i, "INFLOW", surf_Q_in[idx][t]]) + # sfrbndx.append([i, "EVAPORATION", sfr_evaprate]) + + sfr_perioddata.update({t: sfrbndx}) + + # Instantiate SFR observation points + sfr_obs = { + f"{gwfname}.sfr.obs.csv": [ + ("rch1_depth", "depth", 1), + ("rch1_outf", "ext-outflow", 1), + ("rch1_wetwidth", "wet-width", 1), + ], + "digits": 8, + "print_input": True, + "filename": gwfname + ".sfr.obs", + } + + budpth = f"{gwfname}.sfr.cbc" + flopy.mf6.ModflowGwfsfr( + gwf, + save_flows=True, + print_stage=True, + print_flows=True, + print_input=True, + length_conversion=1.0, + time_conversion=1.0, + budget_filerecord=budpth, + mover=False, + nreaches=nreaches, + packagedata=packagedata, + connectiondata=connectiondata, + perioddata=sfr_perioddata, + observations=sfr_obs, + pname="SFR", + filename=f"{gwfname}.sfr", + ) + + # -------------------------------------------------- + # Setup the GWE model for simulating heat transport + # -------------------------------------------------- + gwe = flopy.mf6.ModflowGwe(sim, modelname=gwename) + + # Instantiating solver for GWT + imsgwe = 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"{gwename}.ims", + ) + sim.register_ims_package(imsgwe, [gwename]) + + # Instantiating DIS for GWE + flopy.mf6.ModflowGwedis( + gwe, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + pname="DIS", + filename=f"{gwename}.dis", + ) + + # Instantiate Mobile Storage and Transfer package + flopy.mf6.ModflowGweest( + gwe, + save_flows=True, + porosity=porosity, + heat_capacity_water=Cpw, + density_water=rhow, + latent_heat_vaporization=lhv, + heat_capacity_solid=Cps, + density_solid=rhos, + pname="EST", + filename=f"{gwename}.est", + ) + + # Instantiate Energy Transport Initial Conditions package + flopy.mf6.ModflowGweic(gwe, strt=strt_gw_temp) + + # Instantiate Advection package + flopy.mf6.ModflowGweadv(gwe, scheme="UPSTREAM") + + # Instantiate Dispersion package (also handles conduction) + flopy.mf6.ModflowGwecnd( + gwe, + xt3d_off=True, + ktw=0.5918, + kts=0.2700, + pname="CND", + filename=f"{gwename}.cnd", + ) + + # Instantiating MODFLOW 6 transport source-sink mixing package + # [b/c at least one boundary back is active (SFR), ssm must be on] + sourcerecarray = [[]] + flopy.mf6.ModflowGwessm(gwe, sources=sourcerecarray, filename=f"{gwename}.ssm") + + # Instantiate Streamflow Energy Transport package + sfepackagedata = [] + for irno in range(ncol): + t = (irno, strm_temp, K_therm_strmbed, rbthcnd) + sfepackagedata.append(t) + + sfeperioddata = [] + for irno in range(ncol): + if irno == 0: + sfeperioddata.append((irno, "INFLOW", strm_temp)) + + # Instantiate SFE observation points + abc_obs = [] + if idx == 0: + abc_obs = [ + ("rch1_outftemp", "temperature", 1), + ("rch1_outfener", "ext-outflow", 1), + ("rch1_shf", "shf", 1), + ] + elif idx == 1: + abc_obs = [ + ("rch1_outftemp", "temperature", 1), + ("rch1_outfener", "ext-outflow", 1), + ("rch1_shf", "shf", 1), + ("rch1_lhf", "lhf", 1), + ] + + sfe_obs = { + f"{gwename}.sfe.obs.csv": abc_obs, + "digits": 12, + "print_input": True, + "filename": gwename + ".sfe.obs", + } + + abc_filename = f"{gwename}.sfe.abc" + sfe = flopy.mf6.modflow.ModflowGwesfe( + gwe, + boundnames=False, + save_flows=True, + print_input=False, + print_flows=True, + print_temperature=True, + temperature_filerecord=gwename + ".sfe.bin", + budget_filerecord=gwename + ".sfe.bud", + packagedata=sfepackagedata, + reachperioddata=sfeperioddata, + flow_package_name="SFR", + observations=sfe_obs, + pname="SFE", + filename=f"{gwename}.sfe", + ) + + # Abc utility + abc_spd = {} + for kper in range(len(nstp)): + spd = [] + for irno in range(ncol): + spd.append([irno, "WSPD", wspd]) + spd.append([irno, "TATM", tatm]) + spd.append([irno, "PATM", patm]) + spd.append([irno, "RH", rh]) + abc_spd[kper] = spd + + if idx == 0: + swr_optional_off = True + lwr_optional_off = True + lhf_optional_off = True + shf_optional_off = False + elif idx == 1: + swr_optional_off = True + lwr_optional_off = True + lhf_optional_off = False + shf_optional_off = False + + abc = flopy.mf6.ModflowUtlabc( + sfe, + print_input=True, + swr_off=swr_optional_off, + lwr_off=lwr_optional_off, + lhf_off=lhf_optional_off, + shf_off=shf_optional_off, + density_air=rhoa, + heat_capacity_air=Cpa, + drag_coefficient=c_d, + wind_func_slope=wfslp, + wind_func_int=wfint, + reachperioddata=abc_spd, + filename=abc_filename, + ) + + # Instantiate Output Control package for transport + flopy.mf6.ModflowGweoc( + gwe, + temperature_filerecord=f"{gwename}.ucn", + saverecord=[("TEMPERATURE", "ALL")], + temperatureprintrecord=[("COLUMNS", 3, "WIDTH", 20, "DIGITS", 8, "GENERAL")], + printrecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")], + filename=f"{gwename}.oc", + ) + + # Instantiate Gwf-Gwe Exchange package + flopy.mf6.ModflowGwfgwe( + sim, + exgtype="GWF6-GWE6", + exgmnamea=gwfname, + exgmnameb=gwename, + filename=f"{gwename}.gwfgwe", + ) + + return sim, None + + +def check_output(idx, test): + print("evaluating results...") + msg0 = "Stream channel width less than 1.0, should be 1.0 m" + msg1 = "Latent heat flux does not match externally calculated value." + msg2 = "Sensible heat flux does not match externally calculated value." + + # read flow results from model + name = cases[idx] + gwfname = "gwf-" + name + gwename = "gwe-" + name + + fpth = os.path.join(test.workspace, gwfname + ".sfr.obs.csv") + assert os.path.isfile(fpth) + df = pd.read_csv(fpth) + calc_strm_wid = df.loc[0, "RCH1_WETWIDTH"].copy() + # confirm stream width is 1.0 m + assert np.isclose(calc_strm_wid, 1.0, atol=1e-9), msg0 + + fpth2 = os.path.join(test.workspace, gwename + ".sfe.obs.csv") + assert os.path.isfile(fpth2) + df2 = pd.read_csv(fpth2) + + # calc expected sensible heat exchange + if idx == 0: + tgrad = tatm - (strm_temp + DCTOK) + ener_per_sqm = c_d * rhoa * Cpa * wspd * tgrad + shf = ener_per_sqm * (delr * calc_strm_wid) + + assert np.isclose(shf, df2.loc[0, "RCH1_SHF"].item(), atol=1e-9) + + elif idx > 0: + # read the outflows + fpth2 = os.path.join(test.workspace, gwfname + ".sfr.obs.csv") + assert os.path.isfile(fpth2) + df2 = pd.read_csv(fpth2) + + # read the outflow temperature + fpth3 = os.path.join(test.workspace, gwename + ".sfe.obs.csv") + assert os.path.isfile(fpth3) + df3 = pd.read_csv(fpth3) + + # get output stored in the gwe listing file + fname = gwename + ".lst" + lstfile = os.path.join(test.workspace, fname) + sfe_bud = process_list_file(lstfile) + + # calculate the amount of energy entering the reach + ener_in = surf_Q_in[idx][0] * strm_temp * Cpw * rhow + ener_out = ( + abs(df2.loc[0, "RCH1_OUTF"]) * df3.loc[0, "RCH1_OUTFTEMP"] * Cpw * rhow + ) + net_energy_in_out = ener_out - ener_in + + # calculate the latent and sensible heat fluxes locally and compare + # to the MF6 observations recorded in the corresponding output file + es = calc_e(tatm) + ea = calc_ea(rh, es) + ew = calc_e(df3.loc[0, "RCH1_OUTFTEMP"] + DCTOK) + # calculate latent heat of vaporization + lhv = calc_lhv(df3.loc[0, "RCH1_OUTFTEMP"], rhow) + # calculate evaporation/condensation rate + evap = calc_evap(wfint, wfslp, wspd, ew, ea) + # calculate latent heat flux + lhflx = evap * lhv * rhow + assert np.isclose(lhflx, df3.loc[0, "RCH1_LHF"], atol=1e-9), msg1 + # calculate Bowen ratio + br = calc_bowen(patm, df3.loc[0, "RCH1_OUTFTEMP"] + DCTOK, tatm, ew, ea) + # using the bowen ratio and latent heat flux, calc sensible heat flux + shflx = br * lhflx + assert np.isclose(shflx, df3.loc[0, "RCH1_SHF"], atol=1e-9), msg2 + + +# - No need to change any code below +@pytest.mark.parametrize( + "idx, name", + list(enumerate(cases)), +) +def test_mf6model(idx, name, function_tmpdir, targets): + test = TestFramework( + name=name, + workspace=function_tmpdir, + targets=targets, + build=lambda t: build_models(idx, t), + check=lambda t: check_output(idx, t), + ) + test.run() diff --git a/autotest/test_gwe_sfe_shf01.py b/autotest/test_gwe_sfe_shf01.py new file mode 100644 index 00000000000..df3f8a191b3 --- /dev/null +++ b/autotest/test_gwe_sfe_shf01.py @@ -0,0 +1,663 @@ +# Test the use of the atmospheric boundary condition utility, specifically the +# sensible heat flux calculations, used in conjunction with the +# SFE advanced package. This test checks to make sure that simulated sensible +# heat flux amounts are inline with what would be expected given changes in +# specific input parameters (while holding the others constant). Relative +# changes in stress periods 2 through 4 are compared back to stress period 1 +# or among the 3 reaches within each stress period + +# This test to include: +# - typical ways of specifying input +# - hot gw cell warming an upstream reach +# - thermally hot stream water warming host gw cells +# + + +import os + +import flopy +import numpy as np +import pandas as pd +import pytest +from framework import TestFramework + +cases = ["sfe-shf"] # , "sfe-shf-ts"] + +DCTOK = 273.16 + +# +# The last letter in the names above indicates the following +# n = "no gw/sw exchange" +# i = "gwf into strm" +# o = "strm to gw" +# m = "mixed" (i.e., convection one direction, conductive gradient the other direction?) + + +def get_x_frac(x_coord1, rwid): + x_xsec1 = [val / rwid for val in x_coord1] + return x_xsec1 + + +def get_xy_pts(x, y, rwid): + x_xsec1 = get_x_frac(x, rwid) + x_sec_tab = [[xx, hh] for xx, hh in zip(x_xsec1, y)] + return x_sec_tab + + +# Model units +length_units = "m" +time_units = "days" + +# model domain and grid definition +Lx = 90.0 +Ly = 90.0 +nrow = 3 +ncol = 3 +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 = np.array( + [ + [101.50, 101.25, 101.00], + [101.25, 101.00, 100.75], + [101.50, 101.25, 101.00], + ] +) + +botm = np.array( + [ + [98.5, 98.25, 98.0], + [98.25, 98.0, 97.75], + [98.5, 98.25, 98.0], + ] +) +strthd = 98.75 +chd_on = True + +# Boundary conditions +strt_gw_temp = [4.0, 4.0] +chd_condition = ["n", "i", "o", "m"] + +# NPF parameters +ss = 0.00001 +sy = 0.20 +hani = 1 +laytyp = 1 +k11 = 500.0 +# SFR/SFE +rhk = [0.0, k11] +strm_temp = [18.0, 18.0] +rlen = delr +surf_Q_in = [8.64, 86.4, 8.64, 8.64] # 86400 m^3/d = 1 m^3/s = 35.315 cfs +# SHF +# Remember, from the first to the second stress period, the flow rate increases +# Stress periods 3 and 4 return to the same flow rate as stress period 1 +# For stress period 3, increase wpd (everything else remains as is) +# For stress period 4, increase tatm (everything else remains as is) +wspd = [[1.0, 1.0, 1.0], [1.0, 1.0, 1.0], [6.0, 5.0, 4.0], [1.0, 1.0, 1.0]] +tatm = [ + [10.0, 10.0, 10.0], + [10.0, 10.0, 10.0], + [10.0, 10.0, 10.0], + [30.0, 25.0, 20.0], +] # deg C +tatmK = [[item + DCTOK for item in sublist] for sublist in tatm] + + +# Package boundary conditions +sfr_evaprate = 0.1 +rwid = [9.0, 10.0, 20] +# Channel geometry: trapezoidal +x_sec_tab1 = get_xy_pts( + [0.0, 2.0, 4.0, 5.0, 7.0, 9.0], + [0.66666667, 0.33333333, 0.0, 0.0, 0.33333333, 0.66666667], + rwid[0], +) + +x_sec_tab2 = get_xy_pts( + [0.0, 2.0, 4.0, 6.0, 8.0, 10.0], + [0.5, 0.25, 0.0, 0.0, 0.25, 0.5], + rwid[1], +) + +x_sec_tab3 = get_xy_pts( + [0.0, 4.0, 8.0, 12.0, 16.0, 20.0], + [0.33333333, 0.16666667, 0.0, 0.0, 0.16666667, 0.33333333], + rwid[2], +) +x_sec_tab = [x_sec_tab1, x_sec_tab2, x_sec_tab3] + + +# Transport related parameters +porosity = sy # porosity (unitless) +K_therm = 2.0 # Thermal conductivity # ($W/m/C$) +rhow = 1000 # Density of water ($kg/m^3$) +rhos = 2650 # Density of the aquifer material ($kg/m^3$) +Cpw = 4180 # Heat capacity of water ($J/kg/C$) +Cps = 880 # Heat capacity of the solids ($J/kg/C$) +lhv = 2454000.0 # Latent heat of vaporization ($J/kg$) +# Thermal conductivity of the streambed material ($W/m/C$) +# K_therm_strmbed = [1.5, 1.75, 2.0] +K_therm_strmbed = [0.0, 0.0, 0.0] +rbthcnd = 0.0001 + +# time params +steady = {0: False, 1: False} +transient = {0: True, 1: True} +nstp = [1, 1, 1, 1] +tsmult = [1, 1, 1, 1] +perlen = [1, 1, 1, 1] + +nouter, ninner = 1000, 300 +hclose, rclose, relax = 1e-3, 1e-4, 0.97 + +# +# MODFLOW 6 flopy GWF object +# + + +def build_models(idx, test): + # Base simulation and model name and workspace + ws = test.workspace + name = cases[idx] + + print(f"Building model...{name}") + + # generate names for each model + gwfname = "gwf-" + name + gwename = "gwe-" + 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=f"{gwfname}.ims", + ) + 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 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", "ALL")], + ) + + # Instantiate constant head boundary package + chdelev1 = top[0, 0] - 0.05 + chdelev2 = top[0, -1] - 0.05 + gw_temp = strt_gw_temp[idx] + if chd_on: + chdlist1 = [ + [(0, 0, 0), chdelev1, gw_temp], + [(0, nrow - 1, 0), chdelev1, gw_temp], + [(0, 0, ncol - 1), chdelev2, gw_temp], + [(0, nrow - 1, ncol - 1), chdelev2, gw_temp], + ] + flopy.mf6.ModflowGwfchd( + gwf, + stress_period_data=chdlist1, + print_input=True, + print_flows=True, + save_flows=False, + pname="CHD", + auxiliary="TEMPERATURE", + filename=f"{gwfname}.chd", + ) + + # Instantiate streamflow routing package + # Determine the middle row and store in rMid (account for 0-base) + rMid = 1 + # sfr data + nreaches = ncol + roughness = 0.035 + rbth = 1.0 + strmbd_hk = rhk[idx] + strm_up = 100.25 + strm_dn = 99 + # divide by 10 to further reduce slop + slope = (strm_up - strm_dn) / ((ncol - 1) * delr) / 10 + ustrf = 1.0 + ndv = 0 + strm_incision = 1.0 + + # use trapezoidal cross-section for channel geometry + sfr_xsec_tab_nm1 = f"{gwfname}.xsec.tab1" + sfr_xsec_tab_nm2 = f"{gwfname}.xsec.tab2" + sfr_xsec_tab_nm3 = f"{gwfname}.xsec.tab3" + sfr_xsec_tab_nm = [sfr_xsec_tab_nm1, sfr_xsec_tab_nm2, sfr_xsec_tab_nm3] + crosssections = [] + for n in range(nreaches): + # 3 reaches, 3 cross section types + crosssections.append([n, sfr_xsec_tab_nm[n]]) + + # Setup the tables + for n in range(len(x_sec_tab)): + flopy.mf6.ModflowUtlsfrtab( + gwf, + nrow=len(x_sec_tab[n]), + ncol=2, + table=x_sec_tab[n], + filename=sfr_xsec_tab_nm[n], + pname="sfrxsectable" + str(n + 1), + ) + + init_stgs = [] + packagedata = [] + for irch in range(nreaches): + nconn = 1 + if 0 < irch < nreaches - 1: + nconn += 1 + rp = [ + irch, + (0, rMid, irch), + rlen, + rwid[irch], + slope, + top[rMid, irch] - strm_incision, + rbth, + strmbd_hk, + roughness, + nconn, + ustrf, + ndv, + ] + packagedata.append(rp) + + init_stgs = [ + [0, 100.2533605521219], + [1, 100.0022143308538], + [2, 99.75145344348398], + ] + + 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) + + sfr_perioddata = {} + for t in np.arange(len(surf_Q_in)): + sfrbndx = [] + for i in np.arange(nreaches): + if i == 0: + sfrbndx.append([i, "INFLOW", surf_Q_in[t]]) + # sfrbndx.append([i, "EVAPORATION", sfr_evaprate]) + + sfr_perioddata.update({t: sfrbndx}) + + # Instantiate SFR observation points + sfr_obs = { + f"{gwfname}.sfr.obs.csv": [ + ("rch1_width", "wet-width", 1), + ("rch2_width", "wet-width", 2), + ("rch3_width", "wet-width", 3), + ("rch1_stg", "stage", 1), + ("rch2_stg", "stage", 2), + ("rch3_stg", "stage", 3), + ], + "digits": 20, + "print_input": True, + "filename": gwfname + ".sfr.obs", + } + + budpth = f"{gwfname}.sfr.cbc" + flopy.mf6.ModflowGwfsfr( + gwf, + storage=True, + save_flows=True, + print_stage=True, + print_flows=True, + print_input=True, + length_conversion=1.0, + time_conversion=86400, + budget_filerecord=budpth, + mover=False, + nreaches=nreaches, + packagedata=packagedata, + connectiondata=connectiondata, + crosssections=crosssections, + initialstages=init_stgs, + perioddata=sfr_perioddata, + observations=sfr_obs, + pname="SFR", + filename=f"{gwfname}.sfr", + ) + + # -------------------------------------------------- + # Setup the GWE model for simulating heat transport + # -------------------------------------------------- + gwe = flopy.mf6.ModflowGwe(sim, modelname=gwename) + + # Instantiating solver for GWT + imsgwe = 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"{gwename}.ims", + ) + sim.register_ims_package(imsgwe, [gwename]) + + # Instantiating DIS for GWE + flopy.mf6.ModflowGwedis( + gwe, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + pname="DIS", + filename=f"{gwename}.dis", + ) + + # Instantiate Mobile Storage and Transfer package + flopy.mf6.ModflowGweest( + gwe, + save_flows=True, + porosity=porosity, + heat_capacity_water=Cpw, + density_water=rhow, + latent_heat_vaporization=lhv, + heat_capacity_solid=Cps, + density_solid=rhos, + pname="EST", + filename=f"{gwename}.est", + ) + + # Instantiate Energy Transport Initial Conditions package + flopy.mf6.ModflowGweic(gwe, strt=strt_gw_temp[idx]) + + # Instantiate Advection package + flopy.mf6.ModflowGweadv(gwe, scheme="UPSTREAM") + + # Instantiate Dispersion package (also handles conduction) + flopy.mf6.ModflowGwecnd( + gwe, + xt3d_off=True, + ktw=0.5918, + kts=0.2700, + pname="CND", + filename=f"{gwename}.cnd", + ) + + # Instantiating MODFLOW 6 transport source-sink mixing package + # [b/c at least one boundary back is active (SFR), ssm must be on] + sourcerecarray = [("CHD", "AUX", "TEMPERATURE")] + flopy.mf6.ModflowGwessm(gwe, sources=sourcerecarray, filename=f"{gwename}.ssm") + + # Instantiate Streamflow Energy Transport package + sfepackagedata = [] + for irno in range(ncol): + t = (irno, strm_temp[idx], K_therm_strmbed[irno], rbthcnd) + sfepackagedata.append(t) + + sfeperioddata = [] + for irno in range(ncol): + if irno == 0: + sfeperioddata.append((irno, "INFLOW", strm_temp[idx])) + + # Instantiate SFE observation points + sfe_obs = { + f"{gwename}.sfe.obs.csv": [ + ("rch1_outftemp", "temperature", 1), + ("rch1_outfener", "ext-outflow", 1), + ("rch1_shf", "shf", 1), + ("rch2_outftemp", "temperature", 2), + ("rch2_outfener", "ext-outflow", 2), + ("rch2_shf", "shf", 2), + ("rch3_outftemp", "temperature", 3), + ("rch3_outfener", "ext-outflow", 3), + ("rch3_shf", "shf", 3), + ], + "digits": 8, + "print_input": True, + "filename": gwename + ".sfe.obs", + } + + abc_filename = f"{gwename}.sfe.abc" + sfe = flopy.mf6.modflow.ModflowGwesfe( + gwe, + boundnames=False, + save_flows=True, + print_input=False, + print_flows=False, + print_temperature=True, + temperature_filerecord=gwename + ".sfe.bin", + budget_filerecord=gwename + ".sfe.bud", + packagedata=sfepackagedata, + reachperioddata=sfeperioddata, + observations=sfe_obs, + flow_package_name="SFR", + pname="SFE", + filename=f"{gwename}.sfe", + ) + + # abc (only shf is active) + abc_spd = {} + for kper in range(len(nstp)): + spd = [] + for irno in range(ncol): + spd.append([irno, "WSPD", wspd[kper][irno]]) + spd.append([irno, "TATM", tatmK[kper][irno]]) + abc_spd[kper] = spd + + abc = flopy.mf6.ModflowUtlabc( + sfe, + print_input=True, + density_air=1.225, + heat_capacity_air=717.0, + drag_coefficient=0.002, + reachperioddata=abc_spd, + swr_off=True, + lwr_off=True, + lhf_off=True, + filename=abc_filename, + ) + + # Instantiate Output Control package for transport + flopy.mf6.ModflowGweoc( + gwe, + temperature_filerecord=f"{gwename}.ucn", + saverecord=[("TEMPERATURE", "ALL")], + temperatureprintrecord=[("COLUMNS", 3, "WIDTH", 20, "DIGITS", 8, "GENERAL")], + printrecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")], + filename=f"{gwename}.oc", + ) + + # Instantiate Gwf-Gwe Exchange package + flopy.mf6.ModflowGwfgwe( + sim, + exgtype="GWF6-GWE6", + exgmnamea=gwfname, + exgmnameb=gwename, + filename=f"{gwename}.gwfgwe", + ) + + return sim, None + + +def check_output(idx, test): + print("evaluating results...") + + # read flow results from model + name = cases[idx] + gwfname = "gwf-" + name + gwename = "gwe-" + name + + fname = gwfname + ".sfr.cbc" + fname = os.path.join(test.workspace, fname) + assert os.path.isfile(fname) + + sfrobj = flopy.utils.binaryfile.CellBudgetFile(fname, precision="double") + sfr_wetted_interface_area = sfrobj.get_data(text="gwf") + + # Retrieve simulated top width of each reach + sfr_pth0 = os.path.join(test.workspace, f"{gwfname}.sfr.obs.csv") + assert os.path.isfile(sfr_pth0) + sfroutdf = pd.read_csv(sfr_pth0) + + # assert that each successive reach grows in surface area + msg0 = "surface area not increasing from reach 1 to reach 2 as expected." + msg1 = "surface area not increasing from reach 2 to reach 3 as expected." + for i in np.arange(sfroutdf.shape[0]): + assert sfroutdf.loc[i, "RCH1_WIDTH"] < sfroutdf.loc[i, "RCH2_WIDTH"], msg0 + assert sfroutdf.loc[i, "RCH2_WIDTH"] < sfroutdf.loc[i, "RCH3_WIDTH"], msg1 + + # Retrieve select SFE output + sfe_pth0 = os.path.join(test.workspace, f"{gwename}.sfe.obs.csv") + assert os.path.isfile(sfe_pth0), "missing sfe observation output file" + sfeoutdf = pd.read_csv(sfe_pth0) + + # assert that each successive reach grows in sensible heat flux + # based on the problem setup, each stress period should reflect an increase + # in sensible heat flux relative to the first stress period for different + # reasons. + # From SP1 to SP2: inflow increase, therefore surface area increases + # (trapezoidal cross section) + # From SP1 to SP3: wind speed increases in each successive reach + # From SP1 to SP4: temperature of the atmosphere increases each successive + # reach + for i in np.arange(sfeoutdf.shape[0]): + assert abs(sfeoutdf.loc[i, "RCH1_SHF"]) > abs(sfeoutdf.loc[i, "RCH2_SHF"]), ( + "shf not increasing from reach 1 to 2 as expected in SP " + str(i + 1) + ) + assert abs(sfeoutdf.loc[i, "RCH2_SHF"]) > abs(sfeoutdf.loc[i, "RCH3_SHF"]), ( + "shf not increasing from reach 2 to 3 as expected in SP " + str(i + 1) + ) + + # as a result of the amount of energy leaving from shf, ensure that + # temperatures are changing accordingly in each successive reach + msg2 = "temperatures should decrease moving downstream. Something is amiss in SP " + for i in np.arange(sfeoutdf.shape[0]): + assert sfeoutdf.loc[i, "RCH1_OUTFTEMP"] > sfeoutdf.loc[i, "RCH2_OUTFTEMP"], ( + msg2 + str(i + 1) + ) + assert sfeoutdf.loc[i, "RCH2_OUTFTEMP"] > sfeoutdf.loc[i, "RCH3_OUTFTEMP"], ( + msg2 + str(i + 1) + ) + + # there should be no "external" energy outflow for the first or second + # reaches + msg3 = ( + "should not be external energy outflow from reach 1 in any of the " + "stress periods" + ) + msg4 = ( + "should not be external energy outflow from reach 2 in any of the " + "stress periods" + ) + assert sfeoutdf.loc[:, "RCH1_OUTFENER"].sum() == 0, msg3 + assert sfeoutdf.loc[:, "RCH2_OUTFENER"].sum() == 0, msg4 + + # there should be non-zero energy exported from the final reach in every + # stress period + msg5 = "export of energy from last reach should be negative" + for i in np.arange(sfeoutdf.shape[0]): + assert sfeoutdf.loc[i, "RCH3_OUTFENER"] < 0, msg5 + + +# - No need to change any code below +@pytest.mark.parametrize( + "idx, name", + list(enumerate(cases)), +) +def test_mf6model(idx, name, function_tmpdir, targets): + test = TestFramework( + name=name, + workspace=function_tmpdir, + targets=targets, + build=lambda t: build_models(idx, t), + check=lambda t: check_output(idx, t), + ) + test.run() diff --git a/autotest/test_gwe_sfe_swr00.py b/autotest/test_gwe_sfe_swr00.py new file mode 100644 index 00000000000..4796d019cb5 --- /dev/null +++ b/autotest/test_gwe_sfe_swr00.py @@ -0,0 +1,481 @@ +# Test the use of the shortwave radiation utility used in conjunction with the +# SFE advanced package. This test is a single cell with a single reach. +# Channel flow characteristics are unrealistic: Manning's n is unrealistically +# low and slope is extremely high. These conditions result in an extremely high +# streamflow velocity that results in nearly all of the heat being added to the +# channel exiting at the outlet with very near negligle heat storage increases +# in the channel. The result is a 1 deg C rise in temperature in the +# streamflow - an easy result to confirm in this test. + +import os + +import flopy +import numpy as np +import pandas as pd +import pytest +from framework import TestFramework + +cases = ["sfe-swr"] + +# Model units +length_units = "m" +time_units = "seconds" + +# model domain and grid definition + +nrow = 1 +ncol = 1 +nlay = 1 +delr = 1.0 +delc = 1.0 +xmax = ncol * delr +ymax = nrow * delc + +ibound = 1 +top = 1.0 +botm = 0.0 +strthd = 0.0 +chd_on = False + +# model input parameters +k11 = 500.0 +strt_gw_temp = 99.0 + +ss = 0.00001 +sy = 0.20 +hani = 1 +laytyp = 1 + +# Package boundary conditions +sfr_evaprate = 0.0 +rhk = 0.0 +rwid = 1.0 +strm_temp = 1.0 +surf_Q_in = [ + [10.0], +] +# shortwave radiation parameter values +solr = 47880870.9 # unrealistically high to drive a 1 deg C rise in stream temperature +shd = 0.1 +swrefl = 0.03 + +# Transport related parameters +porosity = sy # porosity (unitless) +K_therm = 2.0 # Thermal conductivity # ($W/m/C$) +rhow = 1000 # Density of water ($kg/m^3$) +rhos = 2650 # Density of the aquifer material ($kg/m^3$) +rhoa = 1.225 # Density of the atmosphere ($kg/m^3$) +Cpw = 4180 # Heat capacity of water ($J/kg/C$) +Cps = 880 # Heat capacity of the solids ($J/kg/C$) +Cpa = 717.0 # Heat capacity of the atmosphere ($J/kg/C$) +lhv = 2454000.0 # Latent heat of vaporization ($J/kg$) +# c_d = 0.002 # Drag coefficient ($unitless$) +# Thermal conductivity of the streambed material ($W/m/C$) +K_therm_strmbed = 0.0 +rbthcnd = 0.0001 + +# time params +steady = {0: True, 1: False} +transient = {0: False, 1: True} +nstp = [1] +tsmult = [1] +perlen = [1] + +nouter, ninner = 1000, 300 +hclose, rclose, relax = 1e-3, 1e-4, 0.97 + +# +# MODFLOW 6 flopy GWF object +# + + +def build_models(idx, test): + # Base simulation and model name and workspace + ws = test.workspace + name = cases[idx] + + print(f"Building model...{name}") + + # generate names for each model + gwfname = "gwf-" + name + gwename = "gwe-" + 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=f"{gwfname}.ims", + ) + 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 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", "ALL")], + ) + + # Instantiate streamflow routing package + # Determine the middle row and store in rMid (account for 0-base) + rMid = 1 + # sfr data + nreaches = ncol + rlen = delr + roughness = 1e-10 + rbth = 0.1 + strmbd_hk = rhk + strm_up = 0.95 + strm_dn = 0.94 + # divide by 10 to further reduce slop + slope = 0.04 + nconn = 0 + ustrf = 1.0 + ndv = 0 + strm_incision = 0.05 + + packagedata = [] + for irch in range(nreaches): + rp = [ + irch, + (0, 0, 0), + rlen, + rwid, + slope, + top - strm_incision, + rbth, + strmbd_hk, + roughness, + nconn, + ustrf, + ndv, + ] + packagedata.append(rp) + + connectiondata = [0] + + sfr_perioddata = {} + for t in np.arange(len(surf_Q_in[idx])): + sfrbndx = [] + for i in np.arange(nreaches): + if i == 0: + sfrbndx.append([i, "INFLOW", surf_Q_in[idx][t]]) + # sfrbndx.append([i, "EVAPORATION", sfr_evaprate]) + + sfr_perioddata.update({t: sfrbndx}) + + # Instantiate SFR observation points + sfr_obs = { + f"{gwfname}.sfr.obs.csv": [ + ("rch1_depth", "depth", 1), + ("rch1_outf", "ext-outflow", 1), + ("rch1_wetwidth", "wet-width", 1), + ], + "digits": 8, + "print_input": True, + "filename": gwfname + ".sfr.obs", + } + + budpth = f"{gwfname}.sfr.cbc" + flopy.mf6.ModflowGwfsfr( + gwf, + save_flows=True, + print_stage=True, + print_flows=True, + print_input=True, + length_conversion=1.0, + time_conversion=1.0, + budget_filerecord=budpth, + mover=False, + nreaches=nreaches, + packagedata=packagedata, + connectiondata=connectiondata, + perioddata=sfr_perioddata, + observations=sfr_obs, + pname="SFR", + filename=f"{gwfname}.sfr", + ) + + # -------------------------------------------------- + # Setup the GWE model for simulating heat transport + # -------------------------------------------------- + gwe = flopy.mf6.ModflowGwe(sim, modelname=gwename) + + # Instantiating solver for GWT + imsgwe = 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"{gwename}.ims", + ) + sim.register_ims_package(imsgwe, [gwename]) + + # Instantiating DIS for GWE + flopy.mf6.ModflowGwedis( + gwe, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + pname="DIS", + filename=f"{gwename}.dis", + ) + + # Instantiate Mobile Storage and Transfer package + flopy.mf6.ModflowGweest( + gwe, + save_flows=True, + porosity=porosity, + heat_capacity_water=Cpw, + density_water=rhow, + latent_heat_vaporization=lhv, + heat_capacity_solid=Cps, + density_solid=rhos, + pname="EST", + filename=f"{gwename}.est", + ) + + # Instantiate Energy Transport Initial Conditions package + flopy.mf6.ModflowGweic(gwe, strt=strt_gw_temp) + + # Instantiate Advection package + flopy.mf6.ModflowGweadv(gwe, scheme="UPSTREAM") + + # Instantiate Dispersion package (also handles conduction) + flopy.mf6.ModflowGwecnd( + gwe, + xt3d_off=True, + ktw=0.5918, + kts=0.2700, + pname="CND", + filename=f"{gwename}.cnd", + ) + + # Instantiating MODFLOW 6 transport source-sink mixing package + # [b/c at least one boundary back is active (SFR), ssm must be on] + sourcerecarray = [[]] + flopy.mf6.ModflowGwessm(gwe, sources=sourcerecarray, filename=f"{gwename}.ssm") + + # Instantiate Streamflow Energy Transport package + sfepackagedata = [] + for irno in range(ncol): + t = (irno, strm_temp, K_therm_strmbed, rbthcnd) + sfepackagedata.append(t) + + sfeperioddata = [] + for irno in range(ncol): + if irno == 0: + sfeperioddata.append((irno, "INFLOW", strm_temp)) + + # Instantiate SFE observation points + sfe_obs = { + f"{gwename}.sfe.obs.csv": [ + ("rch1_outftemp", "temperature", 1), + ("rch1_outfener", "ext-outflow", 1), + ], + "digits": 8, + "print_input": True, + "filename": gwename + ".sfe.obs", + } + + abc_filename = f"{gwename}.sfe.abc" + sfe = flopy.mf6.modflow.ModflowGwesfe( + gwe, + boundnames=False, + save_flows=True, + print_input=False, + print_flows=True, + print_temperature=True, + temperature_filerecord=gwename + ".sfe.bin", + budget_filerecord=gwename + ".sfe.bud", + packagedata=sfepackagedata, + reachperioddata=sfeperioddata, + flow_package_name="SFR", + observations=sfe_obs, + pname="SFE", + filename=f"{gwename}.sfe", + ) + + # Abc utility + abc_spd = {} + for kper in range(len(nstp)): + spd = [] + for irno in range(ncol): + spd.append([irno, "SOLR", solr]) + spd.append([irno, "SHD", shd]) + spd.append([irno, "SWREFL", swrefl]) + abc_spd[kper] = spd + + abc = flopy.mf6.ModflowUtlabc( + sfe, + lwr_off=True, + lhf_off=True, + shf_off=True, + print_input=True, + reachperioddata=abc_spd, + filename=abc_filename, + ) + + # Instantiate Output Control package for transport + flopy.mf6.ModflowGweoc( + gwe, + temperature_filerecord=f"{gwename}.ucn", + saverecord=[("TEMPERATURE", "ALL")], + temperatureprintrecord=[("COLUMNS", 3, "WIDTH", 20, "DIGITS", 8, "GENERAL")], + printrecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")], + filename=f"{gwename}.oc", + ) + + # Instantiate Gwf-Gwe Exchange package + flopy.mf6.ModflowGwfgwe( + sim, + exgtype="GWF6-GWE6", + exgmnamea=gwfname, + exgmnameb=gwename, + filename=f"{gwename}.gwfgwe", + ) + + return sim, None + + +# sim, dum = build_models(0, r"c:\temp\_swr00") +# sim.write_simulation() + + +def check_output(idx, test): + print("evaluating results...") + msg0 = "Stream channel width less than 1.0, should be 1.0 m" + + # read flow results from model + name = cases[idx] + gwfname = "gwf-" + name + gwename = "gwe-" + name + + # calc expected rise in temperature independent of mf6 + + fpth = os.path.join(test.workspace, gwfname + ".sfr.obs.csv") + assert os.path.isfile(fpth) + df = pd.read_csv(fpth) + calc_strm_wid = df.loc[0, "RCH1_WETWIDTH"].copy() + # confirm stream width is 1.0 m + assert np.isclose(calc_strm_wid, 1.0, atol=1e-7), msg0 + + # confirm that the energy added to the stream results in a 1 deg C rise in temp + # temperature gradient + # tgrad = tatm - strm_temp + ener_per_sqm = solr * (1 - shd) * (1 - swrefl) + ener_transfer = ener_per_sqm * (delr * calc_strm_wid) + # calculate expected temperature rise based on energy transfer + temp_rise = ener_transfer / (surf_Q_in[idx][0] * Cpw * rhow) + + fpth = os.path.join(test.workspace, gwename + ".sfe.obs.csv") + assert os.path.isfile(fpth) + df = pd.read_csv(fpth) + + # confirm 1 deg C rise in temp + msg1 = ( + "The MF6 simulated rise in river temperature does not match \ + external calculations. The calculated temperature rise \ + is: " + + str(temp_rise) + ) + + assert np.isclose(df.loc[0, "RCH1_OUTFTEMP"], strm_temp + temp_rise, atol=1e-5), ( + msg1 + " " + str(df.loc[0, "RCH1_OUTFTEMP"]) + ) + + +# - No need to change any code below +@pytest.mark.parametrize( + "idx, name", + list(enumerate(cases)), +) +def test_mf6model(idx, name, function_tmpdir, targets): + test = TestFramework( + name=name, + workspace=function_tmpdir, + targets=targets, + build=lambda t: build_models(idx, t), + check=lambda t: check_output(idx, t), + ) + test.run() diff --git a/doc/Common/gwe-sfeobs.tex b/doc/Common/gwe-sfeobs.tex index b6b602d8049..aad15d6df89 100644 --- a/doc/Common/gwe-sfeobs.tex +++ b/doc/Common/gwe-sfeobs.tex @@ -14,4 +14,9 @@ SFE & runoff & rno or boundname & -- & Runoff rate applied to a reach or group of reaches multiplied by the runoff temperature. \\ SFE & ext-inflow & rno or boundname & -- & Energy inflow into a reach or group of reaches calculated as the external inflow rate multiplied by the inflow temperature. \\ SFE & ext-outflow & rno or boundname & -- & External outflow from a reach or group of reaches to an external boundary. If boundname is not specified for ID, then the external outflow from a specific reach is observed. In this case, ID is the reach rno. \\ -SFE & strmbd-cond & rno or boundname & -- & Amount of heat conductively exchanged with the streambed material. +SFE & strmbd-cond & rno or boundname & -- & Amount of heat conductively exchanged with the streambed material. \\ +SFE & swr & rno or boundname & -- & Amount of energy entering the stream reach from shortwave radiation. \\ +SFE & lwr & rno or boundname & -- & Net amount of longwave radiation energy exchanged with the stream reach. \\ +SFE & lhf & rno or boundname & -- & Amount of energy lost from the stream reach owing to evaporative cooling. \\ +SFE & shf & rno or boundname & -- & Amount of sensible heat exchanged with the atmosphere. \\ +SFE & surfevap & rno & -- & Evaporation rate for a given stream reach. \\ diff --git a/doc/mf6io/gwe/abc.tex b/doc/mf6io/gwe/abc.tex new file mode 100644 index 00000000000..25b255b7e18 --- /dev/null +++ b/doc/mf6io/gwe/abc.tex @@ -0,0 +1,23 @@ +Input to the Atmospheric Boundary Condition (ABC) Package is read from the file that is specified in the ABC6 record of the OPTIONS block in the SFE package. There can only be one ABC utility activated for each instance of SFE. The ABC utility cannot be used without a corresponding SFE Package. + +The ABC utility does not have a dimensions block; instead, dimensions for the ABC utility are set using the dimensions from the corresponding SFE Package. For example, an SFR Package corresponding to an instance of SFE requires specification of the number of reaches (NREACHES). ABC sets the number of reaches equal to NREACHES. Therefore, the PERIOD block below cannot contain boundary information for a reach with an index greater than NREACHES. + +\vspace{5mm} +\subsubsection{Structure of Blocks} +\vspace{5mm} + +\noindent \textit{FOR EACH SIMULATION} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/utl-abc-options.dat} +\vspace{5mm} +\noindent \textit{FOR ANY STRESS PERIOD} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/utl-abc-period.dat} + +\vspace{5mm} +\subsubsection{Explanation of Variables} +\begin{description} +\input{./mf6ivar/tex/utl-abc-desc.tex} +\end{description} + +\vspace{5mm} +\subsubsection{Example Input File} +\lstinputlisting[style=inputfile]{./mf6ivar/examples/utl-abc-example.dat} diff --git a/doc/mf6io/gwe/gwe.tex b/doc/mf6io/gwe/gwe.tex index 9538b3a385e..cd5cb0e2079 100644 --- a/doc/mf6io/gwe/gwe.tex +++ b/doc/mf6io/gwe/gwe.tex @@ -129,6 +129,10 @@ \subsection{Energy Source Loading (ESL) Package} \subsection{Streamflow Energy Transport (SFE) Package} \input{gwe/sfe} +\newpage +\subsection{Atmospheric Boundary Conditions (ABC) Package} +\input{gwe/abc} + \newpage \subsection{Lake Energy Transport (LKE) Package} \input{gwe/lke} diff --git a/doc/mf6io/gwe/shf.tex b/doc/mf6io/gwe/shf.tex new file mode 100644 index 00000000000..bfd17f185b6 --- /dev/null +++ b/doc/mf6io/gwe/shf.tex @@ -0,0 +1,21 @@ +Input to the Sensible Heat Flux (SHF) Package is read from the file that is specified in the SHF6 record of the OPTIONS block in the SFE package. + +\vspace{5mm} +\subsubsection{Structure of Blocks} +\vspace{5mm} + +\noindent \textit{FOR EACH SIMULATION} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/utl-shf-options.dat} +\vspace{5mm} +\noindent \textit{FOR ANY STRESS PERIOD} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/utl-shf-period.dat} + +\vspace{5mm} +\subsubsection{Explanation of Variables} +\begin{description} +\input{./mf6ivar/tex/utl-shf-desc.tex} +\end{description} + +\vspace{5mm} +\subsubsection{Example Input File} +\lstinputlisting[style=inputfile]{./mf6ivar/examples/utl-shf-example.dat} diff --git a/doc/mf6io/gwe/swr.tex b/doc/mf6io/gwe/swr.tex new file mode 100644 index 00000000000..26de26672ab --- /dev/null +++ b/doc/mf6io/gwe/swr.tex @@ -0,0 +1,21 @@ +Input to the Shortwave Radiation (SWR) Package is read from the file that is specified in the SWR6 record of the OPTIONS block in the SFE package. + +\vspace{5mm} +\subsubsection{Structure of Blocks} +\vspace{5mm} + +\noindent \textit{FOR EACH SIMULATION} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/utl-swr-options.dat} +\vspace{5mm} +\noindent \textit{FOR ANY STRESS PERIOD} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/utl-swr-period.dat} + +\vspace{5mm} +\subsubsection{Explanation of Variables} +\begin{description} +\input{./mf6ivar/tex/utl-swr-desc.tex} +\end{description} + +\vspace{5mm} +\subsubsection{Example Input File} +\lstinputlisting[style=inputfile]{./mf6ivar/examples/utl-swr-example.dat} diff --git a/doc/mf6io/mf6ivar/dfn/gwe-sfe.dfn b/doc/mf6io/mf6ivar/dfn/gwe-sfe.dfn index 63027f1f397..f3df96e0064 100644 --- a/doc/mf6io/mf6ivar/dfn/gwe-sfe.dfn +++ b/doc/mf6io/mf6ivar/dfn/gwe-sfe.dfn @@ -255,6 +255,37 @@ optional false longname obs6 input filename description REPLACE obs6_filename {'{#1}': 'SFE', '{#2}': '\\ref{table:gwe-obstypetable}'} +block options +name abc_filerecord +type record abc6 filein abc6_filename +shape +reader urword +tagged true +optional true +longname +description + +block options +name abc6 +type keyword +shape +in_record true +reader urword +tagged true +optional false +longname abc keyword +description keyword to specify that record corresponds to the atmospheric boundary condition (ABC) file. The behavior of the ABC utility and a description of the input file is provided separately. + +block options +name abc6_filename +type string +preserve_case true +in_record true +reader urword +optional false +tagged false +longname file name of ABC information +description defines a sensible heat flux (ABC) input file. Records in the ABC file are used to calculate the flux of thermal energy in or out of a stream reach resulting from sensible heat exchange. # --------------------- gwe sfe packagedata --------------------- diff --git a/doc/mf6io/mf6ivar/dfn/utl-abc.dfn b/doc/mf6io/mf6ivar/dfn/utl-abc.dfn new file mode 100644 index 00000000000..1d92d756582 --- /dev/null +++ b/doc/mf6io/mf6ivar/dfn/utl-abc.dfn @@ -0,0 +1,334 @@ +# --------------------- gwe abc options --------------------- +# flopy subpackage abc_filerecord abc abc_perioddata perioddata +# flopy parent_name_type parent_package MFPackage + +block options +name print_input +type keyword +reader urword +optional true +longname print input to listing file +description REPLACE print_input {'{#1}': 'atmospheric boundary condition'} + +block options +name ts_filerecord +type record ts6 filein ts6_filename +shape +reader urword +tagged true +optional true +longname +description + +block options +name ts6 +type keyword +shape +in_record true +reader urword +tagged true +optional false +longname head keyword +description keyword to specify that record corresponds to a time-series file. + +block options +name filein +type keyword +shape +in_record true +reader urword +tagged true +optional false +longname file keyword +description keyword to specify that an input filename is expected next. + +block options +name ts6_filename +type string +preserve_case true +in_record true +reader urword +optional false +tagged false +longname file name of time series information +description REPLACE timeseriesfile {} + +block options +name density_air +type double precision +reader urword +optional true +longname density of water +description density of the atmospheric air used by calculations related to sensible heat flux. This value is set to 1.225 kg/m3 if no overriding value is specified. A user-specified value should be provided for models that use units other than kilograms and meters or if it is necessary to use a value other than the default. +default_value 1.225 +mf6internal rhoa + +block options +name heat_capacity_air +type double precision +reader urword +optional true +longname heat capacity of air +description mass-based heat capacity of the air used by calculations related to sensible heat flux. This value is set to 717.0 J/kg/C if no overriding value is specified. A user-specified value should be provided for models that use units other than joules, kilograms, and degrees celsius or if it is necessary to use a value other than the default. +default_value 717.0 +mf6internal cpa + +block options +name drag_coefficient +type double precision +reader urword +optional true +longname drag coefficient +description drag coefficient used by calculations related to sensible heat flux. This value is set to 0.002 (unitless) if no overriding value is specified. A user-specified value should be provided for models that require a value other than the default. +default_value 0.0002 +mf6internal cd + +block options +name wind_func_slope +type double precision +reader urword +optional true +longname wind function slope +description wind function slope used by calculations related to the latent heat flux. This value is set to 1.383E-08 (1/mbar) if no overriding value is specified. A user-specified value should be provided for models that require a value other than the default. +default_value 1.383E-08 +mf6internal wf_slope + +block options +name wind_func_int +type double precision +reader urword +optional true +longname wind function intercept +description wind function intercept used by calculations related to the latent heat flux. This value is set to 3.445E-09 (m/s) if no overriding value is specified. A user-specified value should be provided for models that require a value other than the default. +default_value 3.445E-09 +mf6internal wf_int + +block options +name longwave_reflectance +type double precision +reader urword +optional true +longname longwave reflectance +description longwave reflectance used by calculations related to the longwave radiation heat flux. This value is set to 0.03 if no overriding value is specified. A user-specified value should be provided for models that require a value other than the default. +default_value 0.03 +mf6internal lwrefl + +block options +name emissivity_water +type double precision +reader urword +optional true +longname emissivity of water +description emissivity of water used by calculations related to the longwave radiation heat flux. This value is set to 0.95 if no overriding value is specified. A user-specified value should be provided for models that require a value other than the default. +default_value 0.95 +mf6internal emissw + +block options +name emissivity_canopy +type double precision +reader urword +optional true +longname emissivity of riparian canopy +description emissivity of riparian canopy used by calculations related to the longwave radiation heat flux. This value is set to 0.97 if no overriding value is specified. A user-specified value should be provided for models that require a value other than the default. +default_value 0.97 +mf6internal emissr + +block options +name temperature_factor +type double precision +reader urword +optional true +longname temperature adjustment factor +description real value that is used to convert user-specified temperature units into degrees Kelvin, along with the temperature offset amount (TEMPERATURE\_OFFSET). The empirical equations used by the ABC utility require temperatures be in degrees Kelvin. Users can use temperature units of their own choosing, but if specified in units other than Kelvin, for example degrees Celcius or Fahrenheit, then the user must provide a temperature adjustment factor (TEMPERATURE\_FACTOR) and temperature offset (TEMPERATURE\_OFFSET, described below) for converting units. TEMPERATURE\_FACTOR should be set to 1.0 and 5/9 when using temperature units (TEMPERATURE\_UNITS) of Celcius or Fahrenheit in the simulation, respectively. TEMPERATURE\_FACTOR does not need to be specified if TEMPERATURE\_UNITS are Kelvin. +default_value 1.0 +mf6internal tfac + +block options +name temperature_offset +type double precision +reader urword +optional true +longname temperature offset amount +description real value that is used to convert user-specified temperature units into degrees Kelvin, which the ABC utility empirical equations require. If atmospheric temperature is specified in units other than Kelvin, for example degrees Celcius or Fahrenheit, then TEMPERATURE\_OFFSET should be set to 273.15 and 255.3722 when using temperature units (TEMPERATURE\_UNITS) of Celcius or Fahrenheit in the simulation, respectively. TEMPERATURE\_OFFSET does not need to be specified if TEMPERATURE\_UNITS are Kelvin. +default_value 0.0 +mf6internal toff + +block options +name density_factor +type double precision +reader urword +optional true +longname density adjustment factor +description real value that is used to convert user-specified density units into kilograms per cubic meter. The empirical equations used by the ABC utility require density be in kilograms per cubic meter. Users can use density units of their own choosing, but if specified in units other than kilogram per cubic meter, for example pounds per cubic foot, then the user must provide a density adjustment factor (DENSITY\_FACTOR) for converting units. DENSITY\_FACTOR should be set to 1.0 and 5/9 when using density units of kilograms per cubic meter or pounds per cubic foot in the simulation, respectively. DENSITY\_FACTOR will default to 1.0 if not specified which is acceptable when using units of kilograms per cubic meter. +default_value 1.0 +mf6internal dfac + + +block options +name swr_off +tagged true +type keyword +reader urword +optional true +longname deactivate shortwave radiation calculations +description keyword to deactivate shortwave radiation heat exchange with the atmosphere. Shortwave radiation heat exchange between stream reaches and the atmosphere is active by default. + +block options +name lwr_off +tagged true +type keyword +reader urword +optional true +longname deactivate longwave radiation calculations +description keyword to deactivate longwave radiation heat exchange with the atmosphere. Longwave radiation heat exchange between stream reaches and the atmosphere is active by default. + +block options +name lhf_off +tagged true +type keyword +reader urword +optional true +longname deactivate latent heat flux calculations +description keyword to deactivate latent heat exchange with the atmosphere. Latent heat exchange between stream reaches and the atmosphere is active by default. + +block options +name shf_off +tagged true +type keyword +reader urword +optional true +longname deactivate sensible heat flux calculations +description keyword to deactivate sensible heat exchange with the atmosphere. Sensible heat exchange between stream reaches and the atmosphere is active by default. + +# --------------------- gwe abc period --------------------- + +block period +name iper +type integer +block_variable True +in_record true +tagged false +shape +valid +reader urword +optional false +longname stress period number +description REPLACE iper {} + +block period +name reachperioddata +type recarray rno abcsetting +shape +reader urword +longname +description + +block period +name rno +type integer +shape +tagged false +in_record true +reader urword +longname reach number for this entry +description integer value that defines the reach number associated with the specified PERIOD data on the line. RNO must be greater than zero and less than or equal to NREACHES. +numeric_index true + +block period +name abcsetting +type keystring wspd tatm patm solr shd swrefl rh atmc +auxiliaryrecord +shape +tagged false +in_record true +reader urword +longname +description line of information that is parsed into a parameter keyword and values. Property name keywords that can be used to start the ABCSETTING string include: WSPD, TATM, SOLR, SHD, SWREFL, and RH. These settings are used to assign the wind speed, temperature of the atmosphere, pressure of the atmosphere, shortwave solar radiation, channel shade, shortwave reflectance of the water surface, and relative humidity, respectively. + +block period +name wspd +type string +shape +tagged true +in_record true +reader urword +time_series true +longname wind speed (L/T) +description real or character value that defines the wind speed above of each reach. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value. + +block period +name tatm +type double precision +shape +tagged true +in_record true +reader urword +time_series true +longname air temperature +description is the temperature of the atmosphere above the stream reach. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value. + +block period +name patm +type double precision +shape +tagged true +in_record true +reader urword +time_series true +longname atmospheric pressure +description is the atmospheric pressure above the stream reach (mbar). If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value. + +block period +name solr +type string +shape +tagged true +in_record true +reader urword +time_series true +longname shortwave solar radiation (E/T/L^2) +description real or character value that defines the solar radiation above of each reach. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value. + +block period +name shd +type double precision +shape +tagged true +in_record true +reader urword +time_series true +longname riparian shade +description is the shaded fraction of the wetted surface of the stream reach. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value. + +block period +name swrefl +type double precision +shape +tagged true +in_record true +reader urword +time_series true +longname reflectance of shortwave radiation of water surface +description the fraction of solar radiation reflected by water surface. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value. + +block period +name rh +type double precision +shape +tagged true +in_record true +reader urword +time_series true +longname relative humidity of the atmosphere +description the relative humidity of the atmosphere. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value. + +block period +name atmc +type double precision +shape +tagged true +in_record true +reader urword +time_series true +longname atmospheric composition +description empirical adjustment factor for atmospheric composition such as haze, smoke, or aerosols. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value. diff --git a/doc/mf6io/mf6ivar/examples/gwe-sfe-example-obs.dat b/doc/mf6io/mf6ivar/examples/gwe-sfe-example-obs.dat index f3c7b651af5..9004c2f19cb 100644 --- a/doc/mf6io/mf6ivar/examples/gwe-sfe-example-obs.dat +++ b/doc/mf6io/mf6ivar/examples/gwe-sfe-example-obs.dat @@ -22,4 +22,5 @@ BEGIN continuous FILEOUT gwe_sfe02.sfe.obs.csv sfe-5-fjf FLOW-JA-FACE MYREACH1 sfe-6-fjf FLOW-JA-FACE MYREACH2 sfe-7-fjf FLOW-JA-FACE MYREACH3 + sfe-1-shf SHF 1 END continuous diff --git a/doc/mf6io/mf6ivar/examples/utl-abc-example.dat b/doc/mf6io/mf6ivar/examples/utl-abc-example.dat new file mode 100644 index 00000000000..3d288dd22fa --- /dev/null +++ b/doc/mf6io/mf6ivar/examples/utl-abc-example.dat @@ -0,0 +1,29 @@ +BEGIN OPTIONS + PRINT_INPUT + DENSITY_AIR 1.225 + HEAT_CAPACITY_AIR 717.0 + DRAG_COEFFICIENT 0.002 + WIND_FUNC_SLOPE 1.383E-08 + WIND_FUNC_INT 3.445E-09 + LWREFL 0.03 + EMISSW 0.95 + EMISSR 0.97 +END OPTIONS + +BEGIN PERIOD 1 + 1 WSPD 5.00 + 1 TATM 12.97 + 1 SOLR 801.0 + 1 SHD 0.1 + 1 SWREFL 0.03 + 1 RH 30.00 + 1 ATMC 60.00 + 1 ATMC + 2 WSPD 4.95 + 2 TATM 12.95 + 2 SOLR 805.0 + 2 SHD 0.25 + 2 SWREFL 0.03 + 2 RH 31.00 + 2 ATMC 55.00 +END PERIOD 1 diff --git a/msvs/mf6.sln b/msvs/mf6.sln index b70497a5b6a..64b264e94c6 100644 --- a/msvs/mf6.sln +++ b/msvs/mf6.sln @@ -36,10 +36,10 @@ Global {380139B4-29CE-4CF3-BD74-5AA979FD0F43}.Debug|Win32.Build.0 = Debug|Win32 {380139B4-29CE-4CF3-BD74-5AA979FD0F43}.Debug|x64.ActiveCfg = Debug|x64 {380139B4-29CE-4CF3-BD74-5AA979FD0F43}.Debug|x64.Build.0 = Debug|x64 - {380139B4-29CE-4CF3-BD74-5AA979FD0F43}.Profile|Win32.ActiveCfg = Release|Win32 - {380139B4-29CE-4CF3-BD74-5AA979FD0F43}.Profile|Win32.Build.0 = Release|Win32 - {380139B4-29CE-4CF3-BD74-5AA979FD0F43}.Profile|x64.ActiveCfg = Release|x64 - {380139B4-29CE-4CF3-BD74-5AA979FD0F43}.Profile|x64.Build.0 = Release|x64 + {380139B4-29CE-4CF3-BD74-5AA979FD0F43}.Profile|Win32.ActiveCfg = Debug|Win32 + {380139B4-29CE-4CF3-BD74-5AA979FD0F43}.Profile|Win32.Build.0 = Debug|Win32 + {380139B4-29CE-4CF3-BD74-5AA979FD0F43}.Profile|x64.ActiveCfg = Debug|x64 + {380139B4-29CE-4CF3-BD74-5AA979FD0F43}.Profile|x64.Build.0 = Debug|x64 {380139B4-29CE-4CF3-BD74-5AA979FD0F43}.Release|Win32.ActiveCfg = Release|Win32 {380139B4-29CE-4CF3-BD74-5AA979FD0F43}.Release|Win32.Build.0 = Release|Win32 {380139B4-29CE-4CF3-BD74-5AA979FD0F43}.Release|x64.ActiveCfg = Release|x64 diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index d088a615f9b..3bae57cde49 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -2,50 +2,44 @@ - - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -339,15 +333,21 @@ + + + + + + < diff --git a/src/Model/GroundWaterEnergy/PbstBase.f90 b/src/Model/GroundWaterEnergy/PbstBase.f90 new file mode 100644 index 00000000000..293e86901cd --- /dev/null +++ b/src/Model/GroundWaterEnergy/PbstBase.f90 @@ -0,0 +1,265 @@ +!> @brief This module contains common process-based stream temperature functionality +!! +!! This module contains methods for implementing functionality associated with +!! heat fluxes to a stream reach. Four sources of thermal energy commonly +!! accounted for in process-based stream temperature modeling include short- +!! wave radiation, long-wave radiation, sensible heat flux, and latent heat +!! flux. +!< +module PbstBaseModule + use ConstantsModule, only: LINELENGTH, MAXCHARLEN, DZERO, LGP, & + LENPACKAGENAME, TABLEFT, TABCENTER, & + LENVARNAME, DCTOK, DHUNDRED, DONESEVENTH + use KindModule, only: I4B, DP + use NumericalPackageModule, only: NumericalPackageType + use SimModule, only: count_errors, store_error, ustop + use SimVariablesModule, only: errmsg + use TdisModule, only: kper, nper, kstp + use TimeSeriesLinkModule, only: TimeSeriesLinkType + use TimeSeriesManagerModule, only: TimeSeriesManagerType, tsmanager_cr + use TableModule, only: TableType, table_cr + + implicit none + + private + + public :: PbstBaseType + public :: pbstbase_da + + type, extends(NumericalPackageType) :: PbstBaseType + + integer(I4B), pointer :: ncv => null() !< number of control volumes + integer(I4B), dimension(:), pointer, contiguous :: iboundpbst => null() !< package ibound + logical, pointer, public :: active => null() !< logical indicating if a sensible heat flux object is active + character(len=8), dimension(:), pointer, contiguous :: status => null() !< active, inactive, constant + character(len=LENPACKAGENAME) :: text = '' !< text string for package transport term + character(len=LINELENGTH), pointer, public :: inputFilename => null() !< a particular pbst input file name, could be for sensible heat flux or latent heat flux subpackages, for example + type(TimeSeriesManagerType), pointer :: tsmanager => null() + ! + ! -- table objects + type(TableType), pointer :: inputtab => null() !< input table object + + contains + + procedure :: init + procedure :: pbst_ar + procedure :: pbst_options + procedure :: subpck_set_stressperiod + procedure, private :: pbstbase_allocate_scalars + procedure :: da => pbstbase_da + procedure :: pbst_check_valid + procedure, public :: epsa !< function for calculating atmospheric emissivity + procedure, public :: epss !< function for calculating shade-weighted atmospheric emissivity + procedure, public :: evap !< function for calculating evaporation rate using a mass transfer method + + end type PbstBaseType + +contains + + !> @brief Initialize the PbstBaseType object + !! + !! Allocate and initialize data members of the object. + !< + subroutine init(this, name_model, pakname, ftype, inunit, iout, ncv) + ! -- dummy + class(PbstBaseType) :: this + character(len=*), intent(in) :: name_model + character(len=*), intent(in) :: pakname + character(len=*), intent(in) :: ftype + integer(I4B), intent(in) :: inunit + integer(I4B), intent(in) :: iout + integer(I4B), target, intent(in) :: ncv + ! + call this%set_names(1, name_model, pakname, ftype) + call this%pbstbase_allocate_scalars() + this%inunit = inunit + this%iout = iout + this%ncv => ncv + call this%parser%Initialize(this%inunit, this%iout) + end subroutine init + + !> @brief Allocate and read + !! + !! Method to allocate and read static data for the PBST utility packages + !< + subroutine pbst_ar(this) + ! -- dummy + class(PbstBaseType) :: this !< ShfType object + ! + ! -- will be overridden + end subroutine pbst_ar + + !> @brief pbst_set_stressperiod() + !! + !! To be overridden by Pbst sub-packages + !< + subroutine subpck_set_stressperiod(this, itemno, keyword, found) + ! -- dummy + class(PbstBaseType), intent(inout) :: this + integer(I4B), intent(in) :: itemno + character(len=*), intent(in) :: keyword + logical, intent(inout) :: found + ! -- to be overwritten by pbst subpackages (or "utilities") + end subroutine subpck_set_stressperiod + + !> @brief Read additional options for sub-package + !! + !! Read additional options for the SFE boundary package. This method should + !! be overridden by option-processing routine that is in addition to the + !! base options available for all PbstBase packages. + !< + subroutine pbst_options(this, option, found) + ! -- dummy + class(PbstBaseType), intent(inout) :: this !< PbstBaseType object + character(len=*), intent(inout) :: option !< option keyword string + logical(LGP), intent(inout) :: found !< boolean indicating if the option was found + ! + ! -- return with found = .false. + found = .false. + end subroutine pbst_options + + !> @brief Allocate scalar variables + !! + !! Allocate scalar data members of the object. + !< + subroutine pbstbase_allocate_scalars(this) + ! -- modules + use MemoryManagerModule, only: mem_allocate + ! -- dummy + class(PbstBaseType) :: this + ! + allocate (this%active) + allocate (this%inputFilename) + ! + ! -- initialize + this%active = .false. + this%inputFilename = '' + ! + ! -- call standard NumericalPackageType allocate scalars + call this%NumericalPackageType%allocate_scalars() + ! + ! -- allocate + call mem_allocate(this%ncv, 'NCV', this%memoryPath) + ! + ! -- initialize + this%ncv = 0 + ! + ! -- allocate time series manager + allocate (this%tsmanager) + end subroutine pbstbase_allocate_scalars + + !> @brief Deallocate package memory + !! + !! Deallocate package scalars and arrays. + !< + subroutine pbstbase_da(this) + ! -- modules + use MemoryManagerModule, only: mem_deallocate + ! -- dummy + class(PbstBaseType) :: this + ! + deallocate (this%active) + deallocate (this%inputFilename) + ! + ! -- deallocate time series manager + deallocate (this%tsmanager) + ! + ! -- deallocate scalars + call mem_deallocate(this%ncv) + ! + ! -- deallocate arrays + call mem_deallocate(this%iboundpbst) + ! + ! -- deallocate parent + call this%NumericalPackageType%da() + ! + ! -- input table object + if (associated(this%inputtab)) then + call this%inputtab%table_da() + deallocate (this%inputtab) + nullify (this%inputtab) + end if + end subroutine pbstbase_da + + !> @brief Process-based stream temperature transport (or utility) routine + !! + !! Determine if a valid feature number has been specified. + !< + function pbst_check_valid(this, itemno) result(ierr) + ! -- return + integer(I4B) :: ierr + ! -- dummy + class(PbstBaseType), intent(inout) :: this + integer(I4B), intent(in) :: itemno + ! + ! -- initialize + ierr = 0 + ! + if (itemno < 1 .or. itemno > this%ncv) then + write (errmsg, '(a,1x,i6,1x,a,1x,i6)') & + 'Featureno ', itemno, 'must be > 0 and <= ', this%ncv + call store_error(errmsg) + ierr = 1 + end if + end function pbst_check_valid + + !> @brief Calculate atmospheric emissivity + !! + !! Atmospheric emissivity is highly dependent upon ambient vapor + !! concentration of the air (Campbell and Norman 1989). Estimated + !! clear-sky emissivity is calculated using (Brutsaert and Jirka 1984) + !< + function epsa(this, eatm, tatm, atmc) + ! -- dummy + class(PbstBaseType) :: this + real(DP) :: atmc + real(DP) :: eatm + real(DP) :: tatm + ! -- return + real(DP) :: epsa !< atmospheric emissivity + ! + epsa = 1.24_DP * (eatm / tatm)**(DONESEVENTH) * atmc + end function epsa + + !> @brief Calculate shade-altered emissivity + !! + !! Riparian shade can alter the above-channel emissivity which therefore + !! affects the net longwave heat flux. Emissivity above the stream channel + !! is the combination of shade-weighted average of clear sky atmospheric + !! emissivity (epsa) and riparian canopy emissivity + !< + function epss(this, shd, epsa, epsr) + ! -- dummy + class(PbstBaseType) :: this + real(DP) :: shd !< percent shading, expressed as a fraction + real(DP) :: epsa !< atmospheric emissivity + real(DP) :: epsr !< riparian canopy emissivity + ! -- return + real(DP) :: epss !< shade weighted atmospheric emissivity + ! + epss = (1 - shd) * epsa + shd * epsr + end function epss + + !> @brief Calculate Evaporation Rate + !! + !! Used in the latent heat flux calculation and for reporting the calculated + !! evaporation rate as an observation. Units of length per time + !! (Need to circle back to this, but for the time being I believe this + !! calculation has units of mm/day) + !< + function evap(this, wfint, wfslp, wspd, ew, ea) + ! -- dummy + class(PbstBaseType) :: this + real(DP) :: wfint !< wind function intercept + real(DP) :: wfslp !< wind function slope + real(DP) :: wspd !< wind speed + real(DP) :: ew !< saturation vapor pressure at stream temperature + real(DP) :: ea !< ambient vapor pressure at atmospheric temperature + ! -- return + real(DP) :: evap !< evaporation rate + ! + ! -- mass-transfer method for calculating evap rate (A.17) + evap = (wfint + wfslp * wspd) * (ew - ea) + end function evap + +end module PbstBaseModule diff --git a/src/Model/GroundWaterEnergy/gwe-abc.f90 b/src/Model/GroundWaterEnergy/gwe-abc.f90 new file mode 100644 index 00000000000..64be6b7673e --- /dev/null +++ b/src/Model/GroundWaterEnergy/gwe-abc.f90 @@ -0,0 +1,1080 @@ +!> @brief This module contains the atmospheric boundary condition +!! +!! This module contains the methods used to calculate the heat fluxes +!! for surface-water boundaries, like streams and lakes. In its current form, +!! this class acts like a package to a package, similar to the TVK package that +!! can be invoked from the NPF package. Once this package is completed in its +!! prototyped form, it will likely be moved around. +!< + +! SFR flows (sfrbudptr) index var SFE term Equation +! --------------------------------------------------------------------------------- +! -- PBST terms +! SHORTWAVE RADIATION idxbudswr SHORTWAVE (1 - shd) * (1 - swrefl) * solr +! LONGWAVE RADIATION idxbudlwr LONGWAVE longwv_in * (1 - lwrefl) + longwv_out +! LATENT HEAT FLUX idxbudlhf LATENT HEAT evap_rate * latent_heat_vaporization * rho_w +! SENSIBLE HEAT FLUX idxbudshf SENS HEAT cd * rho_a * C_p_a * wspd * (t_air - t_feat) + +module AbcModule + use ConstantsModule, only: LINELENGTH, LENMEMPATH, DZERO, DONE, LENVARNAME, & + LENPACKAGENAME, TABLEFT, TABCENTER, LENMEMTYPE, & + DHUNDRED, DCTOK + use KindModule, only: I4B, DP + use MemoryManagerModule, only: mem_setptr + use MemoryHelperModule, only: create_mem_path + use SimModule, only: store_error, count_errors + use BaseDisModule, only: DisBaseType + use SimVariablesModule, only: errmsg + use PbstBaseModule, only: PbstBaseType, pbstbase_da + use SensHeatModule, only: ShfType, shf_cr + use ShortwaveModule, only: SwrType, swr_cr + use LatHeatModule, only: LhfType, lhf_cr + use LongwaveModule, only: LwrType, lwr_cr + use ObserveModule + use BudgetObjectModule, only: BudgetObjectType, budgetobject_cr + use NumericalPackageModule, only: NumericalPackageType + use TimeArraySeriesManagerModule, only: TimeArraySeriesManagerType + use TimeSeriesLinkModule, only: TimeSeriesLinkType + use TimeSeriesManagerModule, only: TimeSeriesManagerType + use TableModule, only: TableType, table_cr + use BndModule, only: BndType + use GweInputDataModule, only: GweInputDataType + + implicit none + + private + + public :: AbcType + public :: abc_cr + + character(len=LENVARNAME) :: text = ' ABC' + + type, extends(BndType) :: AbcType + + type(GweInputDataType), pointer :: gwecommon => null() !< pointer to shared gwe data used by multiple packages but set in est + + character(len=8), dimension(:), pointer, contiguous :: status => null() !< active, inactive, constant + integer(I4B), pointer :: ncv => null() !< number of control volumes + integer(I4B), dimension(:), pointer, contiguous :: iboundpbst => null() !< package ibound + character(len=LINELENGTH), pointer, public :: inputFilename => null() !< a particular abc input file name, could be for sensible heat flux or latent heat flux subpackages, for example + logical, pointer, public :: active => null() !< logical indicating if a atmospheric boundary condition object is active + ! -- table objects + !type(TableType), pointer :: inputtab => null() !< input table object + ! + logical, pointer, public :: swr_active => null() !< logical indicating if a shortwave radiation heat flux object is active + logical, pointer, public :: lwr_active => null() !< logical indicating if a longwave radiation heat flux object is active + logical, pointer, public :: lhf_active => null() !< logical indicating if a latent heat flux object is active + logical, pointer, public :: shf_active => null() !< logical indicating if a sensible heat flux object is active + ! + type(ShfType), pointer :: shf => null() ! sensible heat flux (shf) object + type(SwrType), pointer :: swr => null() ! shortwave radiation heat flux (swr) object + type(LhfType), pointer :: lhf => null() ! latent heat flux (lhf) object + type(LwrType), pointer :: lwr => null() ! longwave radiation heat flux (lwr) object + ! + ! -- abc budget object + type(BudgetObjectType), pointer :: budobj => null() !< ABC budget object + ! + real(DP), pointer :: rhoa => null() !< desity of air + real(DP), pointer :: cpa => null() !< heat capacity of air + real(DP), pointer :: cd => null() !< drag coefficient + real(DP), pointer :: wfslope => null() !< wind function slope + real(DP), pointer :: wfint => null() !< wind function intercept + real(DP), pointer :: lwrefl => null() !< reflectance of longwave radiation by the water surface + real(DP), pointer :: emissw => null() !< emissivity of water + real(DP), pointer :: emissr => null() !< emissivity of the riparian canopy + real(DP), pointer :: tfac => null() !< temperature adjustment factor + real(DP), pointer :: toff => null() !< temperature units offset + real(DP), pointer :: dfac => null() !< density adjustment factor + real(DP), pointer :: pfac => null() !< atmospheric pressure adjustment factor + ! + real(DP), dimension(:), pointer, contiguous :: wspd => null() !< wind speed + real(DP), dimension(:), pointer, contiguous :: tatm => null() !< temperature of the atmosphere + real(DP), dimension(:), pointer, contiguous :: solr => null() !< solar radiation + real(DP), dimension(:), pointer, contiguous :: shd => null() !< shade fraction + real(DP), dimension(:), pointer, contiguous :: swrefl => null() !< shortwave reflectance of water surface + real(DP), dimension(:), pointer, contiguous :: rh => null() !< relative humidity + real(DP), dimension(:), pointer, contiguous :: atmc => null() !< atmospheric composition adjustment + real(DP), dimension(:), pointer, contiguous :: patm => null() !< atmospheric pressure (mbar) + ! + real(DP), dimension(:), pointer, contiguous :: ea => null() !< ambient vapor pressure of the atmosphere, internally calculated and used by multiple heat calculations + real(DP), dimension(:), pointer, contiguous :: es => null() !< saturation vapor pressure at air temperature, make available for multiple heat flux calculations + real(DP), dimension(:), pointer, contiguous :: ew => null() !< saturation vapor pressure at water temperature, make available for multiple heat flux calculations + + contains + + procedure, public :: abc_df + procedure :: da => abc_da + procedure :: ar + procedure, public :: abc_rp + procedure :: abc_check_valid + procedure :: bnd_options => abc_read_options + procedure :: abc_set_stressperiod + procedure :: abc_allocate_arrays + procedure, private :: abc_allocate_scalars + procedure, public :: abc_cq + procedure, private :: recalc_shared_vars + procedure, private :: calc_eatm !< function for calculating ambient vapor pressure of the atmosphere + procedure, private :: check_for_specified_input !< check for potential omissions in the specified input + procedure, public :: abc_evap !< fetches the evaporation rate for reporting as an observation + + end type AbcType + +contains + + !> @brief Create a new AbcType object + !! + !! Create a new atmospheric boundary condition (AbcType) object. Initially for use with + !! the SFE package. + !< + subroutine abc_cr(abc, name_model, inunit, iout, fname, ncv, gwecommon, dis) + ! -- modules + use TimeSeriesManagerModule, only: tsmanager_cr + use TimeArraySeriesManagerModule, only: tasmanager_cr + ! -- dummy + type(AbcType), pointer, intent(out) :: abc + character(len=*), intent(in) :: name_model + integer(I4B), intent(in) :: inunit + integer(I4B), intent(in) :: iout + character(len=LINELENGTH), intent(in) :: fname + integer(I4B), target, intent(in) :: ncv + type(GweInputDataType), intent(in), target :: gwecommon !< shared data container for use by multiple GWE packages + class(DisBaseType), pointer :: dis !< discretization object + ! + ! -- Create the object + allocate (abc) + ! + ! -- set pointer to dis object for the model + abc%dis => dis + ! + call abc%set_names(1, name_model, 'ABC', 'ABC') + ! + ! -- call parent's define routine + !call abc%bnd_df + ! + ! -- allocate scalars + call abc%abc_allocate_scalars() + ! + abc%inunit = inunit + abc%iout = iout + abc%inputFilename = fname + abc%ncv => ncv + call abc%parser%Initialize(abc%inunit, abc%iout) + ! + ! -- initialize associated abc utilities + call shf_cr(abc%shf, name_model, inunit, iout, ncv) + call swr_cr(abc%swr, name_model, inunit, iout, ncv) + call lhf_cr(abc%lhf, name_model, inunit, iout, ncv) + call lwr_cr(abc%lwr, name_model, inunit, iout, ncv) + ! + ! -- Create time series managers + call tsmanager_cr(abc%tsmanager, abc%iout) + call tasmanager_cr(abc%TasManager, dis, abc%name_model, abc%iout) + ! + ! -- Store pointer to shared data module for accessing cpw, rhow + ! for the heat flux calculations + abc%gwecommon => gwecommon + end subroutine abc_cr + + !> @brief Define routine for ABC + !! + !! Run df routines for the tsmanager + !< + subroutine abc_df(this) + ! -- dummy + class(AbcType), intent(inout) :: this !< AbcType object + ! + ! -- Now that time series will have been read, need to call the df + ! routine to define the manager + call this%tsmanager%tsmanager_df() + call this%tasmanager%tasmanager_df() + end subroutine abc_df + + !> @brief Allocate and read + !! + !! Method to allocate and read static data for the SHF, SWR, and LHF sub-utilities + !< + subroutine ar(this) + ! -- dummy + class(AbcType) :: this !< AbcType object + ! -- formats + character(len=*), parameter :: fmtapt = & + "(1x,/1x,'ABC -- ATMOSPHERIC BOUNDARY CONDITION PACKAGE, VERSION 1, 7/20/2025', & + &' INPUT READ FROM UNIT ', i0, //)" + ! + ! -- print a message identifying the apt package. + write (this%iout, fmtapt) this%inunit + ! + end subroutine ar + + !> @brief ABC read and prepare for setting stress period information + !< + subroutine abc_rp(this) + ! -- module + use TimeSeriesManagerModule, only: read_value_or_time_series_adv + use TdisModule, only: kper, nper + ! -- dummy + class(AbcType) :: this !< AbcType object + ! -- local + integer(I4B) :: ierr + integer(I4B) :: n + logical :: isfound, endOfBlock + character(len=LINELENGTH) :: title + character(len=LINELENGTH) :: line + integer(I4B) :: itemno + ! -- formats + character(len=*), parameter :: fmtblkerr = & + &"('Error. Looking for BEGIN PERIOD iper. Found ', a, ' instead.')" + character(len=*), parameter :: fmtlsp = & + &"(1X,/1X,'REUSING ',A,'S FROM LAST STRESS PERIOD')" + ! + ! -- Set ionper to the stress period number for which a new block of data + ! will be read. + if (this%inunit == 0) return + ! + ! -- get stress period data + if (this%ionper < kper) then + ! + ! -- get period block + call this%parser%GetBlock('PERIOD', isfound, ierr, & + supportOpenClose=.true., & + blockRequired=.false.) + if (isfound) then + ! + ! -- read ionper and check for increasing period numbers + call this%read_check_ionper() + else + ! + ! -- PERIOD block not found + if (ierr < 0) then + ! -- End of file found; data applies for remainder of simulation. + this%ionper = nper + 1 + else + ! -- Found invalid block + call this%parser%GetCurrentLine(line) + write (errmsg, fmtblkerr) adjustl(trim(line)) + call store_error(errmsg) + call this%parser%StoreErrorUnit() + end if + end if + end if + ! + ! -- Read data if ionper == kper + if (this%ionper == kper) then + ! + ! -- setup table for period data + if (this%iprpak /= 0) then + ! + ! -- reset the input table object + title = trim(adjustl(this%text))//' PACKAGE ('// & + trim(adjustl(this%packName))//') DATA FOR PERIOD' + write (title, '(a,1x,i6)') trim(adjustl(title)), kper + call table_cr(this%inputtab, this%packName, title) + call this%inputtab%table_df(1, 4, this%iout, finalize=.FALSE.) + text = 'NUMBER' + call this%inputtab%initialize_column(text, 10, alignment=TABCENTER) + text = 'KEYWORD' + call this%inputtab%initialize_column(text, 20, alignment=TABLEFT) + do n = 1, 2 + write (text, '(a,1x,i6)') 'VALUE', n + call this%inputtab%initialize_column(text, 15, alignment=TABCENTER) + end do + end if + ! + ! -- read data + stressperiod: do + call this%parser%GetNextLine(endOfBlock) + if (endOfBlock) exit + ! + ! -- get feature number + itemno = this%parser%GetInteger() + ! + ! -- read data from the rest of the line + call this%abc_set_stressperiod(itemno) + ! + ! -- write line to table + if (this%iprpak /= 0) then + call this%parser%GetCurrentLine(line) + call this%inputtab%line_to_columns(line) + end if + end do stressperiod + ! + if (this%iprpak /= 0) then + call this%inputtab%finalize_table() + end if + ! + ! -- using stress period data from the previous stress period + else + write (this%iout, fmtlsp) trim(this%filtyp) + end if + ! + ! -- evaluate sufficiency of entered data + call this%check_for_specified_input() + ! + ! -- write summary of stress period error messages + ierr = count_errors() + if (ierr > 0) then + call this%parser%StoreErrorUnit() + end if + end subroutine abc_rp + + !> @brief Print warning messages for input data omissions + !! + !! For the ABC package, some input does not necessarily need to be entered + !! and will work with the default of zero, like relative humidity, for + !! example. However, it may also lead to erroneous results. This routine + !! will print a warning to the listing file where suspected data omissions + !! may exist. + !< + subroutine check_for_specified_input(this) + ! -- module + use SimModule, only: store_warning + use SimVariablesModule, only: warnmsg + ! -- dummy + class(AbcType) :: this !< AbcType object + ! -- local + integer(I4B) :: i + ! + ! -- if + if (this%lwr_active) then + do i = 1, this%ncv + if (this%rh(i) == DZERO) then + write (warnmsg, *) "The longwave heat flux calculations is "// & + "active but relatively humidiy remains at the default value "// & + "of 0.0. This will cause the ambient vapor pressure to be "// & + "equal to zero." + call store_warning(warnmsg) + end if + end do + end if + ! + if (this%lhf_active) then + do i = 1, this%ncv + if (this%rh(i) == DZERO) then + write (warnmsg, *) "The latent heat flux calculations is "// & + "active but relatively humidiy remains at the default value "// & + "of 0.0. This will cause the ambient vapor pressure to be "// & + "equal to zero." + call store_warning(warnmsg) + end if + end do + end if + end subroutine check_for_specified_input + + !> @brief Set options specific to the AbcType + !! + !! This routine overrides TspAptType%gc_options + !< + subroutine abc_read_options(this, option, found) + ! -- modules + use ConstantsModule, only: MAXCHARLEN, LGP + use InputOutputModule, only: urword, getunit, assign_iounit, openfile + ! -- dummy + class(AbcType), intent(inout) :: this + character(len=*), intent(inout) :: option + logical, intent(inout) :: found + ! -- formats + character(len=*), parameter :: fmtaptbin = & + "(4x, a, 1x, a, 1x, ' WILL BE SAVED TO FILE: ', a, & + &/4x, 'OPENED ON UNIT: ', I0)" + ! + found = .true. + select case (option) + case ('DENSITY_AIR') + this%rhoa = this%parser%GetDouble() + if (this%rhoa <= 0.0) then + write (errmsg, '(a)') 'Specified value for the density of & + &the atmosphere must be greater than 0.0.' + call store_error(errmsg) + call this%parser%StoreErrorUnit() + else + write (this%iout, '(4x,a,1pg15.6)') & + "The density of the atmosphere has been set to: ", this%rhoa + end if + case ('HEAT_CAPACITY_AIR') + this%cpa = this%parser%GetDouble() + if (this%cpa <= 0.0) then + write (errmsg, '(a)') 'Specified value for the heat capacity of & + &the atmosphere must be greater than 0.0.' + call store_error(errmsg) + call this%parser%StoreErrorUnit() + else + write (this%iout, '(4x,a,1pg15.6)') & + "The heat capacity of the atmosphere has been set to: ", this%cpa + end if + case ('DRAG_COEFFICIENT') + this%cd = this%parser%GetDouble() + write (this%iout, '(4x,a,1pg15.6)') & + "The surface-atmosphere drag coefficient has been set to: ", this%cd + !end if + case ('WIND_FUNC_SLOPE') + this%wfslope = this%parser%GetDouble() + if (this%wfslope <= 0.0) then + write (errmsg, '(a)') 'Specified value for the wind function slope & + &must be greater than 0.0.' + call store_error(errmsg) + call this%parser%StoreErrorUnit() + else + write (this%iout, '(4x,a,1pg15.6)') & + "The evaporation wind function slope has been set to: ", this%wfslope + end if + case ('WIND_FUNC_INT') + this%wfint = this%parser%GetDouble() + if (this%wfint <= 0.0) then + write (errmsg, '(a)') 'Specified value for the wind function intercept & + &must be greater than 0.0.' + call store_error(errmsg) + call this%parser%StoreErrorUnit() + else + write (this%iout, '(4x,a,1pg15.6)') & + "The evaporation wind function intercept has been set to: ", this%wfint + end if + case ('LONGWAVE_REFLECTANCE') + this%lwrefl = this%parser%GetDouble() + if (this%lwrefl <= 0.0) then + write (errmsg, '(a)') 'Specified value for the reflectance of longwave radiation & + &must be greater than 0.0.' + call store_error(errmsg) + call this%parser%StoreErrorUnit() + else + write (this%iout, '(4x,a,1pg15.6)') & + "The reflectance of longwave radiation has been set to: ", this%lwrefl + end if + case ('EMISSIVITY_WATER') + this%emissw = this%parser%GetDouble() + if (this%emissw <= 0.0) then + write (errmsg, '(a)') 'Specified value for the emissivity of water & + &must be greater than 0.0.' + call store_error(errmsg) + call this%parser%StoreErrorUnit() + else + write (this%iout, '(4x,a,1pg15.6)') & + "The emissivity of water has been set to: ", this%emissw + end if + case ('EMISSIVITY_CANOPY') + this%emissr = this%parser%GetDouble() + if (this%emissr <= 0.0) then + write (errmsg, '(a)') 'Specified value for the emissivity of the riparian canopy & + &must be greater than 0.0.' + call store_error(errmsg) + call this%parser%StoreErrorUnit() + else + write (this%iout, '(4x,a,1pg15.6)') & + "The emissivity of the riparian canopy has been set to: ", this%emissw + end if + case ('TEMPERATURE_FACTOR') + this%tfac = this%parser%GetDouble() + write (this%iout, '(4x,a,1pg15.6)') & + "The temperature adjustment factor has been set to: ", this%tfac + case ('TEMPERATURE_OFFSET') + this%toff = this%parser%GetDouble() + write (this%iout, '(4x,a,1pg15.6)') & + "A temperature offset value has been set to: ", this%toff + case ('DENSITY_FACTOR') + this%dfac = this%parser%GetDouble() + write (this%iout, '(4x,a,1pg15.6)') & + "The density adjustment factor has been set to: ", this%dfac + case ('PRESSURE_FACTOR') + this%pfac = this%parser%GetDouble() + write (this%iout, '(4x,a,1pg15.6)') & + "The atmospheric pressure adjustment factor has been set to: ", this%pfac + case ('SWR_OFF') + this%swr_active = .false. + write (this%iout, '(4x,a)') & + "Shortwave thermal energy exchange between the stream reaches and the & + &atmosphere has been turned off." + case ('LWR_OFF') + this%lwr_active = .false. + write (this%iout, '(4x,a)') & + "Longwave thermal energy exchange between stream reaches and the & + &atmosphere has been turned off." + case ('LHF_OFF') + this%lhf_active = .false. + write (this%iout, '(4x,a)') & + "Latent heat exchange between stream reaches and the atmosphere & + &has been turned off." + case ('SHF_OFF') + this%shf_active = .false. + write (this%iout, '(4x,a)') & + "Sensible heat exchange between the stream reaches and the & + &atmosphere has been turned off." + case default + write (errmsg, '(a,a)') 'Unknown ABC option: ', trim(option) + call store_error(errmsg) + call this%parser%StoreErrorUnit() + end select + end subroutine abc_read_options + + !> @brief Allocate scalars specific to the streamflow energy transport (SFE) + !! package. + !< + subroutine abc_allocate_scalars(this) + ! -- modules + use MemoryManagerModule, only: mem_allocate + ! -- dummy + class(AbcType) :: this + ! + allocate (this%active) + allocate (this%inputFilename) + ! + ! -- initialize + this%active = .false. + this%inputFilename = '' + ! + ! -- allocate + call mem_allocate(this%shf_active, 'SHF_ACTIVE', this%memoryPath) + call mem_allocate(this%swr_active, 'SWR_ACTIVE', this%memoryPath) + call mem_allocate(this%lhf_active, 'LHF_ACTIVE', this%memoryPath) + call mem_allocate(this%lwr_active, 'LWR_ACTIVE', this%memoryPath) + ! + ! -- allocate SHF specific + call mem_allocate(this%rhoa, 'RHOA', this%memoryPath) + call mem_allocate(this%cpa, 'CPA', this%memoryPath) + call mem_allocate(this%cd, 'CD', this%memoryPath) + ! -- allocate LHF specific + call mem_allocate(this%wfslope, 'WFSLOPE', this%memoryPath) + call mem_allocate(this%wfint, 'WFINT', this%memoryPath) + ! -- allocate LWR specific + call mem_allocate(this%lwrefl, 'LWREFL', this%memoryPath) + call mem_allocate(this%emissw, 'EMISSW', this%memoryPath) + call mem_allocate(this%emissr, 'EMISSR', this%memoryPath) + ! -- unit conversions + call mem_allocate(this%tfac, 'TFAC', this%memoryPath) + call mem_allocate(this%toff, 'TOFF', this%memoryPath) + call mem_allocate(this%dfac, 'DFAC', this%memoryPath) + call mem_allocate(this%pfac, 'PFAC', this%memoryPath) + ! + ! -- initialize to default values + this%shf_active = .true. ! Initialize to one for 'on' + this%swr_active = .true. + this%lhf_active = .true. + this%lwr_active = .true. + ! -- initialize to SHF specific default values + this%rhoa = 1.225 ! kg/m3 + this%cpa = 717.0 ! J/kg/C + this%cd = 0.002 ! unitless + ! -- initialize to LHF specific default values + this%wfslope = 1.383e-08 ! 1/mbar Fogg 2023 (change!) + this%wfint = 3.445e-09 ! m/s Fogg 2023 (change!) + ! -- initialize to LWR specific default values + this%lwrefl = 0.03 ! unitless (Anderson, 1954) + this%emissw = 0.95 ! unitless (Dingman, 2015) + this%emissr = 0.97 ! unitless (Sobrino et al, 2005) + ! -- initialize temperature offset such that it assumes atmospheric temperatures are specified in Kelvin + this%tfac = DONE + this%toff = DZERO + this%dfac = DONE + this%pfac = DONE + ! + ! -- call standard NumericalPackageType allocate scalars + call this%BndType%allocate_scalars() + ! + ! -- allocate time series manager + allocate (this%tsmanager) + end subroutine abc_allocate_scalars + + !> @brief Allocate arrays specific to the atmspheric boundary package + !< + subroutine abc_allocate_arrays(this) + ! -- modules + !! use MemoryManagerModule, only: mem_allocate + use MemoryManagerModule, only: mem_setptr, mem_checkin, mem_allocate, & + get_mem_type, mem_reallocate + ! -- dummy + class(AbcType), intent(inout) :: this + ! integer(I4B), dimension(:), pointer, contiguous, optional :: nodelist + !real(DP), dimension(:, :), pointer, contiguous, optional :: auxvar + ! -- local + integer(I4B) :: n + ! + ! -- allocate character array for status + allocate (this%status(this%ncv)) + ! + ! -- initialize arrays + do n = 1, this%ncv + this%status(n) = 'ACTIVE' + end do + ! + ! -- allocate all atmospheric boundary condition vars, initialize to size 0 + call mem_allocate(this%wspd, 0, 'WSPD', this%memoryPath) + call mem_allocate(this%tatm, 0, 'TATM', this%memoryPath) + call mem_allocate(this%solr, 0, 'SOLR', this%memoryPath) + call mem_allocate(this%shd, 0, 'SHD', this%memoryPath) + call mem_allocate(this%swrefl, 0, 'SWREFL', this%memoryPath) + call mem_allocate(this%rh, 0, 'RH', this%memoryPath) + call mem_allocate(this%atmc, 0, 'ATMC', this%memoryPath) + call mem_allocate(this%patm, 0, 'PATM', this%memoryPath) + call mem_allocate(this%ea, 0, 'EA', this%memoryPath) + call mem_allocate(this%es, 0, 'ES', this%memoryPath) + call mem_allocate(this%ew, 0, 'EW', this%memoryPath) + ! + ! -- reallocate abc variables based on which calculations are used + if (this%shf_active .or. this%lhf_active) then + call mem_reallocate(this%wspd, this%ncv, 'WSPD', this%memoryPath) + call mem_reallocate(this%ew, this%ncv, 'ES', this%memoryPath) + do n = 1, this%ncv + this%wspd(n) = DZERO + this%ew(n) = DZERO + end do + end if + ! + if (this%shf_active .or. this%lhf_active .or. this%lwr_active) then + call mem_reallocate(this%tatm, this%ncv, 'TATM', this%memoryPath) + call mem_reallocate(this%ea, this%ncv, 'EA', this%memoryPath) + call mem_reallocate(this%es, this%ncv, 'ES', this%memoryPath) + call mem_reallocate(this%ew, this%ncv, 'EW', this%memoryPath) + do n = 1, this%ncv + this%tatm(n) = DZERO + this%ea(n) = DZERO + this%es(n) = DZERO + this%ew(n) = DZERO + end do + end if + if (this%swr_active) then + call mem_reallocate(this%solr, this%ncv, 'SOLR', this%memoryPath) + call mem_reallocate(this%swrefl, this%ncv, 'SWREFL', this%memoryPath) + do n = 1, this%ncv + this%solr(n) = DZERO + this%swrefl(n) = DZERO + end do + end if + if (this%swr_active .or. this%lwr_active) then + call mem_reallocate(this%shd, this%ncv, 'SHD', this%memoryPath) + do n = 1, this%ncv + this%shd(n) = DZERO + end do + end if + if (this%lhf_active .or. this%lwr_active) then + call mem_reallocate(this%rh, this%ncv, 'RH', this%memoryPath) + do n = 1, this%ncv + this%rh(n) = DZERO + end do + end if + if (this%lwr_active) then + call mem_reallocate(this%atmc, this%ncv, 'ATMC', this%memoryPath) + do n = 1, this%ncv + this%atmc(n) = DONE + end do + end if + if (this%shf_active) then + call mem_reallocate(this%patm, this%ncv, 'PATM', this%memoryPath) + do n = 1, this%ncv + this%patm(n) = DZERO + end do + end if + ! + ! -- call utility ar routines if active + if (this%swr_active) then + call this%swr%pbst_ar() + end if + if (this%lwr_active) then + call this%lwr%pbst_ar() + end if + if (this%lhf_active) then + call this%lhf%pbst_ar() + end if + if (this%shf_active) then + call this%shf%pbst_ar() + end if + end subroutine abc_allocate_arrays + + !> @brief Deallocate memory + !< + subroutine abc_da(this) + ! -- modules + use MemoryManagerModule, only: mem_deallocate + ! -- dummy + class(AbcType) :: this + ! + ! -- SHF (sensible heat flux) + if (this%shf_active) then + call this%shf%da() + deallocate (this%shf) + end if + ! -- SWR (shortwave radiation heat flux) + if (this%swr_active) then + call this%swr%da() + deallocate (this%swr) + end if + ! -- LHF (latent heat flux) + if (this%lhf_active) then + call this%lhf%da() + deallocate (this%lhf) + end if + ! -- LWR (longwave radiation heat flux) + if (this%lwr_active) then + call this%lwr%da() + deallocate (this%lwr) + end if + ! + ! -- Deallocate scalars + call mem_deallocate(this%ncv) + call mem_deallocate(this%shf_active) + call mem_deallocate(this%swr_active) + call mem_deallocate(this%lhf_active) + call mem_deallocate(this%lwr_active) + call mem_deallocate(this%rhoa) + call mem_deallocate(this%cpa) + call mem_deallocate(this%cd) + call mem_deallocate(this%wfslope) + call mem_deallocate(this%wfint) + call mem_deallocate(this%lwrefl) + call mem_deallocate(this%emissw) + call mem_deallocate(this%emissr) + call mem_deallocate(this%tfac) + call mem_deallocate(this%toff) + call mem_deallocate(this%dfac) + call mem_deallocate(this%pfac) + ! + ! -- Deallocate time series manager + deallocate (this%tsmanager) + ! + ! -- Deallocate arrays + call mem_deallocate(this%status) + call mem_deallocate(this%iboundpbst) + call mem_deallocate(this%wspd) + call mem_deallocate(this%tatm) + call mem_deallocate(this%shd) + call mem_deallocate(this%swrefl) + call mem_deallocate(this%solr) + call mem_deallocate(this%rh) + call mem_deallocate(this%atmc) + call mem_deallocate(this%ea) + call mem_deallocate(this%es) + call mem_deallocate(this%ew) + call mem_deallocate(this%patm) + ! + ! -- Deallocate scalars in TspAptType + call this%NumericalPackageType%da() ! this may not work -- revisit and cleanup !!! + end subroutine abc_da + +! !> @brief Calculate observation value and pass it back to APT +! !< +! subroutine abc_bd_obs(this, obstypeid, jj, v, found) +! ! -- dummy +! class(AbcType), intent(inout) :: this +! character(len=*), intent(in) :: obstypeid +! real(DP), intent(inout) :: v +! integer(I4B), intent(in) :: jj +! logical, intent(inout) :: found +! ! -- local +! !integer(I4B) :: n1, n2 +! ! +! found = .true. +! !select case (obstypeid) +! !case ('SHF') +! ! if (this%iboundpak(jj) /= 0) then +! ! call this%sfe_shf_term(jj, n1, n2, v) +! ! end if +! !case ('SWR') +! ! if (this%iboundpak(jj) /= 0) then +! ! call this%sfe_swr_term(jj, n1, n2, v) +! ! end if +! !case default +! ! found = .false. +! !end select +! end subroutine abc_bd_obs + + !> @brief Calculate Atmospheric-Stream Heat Flux + !! + !! Calculate and return the atmospheric heat flux for one reach + !< + subroutine abc_cq(this, ifno, tstrm, abcflx, obstype) + ! -- dummy + class(AbcType), intent(inout) :: this + integer(I4B), intent(in) :: ifno !< stream reach integer id + real(DP), intent(in) :: tstrm !< temperature of the stream reach + real(DP), intent(inout) :: abcflx !< calculated atmospheric boundary flux amount + character(len=*), optional, intent(in) :: obstype !< when present, subroutine will return a specific energy flux for an observation + ! -- local + real(DP) :: swrflx = DZERO + real(DP) :: lwrflx = DZERO + real(DP) :: lhfflx = DZERO + real(DP) :: shfflx = DZERO + ! + ! -- update shared variables + call this%recalc_shared_vars(ifno, tstrm) + ! + ! -- calculate shortwave radiation using HGS equation + if (this%swr_active) then + call this%swr%swr_cq(ifno, swrflx) + end if + ! + ! -- calculate longwave radiation + if (this%lwr_active) then + call this%lwr%lwr_cq(ifno, tstrm, this%tfac, this%toff, lwrflx) + end if + ! + ! -- calculate latent heat flux using Dalton-like mass transfer equation + if (this%lhf_active) then + call this%lhf%lhf_cq(ifno, tstrm, this%gwecommon%gwerhow, this%dfac, this%tfac, this%toff, lhfflx) + end if + ! + ! -- calculate sensible heat flux using HGS equation + if (this%shf_active .and. .not. this%lhf_active) then + call this%shf%shf_cq(ifno, tstrm, this%tfac, this%toff, this%pfac, shfflx) ! default to HGS eqn method ("1") + else if (this%shf_active .and. this%lhf_active) then + call this%shf%shf_cq(ifno, tstrm, this%tfac, this%toff, shfflx, lhfflx) ! use Bowen ratio method ("2") + end if + ! + if (present(obstype)) then + select case (obstype) + case ('swr') + abcflx = swrflx + case ('lwr') + abcflx = lwrflx + case ('lhf') + abcflx = lhfflx + case ('shf') + abcflx = shfflx + case default + errmsg = 'Unrecognized observation type "'// & + trim(obstype)//'" for '// & + trim(adjustl(this%text))//' utility.' + call store_error(errmsg, terminate=.TRUE.) + end select + else + abcflx = swrflx + lwrflx + shfflx - lhfflx + end if + end subroutine abc_cq + + !> @brief Recalculate variables that are used by various heat fluxes + !! + !! For parameters like ambient vapor pressure of the atmosphere that are + !! used in the calculation of multiple heat fluxes, in this case longwave, + !! latent, and sensible heat flux, need to update the value held in memory + !< + subroutine recalc_shared_vars(this, ifno, tstrm) + ! -- dummy + class(AbcType), intent(inout) :: this + integer(I4B), intent(in) :: ifno !< reach id + real(DP), intent(in) :: tstrm !< temperature of the stream reach + ! -- local + real(DP) :: tstrmK + real(DP) :: t_inputK + ! + ! -- calculate saturation vapor pressure at air temperature + t_inputK = this%tfac * this%tatm(ifno) + this%toff + this%es(ifno) = calc_sat_vap_pres(t_inputK) + ! + ! -- calculate ambient vapor pressure at the atmospheric temperature + this%ea(ifno) = this%calc_eatm(ifno) + ! + ! -- calculate saturation vapor pressure at the water temperature + tstrmK = this%tfac * tstrm + this%toff + this%ew(ifno) = calc_sat_vap_pres(tstrmK) + end subroutine recalc_shared_vars + + !> @brief Calculate saturated vapor pressure for a given temperature + !! + !! A function for calculating the saturated vapor pressure given a + !! temperature, commonly either the stream temperature or atmospheric + !! temperature. + !< + function calc_sat_vap_pres(temp) result(e) + ! -- dummy + real(DP), intent(in) :: temp + ! -- return + real(DP) :: e + ! + e = 6.1275_DP * exp(17.2693882_DP * & + ((temp - DCTOK) / (temp - 35.86_DP))) + end function calc_sat_vap_pres + + !> @brief Calculate ambient vapor pressure + !! + !! Calculate ambient vapor pressure of the atmosphere as a function of + !! relative humidity and saturation vapor pressure of the atmosphere + !< + function calc_eatm(this, ifno) result(eatm) + ! -- dummy + class(AbcType) :: this + integer(I4B) :: ifno + ! -- return + real(DP) :: eatm + ! + eatm = this%rh(ifno) / DHUNDRED * this%es(ifno) + end function calc_eatm + + !> @brief Set the stress period attributes based on the keyword + !< + subroutine abc_set_stressperiod(this, itemno) !, keyword, found) MAKING LIKE PBST SET STRESSPERIOD + ! -- module + use TimeSeriesManagerModule, only: read_value_or_time_series_adv + ! -- dummy + class(AbcType), intent(inout) :: this + integer(I4B), intent(in) :: itemno + !logical, intent(inout) :: found + ! -- local + character(len=LINELENGTH) :: text + character(len=LINELENGTH) :: keyword + integer(I4B) :: ierr + integer(I4B) :: jj + real(DP), pointer :: bndElem => null() + ! + ! WIND SPEED + ! TEMPERATURE OF THE ATMOSPHERE + ! SHADE + ! REFLECTANCE OF SHORTWAVE RADIATION OFF WATER SURFACE + ! SOLAR RADIATION + ! RELATIVE HUMIDITY + ! ATMOSPHERIC COMPOSITION + ! ATMOSPHERIC PRESSURE + ! + ! -- read line + call this%parser%GetStringCaps(keyword) + select case (keyword) + case ('STATUS') + ierr = this%abc_check_valid(itemno) + if (ierr /= 0) then + goto 999 + end if + call this%parser%GetStringCaps(text) + this%status(itemno) = text(1:8) + if (text == 'CONSTANT') then + this%iboundpbst(itemno) = -1 + else if (text == 'INACTIVE') then + this%iboundpbst(itemno) = 0 + else if (text == 'ACTIVE') then + this%iboundpbst(itemno) = 1 + else + write (errmsg, '(a,a)') & + 'Unknown '//trim(this%text)//' status keyword: ', text//'.' + call store_error(errmsg) + end if + case ('WSPD') + ierr = this%abc_check_valid(itemno) + if (ierr /= 0) then + goto 999 + end if + call this%parser%GetString(text) + jj = 1 + bndElem => this%wspd(itemno) + call read_value_or_time_series_adv(text, itemno, jj, bndElem, & + this%packName, 'BND', this%tsManager, & + this%iprpak, 'WSPD') + case ('TATM') + ierr = this%abc_check_valid(itemno) + if (ierr /= 0) then + goto 999 + end if + call this%parser%GetString(text) + jj = 1 + bndElem => this%tatm(itemno) + call read_value_or_time_series_adv(text, itemno, jj, bndElem, & + this%packName, 'BND', this%tsManager, & + this%iprpak, 'TATM') + case ('SHD') + ierr = this%abc_check_valid(itemno) + if (ierr /= 0) then + goto 999 + end if + call this%parser%GetString(text) + jj = 1 + bndElem => this%shd(itemno) + call read_value_or_time_series_adv(text, itemno, jj, bndElem, & + this%packName, 'BND', this%tsManager, & + this%iprpak, 'SHD') + case ('SWREFL') + ierr = this%abc_check_valid(itemno) + if (ierr /= 0) then + goto 999 + end if + call this%parser%GetString(text) + jj = 1 + bndElem => this%swrefl(itemno) + call read_value_or_time_series_adv(text, itemno, jj, bndElem, & + this%packName, 'BND', this%tsManager, & + this%iprpak, 'SWREFL') + case ('SOLR') + ierr = this%abc_check_valid(itemno) + if (ierr /= 0) then + goto 999 + end if + call this%parser%GetString(text) + jj = 1 + bndElem => this%solr(itemno) + call read_value_or_time_series_adv(text, itemno, jj, bndElem, & + this%packName, 'BND', this%tsManager, & + this%iprpak, 'SOLR') + case ('RH') + ierr = this%abc_check_valid(itemno) + if (ierr /= 0) then + goto 999 + end if + call this%parser%GetString(text) + jj = 1 + bndElem => this%rh(itemno) + call read_value_or_time_series_adv(text, itemno, jj, bndElem, & + this%packName, 'BND', this%tsManager, & + this%iprpak, 'RH') + case ('ATMC') + ierr = this%abc_check_valid(itemno) + if (ierr /= 0) then + goto 999 + end if + call this%parser%GetString(text) + jj = 1 + bndElem => this%atmc(itemno) + call read_value_or_time_series_adv(text, itemno, jj, bndElem, & + this%packName, 'BND', this%tsManager, & + this%iprpak, 'ATMC') + case ('PATM') + ierr = this%abc_check_valid(itemno) + if (ierr /= 0) then + goto 999 + end if + call this%parser%GetString(text) + jj = 1 + bndElem => this%patm(itemno) + call read_value_or_time_series_adv(text, itemno, jj, bndElem, & + this%packName, 'BND', this%tsManager, & + this%iprpak, 'PATM') + case default + ! + ! -- Keyword not recognized so return to caller with found = .false. + !found = .false. + end select + ! +999 continue + end subroutine abc_set_stressperiod + + !> @brief Determine if a valid feature number has been specified. + !< + function abc_check_valid(this, itemno) result(ierr) + ! -- return + integer(I4B) :: ierr + ! -- dummy + class(AbcType), intent(inout) :: this + integer(I4B), intent(in) :: itemno + ! + ierr = 0 + if (itemno < 1 .or. itemno > this%ncv) then + write (errmsg, '(a,1x,i6,1x,a,1x,i6)') & + 'Featureno ', itemno, 'must be > 0 and <= ', this%ncv + call store_error(errmsg) + ierr = 1 + end if + end function abc_check_valid + + !> @brief Report calculated evaporation rate + !! + !! Allows for evaporation rate from water surface to be reported to an + !! observation output file + !< + subroutine abc_evap(this, ifno, tstrm, evap) + ! -- dummy + class(AbcType) :: this + integer(I4B), intent(in) :: ifno + real(DP), intent(in) :: tstrm !< calculated stream temperature + ! -- return + real(DP), intent(inout) :: evap !< calculated evaporation + ! + ! -- ensure shared variables are updated + call this%recalc_shared_vars(ifno, tstrm) + ! + ! -- calculate the evaporation rate for the latest values + evap = this%lhf%evap(this%wfint, this%wfslope, this%wspd(ifno), & + this%ew(ifno), this%ea(ifno)) + end subroutine abc_evap + +end module AbcModule diff --git a/src/Model/GroundWaterEnergy/gwe-lhf.f90 b/src/Model/GroundWaterEnergy/gwe-lhf.f90 new file mode 100644 index 00000000000..b6e03239655 --- /dev/null +++ b/src/Model/GroundWaterEnergy/gwe-lhf.f90 @@ -0,0 +1,146 @@ +!> @brief This module contains methods for calculating latent heat flux +!! +!! This module contains the methods used to calculate the latent heat flux +!! for surface-water boundaries, like streams and lakes. In its current form, +!! this class acts like a package to a package, similar to the TVK package that +!! can be invoked from the NPF package. Once this package is completed in its +!! prototyped form, it will likely be moved around. +!< + +module LatHeatModule + use ConstantsModule, only: LINELENGTH, LENMEMPATH, DZERO, LENVARNAME, DCTOK + use KindModule, only: I4B, DP + use MemoryManagerModule, only: mem_setptr + use MemoryHelperModule, only: create_mem_path + use TimeSeriesLinkModule, only: TimeSeriesLinkType + use TimeSeriesManagerModule, only: TimeSeriesManagerType, tsmanager_cr + use SimModule, only: store_error + use SimVariablesModule, only: errmsg + use PbstBaseModule, only: PbstBaseType, pbstbase_da + + implicit none + + public :: LhfType + public :: lhf_cr + + character(len=16) :: text = ' LHF' + + type, extends(PbstBaseType) :: LhfType + + real(DP), pointer :: wfslope => null() !< wind function slope + real(DP), pointer :: wfint => null() !< wind function intercept + real(DP), dimension(:), pointer, contiguous :: wspd => null() !< wind speed + real(DP), dimension(:), pointer, contiguous :: tatm => null() !< temperature of the atmosphere + real(DP), dimension(:), pointer, contiguous :: ea => null() !< saturation vapor pressure at air temperature + real(DP), dimension(:), pointer, contiguous :: ew => null() !< saturation vapor pressure at water temperature + + contains + + procedure :: da => lhf_da + procedure :: pbst_ar => lhf_ar_set_pointers + procedure, public :: lhf_cq + + end type LhfType + +contains + + !> @brief Create a new LhfType object + !! + !! Create a new latent heat flux (LhfType) object. Initially for use with + !! the SFE-ABC package. + !< + subroutine lhf_cr(this, name_model, inunit, iout, ncv) + ! -- dummy + type(LhfType), pointer, intent(out) :: this + character(len=*), intent(in) :: name_model + integer(I4B), intent(in) :: inunit + integer(I4B), intent(in) :: iout + integer(I4B), target, intent(in) :: ncv + ! + allocate (this) + call this%init(name_model, 'LHF', 'LHF', inunit, iout, ncv) + this%text = text + ! + end subroutine lhf_cr + + !> @brief Announce package and set pointers to variables + !! + !! Announce package version and set array and variable pointers from the ABC + !! package for access by LHF. + !< + subroutine lhf_ar_set_pointers(this) + ! -- dummy + class(LhfType) :: this + ! -- local + character(len=LENMEMPATH) :: abcMemoryPath + ! + ! -- print a message noting that the SHF utility is active + write (this%iout, '(a)') & + 'LHF -- LATENT HEAT WILL BE INCLUDED IN THE ATMOSPHERIC BOUNDARY '// & + 'CONDITIONS FOR THE STREAMFLOW ENERGY TRANSPORT PACKAGE' + ! + ! -- set pointers to variables hosted in the ABC package + abcMemoryPath = create_mem_path(this%name_model, 'ABC') + call mem_setptr(this%wfslope, 'WFSLOPE', abcMemoryPath) + call mem_setptr(this%wfint, 'WFINT', abcMemoryPath) + call mem_setptr(this%wspd, 'WSPD', abcMemoryPath) + call mem_setptr(this%tatm, 'TATM', abcMemoryPath) + call mem_setptr(this%ea, 'EA', abcMemoryPath) + call mem_setptr(this%ew, 'EW', abcMemoryPath) + ! + ! -- create time series manager + call tsmanager_cr(this%tsmanager, this%iout, & + removeTsLinksOnCompletion=.true., & + extendTsToEndOfSimulation=.true.) + end subroutine lhf_ar_set_pointers + + !> @brief Calculate Latent Heat Flux + !! + !! Calculate and return the latent heat flux for one reach + !< + subroutine lhf_cq(this, ifno, tstrm, rhow, dfac, tfac, toff, lhflx) + ! -- dummy + class(LhfType), intent(inout) :: this + integer(I4B), intent(in) :: ifno !< stream reach integer id + real(DP), intent(in) :: tstrm !< temperature of the stream reach + real(DP), intent(in) :: rhow !< density of water + real(DP), intent(in) :: dfac !< density units adjustment factor + real(DP), intent(in) :: tfac !< temperature units adjustment factor + real(DP), intent(in) :: toff !< temperature units offset + real(DP), intent(inout) :: lhflx !< calculated latent heat flux amount + ! -- local + real(DP) :: l !< latent heat vaporization + real(DP) :: evap + ! + ! -- calculate latent heat of vaporization (water temperature dependent) Eq. A.16 + l = 2499.64_DP - (2.51_DP * (tfac * tstrm + toff) - DCTOK) + ! + ! -- mass-transfer method for calculating evap rate (A.17) + evap = this%evap(this%wfint, this%wfslope, this%wspd(ifno), & + this%ew(ifno), this%ea(ifno)) + ! + ! -- calculate latent heat flux (A.15) + lhflx = evap * l * rhow * dfac + end subroutine lhf_cq + + !> @brief Deallocate package memory + !! + !! Deallocate TVK package scalars and arrays. + !< + subroutine lhf_da(this) + ! -- dummy + class(LhfType) :: this + ! + ! -- nullify pointers to other package variables + nullify (this%wfint) + nullify (this%wfslope) + nullify (this%wspd) + nullify (this%tatm) + nullify (this%ea) + nullify (this%ew) + ! + ! -- deallocate parent + call pbstbase_da(this) + end subroutine lhf_da + +end module LatHeatModule diff --git a/src/Model/GroundWaterEnergy/gwe-lwr.f90 b/src/Model/GroundWaterEnergy/gwe-lwr.f90 new file mode 100644 index 00000000000..d6ea46456f0 --- /dev/null +++ b/src/Model/GroundWaterEnergy/gwe-lwr.f90 @@ -0,0 +1,187 @@ +!> @brief This module contains methods for calculating the longwave radiation heat flux +!! +!! This module contains the methods used to calculate the longwave radiation +!! heat flux for surface-water boundaries, like streams and lakes. In its +!! current form, this class acts like a package to a package, similar to the +!! TVK package (or utiliity) that can be invoked from the NPF package. Once +!! this package is completed in its prototyped form, it will likely be moved +!! around. +!< + +module LongwaveModule + + use ConstantsModule, only: LINELENGTH, LENMEMPATH, DZERO, LENVARNAME, & + DSTEFANBOLTZMANN, DFOUR + use KindModule, only: I4B, DP + use MemoryManagerModule, only: mem_setptr + use MemoryHelperModule, only: create_mem_path + use TimeSeriesLinkModule, only: TimeSeriesLinkType + use TimeSeriesManagerModule, only: TimeSeriesManagerType, tsmanager_cr + use SimModule, only: store_error + use SimVariablesModule, only: errmsg + use PbstBaseModule !, only: PbstBaseType, pbstbase_da + + implicit none + + public :: LwrType + public :: lwr_cr + + character(len=16) :: text = ' LWR' + + type, extends(PbstBaseType) :: LwrType + + real(DP), pointer :: lwrefl => null() !< longwave reflectance of water surface + real(DP), pointer :: emissw => null() !< emissivity of water + real(DP), pointer :: emissr => null() !< emissivity of riparian canopy + real(DP), dimension(:), pointer, contiguous :: shd => null() !< shade fraction + real(DP), dimension(:), pointer, contiguous :: atmc => null() !< atmospheric composition adjustment + real(DP), dimension(:), pointer, contiguous :: tatm => null() !< temperature of the atmosphere + real(DP), dimension(:), pointer, contiguous :: rh => null() !< relative humidity + real(DP), dimension(:), pointer, contiguous :: ea => null() !< ambient vapor pressure of the atmosphere + real(DP), dimension(:), pointer, contiguous :: es => null() !< saturation vapor pressure at air temperature + real(DP), dimension(:), pointer, contiguous :: ew => null() !< saturation vapor pressure at water temperature + + contains + + procedure :: da => lwr_da + procedure :: pbst_ar => lwr_ar_set_pointers + procedure, public :: lwr_cq + procedure, private :: calc_lwr + + end type LwrType + +contains + + !> @brief Create a new LwrType object + !! + !! Create a new longwave radiation flux (LwrType) object. Initially for use with + !! the SFE package. + !< + subroutine lwr_cr(this, name_model, inunit, iout, ncv) + ! -- dummy + type(LwrType), pointer, intent(out) :: this + character(len=*), intent(in) :: name_model + integer(I4B), intent(in) :: inunit + integer(I4B), intent(in) :: iout + integer(I4B), target, intent(in) :: ncv + ! + allocate (this) + call this%init(name_model, 'LWR', 'LWR', inunit, iout, ncv) + this%text = text + end subroutine lwr_cr + + !> @brief Announce package and set pointers to variables + !! + !! Announce package version and set array and variable pointers from the ABC + !! package for access by LWR. + !< + subroutine lwr_ar_set_pointers(this) + ! -- dummy + class(LwrType) :: this + ! -- local + character(len=LENMEMPATH) :: abcMemoryPath + ! + ! -- print a message noting that the LWR utility is active + write (this%iout, '(a)') & + 'LWR -- LONGWAVE RADIATION WILL BE INCLUDED IN THE ATMOSPHERIC '// & + 'BOUNDARY CONDITIONS FOR THE STREAMFLOW ENERGY TRANSPORT PACKAGE' + ! + ! -- set pointers to variables hosted in the ABC package + abcMemoryPath = create_mem_path(this%name_model, 'ABC') + call mem_setptr(this%lwrefl, 'LWREFL', abcMemoryPath) + call mem_setptr(this%shd, 'SHD', abcMemoryPath) + call mem_setptr(this%atmc, 'ATMC', abcMemoryPath) + call mem_setptr(this%tatm, 'TATM', abcMemoryPath) + call mem_setptr(this%emissw, 'EMISSW', abcMemoryPath) + call mem_setptr(this%emissr, 'EMISSR', abcMemoryPath) + call mem_setptr(this%rh, 'RH', abcMemoryPath) + call mem_setptr(this%ea, 'EA', abcMemoryPath) + call mem_setptr(this%es, 'ES', abcMemoryPath) + call mem_setptr(this%ew, 'EW', abcMemoryPath) + ! + ! -- set the riparian canopy emissivity constant + this%emissr = 0.97_DP + ! + ! -- create time series manager + call tsmanager_cr(this%tsmanager, this%iout, & + removeTsLinksOnCompletion=.true., & + extendTsToEndOfSimulation=.true.) + end subroutine lwr_ar_set_pointers + + !> @brief Calculate Longwave Radiation Heat Flux + !! + !! Calculate and return the longwave radiation heat flux for one reach + !< + subroutine lwr_cq(this, ifno, tstrm, tfac, toff, lwrflx) + ! -- dummy + class(LwrType), intent(inout) :: this + integer(I4B), intent(in) :: ifno !< stream reach integer id + real(DP), intent(in) :: tstrm !< temperature of the stream reach + real(DP), intent(in) :: tfac !< temperature units adjustment factor + real(DP), intent(in) :: toff !< temperature units offset + real(DP), intent(inout) :: lwrflx !< calculated longwave radiation heat flux amount + ! -- local + real(DP) :: emissa + real(DP) :: emisss + real(DP) :: lwratm + real(DP) :: lwrstrm + ! + ! -- intermediate calculations + ! + ! -- atmospheric emissivity (A.14) + emissa = this%epsa(this%ea(ifno), tfac * this%tatm(ifno) + toff, this%atmc(ifno)) + ! + ! -- shade-altered above-channel emissivity [Eq. 3, Fogg et al. (2023)] + emisss = this%epss(this%shd(ifno), emissa, this%emissr) + ! + ! -- long wave radiation transmitted from the atmosphere to the water surface (A.12) + lwratm = this%calc_lwr(emisss, tfac * this%tatm(ifno) + toff) + ! + ! -- long wave radiation transmitted from water surface to the atmosphere (A.13) + lwrstrm = this%calc_lwr(this%emissw, tfac * tstrm + toff) + ! + ! -- longwave radiation heat flux + lwrflx = lwratm * (1 - this%lwrefl) - lwrstrm + end subroutine lwr_cq + + !> @brief Deallocate package memory + !! + !! Deallocate TVK package scalars and arrays. + !< + subroutine lwr_da(this) + ! -- dummy + class(LwrType) :: this + ! + ! -- deallocate time series + nullify (this%shd) + nullify (this%lwrefl) + nullify (this%tatm) + nullify (this%emissw) + nullify (this%emissr) + nullify (this%atmc) + nullify (this%rh) + nullify (this%ea) + nullify (this%es) + nullify (this%ew) + ! + ! -- deallocate parent + call pbstbase_da(this) + end subroutine lwr_da + + !> @brief Calculate longwave radiation + !! + !! Function for calculating incoming or outgoing longwave radiation + !< + function calc_lwr(this, eps, temp) result(lwr) + ! -- dummy + class(LwrType) :: this + real(DP) :: eps !< epsilon, representing either emissivity of the atmosphere, the shade-weighted emissivity of the atm, or emissivity of water + real(DP) :: temp !< will be temperature in Kelvin, representing either the temperature of the stream or the atmosphere in degree C + ! -- return + real(DP) :: lwr !< longwave radiation + ! + ! -- generalized equation for calculating longwave radiation + lwr = eps * DSTEFANBOLTZMANN * ((temp)**DFOUR) + end function calc_lwr + +end module LongwaveModule diff --git a/src/Model/GroundWaterEnergy/gwe-sfe.f90 b/src/Model/GroundWaterEnergy/gwe-sfe.f90 index 923ac4db1b4..d95d10e424d 100644 --- a/src/Model/GroundWaterEnergy/gwe-sfe.f90 +++ b/src/Model/GroundWaterEnergy/gwe-sfe.f90 @@ -20,8 +20,11 @@ ! EXT-INFLOW idxbudiflw EXT-INFLOW rhow * cpw * q * tiflw ! EXT-OUTFLOW idxbudoutf EXT-OUTFLOW rhow * cpw * q * tfeat -! -- terms not associated with SFR -! STRMBD-COND none STRMBD-COND ctherm * (tcell - tfeat) (ctherm is a thermal conductance for the streambed) +! -- SFE terms +! STRMBD-COND idxbudsbcd STRMBD-COND ktf * wa / sbthk * (t_cell - t_feat) +! ATM BC idxbudabc ATMOSPHERE swr + lwr + shf - lhf +! SENSIBLE HEAT FLUX idxbudshf SENS HEAT cd * rho_a * C_p_a * wspd * (t_air - t_feat) +! SHORTWAVE RADIATION idxbudswr SHORTWAVE (1 - shd) * (1 - swrefl) * solr ! -- terms from a flow file that should be skipped ! CONSTANT none none none @@ -43,10 +46,13 @@ module GweSfeModule use TspFmiModule, only: TspFmiType use SfrModule, only: SfrType use ObserveModule, only: ObserveType + use BaseDisModule, only: DisBaseType use TspAptModule, only: TspAptType, apt_process_obsID, & apt_process_obsID12 use GweInputDataModule, only: GweInputDataType + use AbcModule use MatrixBaseModule + use InputOutputModule, only: openfile ! implicit none ! @@ -66,6 +72,9 @@ module GweSfeModule integer(I4B), pointer :: idxbudroff => null() !< index of runoff terms in flowbudptr integer(I4B), pointer :: idxbudiflw => null() !< index of inflow terms in flowbudptr integer(I4B), pointer :: idxbudoutf => null() !< index of outflow terms in flowbudptr + integer(I4B), pointer :: idxbudabc => null() !< index of atmospheric boundary condition terms in flowbudptr + + logical, pointer, public :: abc_active => null() !< logical indicating if an atmospheric boundary condition object is active real(DP), dimension(:), pointer, contiguous :: temprain => null() !< rainfall temperature real(DP), dimension(:), pointer, contiguous :: tempevap => null() !< evaporation temperature @@ -76,10 +85,17 @@ module GweSfeModule real(DP), dimension(:), pointer, contiguous :: ktf => null() !< thermal conductivity between the sfe and groundwater cell real(DP), dimension(:), pointer, contiguous :: rfeatthk => null() !< thickness of streambed material through which thermal conduction occurs + type(AbcType), pointer :: abc => null() ! atmospheric boundary condition (abc) object + + integer(I4B), pointer :: inabc => null() ! ABC (atmospheric boundary condition utility) unit number (0 if unused) + contains procedure :: bnd_ad => sfe_ad + !procedure :: bnd_df => sfe_df procedure :: bnd_da => sfe_da + !procedure :: bnd_ar => sfe_ar + procedure :: ancil_rp => utl_rp procedure :: allocate_scalars procedure :: apt_allocate_arrays => sfe_allocate_arrays procedure :: find_apt_package => find_sfe_package @@ -94,12 +110,14 @@ module GweSfeModule procedure :: sfe_iflw_term procedure :: sfe_outf_term procedure, private :: sfe_sbcd_term + procedure, private :: sfe_abc_term procedure :: pak_df_obs => sfe_df_obs procedure :: pak_rp_obs => sfe_rp_obs procedure :: pak_bd_obs => sfe_bd_obs procedure :: pak_set_stressperiod => sfe_set_stressperiod procedure :: apt_get_volumes => sfe_get_volumes procedure :: apt_read_cvs => sfe_read_cvs + procedure :: gc_options => sfe_options end type GweSfeType @@ -109,6 +127,8 @@ module GweSfeModule !< subroutine sfe_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, & fmi, eqnsclfac, gwecommon, dvt, dvu, dvua) + ! -- modules + use InputOutputModule, only: getunit ! -- dummy class(BndType), pointer :: packobj integer(I4B), intent(in) :: id @@ -166,6 +186,69 @@ subroutine sfe_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, & sfeobj%depvarunitabbrev = dvua end subroutine sfe_create + !> @brief Set options specific to the GweSfeType + !! + !! This routine overrides TspAptType%gc_options + !< + subroutine sfe_options(this, option, found) + ! -- modules + use ConstantsModule, only: MAXCHARLEN, LGP + use InputOutputModule, only: urword, getunit, assign_iounit, openfile + ! -- dummy + class(GweSfeType), intent(inout) :: this + character(len=*), intent(inout) :: option + logical, intent(inout) :: found + ! -- local + character(len=LINELENGTH) :: fname + character(len=MAXCHARLEN) :: keyword + ! -- formats + character(len=*), parameter :: fmtaptbin = & + "(4x, a, 1x, a, 1x, ' WILL BE SAVED TO FILE: ', a, & + &/4x, 'OPENED ON UNIT: ', I0)" + ! + found = .true. + select case (option) + case ('ABC6') + ! + call this%parser%GetStringCaps(keyword) + if (trim(adjustl(keyword)) /= 'FILEIN') then + errmsg = 'ABC6 keyword must be followed by "FILEIN" '// & + 'then by filename.' + call store_error(errmsg) + call this%parser%StoreErrorUnit() + end if + if (this%abc_active) then + errmsg = 'Multiple ABC6 keywords detected in OPTIONS block. '// & + 'Only one ABC6 entry allowed for a package.' + call store_error(errmsg) + end if + this%abc_active = .true. + call this%parser%GetString(fname) + ! + ! -- create atmospheric boundary condition object + call openfile(this%inabc, this%iout, fname, 'ABC') + call abc_cr(this%abc, this%name_model, this%inabc, this%iout, fname, & + this%ncv, this%gwecommon, this%numericalpackagetype%dis) !, this%mempath) + call this%abc%read_options() + call this%abc%abc_df() + !this%abc%inputFilename = fname + ! + ! -- call _ar routine for abc sub-package + ! note: this is the best place to call the sub-package (abc) for + ! now because SFE does not override apt_ar(), and since sfe + ! does not run its own customized set of _ar() functionality + ! there is no opportunity to run abc%ar() from a more + ! natural point. In other words, one place to call abc%ar + ! would be from bnd_ar(), however, abc is only callable from + ! sfe and not apt. + !call this%abc%ar() + case default + ! + ! -- No options found + found = .false. + end select + end subroutine sfe_options + !> @brief Find corresponding sfe package !< subroutine find_sfe_package(this) @@ -236,7 +319,9 @@ subroutine find_sfe_package(this) do ip = 1, this%flowbudptr%nbudterm select case (trim(adjustl(this%flowbudptr%budterm(ip)%flowtype))) case ('FLOW-JA-FACE') - this%idxbudfjf = ip + if (this%flowbudptr%budterm(ip)%maxlist /= 0) then + this%idxbudfjf = ip + end if this%idxbudssm(ip) = 0 case ('GWF') this%idxbudgwf = ip @@ -280,6 +365,9 @@ subroutine find_sfe_package(this) ' MAX NO. OF ENTRIES = ', this%flowbudptr%budterm(ip)%maxlist end do write (this%iout, '(a, //)') 'DONE PROCESSING '//ftype//' INFORMATION' + ! + ! -- Experiment + !call this%abc%ar() end subroutine find_sfe_package !> @brief Add matrix terms related to SFE @@ -380,6 +468,17 @@ subroutine sfe_fc_expanded(this, rhs, ia, idxglo, matrix_sln) call matrix_sln%add_value_pos(ipossymoffd, hcofval) end if end do + ! + ! -- Add atmospheric heat flux contribution + if (this%inabc /= 0) then + do j = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist + call this%sfe_abc_term(j, n1, n2, rrate, rhsval, hcofval) + iloc = this%idxlocnode(n1) + iposd = this%idxpakdiag(n1) + call matrix_sln%add_value_pos(iposd, hcofval) + rhs(iloc) = rhs(iloc) + rhsval + end do + end if end subroutine sfe_fc_expanded !> @ brief Add terms specific to sfr to the explicit sfe solve @@ -433,6 +532,7 @@ subroutine sfe_solve(this) end if ! ! Note: explicit streambed conduction terms??? + ! Note: explicit sensible heat flux terms? end subroutine sfe_solve !> @brief Function to return the number of budget terms just for this package. @@ -451,8 +551,14 @@ function sfe_get_nbudterms(this) result(nbudterms) ! 3. runoff ! 4. ext-inflow ! 5. ext-outflow - ! 6. strmbd-cond + ! 6. strmbed-cond nbudterms = 6 + ! -- optional utilities + ! X. atmospheric boundary condition + if (this%inabc /= 0) then + nbudterms = nbudterms + 1 + end if + ! end function sfe_get_nbudterms !> @brief Set up the budget object that stores all the sfe flows @@ -553,6 +659,28 @@ subroutine sfe_setup_budobj(this, idx) n2 = this%flowbudptr%budterm(this%idxbudgwf)%id2(n) call this%budobj%budterm(idx)%update_term(n1, n2, q) end do + ! + ! -- Atmospheric heat flux + if (this%inabc /= 0) then + text = ' ATMOSPHERE' + idx = idx + 1 + maxlist = this%flowbudptr%budterm(this%idxbudgwf)%maxlist + naux = 0 + call this%budobj%budterm(idx)%initialize(text, & + this%name_model, & + this%packName, & + this%name_model, & + this%packName, & + maxlist, .false., .false., & + naux) + call this%budobj%budterm(idx)%reset(maxlist) + q = DZERO + do n = 1, maxlist + n1 = this%flowbudptr%budterm(this%idxbudgwf)%id1(n) + n2 = this%flowbudptr%budterm(this%idxbudgwf)%id2(n) + call this%budobj%budterm(idx)%update_term(n1, n2, q) + end do + end if end subroutine sfe_setup_budobj !> @brief Copy flow terms into this%budobj @@ -642,6 +770,18 @@ subroutine sfe_fill_budobj(this, idx, x, flowja, ccratin, ccratout) flowja(idiag) = flowja(idiag) - q end if end do + ! + ! -- Atmospheric boundary heat exchange + if (this%inabc /= 0) then + idx = idx + 1 + nlist = this%flowbudptr%budterm(this%idxbudgwf)%nlist + call this%budobj%budterm(idx)%reset(nlist) + do j = 1, nlist + call this%sfe_abc_term(j, n1, n2, q) + call this%budobj%budterm(idx)%update_term(n1, n2, q) + call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout) + end do + end if end subroutine sfe_fill_budobj !> @brief Allocate scalars specific to the streamflow energy transport (SFE) @@ -662,6 +802,9 @@ subroutine allocate_scalars(this) call mem_allocate(this%idxbudroff, 'IDXBUDROFF', this%memoryPath) call mem_allocate(this%idxbudiflw, 'IDXBUDIFLW', this%memoryPath) call mem_allocate(this%idxbudoutf, 'IDXBUDOUTF', this%memoryPath) + call mem_allocate(this%idxbudabc, 'IDXBUDABC', this%memoryPath) + call mem_allocate(this%abc_active, 'ABC_ACTIVE', this%memoryPath) + call mem_allocate(this%inabc, 'INABC', this%memoryPath) ! ! -- Initialize this%idxbudrain = 0 @@ -669,6 +812,10 @@ subroutine allocate_scalars(this) this%idxbudroff = 0 this%idxbudiflw = 0 this%idxbudoutf = 0 + this%idxbudabc = 0 + ! + this%abc_active = .false. + this%inabc = 0 end subroutine allocate_scalars !> @brief Allocate arrays specific to the streamflow energy transport (SFE) @@ -703,6 +850,11 @@ subroutine sfe_allocate_arrays(this) this%vnew(n) = DZERO this%vold(n) = DZERO end do + ! + ! -- Call sub-package(s) allocate arrays + if (this%inabc /= 0) then + call this%abc%abc_allocate_arrays() + end if end subroutine sfe_allocate_arrays !> @brief Advance sfe package routine @@ -723,6 +875,22 @@ subroutine sfe_ad(this) end subroutine sfe_ad + !> @brief Call Read and prepare routines for any active pbst subpackages + !! + !! Overrides ancil_rp() subroutine in tsp-apt. The idea being that for a + !! GWE model with add-on packages (like sensible heat flux, for example) + !! they can only be accessed from sfe and not apt. + !< + subroutine utl_rp(this) + ! -- dummy + class(GweSfeType), intent(inout) :: this + ! + ! -- call atmospheric boundary condition sub-package _rp() routine + if (this%inabc /= 0) then + call this%abc%abc_rp() + end if + end subroutine utl_rp + !> @brief Deallocate memory !< subroutine sfe_da(this) @@ -731,12 +899,22 @@ subroutine sfe_da(this) ! -- dummy class(GweSfeType) :: this ! + ! -- ABC (atmospheric boundary condition) + if (this%inabc /= 0) then + call this%abc%da() + deallocate (this%abc) + end if + ! ! -- Deallocate scalars call mem_deallocate(this%idxbudrain) call mem_deallocate(this%idxbudevap) call mem_deallocate(this%idxbudroff) call mem_deallocate(this%idxbudiflw) call mem_deallocate(this%idxbudoutf) + call mem_deallocate(this%idxbudabc) + ! + call mem_deallocate(this%abc_active) + call mem_deallocate(this%inabc) ! ! -- Deallocate time series call mem_deallocate(this%temprain) @@ -798,6 +976,13 @@ subroutine sfe_evap_term(this, ientry, n1, n2, rrate, rhsval, hcofval) n2 = this%flowbudptr%budterm(this%idxbudevap)%id2(ientry) ! -- note that qbnd is negative for evap qbnd = this%flowbudptr%budterm(this%idxbudevap)%flow(ientry) + if (qbnd /= DZERO .and. this%abc_active) then + write (errmsg, '(a,1x,i6)') & + 'ERROR. Specified evaporation in SFR is not compatible ABC latent & + &heat calculations for reach ', ientry + call store_error(errmsg) + call this%parser%StoreErrorUnit() + end if heatlat = this%gwecommon%gwerhow * this%gwecommon%gwelatheatvap if (present(rrate)) rrate = qbnd * heatlat if (present(rhsval)) rhsval = -rrate @@ -922,13 +1107,43 @@ subroutine sfe_sbcd_term(this, ientry, n1, igwfnode, rrate, rhsval, hcofval) if (present(hcofval)) hcofval = ctherm end subroutine sfe_sbcd_term + !> @brief Atmospheric Boundary Condition (ABC) term + !< + subroutine sfe_abc_term(this, ientry, n1, n2, rrate, rhsval, hcofval) + ! -- dummy + class(GweSfeType) :: this + integer(I4B), intent(in) :: ientry + integer(I4B), intent(inout) :: n1 + integer(I4B), intent(inout) :: n2 + real(DP), intent(inout), optional :: rrate + real(DP), intent(inout), optional :: rhsval + real(DP), intent(inout), optional :: hcofval + ! -- local + real(DP) :: atmheat + real(DP) :: strmtemp + integer(I4B) :: auxpos + real(DP) :: sa !< surface area of stream reach, different than wetted area + ! + n1 = this%flowbudptr%budterm(this%idxbudevap)%id1(ientry) + ! -- For now, there is only 1 aux variable under 'EVAPORATION' which + ! is reach surface area + auxpos = this%flowbudptr%budterm(this%idxbudevap)%naux + sa = this%flowbudptr%budterm(this%idxbudevap)%auxvar(auxpos, ientry) + ! + strmtemp = this%xnewpak(n1) + call this%abc%abc_cq(n1, strmtemp, atmheat) + ! + if (present(rrate)) rrate = atmheat * sa + if (present(rhsval)) rhsval = -rrate + if (present(hcofval)) hcofval = DZERO + end subroutine sfe_abc_term + !> @brief Observations !! !! Store the observation type supported by the APT package and override !! BndType%bnd_df_obs !< subroutine sfe_df_obs(this) - ! -- modules ! -- dummy class(GweSfeType) :: this ! -- local @@ -998,6 +1213,37 @@ subroutine sfe_df_obs(this) ! for strmbd-cond observation type. call this%obs%StoreObsType('strmbd-cond', .true., indx) this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID + ! + ! -- Store obs type and assign procedure pointer + ! for net atmospheric boundary condition (abc) flux observation type. + call this%obs%StoreObsType('abc', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID + ! + ! -- Store obs type and assign procedure pointer + ! for shortwave-radiation-flux observation type. + call this%obs%StoreObsType('swr', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID + ! + ! -- Store obs type and assign procedure pointer + ! for longwave radiation flux observation type. + call this%obs%StoreObsType('lwr', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID + ! + ! -- Store obs type and assign procedure pointer + ! for latent heat flux observation type. + call this%obs%StoreObsType('lhf', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID + ! + ! -- Store obs type and assign procedure pointer + ! for sensible heat flux observation type. + call this%obs%StoreObsType('shf', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID + ! + ! -- Store obs type and assign procedure pointer + ! for evaporation rate observation type. + call this%obs%StoreObsType('surfevap', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID + ! end subroutine sfe_df_obs !> @brief Process package specific obs @@ -1027,6 +1273,18 @@ subroutine sfe_rp_obs(this, obsrv, found) call this%rp_obs_byfeature(obsrv) case ('STRMBD-COND') call this%rp_obs_byfeature(obsrv) + case ('ABC') + call this%rp_obs_byfeature(obsrv) + case ('SWR') + call this%rp_obs_byfeature(obsrv) + case ('LWR') + call this%rp_obs_byfeature(obsrv) + case ('LHF') + call this%rp_obs_byfeature(obsrv) + case ('SHF') + call this%rp_obs_byfeature(obsrv) + case ('SURFEVAP') + call this%rp_obs_byfeature(obsrv) case default found = .false. end select @@ -1043,8 +1301,13 @@ subroutine sfe_bd_obs(this, obstypeid, jj, v, found) logical, intent(inout) :: found ! -- local integer(I4B) :: n1, n2 + real(DP) :: strmtemp ! found = .true. + ! + ! -- stream temperature used by observations associated with ABC utility + strmtemp = this%xnewpak(jj) + ! select case (obstypeid) case ('RAINFALL') if (this%iboundpak(jj) /= 0) then @@ -1070,6 +1333,34 @@ subroutine sfe_bd_obs(this, obstypeid, jj, v, found) if (this%iboundpak(jj) /= 0) then call this%sfe_sbcd_term(jj, n1, n2, v) end if + case ('ABC') + if (this%iboundpak(jj) /= 0) then + call this%sfe_abc_term(jj, n1, n2, v) + end if + case ('SWR') + if (this%iboundpak(jj) /= 0) then + !call this%swr_abc_term(jj, n1, n2, v) + call this%abc%abc_cq(jj, strmtemp, v, 'swr') + end if + case ('LWR') + if (this%iboundpak(jj) /= 0) then + !call this%lwr_abc_term(jj, n1, n2, v) + call this%abc%abc_cq(jj, strmtemp, v, 'lwr') + end if + case ('LHF') + if (this%iboundpak(jj) /= 0) then + !call this%lhf_abc_term(jj, n1, n2, v) + call this%abc%abc_cq(jj, strmtemp, v, 'lhf') + end if + case ('SHF') + if (this%iboundpak(jj) /= 0) then + !call this%shf_abc_term(jj, n1, n2, v) + call this%abc%abc_cq(jj, strmtemp, v, 'shf') + end if + case ('SURFEVAP') + if (this%iboundpak(jj) /= 0) then + call this%abc%abc_evap(jj, strmtemp, v) + end if case default found = .false. end select diff --git a/src/Model/GroundWaterEnergy/gwe-shf.f90 b/src/Model/GroundWaterEnergy/gwe-shf.f90 new file mode 100644 index 00000000000..610b7168c41 --- /dev/null +++ b/src/Model/GroundWaterEnergy/gwe-shf.f90 @@ -0,0 +1,154 @@ +!> @brief This module contains methods for calculating sensible heat flux +!! +!! This module contains the methods used to calculate the sensible heat flux +!! for surface-water boundaries, like streams and lakes. In its current form, +!! this class acts like a package to a package, similar to the TVK package that +!! can be invoked from the NPF package. Once this package is completed in its +!! prototyped form, it will likely be moved around. +!< + +module SensHeatModule + use ConstantsModule, only: LINELENGTH, LENMEMPATH, DZERO, LENVARNAME, & + DCTOK + use KindModule, only: I4B, DP + use MemoryManagerModule, only: mem_setptr + use MemoryHelperModule, only: create_mem_path + use TimeSeriesLinkModule, only: TimeSeriesLinkType + use TimeSeriesManagerModule, only: TimeSeriesManagerType, tsmanager_cr + use SimModule, only: store_error + use SimVariablesModule, only: errmsg + use PbstBaseModule, only: PbstBaseType, pbstbase_da + + implicit none + + private + + public :: ShfType + public :: shf_cr + + character(len=16) :: text = ' SHF' + + type, extends(PbstBaseType) :: ShfType + + real(DP), pointer :: rhoa => null() !< ABC density of air + real(DP), pointer :: cpa => null() !< ABC heat capacity of air + real(DP), pointer :: cd => null() !< ABC drag coefficient + real(DP), dimension(:), pointer, contiguous :: wspd => null() !< ABC wind speed + real(DP), dimension(:), pointer, contiguous :: tatm => null() !< ABC temperature of the atmosphere + real(DP), dimension(:), pointer, contiguous :: patm => null() !< ABC atmospheric pressure + real(DP), dimension(:), pointer, contiguous :: ea => null() !< ABC temperature of the atmosphere + real(DP), dimension(:), pointer, contiguous :: ew => null() !< ABC atmospheric pressure + + contains + + procedure :: da => shf_da + procedure :: pbst_ar => shf_ar_set_pointers + procedure, public :: shf_cq + + end type ShfType + +contains + + !> @brief Create a new ShfType object + !! + !! Create a new sensible heat flux (ShfType) object. Initially for use with + !! the SFE-ABC package. + !< + subroutine shf_cr(this, name_model, inunit, iout, ncv) + ! -- dummy + type(ShfType), pointer, intent(out) :: this + character(len=*), intent(in) :: name_model + integer(I4B), intent(in) :: inunit + integer(I4B), intent(in) :: iout + integer(I4B), target, intent(in) :: ncv + ! + allocate (this) + call this%init(name_model, 'SHF', 'SHF', inunit, iout, ncv) + this%text = text + end subroutine shf_cr + + !> @brief Announce package and set pointers to variables + !! + !! Announce package version and set array and variable pointers from the ABC + !! package for access by SHF. + !< + subroutine shf_ar_set_pointers(this) + ! -- dummy + class(ShfType) :: this + ! -- local + character(len=LENMEMPATH) :: abcMemoryPath + ! + ! -- print a message noting that the SHF utility is active + write (this%iout, '(a)') & + 'SHF -- SENSIBLE HEAT WILL BE INCLUDED IN THE ATMOSPHERIC BOUNDARY '// & + 'CONDITIONS FOR THE STREAMFLOW ENERGY TRANSPORT PACKAGE' + ! + ! -- set pointers to variables hosted in the ABC package + abcMemoryPath = create_mem_path(this%name_model, 'ABC') + call mem_setptr(this%rhoa, 'RHOA', abcMemoryPath) + call mem_setptr(this%cpa, 'CPA', abcMemoryPath) + call mem_setptr(this%cd, 'CD', abcMemoryPath) + call mem_setptr(this%wspd, 'WSPD', abcMemoryPath) + call mem_setptr(this%tatm, 'TATM', abcMemoryPath) + call mem_setptr(this%patm, 'PATM', abcMemoryPath) + call mem_setptr(this%ea, 'EA', abcMemoryPath) + call mem_setptr(this%ew, 'EW', abcMemoryPath) + ! + ! -- create time series manager + call tsmanager_cr(this%tsmanager, this%iout, & + removeTsLinksOnCompletion=.true., & + extendTsToEndOfSimulation=.true.) + end subroutine shf_ar_set_pointers + + !> @brief Calculate Sensible Heat Flux + !! + !! Calculate and return the sensible heat flux for one reach + !< + subroutine shf_cq(this, ifno, tstrm, tfac, toff, pfac, shflx, lhflx) + ! -- dummy + class(ShfType), intent(inout) :: this + integer(I4B), intent(in) :: ifno !< stream reach integer id + real(DP), intent(in) :: tstrm !< temperature of the stream reach + real(DP), intent(in) :: tfac !< temperature units adjustment factor + real(DP), intent(in) :: toff !< temperature units offset + real(DP), intent(in) :: pfac !< atmospheric pressure units adjustment factor + real(DP), intent(inout) :: shflx !< calculated sensible heat flux amount + real(DP), optional, intent(in) :: lhflx !< latent heat flux + ! -- local + real(DP) :: shf_const + real(DP) :: br + ! + ! -- calculate sensible heat flux using HGS equation + if (present(lhflx)) then + br = 0.00061_DP * this%patm(ifno) * pfac * & + (((tfac * tstrm + toff) - (tfac * this%tatm(ifno) + toff)) / (this%ew(ifno) - this%ea(ifno))) + shflx = br * lhflx + else + shf_const = this%cd * this%cpa * this%rhoa + shflx = shf_const * this%wspd(ifno) * ((tfac * this%tatm(ifno) + toff) - (tfac * tstrm + toff)) + end if + end subroutine shf_cq + + !> @brief Deallocate package memory + !! + !! Deallocate TVK package scalars and arrays. + !< + subroutine shf_da(this) + ! -- dummy + class(ShfType) :: this + ! + ! -- nullify pointers to other package variables + nullify (this%rhoa) + nullify (this%cpa) + nullify (this%cd) + nullify (this%wspd) + nullify (this%tatm) + nullify (this%patm) + nullify (this%ea) + nullify (this%ew) + ! + ! -- deallocate parent + call pbstbase_da(this) + end subroutine shf_da + +end module SensHeatModule diff --git a/src/Model/GroundWaterEnergy/gwe-swr.f90 b/src/Model/GroundWaterEnergy/gwe-swr.f90 new file mode 100644 index 00000000000..24c027bbe5b --- /dev/null +++ b/src/Model/GroundWaterEnergy/gwe-swr.f90 @@ -0,0 +1,121 @@ +!> @brief This module contains methods for calculating the shortwave radiation heat flux +!! +!! This module contains the methods used to calculate the shortwave radiation heat flux +!! for surface-water boundaries, like streams and lakes. In its current form, +!! this class acts like a package to a package, similar to the TVK package that +!! can be invoked from the NPF package. Once this package is completed in its +!! prototyped form, it will likely be moved around. +!< + +module ShortwaveModule + use ConstantsModule, only: LINELENGTH, LENMEMPATH, DZERO, LENVARNAME + use KindModule, only: I4B, DP + use MemoryManagerModule, only: mem_setptr + use MemoryHelperModule, only: create_mem_path + use TimeSeriesLinkModule, only: TimeSeriesLinkType + use TimeSeriesManagerModule, only: TimeSeriesManagerType, tsmanager_cr + use SimModule, only: store_error + use SimVariablesModule, only: errmsg + use PbstBaseModule, only: PbstBaseType, pbstbase_da + + implicit none + + public :: SwrType + public :: swr_cr + + character(len=16) :: text = ' SWR' + + type, extends(PbstBaseType) :: SwrType + + real(DP), dimension(:), pointer, contiguous :: solr => null() !< solar radiation + real(DP), dimension(:), pointer, contiguous :: shd => null() !< shade fraction + real(DP), dimension(:), pointer, contiguous :: swrefl => null() !< shortwave reflectance of water surface + + contains + + procedure :: da => swr_da + procedure :: pbst_ar => swr_ar_set_pointers + procedure, public :: swr_cq + + end type SwrType + +contains + + !> @brief Create a new SwrType object + !! + !! Create a new shortwave radiation flux (SwrType) object. Initially for use with + !! the SFE package. + !< + subroutine swr_cr(this, name_model, inunit, iout, ncv) + ! -- dummy + type(SwrType), pointer, intent(out) :: this + character(len=*), intent(in) :: name_model + integer(I4B), intent(in) :: inunit + integer(I4B), intent(in) :: iout + integer(I4B), target, intent(in) :: ncv + ! + allocate (this) + call this%init(name_model, 'SWR', 'SWR', inunit, iout, ncv) + this%text = text + end subroutine swr_cr + + !> @brief Announce package and set pointers to variables + !! + !! Announce package version and set array and variable pointers from the ABC + !! package for access by SWR. + !< + subroutine swr_ar_set_pointers(this) + ! -- dummy + class(SwrType) :: this + ! -- local + character(len=LENMEMPATH) :: abcMemoryPath + ! + ! -- print a message noting that the SWR utility is active + write (this%iout, '(a)') & + 'SWR -- SHORTWAVE RADIATION WILL BE INCLUDED IN THE ATMOSPHERIC '// & + 'BOUNDARY CONDITIONS FOR THE STREAMFLOW ENERGY TRANSPORT PACKAGE' + ! + ! -- Set pointers to variables hosted in the ABC package + abcMemoryPath = create_mem_path(this%name_model, 'ABC') + call mem_setptr(this%solr, 'SOLR', abcMemoryPath) + call mem_setptr(this%shd, 'SHD', abcMemoryPath) + call mem_setptr(this%swrefl, 'SWREFL', abcMemoryPath) + ! + ! -- create time series manager + call tsmanager_cr(this%tsmanager, this%iout, & + removeTsLinksOnCompletion=.true., & + extendTsToEndOfSimulation=.true.) + end subroutine swr_ar_set_pointers + + !> @brief Calculate Shortwave Radiation Heat Flux + !! + !! Calculate and return the shortwave radiation heat flux for one reach + !< + subroutine swr_cq(this, ifno, swrflx) + ! -- dummy + class(SwrType), intent(inout) :: this + integer(I4B), intent(in) :: ifno !< stream reach integer id + real(DP), intent(inout) :: swrflx !< calculated shortwave radiation heat flux amount + ! + ! -- calculate shortwave radiation heat flux (version: user input of sol rad data) + swrflx = (1 - this%shd(ifno)) * (1 - this%swrefl(ifno)) * this%solr(ifno) + end subroutine swr_cq + + !> @brief Deallocate package memory + !! + !! Deallocate TVK package scalars and arrays. + !< + subroutine swr_da(this) + ! -- dummy + class(SwrType) :: this + ! + ! -- deallocate time series + nullify (this%shd) + nullify (this%swrefl) + nullify (this%solr) + ! + ! -- deallocate parent + call pbstbase_da(this) + end subroutine swr_da + +end module ShortwaveModule diff --git a/src/Model/GroundWaterFlow/gwf-sfr.f90 b/src/Model/GroundWaterFlow/gwf-sfr.f90 index 320d8bb9d5c..9ae58388fcb 100644 --- a/src/Model/GroundWaterFlow/gwf-sfr.f90 +++ b/src/Model/GroundWaterFlow/gwf-sfr.f90 @@ -2222,7 +2222,7 @@ subroutine sfr_fc(this, rhs, ia, idxglo, matrix_sln) ! -- set initial stage to calculate stage change s0 = this%stage(n) ! - ! -- solve for flow in swr + ! -- solve for flow in sfr if (this%iboundpak(n) /= 0) then call this%sfr_solve(n, hgwf, hhcof, rrhs) else @@ -5121,14 +5121,15 @@ subroutine sfr_setup_budobj(this) text = ' EVAPORATION' idx = idx + 1 maxlist = this%maxbound - naux = 0 + naux = 1 + auxtxt(1) = ' SURFACE-AREA' call this%budobj%budterm(idx)%initialize(text, & this%name_model, & this%packName, & this%name_model, & this%packName, & - maxlist, .false., .false., & - naux) + maxlist, .false., .true., & + naux, auxtxt) ! ! -- text = ' RUNOFF' @@ -5260,7 +5261,7 @@ subroutine sfr_fill_budobj(this) real(DP) :: d real(DP) :: ca real(DP) :: a - real(DP) :: wp + real(DP) :: wp !< wetted perimeter real(DP) :: l ! ! -- initialize counter @@ -5345,11 +5346,18 @@ subroutine sfr_fill_budobj(this) call this%budobj%budterm(idx)%reset(this%maxbound) do n = 1, this%maxbound if (this%iboundpak(n) /= 0) then + if (this%depth(n) > DZERO) then + a = this%calc_surface_area_wet(n, this%depth(n)) + else + a = DZERO + end if + this%qauxcbc(1) = a q = -this%simevap(n) else + this%qauxcbc(1) = DZERO q = DZERO end if - call this%budobj%budterm(idx)%update_term(n, n, q) + call this%budobj%budterm(idx)%update_term(n, n, q, this%qauxcbc) end do ! ! -- RUNOFF diff --git a/src/Model/TransportModel/tsp-apt.f90 b/src/Model/TransportModel/tsp-apt.f90 index 969f1ed6cb9..a14971961ac 100644 --- a/src/Model/TransportModel/tsp-apt.f90 +++ b/src/Model/TransportModel/tsp-apt.f90 @@ -154,6 +154,7 @@ module TspAptModule procedure :: apt_solve procedure :: pak_solve procedure :: bnd_options => apt_options + procedure :: gc_options procedure :: read_dimensions => apt_read_dimensions procedure :: apt_read_cvs procedure :: read_initial_attr => apt_read_initial_attr @@ -174,6 +175,7 @@ module TspAptModule procedure :: pak_setup_budobj procedure :: apt_fill_budobj procedure :: pak_fill_budobj + procedure :: ancil_rp procedure, public :: apt_get_volumes !< made public for sft/sfe procedure, public :: apt_stor_term procedure, public :: apt_tmvr_term @@ -312,7 +314,7 @@ subroutine apt_ar(this) ! -- Get obs setup call this%obs%obs_ar() ! - ! --print a message identifying the apt package. + ! -- print a message identifying the apt package. write (this%iout, fmtapt) this%inunit ! ! -- Allocate arrays @@ -358,8 +360,22 @@ subroutine apt_ar(this) end if end if end if + ! + ! -- The following is not callable from apt_ar, only reachable from sfe_ar + !call this%pbst_shf_ar() end subroutine apt_ar + !> @brief Allocate and read appropriate pbst sub-package + !< + !subroutine pbst_shf_ar(this) + ! ! -- dummy + ! class(TspAptType), intent(inout) :: this + ! ! + ! ! -- This routine to be overridden + ! call store_error('Program error: pbst_shf_ar should be overwritten.', & + ! terminate=.TRUE.) + !end subroutine pbst_shf_ar + !> @brief Advanced package transport read and prepare (rp) routine !! !! This subroutine calls the attached packages' read and prepare routines. @@ -455,7 +471,7 @@ subroutine apt_rp(this) call this%inputtab%line_to_columns(line) end if end do stressperiod - + ! if (this%iprpak /= 0) then call this%inputtab%finalize_table() end if @@ -465,6 +481,10 @@ subroutine apt_rp(this) write (this%iout, fmtlsp) trim(this%filtyp) end if ! + ! -- call abc_rp() routine separately since its stress period info is + ! not strictly tied to whether apt_rp is updated (or not updated) + call this%ancil_rp() + ! ! -- write summary of stress period error messages ierr = count_errors() if (ierr > 0) then @@ -478,10 +498,25 @@ subroutine apt_rp(this) end do end subroutine apt_rp + !> @brief Support for ancillary transport package read and prepare (rp) routine + !! + !! A temporary work-around for calling _rp() routines in the PaBST-family of + !! classes. While specific to GWE at the moment, which cuts against the + !! grain of what TspApt should be all about (i.e., generalized transport code + !! that is not specific to either GWT or GWE), it could be leveraged by GWT + !! as well if any sort of 'package for a package' process was implemented + !! there. At some point, more likely, this is going to have to be lifted out + !!and a more MF6-way of doing things implemented. + subroutine ancil_rp(this) + ! -- dummy + class(TspAptType), intent(inout) :: this + ! -- subroutine available for override by apt-child classes + end subroutine ancil_rp + subroutine apt_ad_chk(this) ! -- dummy class(TspAptType), intent(inout) :: this - ! function available for override by packages + ! -- subroutine available for override by packages end subroutine apt_ad_chk !> @brief Advanced package transport set stress period routine. @@ -592,9 +627,6 @@ subroutine pak_set_stressperiod(this, itemno, keyword, found) integer(I4B), intent(in) :: itemno character(len=*), intent(in) :: keyword logical, intent(inout) :: found - ! -- local - - ! -- formats ! ! -- this routine should never be called found = .false. @@ -616,7 +648,7 @@ function apt_check_valid(this, itemno) result(ierr) ierr = 0 if (itemno < 1 .or. itemno > this%ncv) then write (errmsg, '(a,1x,i6,1x,a,1x,i6)') & - 'Featureno ', itemno, 'must be > 0 and <= ', this%ncv + 'Feature No. ', itemno, 'must be > 0 and <= ', this%ncv call store_error(errmsg) ierr = 1 end if @@ -975,7 +1007,7 @@ subroutine apt_ot_package_flows(this, icbcfl, ibudfl) pertim, totim, this%iout) end if ! - ! -- Print lake flows table + ! -- Print feature flows table if (ibudfl /= 0 .and. this%iprflow /= 0) then call this%budobj%write_flowtable(this%dis, kstp, kper) end if @@ -1175,7 +1207,7 @@ subroutine apt_allocate_arrays(this) ! -- local integer(I4B) :: n ! - ! -- call standard BndType allocate scalars + ! -- call standard BndType allocate arrays call this%BndType%allocate_arrays() ! ! -- Allocate @@ -1320,12 +1352,14 @@ subroutine apt_options(this, option, found) logical, intent(inout) :: found ! -- local character(len=MAXCHARLEN) :: fname, keyword + logical(LGP) :: foundgcclassoption ! -- formats character(len=*), parameter :: fmtaptbin = & "(4x, a, 1x, a, 1x, ' WILL BE SAVED TO FILE: ', a, & &/4x, 'OPENED ON UNIT: ', I0)" ! found = .true. + foundgcclassoption = .false. select case (option) case ('FLOW_PACKAGE_NAME') call this%parser%GetStringCaps(this%flowpackagename) @@ -1394,11 +1428,31 @@ subroutine apt_options(this, option, found) end if case default ! - ! -- No options found - found = .false. + ! -- check for grandchild class options + call this%gc_options(option, foundgcclassoption) + ! + ! -- No grandchild options found + found = foundgcclassoption end select end subroutine apt_options + !> @ brief Read additional grandchild options for package + !! + !! Check whether daughter packages are specified in the options block. + !! This method should be overridden by the (grand)child package to read the + !! options that are in addition to the base options implemented in the + !! boundary package class. + !< + subroutine gc_options(this, option, found) + ! -- dummy + class(TspAptType), intent(inout) :: this !< TspAptType object + character(len=*), intent(inout) :: option !< option keyword string + logical(LGP), intent(inout) :: found !< boolean indicating if the option was found + ! + ! Return with found = .false. + found = .false. + end subroutine gc_options + !> @brief Determine dimensions for this advanced package !< subroutine apt_read_dimensions(this) @@ -1497,7 +1551,7 @@ subroutine apt_read_cvs(this) call mem_allocate(this%lauxvar, this%naux, this%ncv, 'LAUXVAR', & this%memoryPath) ! - ! -- lake boundary and concentrations + ! -- feature boundary and concentrations if (this%imatrows == 0) then call mem_allocate(this%iboundpak, this%ncv, 'IBOUND', this%memoryPath) call mem_allocate(this%xnewpak, this%ncv, 'XNEWPAK', this%memoryPath) @@ -1590,7 +1644,7 @@ subroutine apt_read_cvs(this) nlak = nlak + 1 end do ! - ! -- check for duplicate or missing lakes + ! -- check for duplicate or missing feature ids do n = 1, this%ncv if (nboundchk(n) == 0) then write (errmsg, '(a,1x,i0)') 'No data specified for feature', n @@ -1644,7 +1698,7 @@ subroutine apt_read_initial_attr(this) ! -- todo: read boundname end do ! - ! -- initialize status (iboundpak) of lakes to active + ! -- initialize status (iboundpak) of features do n = 1, this%ncv if (this%status(n) == 'CONSTANT') then this%iboundpak(n) = -1 @@ -1788,14 +1842,14 @@ subroutine apt_accumulate_ccterm(this, ilak, rrate, ccratin, ccratout) q = -rrate this%ccterm(ilak) = this%ccterm(ilak) + q ! - ! -- See if flow is into lake or out of lake. + ! -- See if flow is into feature or out of feature. if (q < DZERO) then ! - ! -- Flow is out of lake subtract rate from ratout. + ! -- Flow is out of feature subtract rate from ratout. ccratout = ccratout - q else ! - ! -- Flow is into lake; add rate to ratin. + ! -- Flow is into feature; add rate to ratin. ccratin = ccratin + q end if end if @@ -2390,7 +2444,6 @@ end subroutine apt_copy2flowp !! - overrides BndType%bnd_obs_supported() !< logical function apt_obs_supported(this) - ! -- modules ! -- dummy class(TspAptType) :: this ! @@ -2405,10 +2458,8 @@ end function apt_obs_supported !! - overrides BndType%bnd_df_obs !< subroutine apt_df_obs(this) - ! -- modules ! -- dummy class(TspAptType) :: this - ! -- local ! ! -- call additional specific observations for lkt, sft, mwt, and uzt call this%pak_df_obs() @@ -2420,10 +2471,8 @@ end subroutine apt_df_obs !! - stores observations supported by the APT package !! - must be overridden by child class subroutine pak_df_obs(this) - ! -- modules ! -- dummy class(TspAptType) :: this - ! -- local ! ! -- this routine should never be called call store_error('Program error: pak_df_obs not implemented.', & @@ -2825,12 +2874,12 @@ end subroutine pak_bd_obs !! and not ID2. !< subroutine apt_process_obsID(obsrv, dis, inunitobs, iout) - ! -- dummy variables + ! -- dummy type(ObserveType), intent(inout) :: obsrv !< Observation object class(DisBaseType), intent(in) :: dis !< Discretization object integer(I4B), intent(in) :: inunitobs !< file unit number for the package observation file integer(I4B), intent(in) :: iout !< model listing file unit number - ! -- local variables + ! -- local integer(I4B) :: nn1 integer(I4B) :: icol integer(I4B) :: istart diff --git a/src/Utilities/Constants.f90 b/src/Utilities/Constants.f90 index 7b2423e6d96..6a6a9b8e9b1 100644 --- a/src/Utilities/Constants.f90 +++ b/src/Utilities/Constants.f90 @@ -65,6 +65,7 @@ module ConstantsModule real(DP), parameter :: DZERO = 0.0_DP !< real constant zero real(DP), parameter :: DQUARTER = 0.25_DP !< real constant 1/3 real(DP), parameter :: DONETHIRD = 1.0_DP / 3.0_DP !< real constant 1/3 + real(DP), parameter :: DONESEVENTH = 1.0_DP / 7.0_DP !< real constant 1/7 real(DP), parameter :: DHALF = 0.5_DP !< real constant 1/2 real(DP), parameter :: DP6 = 0.6_DP !< real constant 3/5 real(DP), parameter :: DTWOTHIRDS = 2.0_DP / 3.0_DP !< real constant 2/3 @@ -132,6 +133,9 @@ module ConstantsModule real(DP), parameter :: DGRAVITY = 9.80665_DP !< real constant gravitational acceleration (m/(s s)) real(DP), parameter :: DCD = 0.61_DP !< real constant weir coefficient in SI units + real(DP), parameter :: DCTOK = 273.16_DP !< real constant conversion of degrees celsius to Kelvin + real(DP), parameter :: DSTEFANBOLTZMANN = 5.670374419e-8_DP !< real constant Stefan-Boltzmann constant + character(len=10), dimension(3, 3), parameter :: & cidxnames = reshape( & [' NODE', ' ', ' ', & diff --git a/src/meson.build b/src/meson.build index fe8cfa02e10..880a1f22009 100644 --- a/src/meson.build +++ b/src/meson.build @@ -197,14 +197,20 @@ modflow_sources = files( 'Model' / 'Geometry' / 'CircularGeometry.f90', 'Model' / 'Geometry' / 'RectangularGeometry.f90', 'Model' / 'GroundWaterEnergy' / 'gwe.f90', + 'Model' / 'GroundWaterEnergy' / 'gwe-abc.f90', 'Model' / 'GroundWaterEnergy' / 'gwe-cnd.f90', 'Model' / 'GroundWaterEnergy' / 'gwe-ctp.f90', 'Model' / 'GroundWaterEnergy' / 'gwe-esl.f90', 'Model' / 'GroundWaterEnergy' / 'gwe-est.f90', + 'Model' / 'GroundWaterEnergy' / 'gwe-lhf.f90', 'Model' / 'GroundWaterEnergy' / 'gwe-lke.f90', + 'Model' / 'GroundWaterEnergy' / 'gwe-lwr.f90', 'Model' / 'GroundWaterEnergy' / 'gwe-mwe.f90', 'Model' / 'GroundWaterEnergy' / 'gwe-sfe.f90', + 'Model' / 'GroundWaterEnergy' / 'gwe-shf.f90', + 'Model' / 'GroundWaterEnergy' / 'gwe-swr.f90', 'Model' / 'GroundWaterEnergy' / 'gwe-uze.f90', + 'Model' / 'GroundWaterEnergy' / 'PbstBase.f90', 'Model' / 'GroundWaterFlow' / 'gwf.f90', 'Model' / 'GroundWaterFlow' / 'gwf-api.f90', 'Model' / 'GroundWaterFlow' / 'gwf-buy.f90',