This page was generated from
ex-gwe-prt.py.
It's also available as a notebook.
Thermal Profile along a Particle Path
Using a Voronoi grid, this example demonstrates the combined use of GWF, GWE, and PRT. A steady flow field is established, followed by the establishment of a steady-state temperature field. Two particles are then routed through the model domain and the temperature of the particle is determined along the flow path.
Initial setup
Import dependencies, define the example name and workspace, and read settings from environment variables.
[1]:
import pathlib as pl
from pprint import pformat
import flopy
import git
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from flopy.discretization import VertexGrid
from flopy.mf6 import MFSimulation
from flopy.utils import GridIntersect
from flopy.utils.triangle import Triangle
from flopy.utils.voronoi import VoronoiGrid
from matplotlib.collections import LineCollection
from modflow_devtools.misc import get_env, timed
from scipy.interpolate import CloughTocher2DInterpolator
from shapely.geometry import LineString, Point
[2]:
# Example name and workspace paths. If this example is running
# in the git repository, use the folder structure described in
# the README. Otherwise just use the current working directory.
sim_name = "ex-gwe-prt"
gwf_name = sim_name + "-gwf"
gwe_name = sim_name + "-gwe"
prt_name = sim_name + "-prt"
try:
root = pl.Path(git.Repo(".", search_parent_directories=True).working_dir)
except:
root = None
data_path = root / "data" / sim_name if root else pl.Path.cwd()
workspace = root / "examples" if root else pl.Path.cwd()
figs_path = root / "figures" if root else pl.Path.cwd()
sim_ws = workspace / sim_name
gwf_ws = sim_ws / "gwf"
gwe_ws = sim_ws / "gwe"
prt_ws = sim_ws / "prt"
gwf_ws.mkdir(exist_ok=True, parents=True)
gwe_ws.mkdir(exist_ok=True, parents=True)
prt_ws.mkdir(exist_ok=True, parents=True)
# Settings from environment variables
write = get_env("WRITE", True)
run = get_env("RUN", True)
plot = get_env("PLOT", True)
plot_show = get_env("PLOT_SHOW", True)
plot_save = get_env("PLOT_SAVE", True)
Define parameters
Define model units, parameters and other settings.
[3]:
# Model units
length_units = "meters"
time_units = "days"
[4]:
# Model parameters
xmin = 0.0 # Left model domain extent ($m$)
xmax = 2000.0 # Right model domain extent ($m$)
ymin = 0.0 # South model domain extent ($m$)
ymax = 1000.0 # North model domain extent ($m$)
top = 1.0 # Top of the model grid ($m$)
botm = 0.0 # Bottom of the model grid ($m$)
angle_min = 30 # Minimum angle of interior angle of a Voronoi grid cell ($^{\circ}$)
area_max = 1000.0 # Maximum area of a Voronoi grid cell ($m^2$)
nlay = 1 # Number of layers ($-$)
porosity = 0.1 # Porosity ($-$)
strt_temp = 10.0 # Starting temperature of model domain ($^{\circ}C$)
scheme = "TVD" # Advection scheme ($-$)
alh = 0.0 # Longitudinal mechanical dispersion term ($m$)
ath1 = 0.0 # Transverse mechanical dispersion term ($m$)
ktw = 0.56 # Thermal conductivity of water ($\frac{W}{m \cdot ^{\circ} C}$)
kts = 2.5 # Thermal conductivity of aquifer material ($\frac{W}{m \cdot ^{\circ} C}$)
rhow = 1000 # Density of water ($kg/m^3$)
cpw = 4180.0 # Heat capacity of water ($\frac{J}{kg \cdot ^{\circ} C}$)
rhos = 2650.0 # Density of dry solid aquifer material ($kg/m^3$)
cps = 900.0 # Heat capacity of dry solid aquifer material ($\frac{J}{kg \cdot ^{\circ} C}$)
# Values that do not show up in the table generated for latex
delr = area_max**0.5
ncol = xmax / delr
nrow = ymax / delr
botm = [botm]
nodes = ncol * nrow
lhv = 2500.0
rpts = [[20, i, 0.5] for i in range(1, 999, 20)]
# Convert kts and ktw units to use days
ktw = ktw * 86400
kts = kts * 86400
[5]:
# ### Solver Parameters
nouter = 1000
ninner = 200
hclose = 1e-6
rclose = 1e-6
relax = 1.0
[6]:
def get_grid(workspace):
workspace.mkdir(exist_ok=True, parents=True)
tri = Triangle(
maximum_area=area_max,
angle=angle_min,
model_ws=workspace,
exe_name="triangle",
)
poly = np.array(((xmin, ymin), (xmax, ymin), (xmax, ymax), (xmin, ymax)))
tri.add_polygon(poly)
tri.build(verbose=False)
return VoronoiGrid(tri)
[7]:
def buildout_grid(ws):
grid = get_grid(ws / "grid")
vgrid = VertexGrid(**grid.get_gridprops_vertexgrid(), nlay=1)
gi = GridIntersect(vgrid)
return grid, vgrid, gi
[8]:
def build_gwf_sim(name):
global cell2d
# create grid
grid, vgrid, gi = buildout_grid(gwf_ws)
# identify cells on left edge
line = LineString([(xmin, ymin), (xmin, ymax)])
cells_left = gi.intersect(line)["cellids"]
cells_left = np.array(list(cells_left))
# identify cells on right edge
line = LineString([(xmax, ymin), (xmax, ymax)])
cells_right = gi.intersect(line)["cellids"]
cells_right = np.array(list(cells_right))
# identify cells on bottom edge
line = LineString([(xmin, ymin), (xmax, ymin)])
cells_bottom = gi.intersect(line)["cellids"]
cells_bottom = np.array(list(cells_bottom))
# identify well cell
points = [Point((1200, 500)), Point((700, 200)), Point((1600, 700))]
well_cells = [vgrid.intersect(p.x, p.y) for p in points]
# create simulation
sim = flopy.mf6.MFSimulation(
sim_name=name, version="mf6", exe_name="mf6", sim_ws=gwf_ws
)
flopy.mf6.ModflowTdis(sim, time_units=time_units, perioddata=[[1.0, 1, 1.0]])
gwf = flopy.mf6.ModflowGwf(sim, modelname=gwf_name, save_flows=True)
flopy.mf6.ModflowIms(
sim,
print_option="SUMMARY",
complexity="complex",
outer_dvclose=1.0e-8,
inner_dvclose=1.0e-8,
filename=f"{gwf_name}.ims",
)
flopy.mf6.ModflowGwfdisv(
gwf, nlay=nlay, **grid.get_disv_gridprops(), top=top, botm=botm
)
# k, j, q
wells = [ # if 5.0 is changed to 0.05 TDIS must be 1e6 250 1.01
[0, c, -0.05, 80.0] for c in well_cells
]
wells[0][0] *= 100
wells[2][-2] *= -100
flopy.mf6.ModflowGwfwel(
gwf,
auxiliary="TEMPERATURE",
maxbound=len(wells),
save_flows=True,
pname="WEL",
stress_period_data={0: wells},
filename=f"{gwf_name}.wel",
)
flopy.mf6.ModflowGwfnpf(
gwf,
xt3doptions=[(True)],
k=10.0,
save_saturation=True,
save_specific_discharge=True,
filename=f"{gwf_name}.npf",
)
flopy.mf6.ModflowGwfsto(
gwf,
ss=0,
sy=0,
steady_state={0: True},
pname="STO",
filename="{}.sto".format(gwf_name),
)
flopy.mf6.ModflowGwfic(gwf, strt=1.0, pname="IC", filename="{}.ic".format(gwf_name))
flopy.mf6.ModflowGwfoc(
gwf,
budget_filerecord=f"{gwf_name}.cbc",
head_filerecord=f"{gwf_name}.hds",
saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
printrecord=[("HEAD", "LAST"), ("BUDGET", "LAST")],
filename="{}.oc".format(gwf_name),
)
chdlist = []
icpl_seen = []
cell2d = grid.get_disv_gridprops()["cell2d"]
for scl, icpl in enumerate(cells_left):
y_coord = cell2d[icpl][2]
chdlist.append([(0, icpl), 2.0, 100.0 * y_coord / ymax])
icpl_seen.append(icpl)
for icpl in cells_right:
chdlist.append([(0, icpl), 1.0, 0.0])
icpl_seen.append(icpl)
for icpl in cells_bottom:
if icpl in icpl_seen:
continue
chdlist.append([(0, icpl), 1.8, 80.0])
flopy.mf6.ModflowGwfchd(
gwf,
auxiliary="TEMPERATURE",
stress_period_data=chdlist,
pname="CHD",
filename="{}.chd".format(gwf_name),
)
return sim
[9]:
def build_gwe_sim(name):
grid, dummy, gi = buildout_grid(gwe_ws)
sim = flopy.mf6.MFSimulation(sim_name=name, sim_ws=gwe_ws, exe_name="mf6")
flopy.mf6.ModflowTdis(sim, time_units=time_units, perioddata=[[1.0e6, 1000, 1.003]])
gwe = flopy.mf6.MFModel(
sim,
model_type="gwe6",
modelname=gwe_name,
model_nam_file="{}.nam".format(gwe_name),
)
imsgwe = flopy.mf6.ModflowIms(
sim,
print_option="SUMMARY",
outer_dvclose=hclose,
outer_maximum=nouter,
under_relaxation="NONE",
inner_maximum=ninner,
inner_dvclose=hclose,
rcloserecord=rclose,
linear_acceleration="BICGSTAB",
scaling_method="NONE",
reordering_method="NONE",
relaxation_factor=relax,
filename=f"{gwe_name}.ims",
)
sim.register_ims_package(imsgwe, [gwe.name])
flopy.mf6.ModflowGwedisv(
gwe,
nlay=nlay,
**grid.get_disv_gridprops(),
top=top,
botm=botm,
pname="DISV",
filename="{}.disv".format(gwe_name),
)
flopy.mf6.ModflowGweic(
gwe, strt=strt_temp, pname="IC", filename="{}.ic".format(gwe_name)
)
flopy.mf6.ModflowGweadv(
gwe, scheme=scheme, pname="ADV", filename="{}.adv".format(gwe_name)
)
if ktw != 0:
flopy.mf6.ModflowGwecnd(
gwe,
alh=alh,
ath1=ath1,
ktw=ktw,
kts=kts,
pname="CND-e",
filename="{}.cnd".format(gwe_name),
)
flopy.mf6.ModflowGweest(
gwe,
porosity=porosity,
cps=cps,
rhos=rhos,
packagedata=[cpw, rhow, lhv],
pname="EST-e",
filename="{}.est".format(gwe_name),
)
sourcerecarray = [
("WEL", "AUX", "TEMPERATURE"),
("CHD", "AUX", "TEMPERATURE"),
]
flopy.mf6.ModflowGwessm(
gwe, sources=sourcerecarray, pname="SSM-e", filename="{}.ssm".format(gwe_name)
)
flopy.mf6.ModflowGweoc(
gwe,
budget_filerecord="{}.cbc".format(gwe_name),
temperature_filerecord="{}.ucn".format(gwe_name),
temperatureprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")],
saverecord={
0: [("TEMPERATURE", "ALL"), ("BUDGET", "ALL")],
},
printrecord=[("TEMPERATURE", "LAST"), ("BUDGET", "LAST")],
)
pd = [
(
"GWFHEAD",
pl.Path(f"../{gwf_ws.name}/{gwf_name}.hds"),
None,
),
(
"GWFBUDGET",
pl.Path(f"../{gwf_ws.name}/{gwf_name}.cbc"),
None,
),
]
flopy.mf6.ModflowGwefmi(gwe, packagedata=pd)
return sim
[10]:
def build_prt_sim(name):
# create grid
grid, vgrid, gi = buildout_grid(prt_ws)
# create simulation
sim = flopy.mf6.MFSimulation(
sim_name=name, version="mf6", exe_name="mf6", sim_ws=prt_ws
)
flopy.mf6.ModflowTdis(sim, time_units=time_units, perioddata=[[1.0, 1, 1.0]])
prt = flopy.mf6.ModflowPrt(sim, modelname=prt_name)
flopy.mf6.ModflowGwfdisv(
prt, nlay=nlay, **grid.get_disv_gridprops(), top=top, botm=botm
)
flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=porosity)
prpdata = [
# index, (layer, cell index), x, y, z
(i, (0, vgrid.intersect(p[0], p[1])), p[0], p[1], p[2])
for i, p in enumerate([rpts[23], rpts[32]]) # first release point crashes
]
prp_track_file = f"{prt_name}.prp.trk"
prp_track_csv_file = f"{prt_name}.prp.trk.csv"
flopy.mf6.ModflowPrtprp(
prt,
pname="prp1",
filename=f"{prt_name}_1.prp",
nreleasepts=len(prpdata),
packagedata=prpdata,
perioddata={0: ["FIRST"]},
track_filerecord=[prp_track_file],
trackcsv_filerecord=[prp_track_csv_file],
boundnames=True,
stop_at_weak_sink=True, # currently required for this problem
exit_solve_tolerance=1e-10,
)
prt_track_file = f"{prt_name}.trk"
prt_track_csv_file = f"{prt_name}.trk.csv"
flopy.mf6.ModflowPrtoc(
prt,
pname="oc",
track_filerecord=[prt_track_file],
trackcsv_filerecord=[prt_track_csv_file],
)
pd = [
(
"GWFHEAD",
pl.Path(f"../{gwf_ws.name}/{gwf_name}.hds"),
None,
),
(
"GWFBUDGET",
pl.Path(f"../{gwf_ws.name}/{gwf_name}.cbc"),
None,
),
]
flopy.mf6.ModflowPrtfmi(
prt,
packagedata=pd,
)
ems = flopy.mf6.ModflowEms(
sim,
pname="ems",
filename=f"{prt_name}.ems",
)
sim.register_solution_package(ems, [prt.name])
return sim
[11]:
def build_models():
gwf_sim = build_gwf_sim(sim_name)
gwe_sim = build_gwe_sim(sim_name)
prt_sim = build_prt_sim(sim_name)
return gwf_sim, gwe_sim, prt_sim
[12]:
def write_models(*sims, silent=False):
for sim in sims:
if isinstance(sim, MFSimulation):
sim.write_simulation(silent=silent)
else:
sim.write_input()
[13]:
@timed
def run_models(*sims, silent=False):
for sim in sims:
if isinstance(sim, MFSimulation):
success, buff = sim.run_simulation(silent=silent, report=True)
else:
success, buff = sim.run_model(silent=silent, report=True)
assert success, pformat(buff)
[14]:
def temp_interp(pls, temperatures):
X = []
Y = []
Z = []
# Match up X, Y locations with the state-steady temperature field
for grdcell in cell2d:
cellid = grdcell[0]
X.append(grdcell[1])
Y.append(grdcell[2])
Z.append(temperatures[cellid])
interp = CloughTocher2DInterpolator(list(zip(X, Y)), Z)
partIDs = pls.name.unique()
part1 = pls[pls["name"] == partIDs[0]].copy()
part2 = pls[pls["name"] == partIDs[1]].copy()
for index, row in part1.iterrows():
x_prt = row["x"]
y_prt = row["y"]
therm = interp(x_prt, y_prt)
part1.loc[index, "therm"] = therm
for index, row in part2.iterrows():
x_prt = row["x"]
y_prt = row["y"]
therm = interp(x_prt, y_prt)
part2.loc[index, "therm"] = therm
thermDF = pd.concat([part1, part2], ignore_index=True)
return thermDF
[15]:
def plot_results(gwf_sim, gwe_sim, prt_sim, silent=True):
# get gwf output
gwf = gwf_sim.get_model()
gwe = gwe_sim.get_model()
temperatures = gwe.output.temperature().get_data()
bdobj = gwf.output.budget()
spdis = bdobj.get_data(text="DATA-SPDIS")[0]
qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf)
# get prt output
prt_track_csv_file = f"{prt_name}.trk.csv"
pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False)
pls = temp_interp(pls, temperatures[0, 0])
plot_2d = True
if plot_2d:
fig, axes = plt.subplots(
3,
1,
figsize=(6, 7),
tight_layout=True,
gridspec_kw={"height_ratios": [3, 1, 1]},
)
ax = axes[0]
ax.set_xlim([0, 2000])
ax.set_ylim([0, 1000])
pmv = flopy.plot.PlotMapView(model=gwf, ax=ax)
ax.set_aspect("equal")
pmv.plot_grid(alpha=0.25)
pmv.plot_ibound(alpha=0.5)
tempmesh = pmv.plot_array(temperatures, alpha=0.7)
cv = pmv.contour_array(temperatures, levels=np.linspace(0, 80, 9))
plt.clabel(cv, colors="k")
plt.colorbar(
tempmesh,
shrink=0.5,
ax=ax,
label="Temperature",
location="bottom",
fraction=0.1,
)
handles = [
mpl.lines.Line2D(
[0],
[0],
marker=">",
linestyle="",
label="Specific discharge",
color="grey",
markerfacecolor="gray",
),
]
handles.append(
mpl.lines.Line2D(
[0],
[0],
marker="o",
linestyle="",
label="Well",
markerfacecolor="red",
),
)
ax.legend(
handles=handles,
loc="lower right",
)
pmv.plot_vector(qx, qy, normalize=True, alpha=0.25)
pmv.plot_bc(ftype="WEL")
mf6_plines = pls.groupby(["iprp", "irpt", "trelease"])
for ipl, ((iprp, irpt, trelease), pl) in enumerate(mf6_plines):
pl.plot(
# title=title,
kind="line",
linestyle="--",
marker="o",
markersize=0,
x="x",
y="y",
xlabel="X",
ylabel="Y",
ax=ax,
legend=False,
color="blue",
)
if ipl == 0:
ax.annotate(
"Particle 1",
xy=(1050, 380),
xycoords="data",
xytext=(30, -20),
textcoords="offset points",
bbox=dict(boxstyle="round", fc="1.0", alpha=0.66),
arrowprops=dict(
arrowstyle="->",
shrinkA=0,
shrinkB=5,
connectionstyle="angle,angleA=0,angleB=135,rad=40",
),
)
else:
ax.annotate(
"Particle 2",
xy=(1050, 610),
xycoords="data",
xytext=(-75, 10),
textcoords="offset points",
bbox=dict(boxstyle="round", fc="1.0", alpha=0.66),
arrowprops=dict(
arrowstyle="->",
shrinkA=0,
shrinkB=5,
connectionstyle="angle,angleA=0,angleB=135,rad=30",
),
)
# Setup 2 additional plots relating the particle's x-position with temperature
# and alternatively the total travel time with temperature.
part1 = pls[pls["name"] == pls.name.unique()[0]][["x", "therm"]]
part2 = pls[pls["name"] == pls.name.unique()[1]][["x", "therm"]]
part1t = pls[pls["name"] == pls.name.unique()[0]][["t", "therm"]]
part2t = pls[pls["name"] == pls.name.unique()[1]][["t", "therm"]]
# Create a set of line segments so that they may be colored individually
# This creates the points as a N x 1 x 2 array so that we can stack points
# together easily to get the segments. The segments array for line collection
# needs to be (numlines) x (points per line) x 2 (for x and y)
points1 = np.array([part1["x"], part1["therm"]]).T.reshape(-1, 1, 2)
points1t = np.array([part1t["t"], part1["therm"]]).T.reshape(-1, 1, 2)
points2 = np.array([part2["x"], part2["therm"]]).T.reshape(-1, 1, 2)
points2t = np.array([part2t["t"], part2["therm"]]).T.reshape(-1, 1, 2)
segments1 = np.concatenate([points1[:-1], points1[1:]], axis=1)
segments1t = np.concatenate([points1t[:-1], points1t[1:]], axis=1)
segments2 = np.concatenate([points2[:-1], points2[1:]], axis=1)
segments2t = np.concatenate([points2t[:-1], points2t[1:]], axis=1)
# Create a continuous norm to map from data points to colors
pad = 5
norm = plt.Normalize(
np.min([part1["therm"].min(), part2["therm"].min()]) - pad,
np.max([part1["therm"].max(), part2["therm"].max()]) + pad,
)
lc1 = LineCollection(segments1, norm=norm)
lc1t = LineCollection(segments1t, norm=norm)
lc2 = LineCollection(segments2, norm=norm)
lc2t = LineCollection(segments2t, norm=norm)
# Set the values used for colormapping
lc1.set_array(part1["therm"])
lc1t.set_array(part1t["therm"])
lc2.set_array(part2["therm"])
lc2t.set_array(part2t["therm"])
lc1.set_linewidth(3)
lc1t.set_linewidth(3)
lc2.set_linewidth(3)
lc2t.set_linewidth(3)
line1 = axes[1].add_collection(lc1)
line2 = axes[1].add_collection(lc2)
axes[1].annotate("Particle 2", xy=(400, 68), xycoords="data")
axes[1].annotate("Particle 1", xy=(400, 50), xycoords="data")
axes[1].set_xlabel("X ($m$)")
axes[1].set_xlim(0, 2000)
xticks = np.arange(0, 2100, 250)
axes[1].set_xticks(xticks)
axes[1].set_ylabel(r"Temperature, $^{\circ}C$")
axes[1].set_ylim(40, 80.0)
yticks = np.arange(40, 81, 10)
axes[1].set_yticks(yticks)
# Third plot - particle travel time vs temperature
line1t = axes[2].add_collection(lc1t)
line2t = axes[2].add_collection(lc2t)
axes[2].annotate("Particle 2", xy=(15000, 68), xycoords="data")
axes[2].annotate("Particle 1", xy=(15000, 50), xycoords="data")
axes[2].set_xlabel("Time ($days$)")
axes[2].set_xlim(0, np.max([part1t["t"].max(), part2t["t"].max()]))
axes[2].set_ylabel(r"Temperature, $^{\circ}C$")
axes[2].set_ylim(40.0, 80.0)
axes[2].set_yticks(yticks)
# show and/or save figure
if plot_show:
plt.show()
if plot_save:
fpth = figs_path / "{}{}".format(sim_name, ".png")
fig.savefig(fpth, dpi=600)
Running the example
Define a function to run the example scenarios and plot results.
[16]:
def scenario(silent=False):
gwf_sim, gwe_sim, prt_sim = build_models()
if write:
write_models(gwf_sim, gwe_sim, prt_sim, silent=silent)
if run:
run_models(gwf_sim, gwe_sim, prt_sim, silent=silent)
if plot:
plot_results(gwf_sim, gwe_sim, prt_sim, silent=silent)
Run the scenario.
[17]:
scenario(silent=True)
run_models took 30024.48 ms
