Automatic calibration of groundwater models is a challenging task due to the complexity of the codes, parameters and workflow involved. Calibration of transient state models is even more complex due to changes in sensibility on the different stages of the aquifer requirement such as pumping and recovery.

Python has awesome packages for machine learning that can be coupled to groundwater models and perform automatic calibration of hydraulic parameters. This tutorial covers the procedure to implement a neural network based on a set of parameters set with corresponding head values; the study case is on a transient pumping test model with an observation point located around 10 meters away from the pump well. The analysis predicts the calibration parameters from the spatially interpolated observed head values.

This is actually the third part of a series of tutorials:

Part 1:

Groundwater modeling of a pumping test over a confined aquifer with MODFLOW-6 and FloPy - Tutorial

Part 2:

Sensibility analysis of transient pumping test with MODFLOW-6, Flopy and SALib - Tutorial

Tutorial

Code

# Import required packages
import os
import json
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import StandardScaler
from sklearn.neural_network import MLPRegressor
np.set_printoptions(precision=2)
from workingTools import paramModel

# Open parameter values
logParamValues = np.load('../Txt/logParamValuesBackup.npy')
logParamValues[:5]

array([[-11.19, -12.92, -12.37, -10.01,   0.16],
       [-11.73, -12.92, -12.37, -10.01,   0.16],
       [-11.19, -12.78, -12.37, -10.01,   0.16],
       [-11.19, -12.92,  -6.04, -10.01,   0.16],
       [-11.19, -12.92, -12.37, -10.84,   0.16]])

# Open head values 
#headValues = np.load('../Txt/headArrayBackup.npy')
headValues = np.load('../Txt/headArrayBackup.npy')[:,14:]
headValues[:5]

array([[25.11, 28.47, 31.19, 33.19, 34.65, 35.6 , 36.07, 36.23, 36.26,
        36.26, 36.26, 36.26, 36.26, 36.26, 36.26],
       [25.09, 28.45, 31.17, 33.17, 34.63, 35.59, 36.06, 36.21, 36.24,
        36.24, 36.25, 36.25, 36.25, 36.25, 36.25],
       [25.13, 28.49, 31.21, 33.2 , 34.67, 35.62, 36.09, 36.24, 36.27,
        36.28, 36.28, 36.28, 36.28, 36.28, 36.28],
       [32.09, 33.25, 34.2 , 34.95, 35.59, 36.1 , 36.43, 36.56, 36.59,
        36.6 , 36.6 , 36.6 , 36.6 , 36.6 , 36.6 ],
       [25.11, 30.1 , 32.96, 34.67, 35.66, 36.11, 36.23, 36.26, 36.26,
        36.26, 36.26, 36.26, 36.26, 36.26, 36.26]])

Working with the neural network

# Define the scaler
headScaler = StandardScaler().fit(headValues)

# Scale the train set
headValuesScaled = headScaler.transform(headValues)
headValuesScaled[:5]

array([[-0.58, -0.09,  0.24,  0.45,  0.62,  0.75,  0.82,  0.83,  0.87,
         0.93,  1.02,  1.15,  1.29,  1.42,  1.53],
       [-0.58, -0.09,  0.24,  0.45,  0.61,  0.74,  0.81,  0.82,  0.85,
         0.9 ,  0.98,  1.09,  1.21,  1.33,  1.42],
       [-0.57, -0.08,  0.25,  0.46,  0.63,  0.76,  0.83,  0.85,  0.9 ,
         0.97,  1.08,  1.22,  1.39,  1.54,  1.66],
       [ 1.93,  1.42,  1.16,  1.02,  0.97,  0.97,  1.03,  1.13,  1.41,
         1.7 ,  2.09,  2.58,  3.12,  3.65,  4.11],
       [-0.58,  0.43,  0.78,  0.94,  0.99,  0.97,  0.91,  0.86,  0.88,
         0.93,  1.02,  1.15,  1.29,  1.42,  1.53]])

# Define the scaler
paramScaler = StandardScaler().fit(logParamValues)

# Scale the train set
paramValuesScaled = paramScaler.transform(logParamValues)
paramValuesScaled[:5]

#paramValuesScaled = logParamValues

array([[ 0.24,  0.68, -1.6 , -0.61, -1.61],
       [-0.16,  0.68, -1.6 , -0.61, -1.61],
       [ 0.24,  0.78, -1.6 , -0.61, -1.61],
       [ 0.24,  0.68,  1.6 , -0.61, -1.61],
       [ 0.24,  0.68, -1.6 , -1.23, -1.61]])

# Check scaler over head values
print(headValuesScaled.mean(axis=0))
print(headValuesScaled.std(axis=0))

[ 6.50e-15 -7.70e-15  5.17e-15 -8.31e-15 -3.65e-15  4.01e-15 -1.63e-14
 -5.13e-15 -2.18e-14  1.55e-14 -4.48e-14 -9.14e-14 -6.70e-14 -2.64e-14
  3.50e-14]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]

# Check scaler over param values
print(paramValuesScaled.mean(axis=0))
print(paramValuesScaled.std(axis=0))

[4.00e-16 5.13e-15 3.27e-15 6.41e-15 4.52e-16]
[1. 1. 1. 1. 1.]

# Split over Train and Test
from sklearn.model_selection import train_test_split

headTrain, headTest, paramTrain, paramTest = train_test_split(headValuesScaled, paramValuesScaled, 
                                                              test_size=0.33, random_state=42)

# Regressor
regr = MLPRegressor(hidden_layer_sizes=(100,100,100),max_iter=50000).fit(headTrain, paramTrain)

# Get predicted values for test array
paramPredict = regr.predict(headTest)
paramPredict[:5]

array([[ 1.3 , -1.11,  0.29,  0.24,  0.27],
       [ 0.24, -0.41, -0.68,  0.57,  1.53],
       [ 0.22, -0.09,  0.92, -1.35,  1.03],
       [ 0.27, -0.48, -0.67,  0.58,  1.57],
       [-0.18, -0.69, -0.3 ,  0.47, -1.11]])

paramNames = ['kLay2','kLay3','kLay4','ss','anistrophy']
fig, [ax0,ax1,ax2,ax3,ax4]= plt.subplots(1,5,figsize=(21,4))
ax0.scatter(paramPredict[:,0],paramTest[:,0])
ax0.set_title(paramNames[0])
ax1.scatter(paramPredict[:,1],paramTest[:,1])
ax1.set_title(paramNames[1])
ax2.scatter(paramPredict[:,2],paramTest[:,2])
ax2.set_title(paramNames[2])
ax3.scatter(paramPredict[:,3],paramTest[:,3])
ax3.set_title(paramNames[3])
ax4.scatter(paramPredict[:,4],paramTest[:,4])
ax4.set_title(paramNames[4])
plt.show()

Get heads for model time

from scipy.interpolate import interp1d

#get calculated values on observations points
obsDf = pd.read_csv('../Model/pumpingTest/pumpingTest.obs.head.csv')
obsDf.head()

	time	11J025	11J029	11J030
0	1.000000	36.052406	36.052393	43.289912
1	339.823529	34.055220	10.303132	43.289912
2	1017.470588	31.723297	6.605409	43.289912
3	2372.764706	29.662453	4.251582	43.289912
4	5083.352941	27.897002	2.393536	43.289912

#total elapsed time
totim = obsDf.time.to_list()

#read observed values
piezoDf = pd.read_csv('../Txt/wtElev_11J025_timeDelta')
piezoDf.head()

	TimeSeconds	WT_Elev_m
0	1.0	36.188904
1	6.0	36.188904
2	11.0	36.188904
3	16.0	36.188904
4	21.0	36.152328

interpFunc = interp1d(piezoDf.TimeSeconds.to_list(), piezoDf.WT_Elev_m.to_list(), kind='cubic')
#headModel = interpFunc(totim)
headModel = interpFunc(totim)[14:]
headModel

array([20.71, 26.65, 28.62, 29.93, 31.05, 32.04, 32.93, 33.71, 34.41,
       34.68, 34.88, 35.03, 35.16, 35.27, 35.36])

#plot two head series
fig, ax = plt.subplots(figsize=(10,5))
ax.plot(piezoDf.TimeSeconds,piezoDf.WT_Elev_m,label='Observed')
#ax.plot(totim,headModel,label='Interpolated', marker='*')
ax.plot(totim[14:],headModel,label='Interpolated', marker='*')
ax.legend()
ax.grid()

Predict model parameters

piezoDataScaled = headScaler.transform([headModel])
piezoDataScaled

array([[-2.15, -0.66, -0.54, -0.6 , -0.7 , -0.85, -1.07, -1.43, -2.07,
        -2.68, -3.34, -4.04, -4.63, -5.11, -5.35]])

predictLogParamValues = regr.predict(piezoDataScaled)
predictLogParamValues

array([[ 1.19e+00, -7.22e+00,  9.59e-02,  1.95e-03,  1.88e+00]])

transPredictParamValues = paramScaler.inverse_transform(predictLogParamValues)
transPredictParamValues

array([[ -9.9 , -23.4 ,  -9.02,  -9.21,   4.79]])

predictParamValues = np.e**(transPredictParamValues)[0]
predictParamValues

array([5.04e-05, 6.87e-11, 1.21e-04, 1.00e-04, 1.20e+02])

np.save('../Txt/predictParamValues.npy',predictParamValues)

def runModelForParams(kLay2Org,kLay2,kLay3Org,kLay3,kLay4Org,kLay4,ss,anisotrophy):
    simName = "mfsim.nam"
    simWs = "../Model/pumpingTest/"
    exeName = "../Bin/mf6.exe"
    destDir = "../workingModel"
    parSim = paramModel(simName=simName,
                          simWs=simWs,
                          exeName=exeName,
                          destDir=destDir)
    parSim.loadSimulation()
    parSim.changeSimPath()
    gwf = parSim.loadModel('pumpingTest')

    npf = gwf.get_package('NPF')
    kArray = np.copy(npf.k.array)
    kArray[kArray == kLay2Org] = kLay2
    kArray[kArray == kLay3Org] = kLay3
    kArray[kArray == kLay4Org] = kLay4
    # npf.k = kArray
    npf.k33 = kArray/anisotrophy

    sto = gwf.get_package('STO')
    ssArray = np.copy(sto.ss.array)
    ssArray = np.ones(ssArray.shape)*ss
    # syArray = np.copy(sto.sy.array)
    # syArray = np.ones(syArray.shape)*sy
    sto.ss = ssArray
    # sto.sy = syArray

    parSim.sim.write_simulation(silent=True)
    parSim.sim.run_simulation(silent=True)

    obs = gwf.get_package('OBS')
    obsFile = list(obs.continuous.data.keys())[0]
    piezoDf = pd.read_csv(os.path.join(destDir,obsFile))
    piezoArray = piezoDf['11J025'].values
    return piezoArray

piezoData = runModelForParams(kLay2Org = 1.e-5,
                              kLay2 = predictParamValues[0],
                              kLay3Org = 1.e-6,
                              kLay3 = predictParamValues[1],
                              kLay4Org = 3.e-5,
                              kLay4 = predictParamValues[2],
                              ss = predictParamValues[3],
                              anisotrophy = predictParamValues[4])
piezoData

Check for files in 'C:\Users\saulm\Documents\webinarSensitivityAnalisysPumpingTestFlopySALib_v1\workingModel'.
Deleted file: C:\Users\saulm\Documents\webinarSensitivityAnalisysPumpingTestFlopySALib_v1\workingModel\mfsim.lst


All files in 'C:\Users\saulm\Documents\webinarSensitivityAnalisysPumpingTestFlopySALib_v1\workingModel' have been deleted.

loading simulation...
  loading simulation name file...
  loading tdis package...
  loading model gwf6...
    loading package dis...
    loading package npf...
    loading package ic...
    loading package sto...
    loading package ghb...
    loading package wel...
    loading package oc...
    loading package obs...
  loading solution package pumpingtest...

 Models in simulation: ['pumpingtest']
Package list: ['DIS', 'NPF', 'IC', 'STO', 'GHB_0', 'WEL_0', 'OC', 'OBS_0']





array([36.03, 34.02, 31.65, 29.49, 27.55, 25.81, 24.44, 23.64, 23.34,
       23.3 , 23.28, 23.27, 23.27, 23.27, 23.27, 25.28, 27.65, 29.81,
       31.75, 33.49, 34.86, 35.66, 35.96, 36.01, 36.03, 36.03, 36.03,
       36.03, 36.03])

#plot two head series
fig, ax = plt.subplots(figsize=(10,5))
ax.plot(piezoDf.TimeSeconds,piezoDf.WT_Elev_m,label='Observed')
ax.plot(totim,piezoData,label='Predicted')
ax.legend()
ax.grid()