This tutorial has the complete process to insert piezometers as a HOB package into a regional groundwater flow in Modflow6. Original location, surface elevation and screen elevation are in csv format that was converted to shapefile in order to be imported by Model Muse 4. The tutorial also covers a Python script to represent the plot of calculated-observed heads with the NRMSE value as a header.

Video

Scripts

This is the Python script to generate the observed - calculated head plot.

#import required packages
import sys, os, re
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.metrics import mean_squared_error
from math import sqrt

#open the modflow6 calculated heads on piezometers
hobSim = pd.read_csv('../Model/Model1.ob_gw_out_head.csv',index_col=0)
hobSim = hobSim.transpose()
hobSim.rename(columns={1.0:'Simulated'}, inplace=True)
hobSim.head()

time	Simulated
OBS_1	3681.363618
OBS_2	4266.586289
OBS_3	4189.450617
OBS_4	3976.768422
OBS_5	4337.856404

#open the piezometric data
hobObs = pd.read_csv('../Txt/Piezometers.csv',index_col=0)
hobObs.head()

	Easting	Northing	surfElev	screenElev	piezometricLevel
Piezometer
OBS_1	618988.927	8355366.926	3794	3616.851	3665.610949
OBS_2	626323.406	8366538.257	4313	4182.432	4236.052783
OBS_3	627046.160	8361952.077	4377	4089.335	4139.527380
OBS_4	622872.221	8360101.290	4144	3901.771	3949.464447
OBS_5	626634.547	8368930.081	4456	4249.271	4308.509850

#create a compiled pandas dataframe and calculate residual
piezoDf = pd.DataFrame()
piezoDf['Observed'] = hobObs['piezometricLevel'] 
piezoDf['Simulated'] = hobSim['Simulated']
piezoDf['Residual'] = piezoDf['Observed'] - piezoDf['Simulated']
piezoDf.head()

	Observed	Simulated	Residual
Piezometer
OBS_1	3665.610949	3681.363618	-15.752669
OBS_2	4236.052783	4266.586289	-30.533506
OBS_3	4139.527380	4189.450617	-49.923237
OBS_4	3949.464447	3976.768422	-27.303975
OBS_5	4308.509850	4337.856404	-29.346554

#get some values about the residual
min=np.min(piezoDf.describe().loc['min'].values[0:2])
max=np.max(piezoDf.describe().loc['max'].values[0:2])
print("Res min: %.2f Res max: %.2f"%(piezoDf['Residual'].min(),piezoDf['Residual'].max()))#

Res min: -49.92 Res max: -3.16

# generate observed - calibrated plot
x = np.linspace(min-5, max+5, 100)
fig = plt.figure(figsize=(10,8))
plt.plot(x, x, linestyle='dashed')
plt.scatter(piezoDf['Observed'],piezoDf['Simulated'], marker='o', c=piezoDf['Residual'])
cbar = plt.colorbar()
cbar.set_label('Residual (m)', fontsize=14)
plt.xlim(min-10,max+10)
plt.ylim(min-10,max+10)
plt.xlabel('Observed Head (m)', fontsize=14)
plt.ylabel('Simulated Head (m)', fontsize=14)
Nrmse = sqrt(mean_squared_error(piezoDf['Observed'].values, piezoDf['Simulated'].values))/\
             (piezoDf['Observed'].max()-piezoDf['Observed'].min())*100
plt.title("NRMSE: "+str(Nrmse))
fig.tight_layout()