examples/Notebooks/flopy3_multi-component_SSM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FloPy\n",
    "\n",
    "## Using FloPy to simplify the use of the MT3DMS ```SSM``` package\n",
    "\n",
    "A multi-component transport demonstration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-15T13:18:58.780064Z",
     "iopub.status.busy": "2022-12-15T13:18:58.779731Z",
     "iopub.status.idle": "2022-12-15T13:19:00.056753Z",
     "shell.execute_reply": "2022-12-15T13:19:00.052245Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.7.15 (default, Oct 18 2022, 09:36:49) \n",
      "[GCC 11.2.0]\n",
      "numpy version: 1.21.6\n",
      "flopy version: 3.3.6\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import sys\n",
    "from tempfile import TemporaryDirectory\n",
    "\n",
    "import numpy as np\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(\"flopy version: {}\".format(flopy.__version__))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we will create a simple model structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-15T13:19:00.108703Z",
     "iopub.status.busy": "2022-12-15T13:19:00.107635Z",
     "iopub.status.idle": "2022-12-15T13:19:00.114311Z",
     "shell.execute_reply": "2022-12-15T13:19:00.113321Z"
    }
   },
   "outputs": [],
   "source": [
    "nlay, nrow, ncol = 10, 10, 10\n",
    "perlen = np.zeros((10), dtype=float) + 10\n",
    "nper = len(perlen)\n",
    "\n",
    "ibound = np.ones((nlay, nrow, ncol), dtype=int)\n",
    "\n",
    "botm = np.arange(-1, -11, -1)\n",
    "top = 0.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create the ```MODFLOW``` packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-15T13:19:00.118768Z",
     "iopub.status.busy": "2022-12-15T13:19:00.118156Z",
     "iopub.status.idle": "2022-12-15T13:19:00.137099Z",
     "shell.execute_reply": "2022-12-15T13:19:00.136084Z"
    }
   },
   "outputs": [],
   "source": [
    "# temporary directory\n",
    "temp_dir = TemporaryDirectory()\n",
    "model_ws = temp_dir.name\n",
    "\n",
    "modelname = \"ssmex\"\n",
    "mf = flopy.modflow.Modflow(modelname, model_ws=model_ws)\n",
    "dis = flopy.modflow.ModflowDis(\n",
    "    mf,\n",
    "    nlay=nlay,\n",
    "    nrow=nrow,\n",
    "    ncol=ncol,\n",
    "    perlen=perlen,\n",
    "    nper=nper,\n",
    "    botm=botm,\n",
    "    top=top,\n",
    "    steady=False,\n",
    ")\n",
    "bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=top)\n",
    "lpf = flopy.modflow.ModflowLpf(mf, hk=100, vka=100, ss=0.00001, sy=0.1)\n",
    "oc = flopy.modflow.ModflowOc(mf)\n",
    "pcg = flopy.modflow.ModflowPcg(mf)\n",
    "rch = flopy.modflow.ModflowRch(mf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll track the cell locations for the ```SSM``` data using the ```MODFLOW``` boundary conditions.\n",
    "\n",
    "\n",
    "Get a dictionary (```dict```) that has the ```SSM``` ```itype``` for each of the boundary types."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-15T13:19:00.141283Z",
     "iopub.status.busy": "2022-12-15T13:19:00.140965Z",
     "iopub.status.idle": "2022-12-15T13:19:00.149011Z",
     "shell.execute_reply": "2022-12-15T13:19:00.148045Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'CHD': 1, 'BAS6': 1, 'PBC': 1, 'WEL': 2, 'DRN': 3, 'RIV': 4, 'GHB': 5, 'MAS': 15, 'CC': -1}\n",
      "[('k', '<i8'), ('i', '<i8'), ('j', '<i8'), ('css', '<f4'), ('itype', '<i8')]\n"
     ]
    }
   ],
   "source": [
    "itype = flopy.mt3d.Mt3dSsm.itype_dict()\n",
    "print(itype)\n",
    "print(flopy.mt3d.Mt3dSsm.get_default_dtype())\n",
    "ssm_data = {}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add a general head boundary (```ghb```).  The general head boundary head (```bhead```) is 0.1 for the first 5 stress periods with a component 1 (comp_1) concentration of 1.0 and a component 2 (comp_2) concentration of 100.0.  Then ```bhead``` is increased to 0.25 and comp_1 concentration is reduced to 0.5 and comp_2 concentration is increased to 200.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-15T13:19:00.153296Z",
     "iopub.status.busy": "2022-12-15T13:19:00.153010Z",
     "iopub.status.idle": "2022-12-15T13:19:00.165754Z",
     "shell.execute_reply": "2022-12-15T13:19:00.164734Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('k', '<i8'), ('i', '<i8'), ('j', '<i8'), ('bhead', '<f4'), ('cond', '<f4')]\n"
     ]
    }
   ],
   "source": [
    "ghb_data = {}\n",
    "print(flopy.modflow.ModflowGhb.get_default_dtype())\n",
    "ghb_data[0] = [(4, 4, 4, 0.1, 1.5)]\n",
    "ssm_data[0] = [(4, 4, 4, 1.0, itype[\"GHB\"], 1.0, 100.0)]\n",
    "ghb_data[5] = [(4, 4, 4, 0.25, 1.5)]\n",
    "ssm_data[5] = [(4, 4, 4, 0.5, itype[\"GHB\"], 0.5, 200.0)]\n",
    "\n",
    "for k in range(nlay):\n",
    "    for i in range(nrow):\n",
    "        ghb_data[0].append((k, i, 0, 0.0, 100.0))\n",
    "        ssm_data[0].append((k, i, 0, 0.0, itype[\"GHB\"], 0.0, 0.0))\n",
    "\n",
    "ghb_data[5] = [(4, 4, 4, 0.25, 1.5)]\n",
    "ssm_data[5] = [(4, 4, 4, 0.5, itype[\"GHB\"], 0.5, 200.0)]\n",
    "for k in range(nlay):\n",
    "    for i in range(nrow):\n",
    "        ghb_data[5].append((k, i, 0, -0.5, 100.0))\n",
    "        ssm_data[5].append((k, i, 0, 0.0, itype[\"GHB\"], 0.0, 0.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add an injection ```well```.  The injection rate (```flux```) is 10.0 with a comp_1 concentration of 10.0 and a comp_2 concentration of 0.0 for all stress periods.  WARNING: since we changed the ```SSM``` data in stress period 6, we need to add the well to the ssm_data for stress period 6."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-15T13:19:00.170130Z",
     "iopub.status.busy": "2022-12-15T13:19:00.169603Z",
     "iopub.status.idle": "2022-12-15T13:19:00.176066Z",
     "shell.execute_reply": "2022-12-15T13:19:00.174958Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('k', '<i8'), ('i', '<i8'), ('j', '<i8'), ('flux', '<f4')]\n"
     ]
    }
   ],
   "source": [
    "wel_data = {}\n",
    "print(flopy.modflow.ModflowWel.get_default_dtype())\n",
    "wel_data[0] = [(0, 4, 8, 10.0)]\n",
    "ssm_data[0].append((0, 4, 8, 10.0, itype[\"WEL\"], 10.0, 0.0))\n",
    "ssm_data[5].append((0, 4, 8, 10.0, itype[\"WEL\"], 10.0, 0.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add the ```GHB``` and ```WEL``` packages to the ```mf``` ```MODFLOW``` object instance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-15T13:19:00.181800Z",
     "iopub.status.busy": "2022-12-15T13:19:00.181243Z",
     "iopub.status.idle": "2022-12-15T13:19:00.186997Z",
     "shell.execute_reply": "2022-12-15T13:19:00.185877Z"
    }
   },
   "outputs": [],
   "source": [
    "ghb = flopy.modflow.ModflowGhb(mf, stress_period_data=ghb_data)\n",
    "wel = flopy.modflow.ModflowWel(mf, stress_period_data=wel_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create the ```MT3DMS``` packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-15T13:19:00.191108Z",
     "iopub.status.busy": "2022-12-15T13:19:00.190631Z",
     "iopub.status.idle": "2022-12-15T13:19:00.204311Z",
     "shell.execute_reply": "2022-12-15T13:19:00.203230Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "found 'rch' in modflow model, resetting crch to 0.0\n",
      "SSM: setting crch for component 2 to zero. kwarg name crch2\n"
     ]
    }
   ],
   "source": [
    "mt = flopy.mt3d.Mt3dms(modflowmodel=mf, modelname=modelname, model_ws=model_ws)\n",
    "btn = flopy.mt3d.Mt3dBtn(mt, sconc=0, ncomp=2, sconc2=50.0)\n",
    "adv = flopy.mt3d.Mt3dAdv(mt)\n",
    "ssm = flopy.mt3d.Mt3dSsm(mt, stress_period_data=ssm_data)\n",
    "gcg = flopy.mt3d.Mt3dGcg(mt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's verify that ```stress_period_data``` has the right ```dtype```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-15T13:19:00.208721Z",
     "iopub.status.busy": "2022-12-15T13:19:00.208217Z",
     "iopub.status.idle": "2022-12-15T13:19:00.212989Z",
     "shell.execute_reply": "2022-12-15T13:19:00.212104Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('k', '<i8'), ('i', '<i8'), ('j', '<i8'), ('css', '<f4'), ('itype', '<i8'), ('cssm(01)', '<f4'), ('cssm(02)', '<f4')]\n"
     ]
    }
   ],
   "source": [
    "print(ssm.stress_period_data.dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create the ```SEAWAT``` packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-15T13:19:00.216812Z",
     "iopub.status.busy": "2022-12-15T13:19:00.216261Z",
     "iopub.status.idle": "2022-12-15T13:19:00.223254Z",
     "shell.execute_reply": "2022-12-15T13:19:00.222082Z"
    }
   },
   "outputs": [],
   "source": [
    "swt = flopy.seawat.Seawat(\n",
    "    modflowmodel=mf,\n",
    "    mt3dmodel=mt,\n",
    "    modelname=modelname,\n",
    "    namefile_ext=\"nam_swt\",\n",
    "    model_ws=model_ws,\n",
    ")\n",
    "vdf = flopy.seawat.SeawatVdf(swt, mtdnconc=0, iwtable=0, indense=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-15T13:19:00.228511Z",
     "iopub.status.busy": "2022-12-15T13:19:00.227626Z",
     "iopub.status.idle": "2022-12-15T13:19:00.268925Z",
     "shell.execute_reply": "2022-12-15T13:19:00.267799Z"
    }
   },
   "outputs": [],
   "source": [
    "mf.write_input()\n",
    "mt.write_input()\n",
    "swt.write_input()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally, modify the ```vdf``` package to fix ```indense```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-15T13:19:00.273306Z",
     "iopub.status.busy": "2022-12-15T13:19:00.272809Z",
     "iopub.status.idle": "2022-12-15T13:19:00.279428Z",
     "shell.execute_reply": "2022-12-15T13:19:00.278550Z"
    }
   },
   "outputs": [],
   "source": [
    "fname = modelname + \".vdf\"\n",
    "f = open(os.path.join(model_ws, fname), \"r\")\n",
    "lines = f.readlines()\n",
    "f.close()\n",
    "f = open(os.path.join(model_ws, fname), \"w\")\n",
    "for line in lines:\n",
    "    f.write(line)\n",
    "for kper in range(nper):\n",
    "    f.write(\"-1\\n\")\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-15T13:19:00.283064Z",
     "iopub.status.busy": "2022-12-15T13:19:00.282594Z",
     "iopub.status.idle": "2022-12-15T13:19:00.287596Z",
     "shell.execute_reply": "2022-12-15T13:19:00.286687Z"
    }
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    # ignore PermissionError on Windows\n",
    "    temp_dir.cleanup()\n",
    "except:\n",
    "    pass"
   ]
  }
 ],
 "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.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}