مدلسازی آب های زیر زمینی منطقه ای با MODFLOW و Flopy :: بیسین - سایت تخصصی مهندسی آب

مدلسازی آب های زیر زمینی منطقه ای با MODFLOW و Flopy


مدل سازی منطقه ای آب های زیرزمینی یک وظیفه مهم در مدیریت آب استراتژیک است که شامل تمام کاربران، فعالیت ها و اکوسیستم های مربوطه و استفاده پایدار از شرایط فعلی و آینده است. برخی از ملاحظات خاص در مورد مدلسازی منطقه ای با توجه به تقسیم بندی اولیه و فضایی وجود دارد؛ یک مدل منطقه ای در نظر گرفته نشده است که پاسخ منطقه آبخوان را برای منطقه مشخص تعیین کند، در عوض شامل ارزیابی جریان آب های زیرزمینی منطقه و اندازه گیری شارژ، تخلیه و سایر فرآیندهای موجود در تعادل آب است.


این آموزش نسخه Flopy از نمونه عددی از حوضه Angascancha با MODFLOW و Model Muse است که شما می توانید در مطالب مشابه پیدا کنید. مثال در حالت پایدار است و با حل کننده NWT حل شده است. نمایندگی مدل خروجی تحت ابزار Flopy / Matplotlib و همچنین بعضی از کد های Python برای ایجاد فایل های VTU و Style در Paraview انجام شده است.


فایل های ورودی

شما می توانید فایل های ورودی این آموزش را در این لینک دانلود کنید.


Python code

This is the described Python code for this tutorial:

Import the required libraries

This tutorial requires some Python core libraries, Scipy libraries (numpy and matplotlib) as well as the Flopy library. For the management and data extraction from raster files we use the GDAL library.

import flopy
import os
import flopy.utils.binaryfile as bf
import numpy as np
from osgeo import gdal
import matplotlib.pyplot as plt

Create a basic MODFLOW model object

As well as any other model generated by Flopy, first we have to setup the model name and the working directory. We strongly recommend to follow the folder/file estructure provided on the "Input Files" of this tutorial. Based on the complex geometry and the differences among layers elevations and layer thickness, the NWT solver has been applied to the model.

modelname = "model1"
modelpath = "../Model/"
#Modflow Nwt
mf1 = flopy.modflow.Modflow(modelname, exe_name="../Exe/MODFLOW-NWT_64.exe", version="mfnwt",model_ws=modelpath)
nwt = flopy.modflow.ModflowNwt(mf1 ,maxiterout=150,linmeth=2)

Open and read raster files

In this part we will use the GDAL library to open the rasters of elevation and boundary conditions and get their geotransformation parameters. The tutorial requires that the raster resolution is the same as the grid resolution. Raster bands are read as Numpy arrays and inserted to the model as model top.

#Raster paths
demPath ="../Rst/DEM_200.tif"
crPath = "../Rst/CR.tif"

#Open files
demDs =gdal.Open(demPath)
crDs = gdal.Open(crPath)
geot = crDs.GetGeoTransform() #Xmin, deltax, ?, ymax, ?, delta y

geot
(613882.5068, 200.0, 0.0, 8375404.9524, 0.0, -200.0)
# Get data as arrays
demData = demDs.GetRasterBand(1).ReadAsArray()
crData = crDs.GetRasterBand(1).ReadAsArray()

# Get statistics
stats = demDs.GetRasterBand(1).GetStatistics(0,1) 
stats
[3396.0, 4873.0, 4243.0958193041, 287.90182492684]

Spatial discretization

Since this tutorial deals with a regional scale, the differences in elevation on the area of study can be more than a hundred meters or more than a thousand meters if we deal with andean basins. The hydrogeological conceptualization for this tutorial is the definition of two upper thin layers (Layer1 and Layer1_2) and thicker layers below.

#Elevation of layers
Layer1 = demData - 7.5
Layer1_2 = demData - 15
Layer2 = (demData - 3000)*0.8 + 3000
Layer3 = (demData - 3000)*0.5 + 3000
Layer4 = (demData>0)*3000

Column and row number and dimension come from the raster attributes

#Boundaries for Dis = Create discretization object, spatial/temporal discretization
ztop = demData
zbot = [Layer1, Layer1_2, Layer2, Layer3, Layer4]
nlay = 5
nrow = demDs.RasterYSize
ncol = demDs.RasterXSize
delr = geot[1]
delc = abs(geot[5])




Definition of the flow packages for the MODFLOW model

In this part we define the packages related to groundwater flow on the MODFLOW model. First we define the DIS package that has the geometry as well as the simulation type (steady / transient). The model run on steady conditions.

dis = flopy.modflow.ModflowDis(mf1, nlay,nrow,ncol,delr=delr,delc=delc,top=ztop,botm=zbot,itmuni=1)

Then we define another MODEL packages. Some boundary conditions are setup from raster valures, while other packages use constant value. These are the required packages for the simulation:

  • the BAS package for setting the constant head cells and inactive cells,
  • the UPW that defines the vertical / horizontal hydraulic conductivity,
  • the RCH package that applies rechage from precipitation
  • the EVT package for the simulation of groundwater consumption from roots
  • the DRN package for the groundwater discharge to the channel network
  • the OC packages for the output record
# Variables for the BAS package
iboundData = np.zeros(demData.shape, dtype=np.int32)
iboundData[crData > 0 ] = 1
bas = flopy.modflow.ModflowBas(mf1,ibound=iboundData,strt=ztop, hnoflo=-2.0E+020)
# Array of hydraulic heads per layer
hk = [1E-4, 1E-5, 1E-7, 1E-8, 1E-9]

# Add UPW package to the MODFLOW model
upw = flopy.modflow.ModflowUpw(mf1, laytyp = [1,1,1,1,0], hk = hk)
#Add the recharge package (RCH) to the MODFLOW model
rch = np.ones((nrow, ncol), dtype=np.float32) * 0.120/365/86400
rch_data = {0: rch}
rch = flopy.modflow.ModflowRch(mf1, nrchop=3, rech =rch_data)
#Add the evapotranspiration package (EVT) to the MODFLOW model
evtr = np.ones((nrow, ncol), dtype=np.float32) * 1/365/86400
evtr_data = {0: evtr}
evt = flopy.modflow.ModflowEvt(mf1,nevtop=1,surf=ztop,evtr=evtr_data, exdp=0.5)
#Add the drain package (DRN) to the MODFLOW model
river = np.zeros(demData.shape, dtype=np.int32)
river[crData == 3 ] = 1
list = []
for i in range(river.shape[0]):
    for q in range(river.shape[1]):
        if river[i,q] == 1:
            list.append([0,i,q,ztop[i,q],0.001]) #layer,row,column,elevation(float),conductanceee
rivDrn = {0:list}

drn = flopy.modflow.ModflowDrn(mf1, stress_period_data=rivDrn)
# Add OC package to the MODFLOW model
oc = flopy.modflow.ModflowOc(mf1) #ihedfm= 1, iddnfm=1

Write files of the MODFLOW model and run simulation

Write the .nam file and the files declared on the .nam file. Then runs the simulation and shows the simulation information on the screen

#Write input files -> write file with extensions
mf1.write_input()

#run model -> gives the solution
mf1.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/02/22 10:19:24

 Solving:  Stress period:     1    Time step:     1    Groundwater-Flow Eqn.
 Run end date and time (yyyy/mm/dd hh:mm:ss): 2019/02/22 10:20:04
 Elapsed run time: 39.475 Seconds

  Normal termination of simulation





(True, [])




Model results post-processing

This tutorial has some model features and output representation done with Flopy and Matplotlib. Another 3D representation on the VTK format is performed on another notebook from this tutorial

#Import model
ml = flopy.modflow.Modflow.load('../Model/model1.nam')

Model grid representation

# First step is to set up the plot
fig = plt.figure(figsize=(15, 15))
ax = fig.add_subplot(1, 1, 1, aspect='equal')

# Next we create an instance of the ModelMap class
modelmap = flopy.plot.ModelMap(sr=ml.dis.sr)

# Then we can use the plot_grid() method to draw the grid
# The return value for this function is a matplotlib LineCollection object,
# which could be manipulated (or used) later if necessary.
linecollection = modelmap.plot_grid(linewidth=0.4)
Cross section of model grid representation
fig = plt.figure(figsize=(15, 6))
ax = fig.add_subplot(1, 1, 1)
# Next we create an instance of the ModelCrossSection class
#modelxsect = flopy.plot.ModelCrossSection(model=ml, line={'Column': 5})
modelxsect = flopy.plot.ModelCrossSection(model=ml, line={'Row': 99})

# Then we can use the plot_grid() method to draw the grid
# The return value for this function is a matplotlib LineCollection object,
# which could be manipulated (or used) later if necessary.
linecollection = modelxsect.plot_grid(linewidth=0.4)
t = ax.set_title('Column 6 Cross-Section - Model Grid')
Active/inactive cells on model extension
fig = plt.figure(figsize=(15, 15))
ax = fig.add_subplot(1, 1, 1, aspect='equal')
modelmap = flopy.plot.ModelMap(model=ml, rotation=0)
quadmesh = modelmap.plot_ibound(color_noflow='cyan')
linecollection = modelmap.plot_grid(linewidth=0.4)
Cross sections of active/inactive cells
fig = plt.figure(figsize=(15, 6))
ax = fig.add_subplot(1, 1, 1)
modelxsect = flopy.plot.ModelCrossSection(model=ml, line={'Column': 5})
patches = modelxsect.plot_ibound(color_noflow='cyan')
linecollection = modelxsect.plot_grid(linewidth=0.4)
t = ax.set_title('Column 6 Cross-Section with IBOUND Boundary Conditions')
Channel network as drain (DRN) package
fig = plt.figure(figsize=(15, 15))
ax = fig.add_subplot(1, 1, 1, aspect='equal')
modelmap = flopy.plot.ModelMap(model=ml, rotation=-14)
quadmesh = modelmap.plot_ibound(color_noflow='cyan')
quadmesh = modelmap.plot_bc('DRN', color='blue')
linecollection = modelmap.plot_grid(linewidth=0.4)
Model grid and heads representation
fname = os.path.join(modelpath, 'model1.hds')
hdobj = flopy.utils.HeadFile(fname)
head = hdobj.get_data()

fig = plt.figure(figsize=(30, 30))
ax = fig.add_subplot(1, 2, 1, aspect='equal')
modelmap = flopy.plot.ModelMap(model=ml, rotation=-14)
#quadmesh = modelmap.plot_ibound()
quadmesh = modelmap.plot_array(head, masked_values=[-2.e+20], alpha=0.8)
linecollection = modelmap.plot_grid(linewidth=0.2)





نظرات (۰)

فرم ارسال نظر

ارسال نظر آزاد است، اما اگر قبلا در بیان ثبت نام کرده اید می توانید ابتدا وارد شوید.
شما میتوانید از این تگهای html استفاده کنید:
<b> یا <strong>، <em> یا <i>، <u>، <strike> یا <s>، <sup>، <sub>، <blockquote>، <code>، <pre>، <hr>، <br>، <p>، <a href="" title="">، <span style="">، <div align="">
تجدید کد امنیتی


درباره بهترين هاي بيسيـــن بدانيد...

Bird

يکي از مهمترين اهداف اين سايت تهيه آموزش هاي روان از ابزارهاي کاربردي علوم آب است.

اهميت مطالعات محيطي با ابزارهاي نوين در چيست؟

امروز با فارغ التحصيلي جمع کثير دانشجويان سالهاي گذشته و حال، با گذر از کمي گرايي ديگر صرف وجود مدارک دانشگاهي حرف اول را در بازار کار نمي زند؛ بلکه سنجش ديگري ملاک؛ و شايسته سالاري به ناچار! باب خواهد شد. يکي از مهم ترين لوازم توسعه علمي در هر کشور و ارائه موضوعات ابتکاري، بهره گيري از ابزار نوين است، بيسين با همکاري مخاطبان مي تواند در حيطه علوم آب به معرفي اين مهم بپردازد.

جستجو در بيسين


بیسین - سایت تخصصی مهندسی آب

سایت مهندسی آب بیسین با معرفی مهم ترین و کاربردی ترین نرم افزارها و مدل های شبیه سازی در حیطه مهندسی آب، تلاش به تهیه خدمات یکپارچه و محلی از محاسبات هیدرولوژیکی و هیدرولیکی می کند

W3Schools


اطلاعات سايت

  • behzadsarhadi@gmail.com
  • بهزاد سرهادي
  • شناسه تلگرام: SubBasin
  • شماره واتساپ: 09190622992-098
  • شماره تماس: 09190622992-098

W3Schools