Coupled Lake-Aquifer Modeling — When Stage Becomes a Solution

15.1 The Conceptual Shift — From BC-Mode to Coupled-Mode

Chapter 14 ended by naming the boundary: SW-as-BC works when stage is something you know and the aquifer is something you solve; it fails when stage itself is something that needs to be solved. This section crystallizes what changes when you make that move.

15.1.1 What Level 2 does and doesn't do

In Level 2 (Ch. 14), a two-way head-dependent lake is configured with:

• A prescribed stage — Const value, TopE from DEM, or transient file (time series of stages from a gauge, say)
• A leakance — controlling the strength of two-way exchange
• A bed bottom — for the connected/disconnected regime logic
• Two-way flux physics — Q = leakance × (h_lake − h_aq), bidirectional as h_aq varies across the lake footprint (the flow-through lake case from Ch. 14 §14.6.4)

What Level 2 does well: gives you regional flow patterns consistent with a known or externally-determined lake stage. Handles flow-through lakes correctly (gaining upstream, losing downstream). Handles seasonal stress patterns when you feed in a transient stage file.

What Level 2 doesn't do: let the lake's stage respond to what the aquifer is doing. If you pump a high-capacity well near a Level-2 lake, the aquifer head drops, the lake starts losing more water (correctly), but the lake stage doesn't change in the model because it was specified as an input. Over a long simulation, this is unphysical — in reality, sustained losing flux would draw the lake down, reducing the stage, reducing the losing flux. Level 2 can't represent this feedback loop.

15.1.2 What coupled (Level 3) adds

In Level 3 coupled modeling, the lake has its own water balance equation. Every time step:

1. The lake's Ins (precipitation, runoff, SW inflow, gaining GW flux) and Outs (evaporation, SW outflow, losing GW flux) are accounted for
2. Net storage change = ΣIns − ΣOuts → change in lake volume → change in lake stage (via bathymetry)
3. The new stage is used for the next time step's aquifer solve

Concretely, this unlocks behaviors that Level 2 cannot represent:

• Lake drawdown from pumping. Pump → aquifer head drops near lake → losing flux increases → lake volume decreases → lake stage drops → new equilibrium at lower stage with reduced pumping impact on aquifer (because the lake is now lower relative to the aquifer). Only possible when stage responds to the coupling.
• Seasonal lake-level response. Wet year: more rainfall + runoff + gaining flux → lake rises. Dry year: evaporation + reduced GW inflow → lake drops. Seasonal dynamics emerge from the water balance.
• Transient equilibration. You start a simulation with an initial stage. Over time, the stage moves toward a dynamic equilibrium consistent with the long-term Ins and Outs. The equilibrium stage is a prediction, not an input.
• Lake augmentation scenarios. Pumping groundwater into the lake (like Barron Lake) raises the stage, affecting everything downstream. The model predicts the resulting stage from the augmentation rate, rather than requiring you to know it.
• Integrated water accounting. Volumes of water moving through the lake over the simulation period are tracked explicitly — useful for water-rights, ecological-flow, and drought analyses.

15.1.3 Why this reframing matters for modeling practice

The move from Level 2 to Level 3 is more than a technical refinement — it shifts what kinds of questions the model can answer:

15.2 The Lake as One SW Cell — Water Balance and Numerical Architecture

The key modeling abstraction that makes coupled lake-aquifer simulation tractable: treat the lake as a single surface-water "cell" with its own water balance. No spatial resolution within the lake; lakes are flat in the small-amplitude limit, so one stage value applies throughout. Only temporal dynamics matter. This section walks through the water balance equation, the Ins and Outs, and the numerical scheme that couples the lake and aquifer solves.

15.2.1 The lake water balance equation

The lake is a single storage reservoir. Its water balance over a time step:

This is a single scalar ordinary differential equation (ODE) in the unknown h_lake(t). No partial differential equation; no spatial discretization within the lake. The entire lake is one computational unit — one stage value, one volume, one water balance. What it has instead of spatial resolution is temporal dynamics: the stage changes through time as the Ins and Outs accumulate.

15.2.2 The Ins — where water enters the lake

15.2.3 The Outs — where water leaves the lake

The gaining-vs-losing distinction is per aquifer cell, not per lake. A flow-through lake (§14.6.4) has some cells gaining (upstream side, where h_aq > h_lake) and some losing (downstream side, where h_aq < h_lake) simultaneously. The GW flux term in the lake's water balance is the net across all cells — it might be slightly positive (lake is net gaining) or slightly negative (lake is net losing) even though individual cells go in opposite directions.

15.2.4 The explicit time scheme — how the coupling works numerically

The crucial design decision: IGW-NET uses an explicit time scheme to couple the lake water balance with the aquifer flow solve. This decouples the two sub-problems within each time step and keeps the numerics simple, fast, and robust.

The algorithm, per time step (from step n to step n+1):

15.2.5 The trade-off — time-step sensitivity

Explicit schemes have stability limits. The "known stage" used in step n's aquifer solve is h_lakeⁿ — but the lake stage is actually changing throughout step n. If the stage changes substantially within a single time step, the value used for the aquifer solve is stale, and errors accumulate. Past a certain time step size, errors compound rather than damp out, and the solution oscillates or diverges.

Symptoms of too-large time step:

• Oscillating lake stage. Stage jumps up and down between successive time steps rather than evolving smoothly.
• Non-physical stage excursions. Lake stage drops below plausible minimums (below the bathymetric floor!) or rises above plausible maximums.
• Mass balance errors growing through the simulation. Each step introduces small errors that accumulate.
• Explicit divergence warnings from the solver.
• "Worked fine until I enabled coupling" — the single most common user report.

15.2.6 What the scheme does NOT handle well

A few situations where the explicit scheme shows its limits, and how to recognize them:

• Very small lakes with very fast dynamics. A kettle pond with high leakance and significant pumping nearby can have stage changing rapidly. You may need very small sub-time-steps to stay stable — so small that the simulation becomes slow. In these cases, question whether Level 3 is really needed or whether a Level-2 transient stage file (from a separate hydrologic model of the pond) is more practical.
• Step changes in stresses. If a pumping rate or SW inflow changes instantaneously (e.g., a pump turns on at day 30), the immediate response can be hard to resolve even with small time steps. Ramp stress changes over a few time steps when possible.
• Lake going dry. When a lake approaches zero volume, the bathymetric stage-volume curve becomes ill-conditioned (small volume changes produce large stage changes). The scheme can become unstable. In practice, most modeling contexts avoid near-dry conditions; if you must model them, use very small time steps and consider a minimum-depth cap.

15.3 The UI — "From Coupled SW/GW Modeling" and What Configures It

The conceptual shift of §15.1 and the numerical architecture of §15.2 map to specific UI elements in IGW-NET. This section walks through the dialogs you actually interact with, using the Barron Lake configuration as the running example.

15.3.1 The setup — from lake polygon to coupled behavior

The workflow starts from a Level-2 two-way lake (Ch. 14 §14.2.2) and adds the coupling via a single checkbox. Steps:

Draw or import the lake polygon as a zone feature

Use ZonePoly (Conceptual Model Tools → Zones → ZonePoly) to draw the lake boundary, or Zone from a shapefile for imported geometry. For the Barron Lake case, the polygon was imported from BarronLake.txt. SaveShape finalizes the polygon and opens the Zone Attributes dialog.

Open Zone Attributes → Sources and Sinks Head Dependent

This is where the Level-2 vs Level-3 decision gets made. The Two-way Head Dependent section has the familiar Level-2 fields (stage, leakance, bed bottom) from Ch. 14 — plus the critical Level-3 checkbox.

Check "Two-way Head Dependent" and configure the lake basics

Lake name; initial stage value (as Const — the initial condition for the coupled simulation); leakance (0.1/day in Barron Lake); River Bed set to a sentinel value like -999999 if bathymetry will handle it; IsBathymetry option if you have a bathymetry file.

The key step — under Stage, check "From coupled SW/GW modeling"

This is the Level-3 primitive. Under the Stage section of the dialog, radio buttons or checkboxes select the stage source: Const, TopE (DEM), Transient (from file), or From coupled SW/GW modeling. Selecting "From coupled SW/GW modeling" tells IGW-NET to compute stage as a solution variable from the lake's water balance rather than taking it as an input. This single selection is what distinguishes a Level-3 coupled lake from a Level-2 prescribed-stage lake.

15.3.2 Bathymetry import — the IsBathymetry option

Without bathymetry, IGW-NET treats the lake as a flat-bottomed prism — constant area × depth. This works fine for small stage changes but produces errors when the lake's stage-area or stage-volume relationship is strongly nonlinear (most real lakes have basins that narrow at the bottom, so area decreases with depth).

The IsBathymetry option in the Zone Attributes dialog lets you import a bathymetry file (WaterDepth.txt in the Barron Lake case) that specifies the lake's depth at a grid of points or a rasterized field. IGW-NET uses this to compute the stage-area and stage-volume relationships accurately:

• At each stage, the submerged area is the integral of "depth greater than zero" across the lake footprint
• At each stage, the volume is the integral of submerged depth across the lake footprint
• These relationships are interpolated efficiently from the bathymetry grid

15.3.3 Surface Water Source and Sink — configuring the Ins and Outs

Once "From coupled SW/GW modeling" is enabled, you need to specify all the Ins and Outs from §15.2.2 and §15.2.3. The Surface Water Source and Sink dialog (accessible from the Zone Attributes Head Dependent tab when coupled is on) exposes these:

You enable the checkboxes for the Ins and Outs that apply to your lake. A simple closed-basin lake might only need Precipitation, Evaporation, and the automatic GW flux (which is always on for coupled lakes). A managed reservoir needs SW Prescribed Source/Sink for controlled inflows, Hydraulic Outlet for controlled discharge, plus the automatic terms. The Barron Lake case uses Precipitation, Evaporation, SW Prescribed Source/Sink (for augmentation pumping), and SW Stage Measurement (for calibration).

15.3.4 Sub-time-step configuration — managing stability

As §15.2.5 established, the explicit coupling scheme has a stability limit set by the time step. The SW sub-time-step is the operational control — a parameter you set in Simulation Settings to govern how finely the lake water balance integrates.

The relationship: the GW simulation has its own time step (typical: 1 day for transient problems). Within each GW step, the SW sub-time-step further subdivides for the coupled lake water balance. For example:

• GW time step = 1 day, SW sub-time-step = 0.1 day → 10 sub-steps of lake balance integration per GW step
• GW time step = 1 day, SW sub-time-step = 0.01 day → 100 sub-steps — much finer, more stable, but slower
• For lakes with slow dynamics, SW sub-time-step = GW time step (no subdivision) may be fine

The sub-time-step is your primary stability control. When coupling breaks convergence, this is the first knob to turn.

15.4 The Barron Lake Case Study — End-to-End Coupled Modeling

Barron Lake, in Niles, Michigan, is a small lake with a real hydrologic challenge: low natural inflow relative to evaporation and groundwater losses means stakeholders rely on pumping augmentation to maintain lake levels during dry periods. The coupled model answers specific questions: how much augmentation is needed? What stage response can we expect? How do 11 private wells near the lake respond?

15.4.1 The site

Barron Lake (41.8479°N, 86.1756°W; Niles, southwest Michigan) is set in glacial terrain with complex stratigraphy — outwash sands, till layers, and bedrock. The aquifer is unconfined and tightly coupled to the lake through sandy lake-bed sediments. Observations cover 2006–2013 (lake stage and nearby private-well water levels); the modeled period is a 135-day transient simulation starting 6/25/2013.

Key physical parameters:

• Lake area: ~1 km² (variable with stage)
• Initial lake stage: 230.0259 m (surveyed reference elevation)
• Lake leakance: 0.1 per day
• Bathymetry imported from WaterDepth.txt (depth values at a fine grid over the lake footprint)
• 11 private monitoring wells within ~2 km of the lake
• 3D geology from T-PROGS (4 materials, K = 53, 5, 0.01, 0.001 ft/day)
• Augmentation pumping rates from Qsw.csv (time-varying, source to the lake)
• Observed lake levels from LakeLevel2013.csv (calibration target)

15.4.2 The coupled configuration

The model integrates all of IGW-NET's major subsystems:

• Coupled lake (this chapter): Level-3 coupled lake with "From coupled SW/GW modeling" enabled; bathymetry imported; SW Prescribed Source/Sink configured with Qsw.csv; SW Stage Measurement with LakeLevel2013.csv; Precipitation and Evaporation from Data Center climate
• T-PROGS 3D geology (Ch. 17): DM Zone with Borehole Simulation selected; BarronLakeTP.tp imported; Kxx values for the 4 materials
• Watershed Solver (Ch. 16): Enabled for runoff calculation; Surface Water Overland Flow Options configured with Land Use, Climate, and Soil Type from the Data Center
• Augmentation well: Pumping well placed at the lake augmentation location; transient pumping rates linked to the SW Prescribed Source/Sink
• Simulation Settings: 135-day transient; SW sub-time-step tuned for stability

15.4.3 What the coupled model produces

The 135-day transient simulation produces a time series of coupled outputs:

• Lake stage over time — compared against LakeLevel2013.csv observations (calibration target)
• Aquifer head fields at each time step — compared against the 11 private well records
• Water balance components — each In and Out tracked explicitly, so you can see how the lake's total water budget resolves: how much came from GW vs runoff vs augmentation, how much went to evaporation vs GW losses vs outflow
• GW-lake flux breakdown per cell — which parts of the lake are gaining, which are losing, how the pattern evolves
• Well capture zones — where does water coming out of the augmentation well originate? With coupled lake, it can come partly from the lake itself (recirculation), which a Level-2 treatment would miss

15.4.4 Calibration workflow

The Barron Lake model is calibrated against both the lake stage observations and the 11 private-well observations using a typical UCODE workflow (Ch. 18):

1. Initial run with default parameters; compare simulated lake stage to LakeLevel2013.csv
2. Adjust lake leakance and aquifer K multipliers to bring simulated stage and well heads closer to observations
3. Iterate — including re-running Watershed Solver runoff calculations as needed
4. Final calibrated model matches observed lake stage and well heads within acceptable tolerances

The value of calibrating against both lake stages and well heads: they constrain different aspects of the model. Lake stages constrain the lake water balance and lake leakance. Well heads constrain the aquifer K field and the lateral flow patterns. Together they pin down a more trustworthy model than either set of observations alone.

15.5 Coupling with Other IGW-NET Subsystems

Coupled lake-aquifer modeling rarely exists in isolation. Lakes receive runoff (Watershed Solver), sit in heterogeneous geology (T-PROGS), respond to calibrated pumping scenarios (UCODE), and deliver results to transport calculations (Ch. 12). This section summarizes how the coupled lake integrates with the rest of the platform.

15.5.1 Watershed Solver — the runoff link

The Watershed Solver (Ch. 16) models overland flow, infiltration, and runoff across the surrounding watershed. When it's enabled and a coupled lake is present, the runoff reaching the lake automatically appears as the "Runoff" In in the lake water balance — no manual CSV import required. Rainfall on the watershed → infiltration or runoff → routed overland → some fraction enters the lake. The coupling is automatic once both systems are active.

For lakes where watershed runoff is a significant fraction of the Ins (most natural lakes, especially in humid climates), the Watershed Solver link is essential. Without it, you're either ignoring runoff entirely (underestimating Ins) or approximating it with a user-specified constant rate (missing temporal dynamics).

15.5.2 T-PROGS 3D geology — realistic aquifer under the lake

T-PROGS (Ch. 17) produces realistic 3D heterogeneous aquifer conductivity fields from borehole data. For a coupled lake, the aquifer cells beneath the lake polygon use the T-PROGS-simulated K field, so the two-way head-dependent flux reflects real geological heterogeneity — not a uniform K assumption. The Barron Lake case uses T-PROGS with 4 material types (K = 53, 5, 0.01, 0.001 ft/day) to capture the glacial outwash-and-till stratigraphy beneath the lake.

This matters because lake-aquifer exchange rates depend critically on the K near the lake-bed interface. A high-K outwash layer adjacent to the lake means strong exchange; a low-K till layer means weak exchange. T-PROGS gives you the geologically-realistic variability rather than forcing a uniform-K assumption.

15.5.3 Calibration with UCODE

Ch. 18 covers UCODE calibration in depth. For coupled lakes, two kinds of calibration targets matter:

• Lake stage observations (via SW Stage Measurement) — constrain the lake water balance directly
• Aquifer head observations at nearby wells — constrain the aquifer K field and the coupling leakance

Calibrating both jointly is what's done at Barron Lake and is generally the right approach — either class of observation alone gives you a single-perspective view; together they constrain the model from multiple angles.

15.5.4 Transient BCs — seasonal and scenario-based runs

Coupled lakes naturally run as transient simulations. Ch. 9 covers transient setup. For coupled lakes specifically:

• Initial stage = starting lake level at t = 0
• Transient SW Prescribed Source/Sink for augmentation, withdrawal, or tributary inflow time series
• Transient climate (precipitation, evaporation) from the Data Center
• Transient pumping from nearby wells
• The lake stage trajectory emerges as a solution variable, responding to all the above

This enables "what if" scenario analysis: what if we pump more? what if the drought continues? what if we install a new well? Each scenario produces a complete predicted time series of lake stage, aquifer heads, and water balance components.

15.5.5 Transport — contaminants through coupled systems

Ch. 12 covers transport. For coupled lakes, transport adds another layer: contaminants in the lake migrate into the aquifer on losing cells and arrive at downstream wells or receptors. The coupled lake model provides the correct flux patterns; transport uses them to route mass. This is especially relevant for lake contamination scenarios (spills, algal blooms producing toxins) where lake water entering the aquifer carries a chemical signal.

15.6 When to Escalate Level 2 → Level 3

Most groundwater modeling never needs Level 3. Ch. 14 Levels 1 and 2 handle the vast majority of questions. But for specific questions, Level 3 is essential — and this section helps you recognize them.

15.6.1 The decision criteria

Escalate from Level 2 to Level 3 when any of the following applies:

• Stage is a thing you need to solve for. The question is "what will the lake stage be under scenario X?" rather than "what will heads be given this stage?"
• Stage-dependent feedback matters. The lake's stage responds to pumping or climate in ways that affect the answer — if you can't let the stage respond, your prediction is wrong.
• Seasonal or multi-year dynamics. The simulation period is long enough that stage changes meaningfully through it; a fixed stage value misrepresents the state.
• Lake water balance itself is a modeling target. You're asked to report how much water came from GW vs runoff vs augmentation — the answer requires explicit water-balance accounting.
• Management scenarios involving lake levels. Augmentation studies, withdrawal scenarios, reservoir operations, habitat management.
• Integrated GW-SW systems where the coupling is the physics. Not just "GW affects lake" or "lake affects GW" but both together, in feedback.

Do NOT escalate to Level 3 when:

• You have observed or externally-modeled stage time series — use Level-2 transient stage input (simpler, same answer for aquifer questions)
• The simulation is steady-state and stage is at a known equilibrium
• Your question is aquifer-only; lake stage is just a BC for the aquifer question
• You don't have the data to parameterize the lake water balance (unknown evaporation, unknown runoff, unknown inflows) — Level 3 with bad inputs is worse than Level 2 with observed stage

15.6.2 The cost-benefit calculation

15.7 Wetlands and Beyond — The Same Framework, Applied Broadly

Lakes are the most common coupled feature, but the framework (one SW cell with water balance, two-way head-dependent aquifer below, stage as solution variable, explicit time scheme) applies equally well to wetlands, reservoirs, and other surface-water bodies. This section briefly covers the variants.

15.7.1 Wetlands

A wetland is, for modeling purposes, a shallow lake with distinctive properties. The same Level-3 configuration works:

• Draw the wetland polygon as a zone
• Check Two-way Head Dependent; set initial stage and leakance (wetland leakance is often lower than open-lake leakance due to organic sediments)
• Check "From coupled SW/GW modeling" under Stage
• Configure the Ins and Outs, with wetland-specific attention to:
  • Evaporation — usually dominant. Wetlands have high ET rates (actual evapotranspiration from wetland vegetation, not just lake-surface evaporation). This often drives the water balance more than any other term.
  • Runoff — often significant. Wetlands occupy topographic lows; their contributing watershed funnels water toward them.
  • GW flux — often bidirectional. Wetlands in flat terrain commonly have GW discharge to the wetland on one side and wetland seepage to the aquifer on the other, like flow-through lakes (Ch. 14 §14.6.4)
  • Bathymetry — simpler. Wetlands are usually close to flat-bottomed; the bathymetric correction is smaller than for lakes

The "From coupled SW/GW modeling" primitive is what makes wetland stage a solution variable. Without it, you'd have a fixed wetland stage that doesn't respond to GW decline or recovery — inadequate for GDE (groundwater-dependent ecosystem) questions, drought impact studies, or climate-change wetland response.

15.7.2 Reservoirs

Engineered reservoirs work the same way, with additions:

• Controlled outflow. Hydraulic Outlet with a rating curve or prescribed operating rule (release schedule, target stage, downstream flow requirement)
• Known operating target stage. Can be useful for validation; the model should produce stages near the operating target under normal conditions
• SW Prescribed Source/Sink for upstream inflow — from tributaries or upstream reservoir discharge
• Withdrawal for water supply — another SW Prescribed Source/Sink (negative rate)

Reservoir-aquifer modeling is often used for water-supply planning, drought-response studies, and environmental-flow compliance — all of which benefit from Level-3 coupling.

15.7.3 Streams — coupled routing

Coupled stream modeling (stage as a solution variable for the stream rather than a fixed input) is qualitatively more complex than coupled lakes because streams have spatial dimension — stage varies along the reach. IGW-NET's approach is to couple stream flow routing with the Watershed Solver (Ch. 16), which handles overland and channel flow simultaneously. For a basic introduction to coupled stream-aquifer modeling, see Ch. 16. Full coupled stream-routing with transient stage is beyond the scope of this chapter.

15.7.4 The unifying pattern