Trouble running transient state transport model

[ougx]

Try to simulate injection to a confined aquifer with MF6 but failed. However it ran with MT3D with the same set up.

I modified autotest/test_gwt_adv04.py and created the model files.

Blow is the error and my python script to create the model (which includes model set up).

forrtl: error (65): floating invalid
Image              PC                Routine            Line        Source             
mf6.exe            00007FF69AE9664C  IMSLINEARMODULE_m        1826  ims8linear.f90
mf6.exe            00007FF69AE93533  IMSLINEARMODULE_m         957  ims8linear.f90
mf6.exe            00007FF69AEF7D63  NUMERICALSOLUTION        2445  NumericalSolution.f90
mf6.exe            00007FF69AEEDADA  NUMERICALSOLUTION        1525  NumericalSolution.f90
mf6.exe            00007FF69AEEA914  NUMERICALSOLUTION        1226  NumericalSolution.f90
mf6.exe            00007FF69AF14D82  SOLUTIONGROUPMODU          90  SolutionGroup.f90
mf6.exe            00007FF69AED50C2  MF6COREMODULE_mp_          36  mf6core.f90
mf6.exe            00007FF69AED4E22  MAIN__                      5  mf6.f90
mf6.exe            00007FF69B01F8EE  Unknown               Unknown  Unknown
mf6.exe            00007FF69B01FCD4  Unknown               Unknown  Unknown
KERNEL32.DLL       00007FFCEBD97034  Unknown               Unknown  Unknown
ntdll.dll          00007FFCECD1D241  Unknown               Unknown  Unknown

import flopy
import os
import numpy as np

def get_model(dir='.'):
    nlay, nrow, ncol = 1, 401, 401
    nper = 1
    perlen = [100.0]
    nstp = [30]
    tsmult = [1.1]
    steady = [False]
    delr = 20.
    delc = 20.
    botm = [0.]
    strt = 20.
    hnoflo = 1e30
    hdry = -1e30
    hk = 2834.8

    ss = 1.e-6
    sy = 0.1

    top = 20.
    laytyp = 0

    # put constant heads all around the box
    chdlist = \
        [(0, i, 0     , strt) for i in range(nrow)] + \
        [(0, i, ncol-1, strt) for i in range(nrow)] + \
        [(0, 0,      i, strt) for i in range(ncol)] + \
        [(0, nrow-1, i, strt) for i in range(ncol)] 
    chdspdict = {0: list(set(chdlist))}

    # injection well with rate (1000 GPM) and concentration of 100.
    q = 192528.0
    c = 100
    w = {0: [[(0, int(nrow / 2), int(ncol / 2)), q, c]]}

    nouter, ninner = 100, 300
    hclose, rclose, relax = 1e-6, 1e-6, 1.

    tdis_rc = []
    for i in range(nper):
        tdis_rc.append((perlen[i], nstp[i], tsmult[i]))

    name = 'SingleWell'

    # build MODFLOW 6 files
    ws = dir
    sim = flopy.mf6.MFSimulation(sim_name=name, version='mf6',
                                 exe_name='mf6',
                                 sim_ws=ws)
    # create tdis package
    tdis = flopy.mf6.ModflowTdis(sim, time_units='DAYS',
                                 nper=nper, perioddata=tdis_rc)

    # create gwf model
    gwfname = 'gwf_' + name
    gwf = flopy.mf6.MFModel(sim, model_type='gwf6', modelname=gwfname,
                            model_nam_file='{}.nam'.format(gwfname))

    # create iterative model solution and register the gwf model with it
    imsgwf = flopy.mf6.ModflowIms(sim, print_option='SUMMARY',
                                  outer_dvclose=hclose,
                                  outer_maximum=nouter,
                                  under_relaxation='NONE',
                                  inner_maximum=ninner,
                                  inner_dvclose=hclose, rcloserecord=rclose,
                                  linear_acceleration='CG',
                                  scaling_method='NONE',
                                  reordering_method='NONE',
                                  relaxation_factor=relax,
                                  filename='{}.ims'.format(gwfname))
    sim.register_ims_package(imsgwf, [gwf.name])

    dis = flopy.mf6.ModflowGwfdis(gwf, nlay=nlay, nrow=nrow, ncol=ncol,
                                  delr=delr, delc=delc,
                                  top=top, botm=botm,
                                  idomain=np.ones((nlay, nrow, ncol), dtype=int),
                                  filename='{}.dis'.format(gwfname))

    # initial conditions
    ic = flopy.mf6.ModflowGwfic(gwf, strt=strt,
                                filename='{}.ic'.format(gwfname))

    # node property flow
    npf = flopy.mf6.ModflowGwfnpf(gwf, save_flows=False,
                                  icelltype=laytyp,
                                  k=hk,
                                  k33=hk)
    # storage
    sto = flopy.mf6.ModflowGwfsto(gwf, save_flows=False,
                                  iconvert=laytyp,
                                  ss=ss, sy=sy,
                                  #steady_state={0: True, 2: True},
                                  transient={0: True})

    # chd files
    chd = flopy.mf6.ModflowGwfchd(gwf,
                                  stress_period_data=chdspdict,
                                  save_flows=False,
                                  pname='CHD-1')

    # wel files
    wel = flopy.mf6.ModflowGwfwel(gwf, print_input=True, print_flows=True,
                                  stress_period_data=w,
                                  save_flows=False,
                                  auxiliary='CONCENTRATION', pname='WEL-1')

    # output control
    oc = flopy.mf6.ModflowGwfoc(gwf,
                                budget_filerecord='{}.cbc'.format(gwfname),
                                head_filerecord='{}.hds'.format(gwfname),
                                headprintrecord=[
                                    ('COLUMNS', 10, 'WIDTH', 15,
                                     'DIGITS', 6, 'GENERAL')],
                                saverecord=[('HEAD', 'LAST')],
                                printrecord=[('HEAD', 'LAST'),
                                             ('BUDGET', 'LAST')])

    # create gwt model
    gwtname = 'gwt_' + name
    gwt = flopy.mf6.MFModel(sim, model_type='gwt6', modelname=gwtname,
                            model_nam_file='{}.nam'.format(gwtname))

    # create iterative model solution and register the gwt model with it
    imsgwt = flopy.mf6.ModflowIms(sim, print_option='SUMMARY',
                                  outer_dvclose=hclose,
                                  outer_maximum=nouter,
                                  under_relaxation='NONE',
                                  inner_maximum=ninner,
                                  inner_dvclose=hclose, rcloserecord=rclose,
                                  linear_acceleration='BICGSTAB',
                                  scaling_method='NONE',
                                  reordering_method='NONE',
                                  relaxation_factor=relax,
                                  filename='{}.ims'.format(gwtname))
    sim.register_ims_package(imsgwt, [gwt.name])

    dis = flopy.mf6.ModflowGwtdis(gwt, nlay=nlay, nrow=nrow, ncol=ncol,
                                  delr=delr, delc=delc,
                                  top=top, botm=botm,
                                  idomain=1,
                                  filename='{}.dis'.format(gwtname))

    # initial conditions
    ic = flopy.mf6.ModflowGwtic(gwt, strt=0.,
                                filename='{}.ic'.format(gwtname))

    # advection
    adv = flopy.mf6.ModflowGwtadv(gwt, scheme='upstream',
                                filename='{}.adv'.format(gwtname))

    # mass storage and transfer
    mst = flopy.mf6.ModflowGwtmst(gwt, porosity=0.1)

    # sources
    sourcerecarray = [('WEL-1', 'AUX', 'CONCENTRATION')]
    ssm = flopy.mf6.ModflowGwtssm(gwt, sources=sourcerecarray,
                                filename='{}.ssm'.format(gwtname))

    # output control
    oc = flopy.mf6.ModflowGwtoc(gwt,
                                budget_filerecord='{}.cbc'.format(gwtname),
                                concentration_filerecord='{}.ucn'.format(gwtname),
                                concentrationprintrecord=[
                                    ('COLUMNS', 10, 'WIDTH', 15,
                                     'DIGITS', 6, 'GENERAL')],
                                saverecord=[('CONCENTRATION', 'LAST')],
                                printrecord=[('CONCENTRATION', 'LAST'),
                                             ('BUDGET', 'LAST')])

    # GWF GWT exchange
    gwfgwt = flopy.mf6.ModflowGwfgwt(sim, exgtype='GWF6-GWT6',
                                     exgmnamea=gwfname, exgmnameb=gwtname,
                                     filename='{}.gwfgwt'.format(name))

    return sim

def main():
    # initialize testing framework
    sim = get_model()

    # build and run the models
    sim.write_simulation()
    sim.run_simulation()

[langevin-usgs]

We've seen this before for these simple types of problems (and your solver settings) where the starting concentration is zero and there is just a right-hand side addition of solute mass. In the first transport pass through the solver there is a divide by zero, which we need to handle better. The very simple answer is to just add a little mass to the initial conditions, such as

    # initial conditions
    cstrt = np.zeros((nlay, nrow, ncol))
    cstrt[0, int(nrow / 2), int(ncol / 2)] = 0.1
    ic = flopy.mf6.ModflowGwtic(gwt, strt=cstrt,
                                filename='{}.ic'.format(gwtname))

With this change the model seems to run fine and gives the answers you'd expect.

But we need to fix this for the next release, so let's leave this issue open.