906 lines
28 KiB
Plaintext
906 lines
28 KiB
Plaintext
{
|
||
"cells": [
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"# Demonstrate use of SFR/LAK/UZF with SFT/LKT/UZT\n",
|
||
"A more comprehensive demonstration of setting up an MT3D-USGS model that uses all of the new packages included in the first release of MT3D-USGS. Also includes RCT.\n",
|
||
"\n",
|
||
"#### Problem Description:\n",
|
||
"* 300 row x 300 col x 3 layer x 2 stress period model\n",
|
||
"* Flow model uses SFR, LAK, and UZF with connections between all three\n",
|
||
"* Transport model simulates streamflow transport (SFT), with connction to a single lake (LKT)\n",
|
||
"* Transport model simulates overland runoff and spring discharge (UZT) to surface water network\n",
|
||
"\n",
|
||
"Start by importing some libraries:"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 1,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"3.8.6 | packaged by conda-forge | (default, Oct 7 2020, 18:42:56) \n",
|
||
"[Clang 10.0.1 ]\n",
|
||
"numpy version: 1.18.5\n",
|
||
"matplotlib version: 3.2.2\n",
|
||
"flopy version: 3.3.3\n"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"import os\n",
|
||
"import sys\n",
|
||
"import platform\n",
|
||
"import numpy as np\n",
|
||
"import pandas as pd\n",
|
||
"import matplotlib as mpl\n",
|
||
"import matplotlib.pyplot as plt\n",
|
||
"\n",
|
||
"\n",
|
||
"# run installed version of flopy or add local path\n",
|
||
"try:\n",
|
||
" import flopy\n",
|
||
"except:\n",
|
||
" fpth = os.path.abspath(os.path.join('..', '..'))\n",
|
||
" sys.path.append(fpth)\n",
|
||
" import flopy\n",
|
||
"\n",
|
||
"print(sys.version)\n",
|
||
"print('numpy version: {}'.format(np.__version__))\n",
|
||
"print('matplotlib version: {}'.format(mpl.__version__))\n",
|
||
"print('flopy version: {}'.format(flopy.__version__))"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"Create a MODFLOW model and store it, in this case in the variable 'mf'. \n",
|
||
"The modelname will be the name given to all MODFLOW files.\n",
|
||
"The exe_name should be the name of the MODFLOW executable. \n",
|
||
"In this case, we want to use version: 'mfnwt' for MODFLOW-NWT"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 2,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"modelpth = os.path.join('.', 'temp', 'no3')\n",
|
||
"modelname = 'no3'\n",
|
||
"mfexe = 'mfnwt'\n",
|
||
"mtexe = 'mt3dusgs'\n",
|
||
"\n",
|
||
"if platform.system() == 'Windows':\n",
|
||
" mfexe += '.exe'\n",
|
||
" mtexe += '.exe'\n",
|
||
"\n",
|
||
"# Make sure modelpth directory exists\n",
|
||
"if not os.path.exists(modelpth):\n",
|
||
" os.makedirs(modelpth)\n",
|
||
"\n",
|
||
"# Instantiate MODFLOW object in flopy\n",
|
||
"mf = flopy.modflow.Modflow(modelname=modelname, exe_name=mfexe, model_ws=modelpth, version='mfnwt')"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Set up model discretization"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 3,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"Lx = 90000.0 \n",
|
||
"Ly = 90000.0\n",
|
||
"nrow = 300\n",
|
||
"ncol = 300\n",
|
||
"nlay = 3\n",
|
||
"\n",
|
||
"delr = Lx / ncol\n",
|
||
"delc = Ly / nrow\n",
|
||
"\n",
|
||
"xmax = ncol * delr\n",
|
||
"ymax = nrow * delc\n",
|
||
"\n",
|
||
"X, Y = np.meshgrid(np.linspace(delr / 2, xmax - delr / 2, ncol), \n",
|
||
" np.linspace(ymax - delc / 2, 0 + delc / 2, nrow))"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Instantiate output control (oc) package for MODFLOW-NWT"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 4,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"oc = flopy.modflow.ModflowOc(mf)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Instantiate solver package for MODFLOW-NWT"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 5,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"# Newton-Raphson Solver: Create a flopy nwt package object\n",
|
||
"\n",
|
||
"headtol = 1.0E-4 \n",
|
||
"fluxtol = 5 \n",
|
||
"maxiterout = 5000\n",
|
||
"thickfact = 1E-06\n",
|
||
"linmeth = 2 \n",
|
||
"iprnwt = 1 \n",
|
||
"ibotav = 1 \n",
|
||
"\n",
|
||
"nwt = flopy.modflow.ModflowNwt(mf, headtol=headtol, fluxtol=fluxtol, maxiterout=maxiterout,\n",
|
||
" thickfact=thickfact, linmeth=linmeth, iprnwt=iprnwt, ibotav=ibotav,\n",
|
||
" options='SIMPLE')"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Instantiate discretization (DIS) package for MODFLOW-NWT"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 6,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"elv_pth = os.path.join('..','data','mt3d_example_sft_lkt_uzt','dis_arrays','grnd_elv.txt')\n",
|
||
"\n",
|
||
"# Top of Layer 1 elevation determined using GW Vistas and stored locally\n",
|
||
"grndElv = np.loadtxt(elv_pth)\n",
|
||
"\n",
|
||
"# Bottom of layer 1 elevation also determined from use of GUI and stored locally\n",
|
||
"bt1_pth = os.path.join('..','data','mt3d_example_sft_lkt_uzt','dis_arrays','bot1.txt')\n",
|
||
"bot1Elv = np.loadtxt(bt1_pth)\n",
|
||
"\n",
|
||
"bot2Elv = np.ones(bot1Elv.shape) * 100\n",
|
||
"bot3Elv = np.zeros(bot2Elv.shape)\n",
|
||
"\n",
|
||
"botm = [bot1Elv, bot2Elv, bot3Elv]\n",
|
||
"botm = np.array(botm)\n",
|
||
"Steady = [False, False]\n",
|
||
"nstp = [1, 1]\n",
|
||
"tsmult = [1., 1.]\n",
|
||
"\n",
|
||
"# Stress periods \n",
|
||
"perlen = [9131.25, 9131.25]\n",
|
||
"\n",
|
||
"# Create the discretization object\n",
|
||
"# itmuni = 4 (days); lenuni = 1 (feet)\n",
|
||
"dis = flopy.modflow.ModflowDis(mf, nlay, nrow, ncol, nper=2, delr=delr, delc=delc, \n",
|
||
" top=grndElv, botm=botm, laycbd=0, itmuni=4, lenuni=1,\n",
|
||
" steady=Steady, nstp=nstp, tsmult=tsmult, perlen=perlen) \n"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Instantiate upstream weighting (UPW) flow package for MODFLOW-NWT\n"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 7,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"# UPW must be instantiated after DIS. Otherwise, during the mf.write_input() procedures, \n",
|
||
"# flopy will crash.\n",
|
||
"\n",
|
||
"# First line of UPW input is: IUPWCB HDRY NPUPW IPHDRY \n",
|
||
"hdry = -1.00e+30\n",
|
||
"iphdry = 0\n",
|
||
"\n",
|
||
"# Next variables are: LAYTYP, LAYAVG, CHANI, LAYVKA, LAYWET\n",
|
||
"laytyp = [1, 3, 3] # >0: convertible\n",
|
||
"layavg = 0 # 0: harmonic mean\n",
|
||
"chani = 1.0 # >0: CHANI is the horizontal anisotropy for the entire layer\n",
|
||
"layvka = 0 # =0: indicates VKA is vertical hydraulic conductivity\n",
|
||
"laywet = 0 # Always set equal to zero in UPW package\n",
|
||
"hk = 20\n",
|
||
"# hani = 1 # Not needed because CHANI > 1\n",
|
||
"vka = 0.5 # Is equal to vert. K b/c LAYVKA = 0\n",
|
||
"ss = 0.00001\n",
|
||
"sy = 0.20\n",
|
||
"\n",
|
||
"upw = flopy.modflow.ModflowUpw(mf, laytyp=laytyp, layavg=layavg, chani=chani, layvka=layvka,\n",
|
||
" laywet=laywet, ipakcb=53, hdry=hdry, iphdry=iphdry, hk=hk, \n",
|
||
" vka=vka, ss=ss, sy=sy)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Instantiate basic (BAS or BA6) package for MODFLOW-NWT"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 8,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"\n",
|
||
"ibnd1_pth = os.path.join('..','data','mt3d_example_sft_lkt_uzt','bas_arrays','ibnd_lay1.txt')\n",
|
||
"ibnd1 = np.loadtxt(ibnd1_pth)\n",
|
||
"ibnd2 = np.ones(ibnd1.shape)\n",
|
||
"ibnd3 = np.ones(ibnd2.shape)\n",
|
||
"\n",
|
||
"ibnd = [ibnd1, ibnd2, ibnd3]\n",
|
||
"ibnd = np.array(ibnd)\n",
|
||
"\n",
|
||
"StHd1_pth = os.path.join('..','data','mt3d_example_sft_lkt_uzt','bas_arrays','strthd1.txt')\n",
|
||
"StHd1 = np.loadtxt(StHd1_pth)\n",
|
||
"\n",
|
||
"StHd2_pth = os.path.join('..','data','mt3d_example_sft_lkt_uzt','bas_arrays','strthd2.txt')\n",
|
||
"StHd2 = np.loadtxt(StHd2_pth)\n",
|
||
"\n",
|
||
"StHd3_pth = os.path.join('..','data','mt3d_example_sft_lkt_uzt','bas_arrays','strthd3.txt')\n",
|
||
"StHd3 = np.loadtxt(StHd3_pth)\n",
|
||
"\n",
|
||
"strtElev = [StHd1, StHd2, StHd3]\n",
|
||
"strtElev = np.array(strtElev)\n",
|
||
"\n",
|
||
"hdry = 999.0\n",
|
||
"\n",
|
||
"bas = flopy.modflow.ModflowBas(mf, ibound=ibnd, hnoflo=hdry, strt=strtElev)\n"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Instantiate general head boundary (GHB) package for MODFLOW-NWT"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 9,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"# GHB boundaries are located along the top (north) and bottom (south) \n",
|
||
"# edges of the domain, all 3 layers.\n",
|
||
"\n",
|
||
"elev_stpt_row1 = 308.82281\n",
|
||
"elev_stpt_row300 = 239.13811\n",
|
||
"elev_slp = (308.82281 - 298.83649) / (ncol - 1)\n",
|
||
"\n",
|
||
"sp = []\n",
|
||
"for k in [0,1,2]: # These indices need to be adjusted for 0-based moronicism\n",
|
||
" for i in [0,299]: # These indices need to be adjusted for 0-based silliness\n",
|
||
" for j in np.arange(0,300,1): # These indices need to be adjusted for 0-based foolishness\n",
|
||
" # Skipping cells not satisfying the conditions below\n",
|
||
" if ((i == 1 and (j < 27 or j > 31)) or (i==299 and (j < 26 or j > 31))):\n",
|
||
" if (i % 2 == 0):\n",
|
||
" sp.append([k, i, j, elev_stpt_row1 - (elev_slp * (j-1)), 11.3636])\n",
|
||
" else:\n",
|
||
" sp.append([k, i, j, elev_stpt_row300 - (elev_slp * (j-1)), 11.3636])\n",
|
||
"\n",
|
||
"\n",
|
||
"for k in [0,1,2]:\n",
|
||
" for j in np.arange(26,32,1):\n",
|
||
" sp.append([k, 299, j, 238.20, 3409.0801])\n",
|
||
"\n",
|
||
"ghb = flopy.modflow.ModflowGhb(mf, stress_period_data = sp)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Instantiate streamflow routing (SFR2) package for MODFLOW-NWT"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 10,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"# Read pre-prepared reach data into numpy recarrays using numpy.genfromtxt()\n",
|
||
"# Remember that the cell indices stored in the pre-prepared NO3_ReachInput.csv file are based on 0-based indexing.\n",
|
||
"# Flopy will convert to 1-based when it writes the files\n",
|
||
"\n",
|
||
"rpth = os.path.join('..','data','mt3d_example_sft_lkt_uzt','sfr_data', 'no3_reachinput.csv')\n",
|
||
"reach_data = np.genfromtxt(rpth, delimiter=',', names=True)\n",
|
||
"reach_data\n",
|
||
"\n",
|
||
"# Read pre-prepared segment data into numpy recarrays using numpy.genfromtxt()\n",
|
||
"\n",
|
||
"spth = os.path.join('..','data','mt3d_example_sft_lkt_uzt','sfr_data', 'no3_segmentdata.csv')\n",
|
||
"ss_segment_data = np.genfromtxt(spth, delimiter=',', names=True)\n",
|
||
"segment_data = {0: ss_segment_data,\n",
|
||
" 1: ss_segment_data}\n",
|
||
"segment_data[0][0:1]['width1']\n",
|
||
"\n",
|
||
"nstrm = len(reach_data)\n",
|
||
"nss = len(segment_data[0])\n",
|
||
"nsfrpar = 0\n",
|
||
"const = 128390.4 # constant for manning's equation, units of cfs\n",
|
||
"dleak = 0.0001\n",
|
||
"ipakcb = 53 # flag for writing SFR output to cell-by-cell budget (on unit 53)\n",
|
||
"istcb2 = 37 # flag for writing SFR output to text file\n",
|
||
"isfropt = 1\n",
|
||
"dataset_5 = {0: [nss, 0, 0], 1: [-1, 0, 0]} # dataset 5 (see online guide) (ITMP, IRDFLG, IPTFLG)\n",
|
||
"\n",
|
||
"# Input arguments generally follow the variable names defined in the Online Guide to MODFLOW\n",
|
||
"sfr = flopy.modflow.ModflowSfr2(mf, nstrm=nstrm, nss=nss, const=const, dleak=dleak, ipakcb=ipakcb,\n",
|
||
" istcb2=istcb2, isfropt=isfropt, reachinput=True,\n",
|
||
" reach_data=reach_data,\n",
|
||
" segment_data=segment_data,\n",
|
||
" dataset_5=dataset_5, unit_number=15)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Instantiate Lake (LAK) package for MODFLOW-NWT"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 11,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"# Read pre-prepared lake arrays\n",
|
||
"LakArr_pth = os.path.join('..','data','mt3d_example_sft_lkt_uzt','lak_arrays','lakarr1.txt')\n",
|
||
"LakArr_lyr1 = np.loadtxt(LakArr_pth)\n",
|
||
"LakArr_lyr2 = np.zeros(LakArr_lyr1.shape)\n",
|
||
"LakArr_lyr3 = np.zeros(LakArr_lyr2.shape)\n",
|
||
"\n",
|
||
"LakArr = [LakArr_lyr1, LakArr_lyr2, LakArr_lyr3]\n",
|
||
"LakArr = np.array(LakArr)\n",
|
||
"\n",
|
||
"nlakes = int(np.max(LakArr))\n",
|
||
"ipakcb = ipakcb # From above\n",
|
||
"theta = -1. # Implicit\n",
|
||
"nssitr = 10 # Maximum number of iterations for Newton’s method\n",
|
||
"sscncr = 1.000e-03 # Convergence criterion for equilibrium lake stage solution\n",
|
||
"surfdep = 2.000e+00 # Height of small topological variations in lake-bottom\n",
|
||
"stages = 268.00 # Initial stage of each lake at the beginning of the run\n",
|
||
"\n",
|
||
"# ITMP > 0, read lake definition data \n",
|
||
"# ITMP1 ≥ 0, read new recharge, evaporation, runoff, and withdrawal data for each lake\n",
|
||
"# LWRT > 0, suppresses printout from the lake package\n",
|
||
"\n",
|
||
"bdlknc_lyr1 = LakArr_lyr1.copy()\n",
|
||
"bdlknc_lyr2 = LakArr_lyr1.copy()\n",
|
||
"bdlknc_lyr3 = np.zeros((LakArr_lyr1.shape))\n",
|
||
"\n",
|
||
"# Need to expand bdlknc_lyr1 non-zero values by 1 in either direction\n",
|
||
"# (left/right and up/down)\n",
|
||
"for i in np.arange(0, LakArr_lyr1.shape[0]):\n",
|
||
" for j in np.arange(0, LakArr_lyr1.shape[1]):\n",
|
||
" im1 = i - 1\n",
|
||
" ip1 = i + 1\n",
|
||
" jm1 = j - 1\n",
|
||
" jp1 = j + 1\n",
|
||
" \n",
|
||
" if(im1 >= 0):\n",
|
||
" if(LakArr_lyr1[i, j] == 1 and LakArr_lyr1[im1, j] == 0):\n",
|
||
" bdlknc_lyr1[im1, j] = 1\n",
|
||
" \n",
|
||
" if(ip1 < LakArr_lyr1.shape[0]):\n",
|
||
" if(LakArr_lyr1[i, j] == 1 and LakArr_lyr1[ip1, j] == 0):\n",
|
||
" bdlknc_lyr1[ip1, j] = 1\n",
|
||
" \n",
|
||
" if(jm1 >= 0):\n",
|
||
" if(LakArr_lyr1[i, j] == 1 and LakArr_lyr1[i, jm1] == 0):\n",
|
||
" bdlknc_lyr1[i, jm1] = 1\n",
|
||
" \n",
|
||
" if(jp1 < LakArr_lyr1.shape[1]):\n",
|
||
" if(LakArr_lyr1[i, j] == 1 and LakArr_lyr1[i, jp1] == 0):\n",
|
||
" bdlknc_lyr1[i, jp1] = 1\n",
|
||
"\n",
|
||
"\n",
|
||
"bdlknc = [bdlknc_lyr1, bdlknc_lyr2, bdlknc_lyr3]\n",
|
||
"bdlknc = np.array(bdlknc)\n",
|
||
"\n",
|
||
"flux_data = {0: [[0.0073, 0.0073, 0.0, 0.0]],\n",
|
||
" 1: [[0.0073, 0.0073, 0.0, 0.0]]}\n",
|
||
"\n",
|
||
"lak = flopy.modflow.ModflowLak(mf, nlakes=nlakes, ipakcb=ipakcb, theta=theta,\n",
|
||
" nssitr=nssitr, sscncr=sscncr, surfdep=surfdep,\n",
|
||
" stages=stages, lakarr=LakArr, bdlknc=bdlknc,\n",
|
||
" flux_data=flux_data, unit_number=16)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Instantiate gage package for use with MODFLOW-NWT package"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 12,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"gages = [[1, 225, 90, 3],\n",
|
||
" [2, 68, 91, 3],\n",
|
||
" [3, 33, 92, 3],\n",
|
||
" [4, 165, 93, 3],\n",
|
||
" [5, 123, 94, 3],\n",
|
||
" [6, 77, 95, 3],\n",
|
||
" [7, 173, 96, 3],\n",
|
||
" [8, 328, 97, 3],\n",
|
||
" [9, 115, 98, 3],\n",
|
||
" [-1, -101, 1]]\n",
|
||
"\n",
|
||
"# gages = [[1,38,61,1],[2,67,62,1], [3,176,63,1], [4,152,64,1], [5,186,65,1], [6,31,66,1]]\n",
|
||
"files = ['no3.gag','seg1_gag.out','seg2_gag.out','seg3_gag.out','seg4_gag.out','seg5_gag.out','seg6_gag.out',\n",
|
||
" 'seg7_gag.out','seg8_gag.out','seg9_gag.out','lak1_gag.out']\n",
|
||
"\n",
|
||
"numgage = len(gages)\n",
|
||
"gage = flopy.modflow.ModflowGage(mf, numgage=numgage, gage_data=gages, filenames = files)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Instantiate Unsaturated-Zone Flow (UZF) package for MODFLOW-NWT"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 13,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"nuztop = 2\n",
|
||
"iuzfopt = 2\n",
|
||
"irunflg = 1\n",
|
||
"ietflg = 0\n",
|
||
"iuzfcb = 52\n",
|
||
"iuzfcb2 = 0\n",
|
||
"ntrail2 = 20\n",
|
||
"nsets2 = 20\n",
|
||
"nuzgag = 2\n",
|
||
"surfdep = 2.\n",
|
||
"\n",
|
||
"eps = 3.\n",
|
||
"thts = 0.30\n",
|
||
"thti = 0.13079\n",
|
||
"\n",
|
||
"fname_uzbnd = os.path.join('..','data','mt3d_example_sft_lkt_uzt','uzf_arrays','iuzbnd.txt')\n",
|
||
"fname_runbnd = os.path.join('..','data','mt3d_example_sft_lkt_uzt','uzf_arrays','irunbnd.txt')\n",
|
||
"\n",
|
||
"iuzfbnd = np.loadtxt(fname_uzbnd)\n",
|
||
"irunbnd = np.loadtxt(fname_runbnd)\n",
|
||
"\n",
|
||
"uzgag = [[106, 160, 121, 3],\n",
|
||
" [ 1, 1, -122, 1]]\n",
|
||
"\n",
|
||
"finf = {0: 1.8250e-03, 1: 1.8250e-03}\n",
|
||
"\n",
|
||
"uzf = flopy.modflow.ModflowUzf1(mf, nuztop=nuztop, iuzfopt=iuzfopt, irunflg=irunflg, ietflg=ietflg, ipakcb=iuzfcb,\n",
|
||
" iuzfcb2=iuzfcb2, ntrail2=ntrail2, nsets=nsets2, surfdep=surfdep, uzgag=uzgag,\n",
|
||
" iuzfbnd=1, irunbnd=0, vks=1.0E-6, eps=3.5, thts=0.35,\n",
|
||
" thtr=0.15, thti=0.20)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Instantiate Drain (DRN) package for MODFLOW-NWT"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 14,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"fname_drnElv = os.path.join('..','data','mt3d_example_sft_lkt_uzt','drn_arrays', 'elv.txt')\n",
|
||
"fname_drnCond = os.path.join('..','data','mt3d_example_sft_lkt_uzt','drn_arrays','cond.txt')\n",
|
||
"\n",
|
||
"drnElv = np.loadtxt(fname_drnElv)\n",
|
||
"drnCond = np.loadtxt(fname_drnCond)\n",
|
||
"\n",
|
||
"drnElv_lst = pd.DataFrame({\n",
|
||
" 'lay': 1,\n",
|
||
" 'row': np.nonzero(drnElv)[0] + 1,\n",
|
||
" 'col': np.nonzero(drnElv)[1] + 1,\n",
|
||
" 'elv': drnElv[np.nonzero(drnElv)],\n",
|
||
" 'cond': drnCond[np.nonzero(drnCond)] }, \n",
|
||
" columns=['lay', 'row', 'col', 'elv', 'cond'])\n",
|
||
"\n",
|
||
"# Convert the DataFrame into a list of lists for the drn constructor\n",
|
||
"stress_period_data = drnElv_lst.values.tolist()\n",
|
||
"\n",
|
||
"# Create a dictionary, 1 entry for each of the two stress periods.\n",
|
||
"stress_period_data = {0: stress_period_data,\n",
|
||
" 1: stress_period_data}\n",
|
||
"\n",
|
||
"drn = flopy.modflow.ModflowDrn(mf, ipakcb=ipakcb, stress_period_data=stress_period_data)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Instantiate linkage with mass transport routing (LMT) package for MODFLOW-NWT (generates linker file)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 15,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"lmt = flopy.modflow.ModflowLmt(mf, output_file_name='NO3.ftl', output_file_header='extended',\n",
|
||
" output_file_format='formatted', package_flows = ['all'])"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"# Now work on MT3D-USGS file creation"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 16,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"# Start by setting up MT3D-USGS class and pass in MODFLOW-NWT object for setting up a number of the BTN arrays\n",
|
||
"\n",
|
||
"mt = flopy.mt3d.Mt3dms(modflowmodel=mf, modelname=modelname, model_ws=modelpth,\n",
|
||
" version='mt3d-usgs', namefile_ext='mtnam', exe_name=mtexe,\n",
|
||
" ftlfilename='no3.ftl', ftlfree=True)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Instantiate basic transport (BTN) package for MT3D-USGS"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 17,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"ncomp = 1\n",
|
||
"lunit = 'FT'\n",
|
||
"sconc = 0.0\n",
|
||
"prsity = 0.3\n",
|
||
"cinact = -1.0\n",
|
||
"thkmin = 0.000001\n",
|
||
"nprs = -2\n",
|
||
"nprobs = 10\n",
|
||
"nprmas = 10\n",
|
||
"dt0 = 0.1\n",
|
||
"nstp = 1\n",
|
||
"mxstrn = 500\n",
|
||
"ttsmult = 1.2\n",
|
||
"ttsmax = 100\n",
|
||
"\n",
|
||
"# These observations need to be entered with 0-based indexing\n",
|
||
"obs = [[0, 104, 158],\n",
|
||
" [1, 104, 158],\n",
|
||
" [2, 104, 158]]\n",
|
||
"\n",
|
||
"btn = flopy.mt3d.Mt3dBtn(mt, MFStyleArr=True, DRYCell=True, lunit=lunit,\n",
|
||
" sconc=sconc, ncomp=ncomp, prsity=prsity, cinact=cinact,\n",
|
||
" obs=obs, thkmin=thkmin, nprs=nprs, nprobs=nprobs, \n",
|
||
" chkmas=True,nprmas=nprmas, dt0=dt0, nstp=nstp, \n",
|
||
" mxstrn=mxstrn, ttsmult=ttsmult, ttsmax=ttsmax)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Instantiate advection (ADV) package for MT3D-USGS"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 18,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"mixelm = 0\n",
|
||
"percel = 1.0000\n",
|
||
"mxpart = 5000\n",
|
||
"nadvfd = 1 # (1 = Upstream weighting)\n",
|
||
"\n",
|
||
"adv = flopy.mt3d.Mt3dAdv(mt, mixelm=mixelm, percel=percel, mxpart=mxpart, \n",
|
||
" nadvfd=nadvfd)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Instantiate generalized conjugate gradient solver (GCG) package for MT3D-USGS"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 19,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"mxiter = 1\n",
|
||
"iter1 = 50\n",
|
||
"isolve = 3\n",
|
||
"ncrs = 0\n",
|
||
"accl = 1.000000\n",
|
||
"cclose = 1.00e-06\n",
|
||
"iprgcg = 5\n",
|
||
"\n",
|
||
"gcg = flopy.mt3d.Mt3dGcg(mt, mxiter=mxiter, iter1=iter1, \n",
|
||
" isolve=isolve, ncrs=ncrs, accl=accl, \n",
|
||
" cclose=cclose, iprgcg=iprgcg)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Instantiate dispersion (DSP) package for MT3D-USGS"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 20,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"al = 0.1 # longitudinal dispersivity\n",
|
||
"trpt = 0.1 # ratio of the horizontal transverse dispersivity to 'AL'\n",
|
||
"trpv = 0.1 # ratio of the vertical transverse dispersitvity to 'AL'\n",
|
||
"dmcoef = 1.0000e-10\n",
|
||
"\n",
|
||
"dsp = flopy.mt3d.Mt3dDsp(mt, al=al, trpt=trpt, trpv=trpv, dmcoef=dmcoef,\n",
|
||
" multiDiff=True)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Instantiate source-sink mixing (SSM) package for MT3D-USGS"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 21,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"# no user-specified concentrations associated with boundary conditions\n",
|
||
"\n",
|
||
"mxss = 11199\n",
|
||
"\n",
|
||
"ssm = flopy.mt3d.Mt3dSsm(mt, mxss=mxss)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Instantiate reaction (RCT) package for MT3D-USGS"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 22,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"isothm = 0\n",
|
||
"ireact = 1\n",
|
||
"irctop = 2\n",
|
||
"igetsc = 0\n",
|
||
"ireaction = 0\n",
|
||
"\n",
|
||
"rc1 = 6.3258e-04 # first-order reaction rate for the dissolved phase\n",
|
||
"rc2 = 0.0 # Decay on Soil Layer\n",
|
||
"\n",
|
||
"rct = flopy.mt3d.Mt3dRct(mt, isothm=isothm, ireact=ireact, igetsc=igetsc,\n",
|
||
" rc1=rc1, rc2=rc2)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Instantiate streamflow transport (SFT) package for MT3D-USGS"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 23,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"nsfinit = len(reach_data)\n",
|
||
"mxsfbc = len(reach_data)\n",
|
||
"icbcsf = 0\n",
|
||
"ioutobs = 92\n",
|
||
"isfsolv = 1\n",
|
||
"wimp = 0.5\n",
|
||
"wups = 1.0\n",
|
||
"cclosesf = 1.0E-6\n",
|
||
"mxitersf = 10\n",
|
||
"crntsf = 1.0\n",
|
||
"iprtxmd = 0\n",
|
||
"coldsf = 0\n",
|
||
"dispsf = 0\n",
|
||
"obs_sf = [225,293,326,491,614,691,864,1192,1307]\n",
|
||
"sf_stress_period_data = {0: [0, 0, 0],\n",
|
||
" 1: [0, 0, 0],\n",
|
||
" 2: [0, 0, 0]}\n",
|
||
"\n",
|
||
"gage_output = [None, None, 'no3.sftobs']\n",
|
||
"\n",
|
||
"sft = flopy.mt3d.Mt3dSft(mt, nsfinit=nsfinit, mxsfbc=mxsfbc, icbcsf=icbcsf,\n",
|
||
" ioutobs=ioutobs, isfsolv=isfsolv, wimp=wimp,\n",
|
||
" wups=wups, cclosesf=cclosesf, mxitersf=mxitersf,\n",
|
||
" crntsf=crntsf, iprtxmd=iprtxmd, coldsf=coldsf,\n",
|
||
" dispsf=dispsf, nobssf=len(obs_sf), obs_sf=obs_sf,\n",
|
||
" sf_stress_period_data = sf_stress_period_data,\n",
|
||
" filenames=gage_output)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Instantiate unsaturated-zone transport (UZT) package for MT3D-USGS"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 24,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"mxuzcon = np.count_nonzero(irunbnd)\n",
|
||
"icbcuz = 45\n",
|
||
"iet = 0\n",
|
||
"wc = np.ones((nlay, nrow, ncol)) * 0.29\n",
|
||
"sdh = np.ones((nlay, nrow, ncol))\n",
|
||
"\n",
|
||
"uzt = flopy.mt3d.Mt3dUzt(mt, mxuzcon=mxuzcon, icbcuz=icbcuz, iet=iet,\n",
|
||
" iuzfbnd=iuzfbnd, sdh=sdh, cuzinf=1.4158e-03,\n",
|
||
" filenames='no3')"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"### Instantiate lake transport (LKT) package for MT3D-USGS"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 25,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": [
|
||
"nlkinit = 1\n",
|
||
"mxlkbc = 720\n",
|
||
"icbclk = 81\n",
|
||
"ietlak = 1\n",
|
||
"coldlak = 1\n",
|
||
"\n",
|
||
"lkt_flux_data = {0: [[0, 1, 0.01667]],\n",
|
||
" 1: [[0, 1, 0.02667]]}\n",
|
||
"\n",
|
||
"lkt = flopy.mt3d.Mt3dLkt(mt, nlkinit=nlkinit, mxlkbc=mxlkbc, icbclk=icbclk,\n",
|
||
" ietlak=ietlak, coldlak=coldlak, \n",
|
||
" lk_stress_period_data=lkt_flux_data)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"#### Write the MT3D-USGS input files for inspecting and running"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 26,
|
||
"metadata": {},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"/Users/jdhughes/Documents/Development/flopy_git/flopy_fork/examples/Notebooks\n"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"pth = os.getcwd()\n",
|
||
"print(pth)\n",
|
||
"\n",
|
||
"mf.write_input()\n",
|
||
"mt.write_input()\n",
|
||
"\n",
|
||
"#mf.run_model()\n",
|
||
"#mt.run_model()"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"metadata": {},
|
||
"outputs": [],
|
||
"source": []
|
||
}
|
||
],
|
||
"metadata": {
|
||
"kernelspec": {
|
||
"display_name": "Python 3",
|
||
"language": "python",
|
||
"name": "python3"
|
||
},
|
||
"language_info": {
|
||
"codemirror_mode": {
|
||
"name": "ipython",
|
||
"version": 3
|
||
},
|
||
"file_extension": ".py",
|
||
"mimetype": "text/x-python",
|
||
"name": "python",
|
||
"nbconvert_exporter": "python",
|
||
"pygments_lexer": "ipython3",
|
||
"version": "3.8.6"
|
||
}
|
||
},
|
||
"nbformat": 4,
|
||
"nbformat_minor": 2
|
||
}
|