# Coupling GoldSim and MODFLOW Using GSPy

This post describes a method for coupling GoldSim with MODFLOW using GSPy (GoldSim's Python interface) and FloPy (a Python package for creating MODFLOW models). The approach builds on existing work like the Dynamic Coupling of GoldSim and MODFLOW/MT3D project but focuses on a simplified implementation suitable for educational purposes.

## Implementation Approach

Traditional model coupling typically requires:

- Writing temporary input files
- Managing external executable calls
- Parsing output files for results
- Handling timing and synchronization
- Managing file system operations

The GSPy + FloPy approach uses:

1. GSPy to handle communication between GoldSim and Python
2. FloPy to create and manage MODFLOW models in memory
3. Python functions that receive GoldSim inputs and return results
4. FloPy's built-in file management for MODFLOW input/output

This implementation uses a Tutorial 2 style MODFLOW setup for clarity and educational value.

## GSPy Interface Function

The core interface function structure:

    def run_modflow_step(*args):
        # Parse GoldSim inputs
        well_rate = float(args[0])
        recharge_rate = float(args[1])
        well_row, well_col = int(args[2]), int(args[3])
        hydraulic_k = float(args[4])

        # Create and run MODFLOW model
        mf = create_modflow_model(well_rate, recharge_rate, well_row, well_col, hydraulic_k)
        mf.write_input()
        success, _ = mf.run_model(silent=True)

        # Extract and return results to GoldSim
        heads = flopy.utils.HeadFile('responsive.hds').get_data()
        return (heads[0, 4, 6], heads[0, 10, 10], heads[0, 13, 13])

The complete implementation is approximately 225 lines of Python code that handles MODFLOW model creation, state persistence, error handling, and simulation detection. FloPy eliminates the need for manual MODFLOW file creation and parsing, while GSPy handles the GoldSim-Python communication that would typically require additional coupling infrastructure.

## MODFLOW Model Configuration

The implementation uses a Tutorial 2 style groundwater model with the following specifications:

### Grid Setup:

- 21×21 cells (441 total cells)
- 100m × 100m cell size (2.1km × 2.1km domain)
- Single confined aquifer layer (50m thick)
- Constant head boundaries (40m) representing surrounding water bodies

### Well and Monitoring Configuration:

- Pumping well at cell (3,4) in northwest quadrant
- Three monitoring points at different distances:
  - Northwest monitor (4,6): 280m from well - maximum response
  - Center monitor (10,10): 850m from well - moderate response
  - Southeast monitor (13,13): 1270m from well - minimum response

Here's the model grid layout:

         0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
     0   C  C  C  C  C  C  C  C  C  C  C  C  C  C  C  C  C  C  C  C  C
     1   C  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  C
     2   C  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  C
     3   C  .  .  .  W  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  C  ← W = Well (3,4)
     4   C  .  .  .  .  .  M  .  .  .  .  .  .  .  .  .  .  .  .  .  C  ← M = NW monitor (4,6)
     5   C  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  C
     6   C  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  C
     7   C  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  C
     8   C  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  C
     9   C  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  C
    10   C  .  .  .  .  .  .  .  .  .  M  .  .  .  .  .  .  .  .  .  C  ← M = Center monitor (10,10)
    11   C  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  C
    12   C  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  C
    13   C  .  .  .  .  .  .  .  .  .  .  .  .  M  .  .  .  .  .  .  C  ← M = SE monitor (13,13)
    14   C  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  C
    15   C  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  C
    16   C  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  C
    17   C  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  C
    18   C  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  C
    19   C  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  C
    20   C  C  C  C  C  C  C  C  C  C  C  C  C  C  C  C  C  C  C  C  C

    Legend: C = Constant head boundary (40m), . = Active aquifer cells, W = Pumping well, M = Monitoring points

## GSPy Interface Configuration

GSPy requires minimal interface configuration. The setup process involves:

1. Write your Python function with the parameters you need
2. Set up GoldSim's External Element to call your function
3. Connect your GoldSim elements to the inputs and outputs

The combined GSPy + FloPy approach handles:

- GSPy: Data type conversion between GoldSim and Python
- GSPy: Error handling and fallback values
- GSPy: Process communication between GoldSim and Python
- FloPy: MODFLOW input file creation and output file parsing
- FloPy: Model object management in memory

The function signature becomes your interface specification:

- 6 inputs: pumping_rate, recharge_rate, well_row, well_col, hydraulic_k, specific_yield
- 3 outputs: northwest_head, center_head, southeast_head

This eliminates the need for XML configuration files, manual MODFLOW file management, or complex setup procedures.

## Implementation Features

Monthly State Persistence: Each MODFLOW run uses groundwater conditions from the previous month as initial conditions, maintaining continuity within multi-year simulations.

Automatic Simulation Detection: The system detects when GoldSim starts a new scenario run (based on time gaps between calls) and resets groundwater conditions to initial state.

Distance-Based Response: Three monitoring points at different distances from the well demonstrate how pumping effects decrease with distance.

Monthly Timing: The implementation accounts for the time lag between pumping decisions and groundwater response.

## Model Operation and Results

### Coupling Mechanics

The model operates on a 5-year simulation period with monthly time steps. The coupling process works as follows:

1. Monthly Execution: At the beginning of each GoldSim month, the External Element calls the Python function with current pumping parameters.
2. State Loading: The Python function loads the previous month's groundwater head distribution from a pickle file (or uses initial conditions for the first month).
3. MODFLOW Setup: FloPy creates a new MODFLOW model using the loaded heads as initial conditions and the current month's pumping rate.
4. 30-Day Simulation: MODFLOW runs a transient 30-day simulation with daily time steps, representing the physical response over the month.
5. Result Extraction: The final heads are extracted from the three monitoring locations and returned to GoldSim.
6. State Persistence: The final head distribution is saved to the pickle file for use in the next month's simulation.

### Simulation Detection and Reset

The system automatically detects when GoldSim starts a new scenario run by monitoring the time gap between function calls. If more than 10 seconds elapse between calls, the system assumes a new simulation has started and resets to initial conditions (uniform 40m heads). This allows multiple scenario runs without manual intervention.

### Cone of Depression Formation

The following image shows the characteristic cone of depression that forms around the pumping well:

The cone of depression demonstrates the radial drawdown pattern typical of confined aquifer pumping. The constant head boundaries (40m) on all edges provide an infinite source of water, allowing steady-state conditions to develop within each 30-day simulation period.

### Monitoring Point Response

The three monitoring points show distance-dependent responses to pumping:

- Northwest Monitor (4,6): Located approximately 280m from the well, shows the strongest response to pumping changes.
- Center Monitor (10,10): Located approximately 850m from the well, shows moderate response.
- Southeast Monitor (13,13): Located approximately 1270m from the well, shows the weakest response.

The following time series demonstrates the cumulative effects of monthly pumping decisions:

### State Persistence Mechanism

Monthly state persistence is achieved through a simple but effective mechanism:

1. State Storage: After each MODFLOW run, the final head distribution is serialized to a pickle file named 'responsive_state.pkl'.
2. State Loading: At the start of each new month, the function checks for the existence of this file and loads the previous heads as initial conditions.
3. Boundary Enforcement: Constant head boundary conditions (40m) are re-applied to ensure boundary consistency.
4. Fallback Handling: If no state file exists (first run or after reset), uniform 40m initial conditions are used.

This approach eliminates the need for complex database management or file synchronization while maintaining the physical continuity required for realistic multi-year groundwater simulations.

### Quantitative Results

The model produces predictable, linear responses to different pumping rates:

- 500 m³/day: 0.64m drawdown at well, 0.07m at NW monitor, 0.15m at SE monitor
- 1000 m³/day: 1.27m drawdown at well, 0.13m at NW monitor, 0.30m at SE monitor
- 2000 m³/day: 2.54m drawdown at well, 0.26m at NW monitor, 0.59m at SE monitor

The linear relationship between pumping rate and drawdown is characteristic of the confined aquifer setup with constant head boundaries, making this model ideal for educational demonstrations of groundwater flow principles.

### Potential Applications

- Educational demonstrations of groundwater concepts
- Training materials for water professionals
- Student research projects
- Proof-of-concept studies
- Learning coupled groundwater-surface water modeling

## Comparison with Traditional Approaches

The complete implementation requires 225 lines of Python code. The GSPy + FloPy combination eliminates the following components typically required:

### Traditional File-Based Approach Would Add:

- 50+ lines for writing MODFLOW input files each month
- 30+ lines for managing external process execution
- 40+ lines for parsing output files and error handling
- 25+ lines for file cleanup and temporary directory management
- 35+ lines for handling file locking, timing, and synchronization
- 20+ lines for GoldSim-specific file I/O protocols

### GSPy + FloPy Approach:

- GSPy: Direct function calls eliminate GoldSim I/O coordination
- FloPy: Handles MODFLOW file creation and parsing automatically
- GSPy: Results return directly as Python objects
- GSPy: Built-in error handling and fallback mechanisms
- FloPy: No manual MODFLOW file management needed

A traditional file-based approach would likely require 400+ lines with additional file management complexity. The GSPy + FloPy approach allows the 225 lines to focus on groundwater modeling logic rather than coupling infrastructure.

State persistence between months uses a simple pickle file operation rather than complex database management or file synchronization.

## Implementation Notes

Key observations from this implementation:

Development Time: The 225 lines of Python code required approximately one day to write and debug. A traditional file-based approach would likely require 3-4 days due to I/O coordination complexity.

Debugging: Issues are typically related to groundwater modeling logic rather than file I/O problems.

Code Structure: Approximately 180 lines focus on MODFLOW model setup and groundwater physics, with 45 lines handling coupling mechanics.

Maintenance: Updates are straightforward since all logic is contained in Python code rather than distributed across multiple file formats.

Error Handling: GSPy provides built-in error handling for GoldSim communication, while FloPy handles MODFLOW file operations, reducing the need for complex file system error management.

## Potential Extensions

Possible enhancements to this implementation:

- Multiple wells with different pumping schedules
- Variable aquifer properties for heterogeneous conditions
- Surface water-groundwater interactions
- Optimization routines for pumping strategies

## Dependencies and Related Work

This implementation relies on:

- FloPy: A Python package for creating, running, and post-processing MODFLOW models, developed and maintained by the USGS. FloPy handles all MODFLOW file operations and model management.
- GSPy: GoldSim's Python interface for communication between GoldSim and Python functions.

For comprehensive applications, refer to GoldSim's Dynamic Coupling framework.

The complete code and documentation are available on GitHub, including the GSPy interface, MODFLOW setup scripts, and example GoldSim models.