Dynamic Coupling of GoldSim and MODFLOW with FloPy

This post describes a method for dynamically coupling GoldSim with MODFLOW using GSPy (GoldSim's Python interface) and FloPy (a Python package for creating MODFLOW models). This approach is ideal for educational purposes and provides a clear, in-memory alternative to traditional file-based model coupling.

Implementation Approach

Traditional model coupling typically requires extensive setup, including:

• 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 simplifies this process by using:

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

This implementation uses a simple, confined aquifer MODFLOW setup for educational purposes.

GSPy Interface Function

The main interface function structure is defined within a GoldSim External Element. This element sends current simulation parameters to Python and receives the calculated groundwater heads back.

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

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()
    # Returns head at monitor locations (layer=0, row, col)
    return (heads[0, 4, 6], heads[0, 10, 10], heads[0, 13, 13])

The complete Python code handles MODFLOW model creation, state persistence, error handling, and scenario reset detection.

MODFLOW Model Configuration

The model uses a simple, transient groundwater flow setup:

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 (W) at cell (3,4) in the northwest quadrant, which can be adjusted using GoldSim Data elements.
• Three monitoring points (M) 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 is the model grid layout, illustrating the locations:

      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

Implementation Features

The combined GSPy + FloPy approach offers sophisticated features for transient modeling:

• Monthly State Persistence: Each MODFLOW run uses groundwater conditions from the previous month as initial conditions, maintaining continuity throughout the simulation using a local pickle file.
• Automatic Simulation Detection: The system detects when GoldSim starts a new scenario run (based on the time gap between function calls) and automatically resets the groundwater conditions to the initial state (uniform 40m heads).
• Distance-Based Response: The three monitoring points clearly demonstrate how pumping effects dissipate with distance.
• Monthly Timing: The MODFLOW simulation is run for a 30-day period during each GoldSim time step, accounting for the physical time lag between pumping and groundwater response.

Model Operation and Results

Coupling Mechanics

The model operates on a multi-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 head distribution (or initial conditions) from a pickle file.
3. MODFLOW Setup: FloPy creates a new MODFLOW model instance, applying 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.
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 back to the pickle file for use in the next month's simulation.

Cone of Depression Formation

This coupling effectively simulates the characteristic cone of depression that forms around the pumping well:

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, demonstrating the radial drawdown pattern typical of confined aquifer pumping.

Monitoring Point Response

The time series demonstrates the cumulative effects of monthly pumping rates at various locations within the grid:

The three monitoring points show the expected distance-dependent responses:

• Northwest Monitor (4,6): Shows the strongest response (closest to the well).
• Center Monitor (10,10): Shows moderate response.
• Southeast Monitor (13,13): Shows the weakest response (furthest from the well).

Contact:

Jason Lillywhite (Lillywhite Water Solutions)

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.

Download the complete example package from here.

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.

Also note that the GSPy dll for this example is called gspy_modflow.json and is located in the gs_modflow folder that is zipped in the below attachment.