How to do a massive Modflow 6 transient model updated from a spreadsheet (*.xlsx) file - Tutorial

Setting up information for several stress periods on normal GUIs tools can be an exhausting since it implies opening/formatting/filling/closing many charts manually. This task becomes more complicated if the model needs to be updated every given time. We wanted to develop a case where a steady state model can be updated to a transient model from data of a spreadsheet file. 

This tutorial cover the implementation of a transient model of a river/well/aquifer interaction where the number of stress periods, the length of stress periods, time steps, multiplier, steady/transient mode, pumping rate, river stage and output mode are updated from a excel (.xlsx) file. The tutorial is divided into 3 steps where each part shows sequential parts of this model transformation.

Tutorial

Code

First part: Inspect model and create model copy

import flopy, os

#import simulation
simPath = os.path.abspath('../out/model')
sim = flopy.mf6.MFSimulation.load('sim',sim_ws=simPath)

#check if model is transient or steady
sim.get_package('tdis')

#we cant say that is, but we are going to redefine the simulation steps, but first
#move the model to another folder because we dont want to mess the original data
sim.set_sim_path('../out/modelTrans/')
sim.write_simulation()

Second part: Modify Tdis and Sto

#import required packages
import flopy, os
import pandas as pd

#import simulation
simPath = os.path.abspath('../out/modelTrans')
exePath = '../exe/mf6.exe'
sim = flopy.mf6.MFSimulation.load('sim',sim_ws=simPath, exe_name=exePath)

#open model
model = sim.get_model('geo_model')

#set the temporal distribution 
timeDf = pd.read_excel('../in/xls/riverWellSchedule.xlsx', parse_dates=True, index_col=0)
#create timedelta column
timeDf['Timedelta'] = 0
timeDf.head()

#creating a timedelta list with one second for the first steady stress period
timeDeltaList = [1]
timeDeltaIndex = timeDf.iloc[1:].index - timeDf.iloc[:-1].index
timeDeltaList += list(timeDeltaIndex.total_seconds())
timeDeltaList[:5]

#add the time list in seconds to column
timeDf['Timedelta'] = timeDeltaList
timeDf.head()

#open tdis object
tdis = sim.get_package('tdis')

#show former perioddata
tdis.perioddata

#new number of stress periods
nper = timeDf.shape[0]
nper

#list of stress periods length, time steps and multiplier
tdisPeriodData = []
for index, row in timeDf.iterrows():
    tdisPeriodData.append([row.Timedelta, row['Time Steps'], row.Multiplier])
tdisPeriodData[:5]

#assign the number of stress periods and stress period data
tdis.nper= nper
tdis.perioddata = tdisPeriodData
tdis

#write simulation and check that the tdis has change
#sim.write_simulation()

#open sto object
sto = model.get_package('sto')

steadyDict = {}
transientDict = {}

i=0
for index, row in timeDf.iterrows():
    if row.Mode == 'Steady':
        steadyDict[i] = True
        transientDict[i] = False
    else:
        steadyDict[i] = False
        transientDict[i] = True
    #print preview
    i+=1
    if i == 5:
        print(steadyDict)
        print(transientDict)

#define the stress period types
sto.steady_state.set_data(steadyDict)
sto.transient.set_data(transientDict)

#write simulation and check that sto has changed correctly
sim.write_simulation()
#sim.run_simulation()

timeDf.to_excel('../In/xls/riverWellScheduleTimedelta.xlsx')

Third part: Modify boundary conditions (WEL and RIV) as setting output control

import flopy, os
import pandas as pd

#import simulation
simPath = os.path.abspath('../out/modelTrans')
exePath = '../exe/mf6.exe'
sim = flopy.mf6.MFSimulation.load('sim',sim_ws=simPath, exe_name=exePath)

#open model
model = sim.get_model('geo_model')

#open the temporal data 
timeDf = pd.read_excel('../in/xls/riverWellScheduleTimedelta.xlsx', parse_dates=True, index_col=0)
timeDf.head()

#working with the WEL package
wel = model.get_package('wel')
wel.stress_period_data.data

cellids = wel.stress_period_data.data[1].cellid

#extracting cellid for wells
wellDict={}
for index,cell in enumerate(cellids):
    wellDict[index]=cell[1]
wellDict

i = 0
#create empty stress period data dict
welSpd = {}
#looping over each stress period
for index, row in timeDf.iterrows():
    #for well 1, 2 and 3, with sign correction due to pumping
    wel1 = [(1,wellDict[0]), -row['Well 1']]
    wel2 = [(1,wellDict[1]), -row['Well 2']]
    wel3 = [(1,wellDict[2]), -row['Well 3']]
    welSpd[i] = [wel1, wel2, wel3]
    i+=1

    if i == 5:
        print(welSpd)
#assign to the wel package    
wel.stress_period_data = welSpd

#working with the RIVER package

import numpy as np

riv = model.get_package('riv')
tempRecarray = np.recarray.copy(riv.stress_period_data.data[0])
tempRecarray

i = 0
#create empty stress period data dict
rivSpd = {}

tempRecarray
for index, row in timeDf.iterrows():
    tempRecArray = np.recarray.copy(riv.stress_period_data.data[0])
    tempRecArray.stage = tempRecarray.rbot + row['River Stage']
    rivSpd[i] = tempRecArray
    i+=1

    if i==2:
        print(rivSpd)

riv.stress_period_data = rivSpd

#riv.stress_period_data.data

# working with the output control
oc = model.get_package('oc')

oc.saverecord.data

i = 0
#create empty stress period data dict
saveRecordSpd = {}
printRecordSpd = {}
#looping over each stress period
for index, row in timeDf.iterrows():
    if row['Output Control'] == 'Yes':
        saveRecordSpd[i] = [('HEAD', 'LAST'),('BUDGET', 'LAST')]
        printRecordSpd[i] = [('BUDGET', 'LAST')]
    else:
        saveRecordSpd[i] = []
        printRecordSpd[i] = []
    i+=1
#assign to the wel package    
oc.saverecord = saveRecordSpd
oc.printrecord = printRecordSpd

sim.write_simulation()
sim.run_simulation()

Input data

You can download the input data from this link.