Pastas - Groundwater Time Series Analysis

Quick Reference

import pastas as ps
import pandas as pd

# Load data
head = pd.read_csv('well.csv', index_col=0, parse_dates=True).squeeze()
precip = pd.read_csv('precip.csv', index_col=0, parse_dates=True).squeeze()
evap = pd.read_csv('evap.csv', index_col=0, parse_dates=True).squeeze()

# Create model
ml = ps.Model(head, name='Well_001')

# Add recharge stress
sm = ps.RechargeModel(precip, evap, rfunc=ps.Gamma(), name='recharge')
ml.add_stressmodel(sm)

# Solve and plot
ml.solve()
ml.plot()

Key Classes

Essential Operations

Create and Solve Model

ml = ps.Model(head, name='well')
ml.add_stressmodel(ps.RechargeModel(precip, evap, rfunc=ps.Gamma(), name='recharge'))
ml.solve()

Add Pumping Well

pumping = pd.read_csv('pumping.csv', index_col=0, parse_dates=True).squeeze()
ml.add_stressmodel(ps.StressModel(pumping, rfunc=ps.Hantush(),
                                   name='pumping', up=False))  # up=False for drawdown

Model Diagnostics

print(f"EVP: {ml.stats.evp():.1f}%")      # Explained variance
print(f"RMSE: {ml.stats.rmse():.3f} m")   # Root mean square error
print(f"AIC: {ml.stats.aic():.1f}")       # Model selection criterion

ml.plots.diagnostics()                     # Diagnostic plots
ml.plots.acf()                            # Autocorrelation

Get Contributions

contributions = ml.get_contributions()
for name, contrib in contributions.items():
    print(f"{name}: mean={contrib.mean():.2f}")

Step and Impulse Response

step = ml.get_step_response('recharge')    # Step response
block = ml.get_block_response('recharge')  # Impulse response

Export and Load

ml.to_json('model.pas')                    # Save model
ml_loaded = ps.io.load('model.pas')        # Load model

sim = ml.simulate()
sim.to_csv('simulation.csv')               # Export results

Model Statistics

Common Patterns

Compare Response Functions

for rfunc in [ps.Gamma(), ps.Exponential(), ps.Hantush()]:
    ml = ps.Model(head)
    ml.add_stressmodel(ps.RechargeModel(precip, evap, rfunc=rfunc, name='r'))
    ml.solve(report=False)
    print(f"{rfunc.name}: EVP={ml.stats.evp():.1f}%, AIC={ml.stats.aic():.1f}")

Forecast Future Levels

ml.solve()
forecast = ml.simulate(tmin='2024-01-01', tmax='2025-12-31')
ml.plot(tmax='2025-12-31')

River or Custom Stress

river = pd.read_csv('river_stage.csv', index_col=0, parse_dates=True).squeeze()
sm = ps.StressModel(river, rfunc=ps.Exponential(), name='river',
                    settings='waterlevel')
ml.add_stressmodel(sm)

When to Use vs Alternatives

Choose Pastas when: You have groundwater level time series and want to model responses to precipitation, evaporation, or pumping using transfer function noise models. Excellent for rapid model building with diagnostics.

Choose FloPy when: You need spatially distributed groundwater flow modelling with MODFLOW, including multiple layers, boundary conditions, and transport.

Choose custom scipy when: You only need to fit a simple analytical model (e.g., Theis equation) to pumping test data without time series decomposition.

Common Workflows

Groundwater Response Model with Diagnostics

• Load head time series and stress data (precipitation, evaporation, pumping)
• Inspect data: check for gaps, outliers, and time coverage
• Create ps.Model(head) with observation data
• Add recharge stress with ps.RechargeModel(precip, evap, rfunc=ps.Gamma())
• Add pumping or river stresses if applicable
• Solve model with ml.solve()
• Check EVP (>70%), RMSE, and AIC
• Run ml.plots.diagnostics() to inspect residuals
• Check residual autocorrelation; enable noise model if needed: ml.solve(noise=True)
• Compare response functions (Gamma vs Exponential vs Hantush) using AIC
• Extract step/block responses to interpret aquifer behavior
• Decompose signal into individual stress contributions
• Export model to JSON and simulation results to CSV

Tips

1. Start simple - Add stresses incrementally
2. Check residuals - Should be white noise (use ml.plots.diagnostics())
3. Compare response functions - Use AIC/BIC to select best model
4. Use daily data - Pastas works best with daily time series
5. Normalize units - Precipitation in mm/day, head in meters

Common Issues

References

• Stress Models - Available stress model types
• Response Functions - Response function selection

Scripts

• scripts/groundwater_model.py - Complete groundwater modeling workflow