Lake Rotorua Integrated Modelling — Project Update

Project Overview

• Overarching goal: Develop an integrated catchment–lake modelling framework to support PC10 review and implementation
• Two workstreams running in parallel:
  • Workstream 1 — Catchment (SWAT+gwflow): Pilot application in three sub-catchments (Puarenga, Hamurana, Waingaehe) to simulate nutrient transport via surface runoff and groundwater pathways; running November 2025 – April 2026
  • Workstream 2 — Lake (GLM-AED): Upgrade of the lake model to GLM-AED via the AEME framework, enabling improved simulation of vertical mixing, nutrient cycling, sediment feedbacks, and climate scenarios; calibration underway with initial climate scenarios targeted by September 2026
• Two-phase programme:
  • Phase 1 (current): Pilot catchment models + lake model upgrade and calibration
  • Phase 2 (future): Full catchment-scale SWAT+gwflow coupled with lake model for scenario simulation at scale — contingent on Phase 1 outcomes

Current Model Setup — Puarenga (& Pilot Catchments)

Data Inputs Currently Incorporated

Data Available but Not Yet Incorporated

SWAT+ Coupling: MODFLOW vs gwflow

Goal 2: Simulating Multi-Decade Concentration Lags

• Groundwater residence times in the Rotorua catchment can extend decades, particularly through the deep pumice aquifer system
• Delayed nutrient delivery means that current lake water quality reflects land use decisions made generations ago — and that future improvements will similarly lag management changes
• Understanding this lag is central to setting realistic expectations for PC10 outcomes

Approach in SWAT+gwflow

• gwflow simulates spatially distributed groundwater head across the catchment using a cell-by-cell water balance solved at each timestep
• Nutrient concentrations (e.g. NO₃⁻) are transported through the gwflow grid via advective flux — mass is carried with groundwater flow between cells and ultimately discharged to stream reaches
• Lag times emerge naturally from the physics of the system: slow-moving groundwater in low-conductivity or thick aquifer cells retains solutes for longer before discharging to surface water
• This is preferable to conceptual lag approaches (e.g. ROTAN’s exponential delay functions) which impose a lag distribution rather than simulating the underlying flow dynamics

Goal 3: Surface Water–Groundwater Interaction

3a. GW Flow Mapped onto HRUs

• gwflow uses a regular grid overlaid on the catchment, with each cell linked to SWAT+ HRUs for exchange of recharge and groundwater fluxes
• Grid resolution can be tuned independently of HRU delineation, allowing finer GW spatial representation where needed
• Fluxes (recharge, baseflow, SW–GW exchange) are aggregated from HRUs to gwflow cells at each timestep

3b. Layer Configuration

• Current implementation uses a single unconfined layer via gwflow (Dupuit–Forchheimer assumption)
• Appropriate for the Puarenga pilot given its shallow GW system
• Multi-layer representation (e.g. for deeper pumice aquifer systems like Hamurana) would require moving to SWAT+MODFLOW

3c. Residence Time & Model Stability

• Successfully simulated 40 years (1980–2020) in Puarenga catchment without numerical instability
• This demonstrates the model can sustain long-duration runs — a prerequisite for capturing multi-decade GW lag times
• Groundwater residence times in the Rotorua catchment can extend decades; the stable 40-year run is an encouraging proof of concept

Challenges: Dynamic Land Cover in SWAT+

• SWAT+ supports land use change, but implementing it is non-trivial
• The model requires discrete land use snapshots rather than continuous change — transitions must be manually defined and scheduled
• Key challenges include:
  • Ensuring HRU definitions remain consistent across land use scenarios
  • Validating the hydrological and nutrient impacts of land use change at each transition point

Model Configuration

• Land use change in SWAT+ is implemented via management schedules — each transition must be explicitly defined in the input files
• HRU definitions are static at setup; a change in land cover fraction within a sub-basin requires re-mapping HRUs or use of the land use update functionality, which has limited documentation
• Interactions between changing land use and groundwater recharge, evapotranspiration, and nutrient loading need to be carefully validated at each transition point

Workflow Implications

• In a GUI-based setup, each land cover scenario requires substantial manual re-configuration — compounding the bottleneck described previously
• A scripted R-based workflow could make scenario switching significantly more tractable

Bottleneck: GUI-Dependent Workflow in QSWAT+

• QSWAT+ relies on a QGIS-based graphical interface for model setup
• Key steps — watershed delineation, HRU definition, input file generation — are currently performed manually through the GUI
• This creates significant bottlenecks for:
  • Reproducibility and version control
  • Iterative scenario testing and sensitivity analysis
  • Scaling to full catchment or multiple FMUs
  • Collaboration and code review

Moving Towards a Programmable Workflow

• Work is underway to develop R packages that replicate QSWAT+ GUI steps programmatically
• Key functions being implemented include:
  • Watershed and sub-basin delineation
  • HRU generation from land use, soil, and slope inputs
  • SWAT+ input file construction and editing
• This would allow model setup to be fully scripted, version-controlled, and integrated into reproducible workflows (e.g. via renv, Git)

Progress Summary

Key Risks & Issues

• GUI workflow bottleneck — QSWAT+ GUI dependency slows iteration, reproducibility, and scenario testing; R-based automation is in development but not yet production-ready, creating a near-term constraint on model development speed

• Global datasets as placeholders — soils (SWAT+ Global), aquifer thickness (global bedrock), and permeability (GLHYMPS) are coarse proxies; NZ-specific alternatives (LRIS, hydraulic head, aquifer potential maps) exist but require conversion effort before they can improve model realism

• Cross-boundary groundwater flow — particularly in Hamurana, where regional pumice aquifer flow crosses surface catchment divides; incongruous GW–SW boundaries not yet configured and will require hydraulic head data to constrain

• Dynamic land cover complexity — implementing time-varying land use in SWAT+ requires manual scheduling of HRU transitions; compounded by differences in classification schemes across LCDB versions

• Nutrient transport calibration — advective transport of NO₃⁻ through gwflow not yet validated against observed lake inflow concentrations; this is the critical remaining step before lag times can be formally quantified

Risk: gwflow Limitations & Alternative Modelling Platforms

The gwflow implementation makes pragmatic trade-offs that may limit scientific defensibility for some aspects of the Rotorua problem:

• Single unconfined layer — this assumption is a significant simplification for Rotorua’s multi-layer volcanic aquifer geology (interbedded ignimbrites, rhyolite domes, semi-confined sequences, etc.)

• No native solute transport — nitrogen is routed through SWAT+’s HRU-level nutrient routines rather than tracked physically through the aquifer; the decades-long subsurface nitrogen transit central to the PC10 management question cannot be explicitly represented in the current framework

• Surface catchment boundary constraint — gwflow cannot readily extend beyond topographic divides, yet an estimated 45–60 km² of additional groundwater catchment area outside surface boundaries contributes flow to the lake (as identified in ROTAN)

Alternative Platforms

• MODFLOW 6 + GWT + MODPATH — potentially addresses all three limitations; 3D unstructured grid can represent volcanic stratigraphy, GWT handles solute transport, and MODPATH particle tracking produces physically-based MRT distributions directly comparable to tritium data; precedent exists in the same volcanic geology in the Upper Waikato catchment

• However, setup demands are substantially higher and would require a well-constrained 3D geological grid

• INCA-N / INCA-P — lighter-weight alternative; faster to set up than MODFLOW-GWT; adds P transport capability absent from the current framework; better suited to scenario screening and policy communication than gwflow; worth considering if gwflow transport limitations prove critical for Phase 2

GLM-AED Lake Model — Overview

• Lake model upgrade from previous generation to GLM-AED via the AEME (Aquatic Ecosystem Model Ensemble) R package framework
• Pipeline built using targets — a reproducible, dependency-aware workflow manager ensuring all steps are tracked, versioned, and re-run only when inputs change
• Parallel processing enabled via crew for computationally intensive steps (e.g. CMIP6 processing, GCM scenario grid)
• Model configured with the following AED biogeochemical modules active:
  • Sediment flux, oxygen, silica, nitrogen, phosphorus, organic matter, phytoplankton, totals
  • Zooplankton and macrophyte modules defined but currently disabled

Lake Model — Data Inputs & Preparation

Observational Data (Completed)

• Lake level — Bay of Plenty (BOP) bulk export data loaded, QC-filtered (qc_code > 200), and aggregated to daily means from 2014 onward
• CTD profiles — full record extracted from BOP Excel files; depth-resolved temperature and water quality profiles
• Lake chemistry — BOP water quality data 1989–2025 loaded (Lakes Data 1989_2025.xlsx)
• Buoy data — profiler and meteorological data retrieved from LERNZMP API; formatted for AEME input
• NIWA met data — daily and hourly station data (rain, temperature, wind, radiation, evaporation) read and converted to AEME format

Hypsograph

• Hypsograph retrieved from LERNZMP and corrected to align with median observed lake level (masl)
• Sediment zones and parameters estimated from hypsograph for AED sediment flux module

Lake Model — Climate Scenario Setup

CMIP6 Climate Forcing (Completed)

• 7 GCMs configured: ACCESS-CM2, AWI-CM-1-1-MR, CNRM-CM6-1, EC-Earth3, GFDL-ESM4, NZESM, NorESM2-MM
• 4 SSP scenarios: ssp126, ssp245, ssp370, ssp585 (NZESM excluded from ssp585 — data unavailable)
• 5 variables: air temperature, humidity, precipitation, shortwave radiation, wind speed (tas, hurs, pr, rsds, sfcWind)
• CMIP6 data clipped to Rotorua catchment bbox and aligned to VCSN grid points; standardised to Gregorian calendar
• Point-scale GCM data extracted at lake centroid and written to CSV

Simulation Periods Defined

Simulation Grid

• Full crossing of 7 GCMs × 5 scenarios = up to 35 simulation runs per model configuration

Lake Model — Current Status & Next Steps

Completed ✅

• Reproducible targets pipeline built and operational
• All observational data loaded, QC’d, and formatted for AEME
• Hypsograph corrected and sediment zones parameterised
• CMIP6 data processed, standardised, and extracted for all GCMs and scenarios
• GCM time series summaries and plots generated

In Progress 🟡

• Baseline model calibration — GLM-AED model build and run targets are drafted but currently commented out pending resolution of water balance and hypsograph configuration
• Inflow configuration — precipitation-as-inflow approach identified; inflow data integration not yet finalised
• BGC calibration — AED biogeochemical modules configured but not yet calibrated against observed CTD and chemistry data

Pending 🔴

• Scenario simulations across full GCM × SSP grid (35 runs)
• Alum–phosphorus mechanistic module (contingent on UWA collaboration)
• Integration of catchment model nutrient loads as lake inflow boundary conditions (contingent on Phase 1 catchment outcomes)
• Reporting and visualisation targets

Next Steps

Near-term (to March 2026)

• Incorporate NZ-specific soils (LRIS) and hydraulic head data into gwflow calibration for Puarenga
• Configure and test incongruous GW–SW boundary conditions in gwflow for Hamurana
• Calibrate nutrient transport against observed lake inflow concentrations across all three pilot catchments

Mid-term (to September 2026)

• Evaluate pilot results and prepare Phase 2 scaling decision (early 2026)
• Replace ERA5 with VCSN for historical calibration period
• Implement CMIP6 climate scenarios for lake model projections
• Deliver initial climate change scenarios for lake model
• Progress R-based programmable workflow to replace GUI setup steps
• Develop and validate dynamic land cover change workflow using LCDB versions

Longer-term (to February 2027)

• Preliminary advice on Lake Rotoiti / Ōhau Wall downstream impacts (Sep–Oct 2026)
• Finalised lake model outputs including full climate and management scenario set
• Full documentation and reporting