Cómo insertar Piezómetros (HOB) a un Modelo MODFLOW6 y mostrar los Niveles Observados - Calculados

Este tutorial muestra el proceso completo para insertar piezómetros como paquete HOB en un modelo de flujo regional de aguas subterráneas en Modflow6. La ubicación original, la elevación de la superficie y la elevación de la rejilla del piezometro están en formato csv que se convierte a shapefile para ser importado por Model Muse 4. El tutorial también contiene un script de Python para representar el gráfico de cargas hidráulicas calculadas - observadas con el valor NRMSE como encabezado.

Video

Código

Este es el código en Python para generar el gráfico de cargas hidráulicas observadas - calculadas

#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()