Coupling GoldSim and MODFLOW using a Python Interface

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 educational purposes.

GSPy Interface Function

The main interface function structure (defined using a GoldSim External element):

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 simple 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 some 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 simulation 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 rates at various locations within the grid:

State Persistence

Monthly state persistence is achieved through:

1. State Storage: After each MODFLOW run, the final head distribution is serialized to a pickle file named 'responsive_state.pkl'. This format is used to preserve the Numpy structure of the data.
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 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 useful for educational demonstrations of groundwater flow principles.

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 Downloads

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.

Click here to download the complete code and documentation, including the GSPy interface, MODFLOW setup scripts, and example GoldSim models. Here is a quick start guide. Here is the GoldSim model used to interact with MODFLOW.