{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n# 1D Solute Transport Benchmarks\n\nThis example is taken from the MODFLOW6 Examples, number 35.\n\nAs explained there, the setup is a simple 1d homogeneous aquifer with a steady\nstate flow field of constant velocity. The benchmark consists of four transport\nproblems that are modeled using this flow field. Here we have modeled these\nfour transport problems as a single simulation with multiple species. In all\ncases the initial concentration in the domain is zero, but water entering the\ndomain has a concentration of one:\n\n* species a is transported with zero diffusion or dispersion and the\n concentration distribution should show a sharp front, but due to the\n numerical method we see some smearing, which is expected.\n* species b has a sizeable dispersivity and hence shows more smearing than\n species a but the same centre of mass.\n* Species c has linear sorption and therefore the concentration doesn't enter\n the domain as far as species a or species b, but the front of the solute\n plume has the same overall shape as for species a or species b.\n* Species d has linear sorption and first order decay, and this changes the\n shape of the front of the solute plume.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\nimport pandas as pd\nimport xarray as xr\n\nimport imod\n\n\ndef create_transport_model(flowmodel, speciesname, dispersivity, retardation, decay):\n \"\"\"\n Function to create a transport model, as we intend to create four similar\n transport models.\n\n Parameters\n ----------\n flowmodel: GroundwaterFlowModel\n speciesname: str\n dispersivity: float\n retardation: float\n decay: float\n\n Returns\n -------\n transportmodel: GroundwaterTransportModel\n \"\"\"\n\n rhobulk = 1150.0\n porosity = 0.25\n\n tpt_model = imod.mf6.GroundwaterTransportModel()\n tpt_model[\"ssm\"] = imod.mf6.SourceSinkMixing.from_flow_model(flowmodel, speciesname)\n tpt_model[\"adv\"] = imod.mf6.AdvectionUpstream()\n tpt_model[\"dsp\"] = imod.mf6.Dispersion(\n diffusion_coefficient=0.0,\n longitudinal_horizontal=dispersivity,\n transversal_horizontal1=0.0,\n xt3d_off=False,\n xt3d_rhs=False,\n )\n\n # Compute the sorption coefficient based on the desired retardation factor\n # and the bulk density. Because of this, the exact value of bulk density\n # does not matter for the solution.\n if retardation != 1.0:\n sorption = \"linear\"\n kd = (retardation - 1.0) * porosity / rhobulk\n else:\n sorption = None\n kd = 1.0\n\n tpt_model[\"mst\"] = imod.mf6.MobileStorageTransfer(\n porosity=porosity,\n decay=decay,\n decay_sorbed=decay,\n bulk_density=rhobulk,\n distcoef=kd,\n first_order_decay=True,\n sorption=sorption,\n )\n\n tpt_model[\"ic\"] = imod.mf6.InitialConditions(start=0.0)\n tpt_model[\"oc\"] = imod.mf6.OutputControl(\n save_concentration=\"all\", save_budget=\"last\"\n )\n tpt_model[\"dis\"] = flowmodel[\"dis\"]\n return tpt_model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create the spatial discretization.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "nlay = 1\nnrow = 2\nncol = 101\ndx = 10.0\nxmin = 0.0\nxmax = dx * ncol\nlayer = [1]\ny = [1.5, 0.5]\nx = np.arange(xmin, xmax, dx) + 0.5 * dx\n\ngrid_dims = (\"layer\", \"y\", \"x\")\ngrid_coords = {\"layer\": layer, \"y\": y, \"x\": x}\ngrid_shape = (nlay, nrow, ncol)\ngrid = xr.DataArray(np.ones(grid_shape, dtype=int), coords=grid_coords, dims=grid_dims)\nbottom = xr.full_like(grid, -1.0, dtype=float)\n\ngwf_model = imod.mf6.GroundwaterFlowModel()\ngwf_model[\"ic\"] = imod.mf6.InitialConditions(0.0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create the input for a constant head boundary and its associated concentration.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "constant_head = xr.full_like(grid, np.nan, dtype=float)\nconstant_head[..., 0] = 60.0\nconstant_head[..., 100] = 0.0\n\nconstant_conc = xr.full_like(grid, np.nan, dtype=float)\nconstant_conc[..., 0] = 1.0\nconstant_conc[..., 100] = 0.0\nconstant_conc = constant_conc.expand_dims(\n species=[\"species_a\", \"species_b\", \"species_c\", \"species_d\"]\n)\n\ngwf_model[\"chd\"] = imod.mf6.ConstantHead(constant_head, constant_conc)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Add other flow packages.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "gwf_model[\"npf\"] = imod.mf6.NodePropertyFlow(\n icelltype=1,\n k=xr.full_like(grid, 1.0, dtype=float),\n variable_vertical_conductance=True,\n dewatered=True,\n perched=True,\n)\ngwf_model[\"dis\"] = imod.mf6.StructuredDiscretization(\n top=0.0,\n bottom=bottom,\n idomain=grid,\n)\ngwf_model[\"oc\"] = imod.mf6.OutputControl(save_head=\"all\", save_budget=\"all\")\ngwf_model[\"sto\"] = imod.mf6.SpecificStorage(\n specific_storage=1.0e-5,\n specific_yield=0.15,\n transient=False,\n convertible=0,\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create the simulation.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "simulation = imod.mf6.Modflow6Simulation(\"1d_tpt_benchmark\")\nsimulation[\"flow\"] = gwf_model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Add four transport simulations, and setup the solver flow and transport.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "simulation[\"tpt_a\"] = create_transport_model(gwf_model, \"species_a\", 0.0, 1.0, 0.0)\nsimulation[\"tpt_b\"] = create_transport_model(gwf_model, \"species_b\", 10.0, 1.0, 0.0)\nsimulation[\"tpt_c\"] = create_transport_model(gwf_model, \"species_c\", 10.0, 5.0, 0.0)\nsimulation[\"tpt_d\"] = create_transport_model(gwf_model, \"species_d\", 10.0, 5.0, 0.002)\n\nsimulation[\"flow_solver\"] = imod.mf6.Solution(\n modelnames=[\"flow\"],\n print_option=\"summary\",\n csv_output=False,\n no_ptc=True,\n outer_dvclose=1.0e-4,\n outer_maximum=500,\n under_relaxation=None,\n inner_dvclose=1.0e-4,\n inner_rclose=0.001,\n inner_maximum=100,\n linear_acceleration=\"bicgstab\",\n scaling_method=None,\n reordering_method=None,\n relaxation_factor=0.97,\n)\nsimulation[\"transport_solver\"] = imod.mf6.Solution(\n modelnames=[\"tpt_a\", \"tpt_b\", \"tpt_c\", \"tpt_d\"],\n print_option=\"summary\",\n csv_output=False,\n no_ptc=True,\n outer_dvclose=1.0e-4,\n outer_maximum=500,\n under_relaxation=None,\n inner_dvclose=1.0e-4,\n inner_rclose=0.001,\n inner_maximum=100,\n linear_acceleration=\"bicgstab\",\n scaling_method=None,\n reordering_method=None,\n relaxation_factor=0.97,\n)\n\nduration = pd.to_timedelta(\"2000d\")\nstart = pd.to_datetime(\"2000-01-01\")\nsimulation.create_time_discretization(additional_times=[start, start + duration])\nsimulation[\"time_discretization\"][\"n_timesteps\"] = 100" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Run the simulation.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "modeldir = imod.util.temporary_directory()\nsimulation.write(modeldir, binary=False)\nsimulation.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Open the concentration results and store them in a single DataArray.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "concentrations = []\nfor species in [\"a\", \"b\", \"c\", \"d\"]:\n conc = imod.mf6.out.open_conc(\n modeldir / f\"tpt_{species}/tpt_{species}.ucn\",\n modeldir / \"flow/dis.dis.grb\",\n ).assign_coords(species=species)\n concentrations.append(conc)\nconcentration = xr.concat(concentrations, dim=\"species\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Visualize the last concentration profiles of the model run for the different\nspecies.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "concentration.isel(time=-1, y=0).plot(x=\"x\", hue=\"species\")" ] } ], "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.10.11" } }, "nbformat": 4, "nbformat_minor": 0 }