Tutorial 31: Variable-Density Flow & Realtime Multi-Species Transport (SEAWAT / MODFLOW-6)

1Overview — Why Density Flow Matters

Most groundwater models assume water has constant density — one kg/m³ worth of freshwater anywhere. But in coastal aquifers, brine pockets, or around dense contaminant plumes, density differences drive flow on their own. Freshwater floats on saltwater. Dense plumes sink. Pumping near a coast can pull a saltwater wedge inland.

What this tutorial covers

You'll learn how to:

• Enable Density Flow in IGW-NET via the Solver Options dialog
• Choose between SEAWAT V4 (older USGS, 2012) and MODFLOW-6 density solvers
• Run the classic Henry Problem benchmark for saltwater intrusion
• Work with a real-world case (Salinas Valley, California) showing pumping-induced intrusion
• Use the Realtime Multi-Species Transport interface (coupled with transient flow) — contrasted with the post-analysis alternative in Tutorial 30
• Assign species properties per zone via ZoneAttr

How this fits with other tutorials

This tutorial sits at the top of the transport stack. If you're new to contaminant transport in IGW-NET, walk this path in order:

1. Tutorial 1 — Steady 2D Flow for basic model setup
2. Tutorial 5 — Contaminant Transport for single-species transport
3. Tutorial 29 — MT3D Advection Solvers for engine choice
4. Tutorial 30 — MT3DMS Multi-Species (Post-Analysis) for steady-flow reactive chains
5. This tutorial — for transient-flow, density-coupled, realtime reactive transport

You'll also want Tutorial 17 — Monte Carlo Transport if uncertainty quantification matters for your problem.

2Concept — Coupled Flow and Transport

Before touching the interface, understand what the solver actually does.

Variable-density flow is a two-way coupling

At every time step, the simulator alternates between solving flow and solving transport:

1. Solve flow for a given density field → produces heads and velocities
2. Solve transport using those velocities → produces updated concentrations
3. Update density from the new concentration field (typically via a linear density–concentration law, e.g., +0.000714 kg/m³ per ppm TDS)
4. Go to step 1 with the new density field

Two coupling strategies

• Explicit coupling (one-step lag) — flow and transport are solved once per time step with the density field from the previous step. Fast but can drift if density gradients are strong.
• Iterative coupling (Picard) — the flow-transport cycle is repeated within each time step until density stops changing. Slower but stable for strong gradients.

Two implementations in IGW-NET

3Enable Density Flow in Solver Options

Variable-density flow is not a default. You must explicitly turn it on, and you must make sure the flow and transport engines are compatible.

Step 1 — Open the Solver Options dialog

• Expand Conceptual Model Tools and click DomainAttr
• Switch to the Simulation Settings tab
• Click Solver Options

Step 2 — Select the compatible engines

In the Solver Options dialog, check the following:

1. Flow solver: select ModFlow (not the internal Magnet Flow solver — variable-density coupling requires MODFLOW).
2. Transport solver: select MT3D. For background on why and when to switch from Magnet Transport to MT3D, see Tutorial 29 — MT3D Advection Solvers.
3. Density Flow: check the box. Then click Edit SWT Solver to configure parameters.

4SEAWAT Solver Parameters

Click Edit SWT Solver in Solver Options to open the density-flow parameter dialog. At the top, select Mode: SEAWAT (the other option is MODFLOW-6, covered in the next section).

Key parameters explained

• Flow and transport coupling — "explicit (one-substep lag)" is the default and fastest. Switch to iterative if you see instability or strong density gradients.
• Initial transport substep (days) — first substep MT3D tries. Leave at 0.01 unless you have stability problems.
• Conductance scheme for density calcs — upstream (default) vs. central-in-space. Upstream is more stable; central-in-space is slightly more accurate for slow flows.
• Apply variable-density water table corrections — enable this if your model has a free water table where density varies across the interface. Corrects for the hydrostatic pressure effect.
• Density-concentration relationship — linear, defined by slope. 0.000714 kg/m³ per ppm is the standard value for NaCl solutions.
• Reference density / concentration — the baseline from which density is computed. Keep at 1000 kg/m³ / 0 ppm for freshwater-at-0 as the zero point.

5MODFLOW-6 Solver Parameters

MODFLOW-6 is the modern USGS replacement for SEAWAT. The interface offers cleaner density coupling and better numerical methods for dispersion.

What differs from SEAWAT

• Advective transport solver — MF6 ships with FDM (standard finite-difference) by default. The MOC/MMOC/HMOC/ULTIMATE schemes available in classic MT3D are not in MF6's transport model.
• Dispersion options — "Complete XT3D formulation" gives you proper anisotropic dispersion on unstructured grids. If you're doing 2D structured-grid work, this is overkill but harmless.
• Number of substeps — controls how many flow-transport iterations happen per outer step. 10 is a reasonable default.
• Off-diagonal terms to right-hand side — a numerical trick for reducing matrix asymmetry. Don't enable unless you have a specific convergence problem it solves.
• IMS6 solver — MODFLOW-6's native iterative solver. Configure via ... button if defaults don't converge.

6Benchmark — The Henry Problem

The Henry Problem (Henry, 1964) is the standard validation case for saltwater intrusion simulators. Every variable-density code ever written has been run against it. It's the "Hello World" of coupled flow–transport.

Setup

• A vertical slice of aquifer (2D cross-section)
• Isotropic, homogeneous, confined
• Constant freshwater flux on the inland (left) boundary
• Seaward (right) boundary exposed to saltwater (35,000 ppm, density 1025 kg/m³)
• Top and bottom are no-flow

The expected result: a classical saltwater wedge intrudes from the seaward boundary, with its toe position determined by the freshwater flux rate. Higher flux → wedge pushed back out to sea. Lower flux or pumping → wedge migrates inland.

7Real-World Case Study — Salinas Valley

Henry's toy problem is nice, but the point is modeling real aquifers. Salinas Valley, California is a classic saltwater intrusion site: a coastal agricultural valley where decades of pumping have drawn seawater inland from Monterey Bay.

What the model shows

The Salinas Valley model is built the same way as any IGW-NET model — import a DEM, build a conceptual model, define zones, wells, and river boundaries. The difference is that it's run with Density Flow enabled so the saltwater wedge migrates realistically with pumping scenarios.

Scenarios to explore:

• Baseline — current pumping; watch where the wedge equilibrates over decades
• Reduced pumping — cut extraction by X%; see how fast the wedge retreats
• Managed aquifer recharge — add injection wells inland of the coast; see the hydraulic barrier effect

8Realtime Multi-Species Transport

Now the other major feature of this tutorial: Realtime Multi-Species Transport — reactive transport with multiple species, fully coupled with transient flow, inside the main simulation loop.

Realtime vs. Post-Analysis — which to use?

IGW-NET offers two ways to run multi-species reactive transport. Same chemistry options, same species definitions, but different execution contexts:

What the two share

Both realtime and post-analysis use the same chemistry infrastructure:

• Same processes: linear sorption, first-order sorption, dual-domain mass transfer
• Same chemical reactions: first-order rate, zero-order, Monod kinetic, first-order chain
• Same reaction model templates: BTEX Instant Aerobic, BTEX Kinetic with Multiple EAs, PCE Sequential Decay, PCE/TCE aerobic/anaerobic, Single-pair Instant EA/ED

9Configure Realtime Multi-Species

The interface is structurally the same as Tutorial 30's post-analysis dialog, but you access it from Solver Options, not from the Analysis menu.

Step 3 — Enable MT3D in Solver Options

You've already done this if you followed Section 3 above. If not: DomainAttr → Simulation Settings → Solver Options → enable MT3D. Then click Edit MT3D Solver and configure the multi-species options within.

Step 4 — Choose Process / Reaction / Template

Exactly as in Tutorial 30, three dropdowns govern the chemistry:

Step 5 — Build/Edit defaults

Click Build/Edit to define each species' default properties (initial concentration, retardation, decay rate). These defaults apply wherever you haven't explicitly overridden them per-zone.

10Assign Species Properties to Zones

Here's a feature unique to realtime multi-species: per-zone species overrides. You can specify different initial concentrations, reaction rates, or retardation factors in different parts of the model — all within the main simulation loop.

Step 6 — Open ZoneAttr

Click a zone on the map (or right-click and choose Edit Zones), then switch to the Biochemical tab in the Zone Properties dialog.

Step 7 — Check "multiple Species" and edit per-species

At the bottom of the Biochemical tab, check multiple Species and click Edit multiple Species Property to open the per-species editor for this zone.

Typical per-zone workflows

• Source zone — set initial PCE concentration to 100 ppm, leave TCE/DCE/VC at 0 (they'll appear only as daughter products of the decay chain).
• Treatment zone — specify an accelerated decay rate to represent soil amendments or in-situ chemical oxidation.
• Tight clay lens — specify high retardation factor to represent matrix diffusion that the overall dual-domain model doesn't capture.
• Seaward coast (in combination with density flow) — set constant saltwater concentration as a boundary condition.

11Run the Simulation and View Results

Step 8 — Simulate

Click Simulate. With density flow + multi-species + transient flow all active, this is the heaviest IGW-NET simulation you can run. Expect it to be slower than a standalone transport run — SEAWAT is solving flow, MT3D is solving transport, and the coupling runs multiple times per step. For a modest 2D model, expect a few minutes; for large 3D models with chain reactions, expect tens of minutes.

Step 9 — Visualize results

While the simulation runs and after it completes, use any of these tools to visualize:

• Map view — default; shows head contours + concentration colormap. Select species and time via top-bar dropdowns.
• Cross-section — click Analysis → X-Section, draw a line on the map. Essential for seeing the saltwater wedge in vertical section.
• Display Charts — breakthrough curves at monitoring wells, plotted as time series for each species.
• 3D visualization — see Tutorial 25 — 3D Flow Visualization to render the plume in 3D with head contours.

Uncertainty quantification

If your rate coefficients, boundary concentrations, or recharge rates are uncertain, wrap the whole workflow in Monte Carlo simulation — see Tutorial 17 — Monte Carlo Transport. The realtime multi-species setup works inside the MC loop just like any other simulation.

12Challenges and Best Practices

Density-flow + reactive-transport is the most numerically demanding workflow in IGW-NET. Here's what can go wrong and how to handle it.

SEAWAT instability

The USGS SEAWAT V4 executable has known stability issues on certain configurations — convergence stalls, NaN propagation, or runaway iterations. Symptoms include:

• Simulation hangs on one timestep with high iteration count
• Concentrations explode to unphysical values (10⁹ ppm or similar)
• Heads go NaN and the whole solve fails

Workarounds:

1. Switch to MODFLOW-6 in the SWT Solver dialog. Its density coupling is better-behaved and it's the path USGS is actively maintaining.
2. Switch flow-transport coupling from explicit to iterative (Picard) in the SWT dialog.
3. Lower the initial transport substep (try 0.001 instead of 0.01).
4. Increase Max Iteration in the MT3D solver settings.

Timestep sensitivity

Variable-density flow requires timesteps that are just right:

• Too long → density changes too much between flow updates; the one-step lag becomes inaccurate; solver may diverge
• Too short → each step shows no meaningful density-driven flow change; simulation is wastefully slow

Start with a timestep around Δt ≈ (cell size) / (10 × flow velocity). If it converges cleanly, try larger steps. If it's unstable, halve the timestep and try again.

Initial conditions matter

Because density coupling is sensitive, a bad initial condition can make the first few timesteps unstable. For saltwater intrusion problems:

• Start with a known-steady-state density distribution (run a long initialization run with density coupling enabled, save, reload)
• Or start with uniformly fresh water and let the saltwater boundary propagate the wedge naturally — this is slower but more stable

13Where to Go Next

You've now covered the full transport stack in IGW-NET. A few cross-reference directions depending on what you want to build:

Other IGW-NET Quick Tutorials

• Tutorial 5 — Contaminant Transport — single-species starting point
• Tutorial 17 — Monte Carlo Transport — uncertainty quantification on top of any transport run
• Tutorial 27 — DataNET Groundwater Model — data preprocessing for real-world sites
• Tutorial 29 — MT3D Advection Solvers — choosing MMOC/FDM/MOC/HMOC/ULTIMATE
• Tutorial 30 — Multi-Species (Post-Analysis) — steady-flow alternative to this tutorial

Reference material

• IGW-NET Users' Reference Manual — detailed conceptual chapters
• Pinder Beginner's Manual — gentler mathematical background on transport and density flow
• IGW-NET Realtime Help — in-context help for every UI dialog
• — short procedures for common tasks

Other platforms

• DataNET — geospatial data preparation (DEMs, boreholes, interpolation)
• SwaNET — watershed-scale surface-water modeling (SWAT-based)
• StormNET — stormwater modeling (SWMM-based)
• ConduitNET — water distribution network modeling (EPANET-based)

External documentation

• USGS SEAWAT — the original solver reference
• USGS MODFLOW-6 — current USGS model with Groundwater Transport (GWT)
• USGS MT3D-USGS — full parameter reference