mf6adj demonstration using the San Pedro model
In this notebook, we show how mf6adj can be used with a MODFLOW 6 version of the well-known San Pedro model of Leake et al. (2010).
Import Packages
[1]:
import os
import pathlib as pl
import platform
import shutil
import sys
from datetime import datetime
import flopy
import h5py
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pyemu
[2]:
try:
import mf6adj
except ImportError:
sys.path.insert(0, str(pl.Path("../").resolve()))
import mf6adj
[3]:
mf6_bin, lib_name = mf6adj.get_conda_mf6_paths()
print(f"Using MF6 binary: {mf6_bin}")
print(f"Using MF6 library: {lib_name}")
Using MF6 binary: /home/runner/work/mf6adj/mf6adj/.pixi/envs/default/bin/mf6
Using MF6 library: /home/runner/work/mf6adj/mf6adj/.pixi/envs/default/bin/libmf6.so
Copy the MODFLOW Model
[4]:
org_ws = pl.Path("..") / "autotest" / "sanpedro" / "mf6_transient_ghb"
assert pl.Path(org_ws).exists()
Set up a local copy of the model files.
[5]:
ws = pl.Path("sanpedro")
if pl.Path(ws).exists():
shutil.rmtree(ws)
shutil.copytree(org_ws, ws)
[5]:
PosixPath('sanpedro')
[6]:
sim = flopy.mf6.MFSimulation.load(sim_ws=ws)
m = sim.get_model()
loading simulation...
loading simulation name file...
loading tdis package...
loading model gwf6...
loading package dis...
loading package rch...
loading package ghb...
loading package evt...
loading package drn...
loading package oc...
loading package npf...
loading package sto...
loading package ic...
loading package sfr...
loading solution package sp_mf6...
[7]:
ib = m.dis.idomain.array.astype(float)
ib[ib > 0] = np.nan
ib[ib <= 0] = 0.5
ib_cmap = plt.get_cmap("Greys_r")
ib_cmap.set_bad(alpha=0.0)
X, Y = m.modelgrid.xcellcenters, m.modelgrid.ycellcenters
rdf = pd.DataFrame.from_records(m.sfr.packagedata.array)
rdf["i"] = rdf.cellid.apply(lambda x: int(x[1]))
rdf["j"] = rdf.cellid.apply(lambda x: int(x[2]))
sfr_arr = np.zeros((m.dis.nrow.data, m.dis.ncol.data))
sfr_arr[rdf.i, rdf.j] = 1
sfr_arr[sfr_arr == 0] = np.nan
sfr_cmap = plt.get_cmap("Blues")
sfr_cmap.set_bad(alpha=0.0)
def plot_model(k, arr, units=None, cmap="viridis"):
arr[~np.isnan(ib[k, :, :])] = np.nan
fig, ax = plt.subplots(1, 1, figsize=(6, 6))
ax.set_aspect("equal")
cb = ax.pcolormesh(X, Y, arr, cmap=cmap)
plt.colorbar(cb, ax=ax, label=units)
ax.pcolormesh(X, Y, ib[k, :, :], cmap=ib_cmap, alpha=0.5)
ax.pcolormesh(X, Y, sfr_arr)
ax.set_xticks([])
ax.set_yticks([])
return fig, ax
[8]:
# fig, ax = plot_model(4, m.dis.botm.array[4, :, :])
fig, ax = plot_model(4, m.dis.top.array)
_ = ax.set_title("top")
Run the existing model in our local workspace.
[9]:
pyemu.os_utils.run(mf6_bin.name, cwd=ws)
mf6
MODFLOW 6
U.S. GEOLOGICAL SURVEY MODULAR HYDROLOGIC MODEL
VERSION 6.8.0.dev0+a6a4984.dirty
***DEVELOP MODE***
MODFLOW 6 compiled Jun 01 2026 18:46:14 with Intel(R) Fortran Intel(R) 64
Compiler Classic for applications running on Intel(R) 64, Version 2021.6.0
Build 20220226_000000
This software is preliminary or provisional and is subject to
revision. It is being provided to meet the need for timely best
science. The software has not received final approval by the U.S.
Geological Survey (USGS). No warranty, expressed or implied, is made
by the USGS or the U.S. Government as to the functionality of the
software and related material nor shall the fact of release
constitute any such warranty. The software is provided on the
condition that neither the USGS nor the U.S. Government shall be held
liable for any damages resulting from the authorized or unauthorized
use of the software.
MODFLOW runs in SEQUENTIAL mode
Run start date and time (yyyy/mm/dd hh:mm:ss): 2026/06/02 16:38:53
Writing simulation list file: mfsim.lst
Using Simulation name file: mfsim.nam
Solving: Stress period: 1 Time step: 1
Solving: Stress period: 1 Time step: 2
Solving: Stress period: 1 Time step: 3
Solving: Stress period: 1 Time step: 4
Solving: Stress period: 1 Time step: 5
Solving: Stress period: 1 Time step: 6
Solving: Stress period: 1 Time step: 7
Solving: Stress period: 1 Time step: 8
Solving: Stress period: 1 Time step: 9
Solving: Stress period: 1 Time step: 10
Run end date and time (yyyy/mm/dd hh:mm:ss): 2026/06/02 16:39:03
Elapsed run time: 10.240 Seconds
Normal termination of simulation.
Now plot some heads.
[10]:
hds = flopy.utils.HeadFile(pl.Path(ws) / "sp_mf6.hds")
final_arr = hds.get_data()
fig, ax = plot_model(4, final_arr[4, :, :], units="meters")
ax.set_title("layer 5 final heads")
[10]:
Text(0.5, 1.0, 'layer 5 final heads')
Nice.
The main requirement for using mf6adj is an input file that describes the performance measures. Fortunately, this file follows a modern format similar to other MF6 input files. Here we will create one programmatically. mf6adj supports so-called flux-based performance measures, which yield the sensitivity of a simulated flux to the model inputs. These measures can be described just as granularly as head-based performance measures, so let’s look at the sensitivity of simulated surface-water/groundwater flux between the groundwater system and SFR across all output times.
Run mf6adj
[11]:
pm_fname = "sfr_perfmeas.dat"
with open(pl.Path(ws) / pm_fname, "w") as fpm:
sfr_data = pd.DataFrame.from_records(m.sfr.packagedata.array)
fpm.write("begin performance_measure swgw\n")
for kper in range(sim.tdis.nper.data):
for kij in sfr_data.cellid.values:
fpm.write(
f"{kper + 1} 1 {kij[0] + 1} {kij[1] + 1} {kij[2] + 1} "
+ "sfr-1 direct 1.0 -1.0e+30\n"
)
fpm.write("end performance_measure\n\n")
Now we should be ready to go. The adjoint solution process requires one forward model run and then a solve for the adjoint state, which uses the forward solution components (for example, the conductance matrix, RHS, heads, and saturation). The adjoint solve has two important characteristics: it is linear, regardless of the forward model’s linearity, and it proceeds backward in time, starting with the last stress period.
The adjoint solve is slower than the forward run because most of the time is spent in the NumPy sparse linear solve.
[12]:
forward_hdf5_name = "forward.hdf5"
start = datetime.now()
adj = mf6adj.Mf6Adj(
pm_fname,
lib_name,
logging_level="INFO",
working_directory=ws,
)
adj.solve_forward_model(
hdf5_name=forward_hdf5_name
) # solve the standard forward solution
dfsum = adj.solve_adjoint() # solve the adjoint state for each performance measure
adj.finalize() # release components
duration = (datetime.now() - start).total_seconds()
print("took:", duration)
2026-06-02 16:39:04,229 - Logger instance 'Mf6Adj-sfr_perfmeas-11247371' created.
2026-06-02 16:39:04,230 - Running from /home/runner/work/mf6adj/mf6adj/examples/sanpedro
2026-06-02 16:39:04,231 - Structured grid found
2026-06-02 16:39:05,335 - MODFLOW 6 version: 6.8.0.dev0
2026-06-02 16:39:05,337 - Processing adjoint file: sfr_perfmeas.dat
2026-06-02 16:39:05,373 - Starting flow solution
2026-06-02 16:39:08,002 - Flow (stress period,time step) (1,1) converged in 4 iters, took 0.043785 mins
2026-06-02 16:39:08,876 - Flow (stress period,time step) (1,2) converged in 3 iters, took 0.0066188 mins
2026-06-02 16:39:09,747 - Flow (stress period,time step) (1,3) converged in 3 iters, took 0.0066478 mins
2026-06-02 16:39:10,605 - Flow (stress period,time step) (1,4) converged in 3 iters, took 0.0065127 mins
2026-06-02 16:39:11,460 - Flow (stress period,time step) (1,5) converged in 3 iters, took 0.0064692 mins
2026-06-02 16:39:12,408 - Flow (stress period,time step) (1,6) converged in 4 iters, took 0.0079772 mins
2026-06-02 16:39:13,274 - Flow (stress period,time step) (1,7) converged in 3 iters, took 0.006611 mins
2026-06-02 16:39:14,227 - Flow (stress period,time step) (1,8) converged in 4 iters, took 0.0080679 mins
2026-06-02 16:39:15,081 - Flow (stress period,time step) (1,9) converged in 3 iters, took 0.0064179 mins
2026-06-02 16:39:15,844 - Flow (stress period,time step) (1,10) converged in 2 iters, took 0.0048983 mins
2026-06-02 16:39:16,311 - Flow solution finished and took 0.1823 minutes
2026-06-02 16:39:16,356 - Starting solve_adjoint at 2026-06-02 16:39:16.356496
2026-06-02 16:39:16,366 - Structured grid found, shape: (5, 440, 320)
2026-06-02 16:39:16,371 - Starting adjoint solve for PerfMeas: swgw (kper, kstp) (1, 10)
2026-06-02 16:39:16,379 - Solving for lambda
2026-06-02 16:39:17,147 - Solving with bicgstab with solver options: {'rtol': 0.0, 'atol': 0.0, 'maxiter': 200, 'M': <115531x115531 _CustomLinearOperator with dtype=float64>}
2026-06-02 16:39:17,148 - Solver dvclose value: 1e-06
2026-06-02 16:39:17,148 - Solver rclose value: 0.001
2026-06-02 16:39:17,151 - Solver return code: 0 iterations: 0 solver norms: L2 (0.0) infinity (0.0)
2026-06-02 16:39:17,152 - Solving for lambda took: 0.773086 seconds
2026-06-02 16:39:19,032 - Adjoint solve took: 2.660634 seconds to solve adjoint solution for PerfMeas: swgw (kper, kstp) (1, 10)
2026-06-02 16:39:19,032 - Write group to hdf file
2026-06-02 16:39:19,405 - Starting adjoint solve for PerfMeas: swgw (kper, kstp) (1, 9)
2026-06-02 16:39:19,410 - Solving for lambda
2026-06-02 16:39:20,183 - Solving with bicgstab with solver options: {'rtol': 0.0, 'atol': 0.0, 'maxiter': 200, 'M': <115531x115531 _CustomLinearOperator with dtype=float64>}
2026-06-02 16:39:20,184 - Solver dvclose value: 1e-06
2026-06-02 16:39:20,185 - Solver rclose value: 0.001
2026-06-02 16:39:20,190 - Solver return code: 0 iterations: 0 solver norms: L2 (0.0) infinity (0.0)
2026-06-02 16:39:20,190 - Solving for lambda took: 0.780335 seconds
2026-06-02 16:39:22,021 - Adjoint solve took: 2.615849 seconds to solve adjoint solution for PerfMeas: swgw (kper, kstp) (1, 9)
2026-06-02 16:39:22,022 - Write group to hdf file
2026-06-02 16:39:22,394 - Starting adjoint solve for PerfMeas: swgw (kper, kstp) (1, 8)
2026-06-02 16:39:22,398 - Solving for lambda
2026-06-02 16:39:23,166 - Solving with bicgstab with solver options: {'rtol': 0.0, 'atol': 0.0, 'maxiter': 200, 'M': <115531x115531 _CustomLinearOperator with dtype=float64>}
2026-06-02 16:39:23,167 - Solver dvclose value: 1e-06
2026-06-02 16:39:23,167 - Solver rclose value: 0.001
2026-06-02 16:39:23,171 - Solver return code: 0 iterations: 0 solver norms: L2 (0.0) infinity (0.0)
2026-06-02 16:39:23,172 - Solving for lambda took: 0.774286 seconds
2026-06-02 16:39:25,024 - Adjoint solve took: 2.63008 seconds to solve adjoint solution for PerfMeas: swgw (kper, kstp) (1, 8)
2026-06-02 16:39:25,025 - Write group to hdf file
2026-06-02 16:39:25,394 - Starting adjoint solve for PerfMeas: swgw (kper, kstp) (1, 7)
2026-06-02 16:39:25,398 - Solving for lambda
2026-06-02 16:39:26,172 - Solving with bicgstab with solver options: {'rtol': 0.0, 'atol': 0.0, 'maxiter': 200, 'M': <115531x115531 _CustomLinearOperator with dtype=float64>}
2026-06-02 16:39:26,173 - Solver dvclose value: 1e-06
2026-06-02 16:39:26,173 - Solver rclose value: 0.001
2026-06-02 16:39:26,186 - Solver return code: 0 iterations: 0 solver norms: L2 (0.0) infinity (0.0)
2026-06-02 16:39:26,186 - Solving for lambda took: 0.78807 seconds
2026-06-02 16:39:28,043 - Adjoint solve took: 2.648763 seconds to solve adjoint solution for PerfMeas: swgw (kper, kstp) (1, 7)
2026-06-02 16:39:28,043 - Write group to hdf file
2026-06-02 16:39:28,407 - Starting adjoint solve for PerfMeas: swgw (kper, kstp) (1, 6)
2026-06-02 16:39:28,411 - Solving for lambda
2026-06-02 16:39:29,172 - Solving with bicgstab with solver options: {'rtol': 0.0, 'atol': 0.0, 'maxiter': 200, 'M': <115531x115531 _CustomLinearOperator with dtype=float64>}
2026-06-02 16:39:29,173 - Solver dvclose value: 1e-06
2026-06-02 16:39:29,174 - Solver rclose value: 0.001
2026-06-02 16:39:29,180 - Solver return code: 0 iterations: 0 solver norms: L2 (0.0) infinity (0.0)
2026-06-02 16:39:29,185 - Solving for lambda took: 0.774617 seconds
2026-06-02 16:39:31,038 - Adjoint solve took: 2.631226 seconds to solve adjoint solution for PerfMeas: swgw (kper, kstp) (1, 6)
2026-06-02 16:39:31,038 - Write group to hdf file
2026-06-02 16:39:31,401 - Starting adjoint solve for PerfMeas: swgw (kper, kstp) (1, 5)
2026-06-02 16:39:31,406 - Solving for lambda
2026-06-02 16:39:32,172 - Solving with bicgstab with solver options: {'rtol': 0.0, 'atol': 0.0, 'maxiter': 200, 'M': <115531x115531 _CustomLinearOperator with dtype=float64>}
2026-06-02 16:39:32,172 - Solver dvclose value: 1e-06
2026-06-02 16:39:32,173 - Solver rclose value: 0.001
2026-06-02 16:39:32,177 - Solver return code: 0 iterations: 0 solver norms: L2 (0.0) infinity (0.0)
2026-06-02 16:39:32,177 - Solving for lambda took: 0.771512 seconds
2026-06-02 16:39:34,017 - Adjoint solve took: 2.615502 seconds to solve adjoint solution for PerfMeas: swgw (kper, kstp) (1, 5)
2026-06-02 16:39:34,018 - Write group to hdf file
2026-06-02 16:39:34,382 - Starting adjoint solve for PerfMeas: swgw (kper, kstp) (1, 4)
2026-06-02 16:39:34,386 - Solving for lambda
2026-06-02 16:39:35,155 - Solving with bicgstab with solver options: {'rtol': 0.0, 'atol': 0.0, 'maxiter': 200, 'M': <115531x115531 _CustomLinearOperator with dtype=float64>}
2026-06-02 16:39:35,156 - Solver dvclose value: 1e-06
2026-06-02 16:39:35,156 - Solver rclose value: 0.001
2026-06-02 16:39:35,158 - Solver return code: 0 iterations: 0 solver norms: L2 (0.0) infinity (0.0)
2026-06-02 16:39:35,159 - Solving for lambda took: 0.773012 seconds
2026-06-02 16:39:36,989 - Adjoint solve took: 2.606825 seconds to solve adjoint solution for PerfMeas: swgw (kper, kstp) (1, 4)
2026-06-02 16:39:36,989 - Write group to hdf file
2026-06-02 16:39:37,350 - Starting adjoint solve for PerfMeas: swgw (kper, kstp) (1, 3)
2026-06-02 16:39:37,355 - Solving for lambda
2026-06-02 16:39:38,115 - Solving with bicgstab with solver options: {'rtol': 0.0, 'atol': 0.0, 'maxiter': 200, 'M': <115531x115531 _CustomLinearOperator with dtype=float64>}
2026-06-02 16:39:38,116 - Solver dvclose value: 1e-06
2026-06-02 16:39:38,116 - Solver rclose value: 0.001
2026-06-02 16:39:38,119 - Solver return code: 0 iterations: 0 solver norms: L2 (0.0) infinity (0.0)
2026-06-02 16:39:38,120 - Solving for lambda took: 0.76571 seconds
2026-06-02 16:39:39,957 - Adjoint solve took: 2.607225 seconds to solve adjoint solution for PerfMeas: swgw (kper, kstp) (1, 3)
2026-06-02 16:39:39,958 - Write group to hdf file
2026-06-02 16:39:40,320 - Starting adjoint solve for PerfMeas: swgw (kper, kstp) (1, 2)
2026-06-02 16:39:40,325 - Solving for lambda
2026-06-02 16:39:41,084 - Solving with bicgstab with solver options: {'rtol': 0.0, 'atol': 0.0, 'maxiter': 200, 'M': <115531x115531 _CustomLinearOperator with dtype=float64>}
2026-06-02 16:39:41,085 - Solver dvclose value: 1e-06
2026-06-02 16:39:41,085 - Solver rclose value: 0.001
2026-06-02 16:39:41,092 - Solver return code: 0 iterations: 0 solver norms: L2 (0.0) infinity (0.0)
2026-06-02 16:39:41,092 - Solving for lambda took: 0.767647 seconds
2026-06-02 16:39:42,961 - Adjoint solve took: 2.640921 seconds to solve adjoint solution for PerfMeas: swgw (kper, kstp) (1, 2)
2026-06-02 16:39:42,962 - Write group to hdf file
2026-06-02 16:39:43,319 - Starting adjoint solve for PerfMeas: swgw (kper, kstp) (1, 1)
2026-06-02 16:39:43,324 - Solving for lambda
2026-06-02 16:39:44,088 - Solving with bicgstab with solver options: {'rtol': 0.0, 'atol': 0.0, 'maxiter': 200, 'M': <115531x115531 _CustomLinearOperator with dtype=float64>}
2026-06-02 16:39:44,089 - Solver dvclose value: 1e-06
2026-06-02 16:39:44,089 - Solver rclose value: 0.001
2026-06-02 16:39:45,505 - Custom convergence reached: dvmax = 5.672088315032339e-07 < 1e-06 rmax = 4.751487018950229e-05 < 0.001
2026-06-02 16:39:45,506 - Solver return code: 0 iterations: 38 solver norms: L2 (0.0005934208908827322) infinity (4.751487018950229e-05)
2026-06-02 16:39:45,508 - Solving for lambda took: 2.183584 seconds
2026-06-02 16:39:47,380 - Adjoint solve took: 4.061453 seconds to solve adjoint solution for PerfMeas: swgw (kper, kstp) (1, 1)
2026-06-02 16:39:47,381 - Write group to hdf file
2026-06-02 16:39:47,747 - Formulate composite sensitivities
2026-06-02 16:39:47,748 - Writing composite sensitivities
2026-06-02 16:39:48,032 - Adjoint solve took: 31.676 seconds for pm 'swgw' at (kper,kstp) (1, 1)
2026-06-02 16:39:48,035 - Finalizing Mf6Adj
took: 43.963097
Plot the Results
Done. Let’s see what happened.
[13]:
[f for f in os.listdir(ws) if f.endswith("hdf5")]
[13]:
['forward.hdf5', 'adjoint_solution_swgw.hdf5']
mf6adj uses the widely available HDF5 format to store information. These files hold low-level details about the adjoint solution. However, mf6adj.solve_adjoint() also returns a higher-level summary of the results. Let’s look at that first:
[14]:
type(dfsum)
[14]:
dict
[15]:
list(dfsum.keys())
[15]:
['swgw']
[16]:
dfhw = dfsum["swgw"]
dfhw
[16]:
| k11 | k33 | wel6_q | rch6_recharge | drn-1_elev | drn-1_cond | ghb-1_bhead | ghb-1_cond | sfr-1_stage | sfr-1_cond | ss | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| node | |||||||||||
| 26747 | -1.313218e+00 | -0.154118 | -9.281390e-01 | -9.281390e-01 | 0.0 | 0.0 | 0.0 | 0.0 | 9951.831461 | -0.521997 | -2.402062e-05 |
| 26748 | -1.418439e-02 | -0.119453 | -9.366011e-01 | -9.366011e-01 | 0.0 | 0.0 | 0.0 | 0.0 | 8938.324159 | 0.111638 | -1.631893e-05 |
| 27067 | -1.716591e+00 | -0.033958 | -4.548739e-01 | -4.548739e-01 | 0.0 | 0.0 | 0.0 | 0.0 | 164.653497 | 0.613565 | -6.044541e-06 |
| 27068 | 1.089327e+00 | -0.056021 | -5.597631e-01 | -5.597631e-01 | 0.0 | 0.0 | 0.0 | 0.0 | 7163.251654 | -0.241463 | -1.072195e-05 |
| 27387 | 2.995759e-02 | -0.008328 | -3.346814e-01 | -3.346814e-01 | 0.0 | 0.0 | 0.0 | 0.0 | 0.000000 | 0.000000 | -1.957426e-06 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 694854 | -5.237353e-18 | 0.000000 | -3.019869e-20 | -3.019869e-20 | 0.0 | 0.0 | 0.0 | 0.0 | 0.000000 | 0.000000 | -1.413234e-24 |
| 694855 | -1.294803e-17 | 0.000000 | -5.065685e-20 | -5.065685e-20 | 0.0 | 0.0 | 0.0 | 0.0 | 0.000000 | 0.000000 | 7.057297e-25 |
| 694856 | -1.471397e-17 | 0.000000 | -7.949939e-20 | -7.949939e-20 | 0.0 | 0.0 | 0.0 | 0.0 | 0.000000 | 0.000000 | -1.728057e-23 |
| 695174 | -8.098884e-19 | 0.000000 | -1.564736e-20 | -1.564736e-20 | 0.0 | 0.0 | 0.0 | 0.0 | 0.000000 | 0.000000 | -1.470258e-24 |
| 695175 | -2.228106e-18 | 0.000000 | -2.045437e-20 | -2.045437e-20 | 0.0 | 0.0 | 0.0 | 0.0 | 0.000000 | 0.000000 | 7.405923e-24 |
115531 rows × 11 columns
Those are the node-scale sensitivities for the SFR flux-based performance measure. Plotting them is easiest using the HDF5 file itself.
[17]:
result_hdf = "adjoint_solution_swgw.hdf5"
hdf = h5py.File(pl.Path(ws) / result_hdf, "r")
keys = list(hdf.keys())
keys.sort()
print(keys)
['composite', 'solution_kper:00000_kstp:00000', 'solution_kper:00000_kstp:00001', 'solution_kper:00000_kstp:00002', 'solution_kper:00000_kstp:00003', 'solution_kper:00000_kstp:00004', 'solution_kper:00000_kstp:00005', 'solution_kper:00000_kstp:00006', 'solution_kper:00000_kstp:00007', 'solution_kper:00000_kstp:00008', 'solution_kper:00000_kstp:00009']
The composite group contains sensitivities of the performance measure to the model inputs, summed across all adjoint solutions.
[18]:
grp = hdf["composite"]
plot_keys = [
i
for i in grp.keys()
if len(grp[i].shape) == 3 and ("k11" in i or "wel" in i or "ss" in i)
]
plot_keys
[18]:
['k11', 'ss', 'wel6_q']
Here is a simple routine to plot these sensitivities.
[19]:
for pkey in plot_keys:
arr = grp[pkey][:]
for k, karr in enumerate(arr):
karr[karr == 0.0] = np.nan
fig, ax = plot_model(k, karr)
ax.set_title(pkey + f", layer:{k + 1}", loc="left")
One especially well-known plot in this context is the so-called capture fraction map:
[20]:
arr = grp["wel6_q"][3, :, :]
arr[arr == 0.0] = np.nan
fig, ax = plot_model(3, np.abs(arr), cmap="viridis_r")
ax.set_title("capture fraction layer 4", loc="left");
[21]:
arr = grp["wel6_q"][4, :, :]
arr[arr == 0.0] = np.nan
fig, ax = plot_model(4, np.abs(arr), cmap="viridis_r")
ax.set_title("capture fraction layer 5", loc="left");
This plot shows capture fraction: the proportion of groundwater captured from the simulated surface-water/groundwater flux if a groundwater well were added in a given model cell. Traditionally, this would be calculated by adding a well or other specified-flux boundary in each model cell, rerunning the model, recording how the surface water/groundwater flux changed, and normalizing that change by the applied rate for every active cell. That can take a long time. With the adjoint approach, the so-called adjoint state for this performance measure is simply the negative of the capture fraction: it shows how the simulated gw-sw flux changes in response to a unit injection of water in each active model cell.