Cómo crear un Modelo de Aguas Subterráneas totalmente geoespacial con MODFLOW y Flopy - Tutorial

infahatari.jpg

La naturaleza y todos los procesos físicos relacionados con el agua subterránea y el régimen de transporte están espacialmente distribuidos. Los modelos de agua subterránea se basan en una estructura grillada que se discretizan en celdas con arreglos de filas y columnas; el nivel de desconexión espacial de una parte del medio poroso con las filas y columnas de las celdas son algunos desafíos para la gestión sostenible de los recursos de aguas subterráneas.

Tenemos que crear o recrear la dualidad entre la grilla geoespacial y la grilla del modelo, que sería similar a la dualidad de un objeto vectorial GIS y sus metadatos, pero con una cierta complejidad. Los usuarios principiantes no enfrentarán tales dificultades cuando trabajen con un software de pre-procesamiento y post-procesamiento como Model Muse, sin embargo, para investigadores y profesionales relacionados con conocimientos intermedios en modelamiento de aguas subterráneas como el análisis de sensibilidad de la discretización, automatización de la calibración y el acoplamiento de modelos de aguas subterráneas con ecosistemas, sería mucho más eficiente trabajar los modelos de aguas subterráneas referenciados geoespacialmente.

Afortunadamente, Flopy y las librerías de Python permiten construir y simular modelos de MODFLOW, y tienen herramientas para georreferenciar las grillas del modelo incluso con opciones de rotación; sin embargo, el proceso es algo explícito que significa que el modelador tenga conocimientos de Flopy y Python. Este tutorial muestra todo el procedimiento para crear un modelo de agua subterránea totalmente geoespacial con MODFLOW y Flopy.

Problemas sobre un “modelo totalmente geoespacial”

A pesar de que hemos realizado algunos tutoriales que representan los resultados de MODFLOW en QGIS, los ejemplos no fueron completamente geoespaciales. Para crear un modelo numérico “totalmente geoespacial", tenemos que crear la extensión del modelo desde shapefiles, la discretización del modelo desde shapefiles, importar la elevación de la superficie desde rásters y exportar los archivos de salida del modelo a plataformas GIS como QGIS. Se necesita un acoplamiento completo de la naturaleza espacial del régimen de flujo de agua subterránea en nuestro modelo.


Tutorial

Código de Python

Este es el código de Python usado:

Import the required libraries

import os
import numpy as np
import matplotlib.pyplot as plt
import flopy
import shapefile as sf #in case you dont have it, form anaconda prompt: pip install pyshp
from flopy.utils.gridgen import Gridgen
from flopy.utils.reference import SpatialReference
import pandas as pd
from scipy.interpolate import griddata
flopy is installed in E:\Software\Anaconda37\lib\site-packages\flopy

Creation of model object and application of MODFLOW NWT

modelname='Model1'
model_ws= '../Model'
mf = flopy.modflow.Modflow(modelname, exe_name="../Exe/MODFLOW-NWT_64.exe", version="mfnwt",model_ws=model_ws)

#Definition of MODFLOW NWT Solver
nwt = flopy.modflow.ModflowNwt(mf ,maxiterout=200,linmeth=2,headtol=0.01)

Refinement areas and spatial / temporal discretization

# Open the shapefiles from the model limit and the refiment area around the wells
modelLimitShp = sf.Reader('../Shp/ModelLimit1.shp')
modelWelShp = sf.Reader('../Shp/ModelWell2.shp')
#We create a preview figure of the model limits and refinement area 

#Numpy array of model limit and refinement area (well area)
limitArray = np.array(modelLimitShp.shapeRecords()[0].shape.points)
wellArray = np.array([point.shape.points[0] for point in modelWelShp.shapeRecords()]) # looks complicated but it is just a list comprehension converted to numpy array
WBB = modelWelShp.bbox                                                                # WBB comes from Well Bounding Box

fig = plt.figure()
plt.plot(limitArray[:,0],limitArray[:,1])
plt.plot([WBB[0],WBB[0],WBB[2],WBB[2],WBB[0]],[WBB[1],WBB[3],WBB[3],WBB[1],WBB[1]])
plt.scatter(wellArray[:,0],wellArray[:,1])
plt.show()
output_6_0.png
#Since we have geospatial and georeferenced data, we can add a satellite image in the background with the Mplleaflet package

import mplleaflet

crs = {'init' :'epsg:32718'}
mplleaflet.display(fig,crs=crs)
E:\Software\Anaconda37\lib\site-packages\IPython\core\display.py:689: UserWarning: Consider using IPython.display.IFrame instead
  warnings.warn("Consider using IPython.display.IFrame instead")
# Coordinates of the global and local discretization
GloRefBox = modelLimitShp.bbox 
LocRefBox = modelWelShp.bbox 

print(GloRefBox)
print(LocRefBox)
[353300.0, 8547900.0, 356400.0, 8549600.0]
[354369.8840337161, 8548609.34745342, 354904.6274241918, 8549128.04875159]
#Calculating Global Model (Glo) and Local Refinement (Loc) dimensions
GloLx = GloRefBox[2] - GloRefBox[0] #x_max - x_min
GloLy = GloRefBox[3] - GloRefBox[1]
print('Global Refinement Dimension. Easting Dimension: %8.1f, Northing Dimension: %8.1f' % (GloLx,GloLy))

LocLx = LocRefBox[2] - LocRefBox[0] #x_max - x_min
LocLy = LocRefBox[3] - LocRefBox[1]
print('Local Refinement Dimension. Easting Dimension: %8.1f, Northing Dimension: %8.1f' % (LocLx,LocLy))
Global Refinement Dimension. Easting Dimension:   3100.0, Northing Dimension:   1700.0
Local Refinement Dimension. Easting Dimension:    534.7, Northing Dimension:    518.7
#Defining Global and Local Refinements, for purpose of simplicity cell x and y dimension will be the same
celGlo = 100
celRef = 50
def arrayGeneratorCol(gloRef, locRef, gloSize, locSize):

    cellArray = np.array([])

    while cellArray.sum() + gloRef[0] < locRef[0] - celGlo:
        cellArray = np.append(cellArray,[gloSize])
    while cellArray.sum() + gloRef[0] > locRef[0] - celGlo and cellArray.sum() + gloRef[0] < locRef[2] + celGlo:
        cellArray = np.append(cellArray,[locSize])
    while cellArray.sum() + gloRef[0] > locRef[2] + celGlo and cellArray.sum() + gloRef[0] < gloRef[2]:
        cellArray = np.append(cellArray,[gloSize])

    return cellArray
def arrayGeneratorRow(gloRef, locRef, gloSize, locSize):

    cellArray = np.array([])
    accumCoordinate =  gloRef[3] - cellArray.sum()

    while gloRef[3] - cellArray.sum() > locRef[3] + celGlo:
        cellArray = np.append(cellArray,[gloSize])
    while gloRef[3] - cellArray.sum() < locRef[3] + celGlo and gloRef[3] - cellArray.sum() > locRef[1] - celGlo:
        cellArray = np.append(cellArray,[locSize])
    while gloRef[3] - cellArray.sum() < locRef[1] - celGlo and gloRef[3] - cellArray.sum() > gloRef[1]:
        cellArray = np.append(cellArray,[gloSize])

    return cellArray
#Remember that DELR is space between rows, so it is the column dimension array
delRArray = arrayGeneratorCol(GloRefBox, LocRefBox, celGlo, celRef)
delRArray
array([100., 100., 100., 100., 100., 100., 100., 100., 100., 100.,  50.,
        50.,  50.,  50.,  50.,  50.,  50.,  50.,  50.,  50.,  50.,  50.,
        50.,  50.,  50., 100., 100., 100., 100., 100., 100., 100., 100.,
       100., 100., 100., 100., 100., 100.])
#And DELC is the space between columns, so its the row dimension array
delCArray = arrayGeneratorRow(GloRefBox, LocRefBox, celGlo, celRef)
delCArray
array([100., 100., 100., 100.,  50.,  50.,  50.,  50.,  50.,  50.,  50.,
        50.,  50.,  50.,  50.,  50.,  50.,  50., 100., 100., 100., 100.,
       100., 100.])
#Calculating number or rows and cols since they are dependant from the discretization
nrows = delCArray.shape[0]
ncols = delRArray.shape[0]
print('Number of rows: %d and number of cols: %d' % (nrows,ncols))
Number of rows: 24 and number of cols: 39
#Define some parametros and values for the spatial and temporal discretization (DIS package)

#Number of layers and layer elevations
nlays = 3
mtop = 0
botm = [-10,-20,-30]

#Number of stress periods and time steps
#In this case we will simulate the effect on the aquifer of 2000 days divided in 10 stress periods of 200 days,
#each stress period is divided in 4 time steps
#there is a steady state stress period at the beginning of the simulation
nper = 11
perlen = np.ones(nper)
perlen[0] = 1
perlen[1:] = 200 * 86400

#Definition of time steps
nstp = np.ones(11) 
nstp[0] = 1
nstp[1:] = 4

#Definition of stress period type, transient or steady state
periodType = np.zeros(11, dtype=bool)
periodType[0] = True
print('This is the lenght of stress periods: ', perlen)
print('This is the number of time steps: ', nstp)
print('This is the stress period type: ', periodType)
This is the lenght of stress periods:  [1.000e+00 1.728e+07 1.728e+07 1.728e+07 1.728e+07 1.728e+07 1.728e+07
 1.728e+07 1.728e+07 1.728e+07 1.728e+07]
This is the number of time steps:  [1. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.]
This is the stress period type:  [ True False False False False False False False False False False]
# Apply the spatial and temporal discretization parameters to the DIS package
dis = flopy.modflow.ModflowDis(mf, nlay=nlays,
                               nrow=delCArray.shape[0],ncol=delRArray.shape[0],
                               delr=delRArray,delc=delCArray,
                               top=mtop,botm=botm,
                               itmuni=1.,nper=nper,perlen=perlen,nstp=nstp,steady=periodType)
#Assign spatial reference
mf.dis.sr = SpatialReference(delr=delRArray,delc=delCArray, xul=GloRefBox[0], yul= GloRefBox[3],epsg=32718)
mf.dis.sr.epsg
32718
mf.modelgrid.set_coord_info(xoff=GloRefBox[0], yoff=GloRefBox[3]-delCArray.sum(), angrot=0,epsg=32717)
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(1, 1, 1, aspect='equal')
modelmap = flopy.plot.PlotMapView(model=mf)
linecollection = modelmap.plot_grid(linewidth=0.5, color='royalblue')
output_22_0.png
#Representation of model geometry with all the boundary conditions
fig = plt.figure(figsize=(8,8))

modelmap = flopy.plot.PlotMapView(model=mf)
linecollection = modelmap.plot_grid(linewidth=0.5, color='cyan')
shpRiver = flopy.plot.plot_shapefile('../Shp/ModelRiver2', facecolor='none') #RIV boundary condition
shpGHB = flopy.plot.plot_shapefile('../Shp/ModelGHB1', facecolor='none')     #GHB boundary condition

plt.plot(limitArray[:,0],limitArray[:,1])
plt.scatter(wellArray[:,0],wellArray[:,1])

crs = {'init' :'epsg:32718'}
mplleaflet.display(fig,crs=crs)

Definition of Model Top and Layer Bottom Elevation for the UPW package

dem = pd.read_csv('../Rst/ModeloDem2.csv')
dem.head()
Easting Northing Elevation
0 353214.799 8549716.14 58.925
1 353245.085 8549716.14 59.267
2 353275.372 8549716.14 59.609
3 353305.658 8549716.14 59.951
4 353335.945 8549716.14 60.293
points = dem[['Easting','Northing']].values
values = dem['Elevation'].values
grid_x = mf.modelgrid.xcellcenters
grid_y = mf.modelgrid.ycellcenters
mtop = griddata(points, values, (grid_x, grid_y), method='nearest')
mtop[:5,:5]
array([[60.293, 61.661, 62.686, 63.712, 65.08 ],
       [60.293, 61.661, 62.686, 63.712, 65.08 ],
       [60.293, 61.661, 62.686, 63.712, 65.08 ],
       [60.293, 61.661, 62.686, 63.712, 65.08 ],
       [60.293, 61.661, 62.686, 63.712, 65.08 ]])
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect('equal')
modelmap = flopy.plot.PlotMapView(model=mf,ax=ax)
quadmesh = modelmap.plot_array(mtop)
output_27_0.png
#Asign surface elevations
mf.dis.top  = mtop

AcuifInf_Bottom = -120

AcuifMed_Bottom = AcuifInf_Bottom + (0.5 * (mtop - AcuifInf_Bottom))

AcuifSup_Bottom = AcuifInf_Bottom + (0.75 * (mtop - AcuifInf_Bottom))

zbot = np.zeros((nlays,nrows,ncols))

zbot[0,:,:] = AcuifSup_Bottom
zbot[1,:,:] = AcuifMed_Bottom
zbot[2,:,:] = AcuifInf_Bottom

#Asign layer bottom elevations
mf.dis.botm = zbot
fig = plt.figure(figsize=(20, 3))
ax = fig.add_subplot(1, 1, 1)
modelxsect = flopy.plot.PlotCrossSection(model=mf, line={'Row': 20})
linecollection = modelxsect.plot_grid()
output_29_0.png
Kx = np.zeros((nlays,nrows,ncols))
Kx[0,:,:] = 1E-5 #first layer of lime
Kx[1,:,:] = 5E-4 #second layer of sand
Kx[2,:,:] = 2E-4 #third layer of sandy gravel
fig = plt.figure(figsize=(20, 3))
ax = fig.add_subplot(1, 1, 1)
modelxsect = flopy.plot.PlotCrossSection(model=mf, line={'Row': 20})
linecollection = modelxsect.plot_grid()
modelxsect.plot_array(Kx)
print(Kx[:,10,10]) #print values of K for all layers in cell row / col = 10 , 10
[1.e-05 5.e-04 2.e-04]
output_31_1.png
# Definition of layer type, the first two layers are convertible
laytyp = [1,1,0]
#lpf = flopy.modflow.ModflowLpf(mf, laytyp = laytyp, hk = Kx)
upw = flopy.modflow.ModflowUpw(mf, laytyp = laytyp, hk = Kx,
                              ss=1e-05, sy=0.15)

Configuration of active zones and initial heads: BAS package

ibound  = np.ones([nlays,nrows,ncols])
bas = flopy.modflow.ModflowBas(mf, ibound=ibound,strt=mtop)

Definition of the GRIDGEN object

For the manipulation of spatial data to determine hydraulic parameters or boundary conditions

g = Gridgen(mf.dis, model_ws=model_ws,exe_name='../Exe/gridgen_x64.exe',surface_interpolation='replicate')
g.build(verbose=False)

Recharge and Evapotranspiration as RCH and EVT boundary condition

# We apply recharge to the extension of the irrigated areas and evapotranspiration to the whole extent,
# The recharge rate is estimated in 200 mm /yr and the potential evapotraspiration is 1200 mm/yr

# Evapotranspiration rate
evtr = 1.2 / 365. / 86400. # 1200 mm/yr in m/s
evt = flopy.modflow.ModflowEvt(mf,evtr=evtr, surf=mtop, exdp=0.5, nevtop=1)
# Recharge rate
rechRate = 0.2 / 365. / 86400. # 200 mm/yr in m/s
rech_intersect = g.intersect('../Shp/ModelRechargeZone1','polygon',0)
rech_spd = {}                                      # empty dictionary for all stress periods, actually spd comes from stress period data
rech_spd[0] = np.zeros([nrows,ncols])              # empty list as the first item of dictionary for the first stress period
rech_unique = np.unique(rech_intersect.nodenumber) # unique cells where recharge is applied, the np.unique also sorts the cell number
print(rech_unique[:5])  # a sample of the sorted list
[0 1 2 3 4]
for i in np.arange(rech_unique.shape[0]):
    x,y = g.get_center(rech_unique[i]) #get the cell centroid x and y
    i,j = mf.sr.get_ij(x,y)            #get the cell row and column
    rech_spd[0][i,j] =rechRate         #add the recharge in the cell to the array for the first stress period

# Recharge is applied to the highest active cell nrchop = 3 (default)
rec = flopy.modflow.ModflowRch(mf, nrchop=3, rech=rech_spd)
# Representation of applied recharge values for the last stress period. 
# For flopy, and modflow if no recharge for other stress period is specified, the first recharge will remain for the entire simulation. 
fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(1, 1, 1, aspect='equal')
modelmap = flopy.plot.PlotMapView(model=mf)
quadmesh = modelmap.plot_array(mf.rch.rech.array[-1,0,:,:], alpha=0.3)
output_41_0.png
wel_intersect = g.intersect('../Shp/ModelWell2','point',0)
wel_spd = {}
wel_spd[0] = [0,0,0,0] #wells are not pumping on steady state (stress period 1) but we have to insert zero pumping rate to at least one cells
wel_spd[1] = []
wel_unique = np.unique(wel_intersect.nodenumber) #For the well case this is kind of redundant because wells are located apart form each other
#It could be good if the shapefile have overlying points
for i in np.arange(wel_unique.shape[0]):
    x,y = g.get_center(wel_unique[i])
    i,j = mf.sr.get_ij(x,y)
    wel_spd[1].append([1,i,j,-0.03]) #well are located on the second layer and each of them pump 10 l/s (-0.001 m3/s)
wel = flopy.modflow.ModflowWel(mf,stress_period_data=wel_spd)
# Plot the well shapefile and the cells conceptualized as wells
fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(1, 1, 1, aspect='equal')
modelmap = flopy.plot.PlotMapView(model=mf)
g.plot(ax, linewidth=0.5)
quadmesh = modelmap.plot_array(mf.wel.stress_period_data.array['flux'][-1,1,:,:], cmap='Spectral',ax=ax)
shp = flopy.plot.plot_shapefile('../Shp/ModelWell2', ax=ax, radius=10)
output_44_0.png

River path as RIV boundary condition

river_intersect = g.intersect('../Shp/ModelRiver2', 'polygon', 0)
river_spd = {}
river_spd[0] = []
river_unique = np.unique(river_intersect.nodenumber)
for i in np.arange(river_unique.shape[0]):
    x,y = g.get_center(river_unique[i])                           # river_intersect.nodenumber[i]
    i,j = mf.sr.get_ij(x, y)
    river_spd[0].append([0,i,j,mtop[i,j],0.01,mtop[i,j]-1])         # the array follow this order: [lay, row, col, stage, cond, rbot]
riv = flopy.modflow.ModflowRiv(mf, stress_period_data=river_spd)
# Plot the river shapefile and the cells conceptualized as river
fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(1, 1, 1, aspect='equal')
modelmap = flopy.plot.PlotMapView(model=mf)
g.plot(ax, linewidth=0.5)
quadmesh = modelmap.plot_array(mf.riv.stress_period_data.array['cond'][-1,0,:,:], cmap='Spectral',ax=ax, alpha=0.5)
shp = flopy.plot.plot_shapefile('../Shp/ModelRiver2', ax=ax,facecolor='none')
output_48_0.png

Regional flow as GHB boundary condition

ghb_intersect = g.intersect('../Shp/ModelGHB1', 'line', 0)
ghb_spd = {}
ghb_spd[0] = []
ghb_unique = np.unique(ghb_intersect.nodenumber)
ghb_unique
array([  0,  38,  39,  77,  78, 116, 117, 155, 156, 194, 195, 233, 234,
       272, 273, 312, 351, 390, 429, 506, 545, 584, 623, 662, 663, 701,
       702, 740, 741, 779, 780, 818, 819, 857, 858, 896, 897, 935])
for i in np.arange(ghb_unique.shape[0]):
    x,y = g.get_center(ghb_unique[i])                           # river_intersect.nodenumber[i]
    i,j = mf.sr.get_ij(x, y)
    if j < 5:
        ghb_spd[0].append([0,i,j,55,0.01]) # the array follow this order: [lay, row, col, elev, cond]
    elif j > ncols -5:
        ghb_spd[0].append([0,i,j,90,0.01])
    else: print("GHB on the wrong position")
ghb = flopy.modflow.ModflowGhb(mf, stress_period_data=ghb_spd)
# Plot the GHB shapefile and the cells conceptualized as GHB
fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(1, 1, 1, aspect='equal')
modelmap = flopy.plot.PlotMapView(model=mf)
g.plot(ax, linewidth=0.5)
quadmesh = modelmap.plot_array(mf.ghb.stress_period_data.array['cond'][-1,0,:,:], cmap='RdGy',ax=ax, alpha=0.2)
shp = flopy.plot.plot_shapefile('../Shp/ModelGHB1', ax=ax,linewidth=8)
output_53_0.png
### Set the Output Control and run simulation
oc_spd = {(0, 0): ['save head']}

for i in range(mf.dis.nstp.shape[0]):
    oc_spd[(i,mf.dis.nstp.array[3]-1)] = ['save head']
oc_spd
{(0, 0): ['save head'],
 (0, 3): ['save head'],
 (1, 3): ['save head'],
 (2, 3): ['save head'],
 (3, 3): ['save head'],
 (4, 3): ['save head'],
 (5, 3): ['save head'],
 (6, 3): ['save head'],
 (7, 3): ['save head'],
 (8, 3): ['save head'],
 (9, 3): ['save head'],
 (10, 3): ['save head']}
oc = flopy.modflow.ModflowOc(mf, stress_period_data=oc_spd)
mf.write_input()
mf.run_model()
FloPy is using the following  executable to run the model: ../Exe/MODFLOW-NWT_64.exe

                                  MODFLOW-NWT-SWR1 
    U.S. GEOLOGICAL SURVEY MODULAR FINITE-DIFFERENCE GROUNDWATER-FLOW MODEL
                             WITH NEWTON FORMULATION
                             Version 1.1.4 4/01/2018                         
                    BASED ON MODFLOW-2005 Version 1.12.0 02/03/2017                       

                    SWR1 Version 1.04.0 09/15/2016                       

 Using NAME file: Model1.nam 
 Run start date and time (yyyy/mm/dd hh:mm:ss): 2019/06/10 10:57:55

 Solving:  Stress period:     1    Time step:     1    Groundwater-Flow Eqn.
 Solving:  Stress period:     2    Time step:     1    Groundwater-Flow Eqn.
 Solving:  Stress period:     2    Time step:     2    Groundwater-Flow Eqn.
 Solving:  Stress period:     2    Time step:     3    Groundwater-Flow Eqn.
 Solving:  Stress period:     2    Time step:     4    Groundwater-Flow Eqn.
 Solving:  Stress period:     3    Time step:     1    Groundwater-Flow Eqn.
 Solving:  Stress period:     3    Time step:     2    Groundwater-Flow Eqn.
 Solving:  Stress period:     3    Time step:     3    Groundwater-Flow Eqn.
 Solving:  Stress period:     3    Time step:     4    Groundwater-Flow Eqn.
 Solving:  Stress period:     4    Time step:     1    Groundwater-Flow Eqn.
 Solving:  Stress period:     4    Time step:     2    Groundwater-Flow Eqn.
 Solving:  Stress period:     4    Time step:     3    Groundwater-Flow Eqn.
 Solving:  Stress period:     4    Time step:     4    Groundwater-Flow Eqn.
 Solving:  Stress period:     5    Time step:     1    Groundwater-Flow Eqn.
 Solving:  Stress period:     5    Time step:     2    Groundwater-Flow Eqn.
 Solving:  Stress period:     5    Time step:     3    Groundwater-Flow Eqn.
 Solving:  Stress period:     5    Time step:     4    Groundwater-Flow Eqn.
 Solving:  Stress period:     6    Time step:     1    Groundwater-Flow Eqn.
 Solving:  Stress period:     6    Time step:     2    Groundwater-Flow Eqn.
 Solving:  Stress period:     6    Time step:     3    Groundwater-Flow Eqn.
 Solving:  Stress period:     6    Time step:     4    Groundwater-Flow Eqn.
 Solving:  Stress period:     7    Time step:     1    Groundwater-Flow Eqn.
 Solving:  Stress period:     7    Time step:     2    Groundwater-Flow Eqn.
 Solving:  Stress period:     7    Time step:     3    Groundwater-Flow Eqn.
 Solving:  Stress period:     7    Time step:     4    Groundwater-Flow Eqn.
 Solving:  Stress period:     8    Time step:     1    Groundwater-Flow Eqn.
 Solving:  Stress period:     8    Time step:     2    Groundwater-Flow Eqn.
 Solving:  Stress period:     8    Time step:     3    Groundwater-Flow Eqn.
 Solving:  Stress period:     8    Time step:     4    Groundwater-Flow Eqn.
 Solving:  Stress period:     9    Time step:     1    Groundwater-Flow Eqn.
 Solving:  Stress period:     9    Time step:     2    Groundwater-Flow Eqn.
 Solving:  Stress period:     9    Time step:     3    Groundwater-Flow Eqn.
 Solving:  Stress period:     9    Time step:     4    Groundwater-Flow Eqn.
 Solving:  Stress period:    10    Time step:     1    Groundwater-Flow Eqn.
 Solving:  Stress period:    10    Time step:     2    Groundwater-Flow Eqn.
 Solving:  Stress period:    10    Time step:     3    Groundwater-Flow Eqn.
 Solving:  Stress period:    10    Time step:     4    Groundwater-Flow Eqn.
 Solving:  Stress period:    11    Time step:     1    Groundwater-Flow Eqn.
 Solving:  Stress period:    11    Time step:     2    Groundwater-Flow Eqn.
 Solving:  Stress period:    11    Time step:     3    Groundwater-Flow Eqn.
 Solving:  Stress period:    11    Time step:     4    Groundwater-Flow Eqn.
 Run end date and time (yyyy/mm/dd hh:mm:ss): 2019/06/10 10:57:55
 Elapsed run time:  0.242 Seconds

  Normal termination of simulation





(True, [])

Model output visualization

mfheads = flopy.utils.HeadFile('../Model/Model1.hds')
mfheads.get_times()
[1.0,
 17280000.0,
 34560000.0,
 51840000.0,
 69120000.0,
 86400000.0,
 103680000.0,
 120960000.0,
 138240000.0,
 155520000.0,
 172800000.0]
head = mfheads.get_data()
head.shape
(3, 24, 39)
# Plot the heads for a defined layer and boundary conditions
fig = plt.figure(figsize=(36, 24))
ax = fig.add_subplot(1, 1, 1, aspect='equal')
modelmap = flopy.plot.PlotMapView(model=mf)
g.plot(ax, linewidth=0.5, alpha=0.4)
contour = modelmap.contour_array(head[1],ax=ax)

quadmesh = modelmap.plot_array(mf.riv.stress_period_data.array['cond'][-1,0,:,:], cmap='Spectral',ax=ax, alpha=0.5)

cellhead = modelmap.plot_array(head[1],ax=ax, cmap='Blues', alpha=0.2)
ax.clabel(contour)

plt.tight_layout()

# Boundary conditions
shpGHB = flopy.plot.plot_shapefile('../Shp/ModelGHB1', ax=ax,linewidth=8)
shpWel = flopy.plot.plot_shapefile('../Shp/ModelWell2', ax=ax,radius=20)
shpRiv = flopy.plot.plot_shapefile('../Shp/ModelRiver2', ax=ax)

# Export heads on layer two at the beggining of simulation
contourShp = mfheads.to_shapefile(filename='../Shp/ModelHeadsLay2',kstpkper=(0,0),mflay=1)
fig.savefig('infahatari.png')
wrote ../Shp/ModelHeadsLay2
output_61_1.png
#Representation of model geometry with all the boundary conditions
fig = plt.figure(figsize=(18,8))

modelmap = flopy.plot.PlotMapView(model=mf)
linecollection = modelmap.plot_grid(linewidth=0.5, color='cyan')
shpRiver = flopy.plot.plot_shapefile('../Shp/ModelRiver2', facecolor='none') #RIV boundary condition
shpGHB = flopy.plot.plot_shapefile('../Shp/ModelGHB1', facecolor='none')     #GHB boundary condition

contour = modelmap.contour_array(head[1])
plt.plot(limitArray[:,0],limitArray[:,1])
plt.scatter(wellArray[:,0],wellArray[:,1])

tiles=('http://mt0.google.com/vt/lyrs=s&hl=en&x={x}&y={y}&z={z}', 
       '<a href=https://www.openstreetmap.org/about">© OpenStreetMap</a>')
crs = {'init' :'epsg:32718'}
mplleaflet.display(fig, crs=crs, tiles=tiles)
fig = plt.figure(figsize=(20, 3))
ax = fig.add_subplot(1, 1, 1)
modelxsect = flopy.plot.PlotCrossSection(model=mf, line={'Row': 10})
linecollection = modelxsect.plot_grid()
modelxsect.contour_array(head)
output_62_1.png


Datos de entrada

Usted puede descargar los datos de entrada para este tutorial de este enlace.

Smiley face

Suscríbase a nuestro boletín gratuito para recibir noticias, datos interesantes y fechas de nuestros cursos en recursos hídricos.

Posted on June 13, 2019 and filed under TutorialModflow.