modflowpy / flopy

A Python package to create, run, and post-process MODFLOW-based models.
https://flopy.readthedocs.io
Other
520 stars 314 forks source link

UZF Cells Package Data not processed in simulation #2277

Closed iaped closed 3 months ago

iaped commented 3 months ago

I have built my UZF file in flopy and it is in the exact format as it should be. However, when I run the simulation, in the list file it does not show it being completed. This additionally messes up the UZT as it cannot find the corresponding UZF file. I even reduced my UZF file to just 1 UZF cell, however it still wont process or even produce an error. Does anyone know how to fix this? I don't want to give up on the UZF package.

This is what the list file says:

UZF -- UZF CELLS PACKAGE, VERSION 8, 2/22/2014 INPUT READ FROM UNIT 1022

File generated by Flopy version 3.3.6 on 07/30/2024 at 20:37:58.

PROCESSING UZF CELLS OPTIONS END OF UZF CELLS OPTIONS

PROCESSING UZF CELLS DIMENSIONS NUZFCELLS = 1 NTRAILWAVES = 7 NTRAILSETS = 40 END OF UZF CELLS DIMENSIONS

PROCESSING UZF CELLS PACKAGEDATA STRT set from input file

STO -- STORAGE PACKAGE, VERSION 1, 5/19/2014 INPUT READ FROM UNIT 1011

emorway-usgs commented 3 months ago

Hello @iaped, would you be willing to share your script? There are a number of UZF/UZT and UZF/UZE test models - all built with FloPy scripts - so the work flow you describe has been tested. It would be helpful to see a bit more of the details.

iaped commented 3 months ago

Create the UZF File

UZF options

none, not modeling ET

Define UZF dimensions

nuzfcells = 2200 ntrailwaves= 7 nwavesets= 40

Define the UZF package data

landflag= 1 #indicates if the surface is land that can be infiltrated ivertcon=-1 #indicates the UZF cell is not connected to an underlying UZF cell surfdep=0.33 #surface depression of ~1 ft to allow for ponding uzf_thts = 0.41 # saturated water content (Qs) uzf_thtr = 0.057 # residual water content (Qr) uzf_thti = 0.2 # initial water content uzf_vks = 3.5 # saturated vertical hydraulic conductivity (m/d) eps= 4.38 #Brooks Corey Exponent from Dingman physical hydrology for loamy sand finf= 0.8 #general steady infiltration rate from Cornell

get source row and columns from df

uzf_packagedata = []

Iterate through each row and column

for i, row in sources.iterrows(): j, k = row['row']-1, row['column']-1 cellid = (0, j, k) uzf_packagedata.append([i, cellid, landflag, ivertcon, surfdep, uzf_vks, uzf_thtr, uzf_thts, uzf_thti, eps])

Define the UZF period data (e.g., infiltration rate)

uzf_perioddata = {} for iper in range(446): periodstep = [] for i, row in sources.iterrows(): j, k = row['row'], row['column'] cellid = (1, j, k) periodstep.append([i, finf, 0, 0, 0, 0, 0, 0]) uzf_perioddata[iper] = periodstep

uzf = flopy.mf6.ModflowGwfuzf( gwf, nuzfcells=nuzfcells, ntrailwaves=ntrailwaves, nwavesets=nwavesets, packagedata=uzf_packagedata, perioddata=uzf_perioddata, pname= 'VADOSE' )

iaped commented 3 months ago

This is how the UZF file is built in Flopy, and it does not throw an error in MODFLOW however MODFLOW does not say "END OF UZF PACKAGE DATA" so it is failing to read the package data in the input file. I have manipulated it to be indentical to the UZF file in the input/output documentation for mf6 but even this does not work

emorway-usgs commented 3 months ago

Hello @iaped, I can confirm that UZF was not writing a msg to the standard GWF listing file denoting when it finished processing the UZF input. I opened up a PR to fix this oversight.

When I asked if you would be willing to share your script, I was wondering if you'd be willing to share the script in its entirety? This will allow me to try generating the input locally and then running with MF6 to more easily identify what might be happening. It's difficult to tell from the snippet of script you sent what the trouble is. If you prefer not to share the script, would be willing to post the UZF input file?

iaped commented 3 months ago
#Import necessary Packages
import os
import matplotlib.pyplot as plt
import numpy as np
import flopy
import pandas as pd

Load Necessary Data files for model

#Load data files
dsp_dat = np.loadtxt('baseModelTran-disp.dat') #dispersivities
por_dat=  np.loadtxt('baseModelTran-porosity-mid.dat') #porosity
kd_dat= np.loadtxt('Kd-pfos.dat') # Kd for PFOS
rb_dat= np.loadtxt('baseModelTran-rho_b.dat') # bulk density
sr_dat= np.loadtxt('baseModelTran-sr.dat') #immobile porosity
#Load vadose zone concentration file
conc_in = pd.read_excel('input_concentrations.xlsx', skiprows=[1,2])
conc_in.rename(columns={'Time.1':'time', 'PFHxS.1':'PFHxS', 'PFOA.1':'PFOA', 'PFOS.1':'PFOS'}, inplace=True)
#Load in source cells
sources = pd.read_excel('Sources.xlsx', sheet_name='Source Cells')
#Create the Groundwater Flow Model
sim = flopy.mf6.MFSimulation.load(sim_name='baseModel', sim_ws='.')
gwf = sim.get_model()
nlay, nrow, ncol = 1, 972, 700
delr = delc = 20.0
loading simulation...
  loading simulation name file...
  loading tdis package...
  loading model gwf6...
    loading package dis...
    loading package npf...
    loading package ic...
    loading package oc...
    loading package obs...
    loading package chd...
    loading package rch...
    loading package riv...
    loading package riv...
    loading package sto...
    loading package wel...
  loading ims package baseflowmod...
#Create the UZF File

#UZF options
#none, not modeling ET

#Define UZF dimensions
nuzfcells = 2200
ntrailwaves= 7
nwavesets= 40

#Define the UZF package data 
landflag= 1 #indicates if the surface is land that can be infiltrated
ivertcon=-1 #indicates the UZF cell is not connected to an underlying UZF cell
surfdep=0.33 #surface depression of ~1 ft to allow for ponding
uzf_thts = 0.41 # saturated water content (Qs)
uzf_thtr = 0.057 # residual water content (Qr)
uzf_thti = 0.2  # initial water content
uzf_vks = 3.5   # saturated vertical hydraulic conductivity (m/d)
eps= 4.38 #Brooks Corey Exponent from Dingman physical hydrology for loamy sand
finf= 0.8 #general steady infiltration rate from Cornell 

#get source row and columns from df
uzf_packagedata = []
# Iterate through each row and column
for i, row in sources.iterrows():
    j, k = row['row']-1, row['column']-1
    cellid = (0, j, k)
    uzf_packagedata.append([i, cellid, landflag, ivertcon, surfdep, uzf_vks, uzf_thtr, uzf_thts, uzf_thti, eps])

# Define the UZF period data (e.g., infiltration rate)
uzf_perioddata = {}
for iper in range(446):
    periodstep = []
    for i, row in sources.iterrows():
        j, k = row['row'], row['column']
        cellid = (1, j, k)

        periodstep.append([i, finf, 0, 0, 0, 0, 0, 0])

    uzf_perioddata[iper] = periodstep
uzf = flopy.mf6.ModflowGwfuzf(
    gwf,
    nuzfcells=nuzfcells,
    ntrailwaves=ntrailwaves,
    nwavesets=nwavesets,
    packagedata=uzf_packagedata,
    perioddata=uzf_perioddata,
    pname= 'VADOSE'
)

sim.write_simulation()# Transport Model

#create the tranport model
gwt= flopy.mf6.ModflowGwt(sim, modelname='baseModelTran')
#Create the DIS package
nlay, nrow, ncol = 1, 972, 700
delr, delc = 20.0, 20.0
top = "baseModelGrid-topElev.dat"
botm = "baseModelGrid-bedrockElev.dat"
idomain= "baseModelGrid-idomain.dat"
dis = flopy.mf6.ModflowGwtdis(gwt, length_units= 'meters', nlay=nlay, nrow=nrow, ncol=ncol, delr=delr, delc=delc, top=top, botm=botm, idomain=idomain)
#Create the IC package
strt = 0.0
ic = flopy.mf6.ModflowGwtic(gwt, strt=strt)
#Create the ADV package
adv = flopy.mf6.ModflowGwtadv(gwt, scheme='TVD')
#Create the DSP package
dsp = flopy.mf6.ModflowGwtdsp(
    gwt, 
    alh=dsp_dat, 
    ath1= dsp_dat * .1, 
    atv= dsp_dat*0.01
)
#Create the UZT Pacakge

#Define Package Data
strt= 0 #starting concentration
packagedata=[(i, strt) for i in range(2200)]

#Define period data
uztperioddata ={}
for i in range(446):
    perioddata = []
    for j in range(2200):
        perioddata.append([j, 'CONCENTRATION', 'PFOS']) 
    uztperioddata[i] = perioddata
#Create UZT File
#Manually write .ts FILEIN in the input file

uzt= flopy.mf6.ModflowGwtuzt(
    gwt, 
    # timeseries= 'input-conc.ts', 
    packagedata=packagedata, 
    uztperioddata= uztperioddata,
    pname= 'VADOSE'
)
#Create the SSM package
ssm = flopy.mf6.ModflowGwtssm(gwt, sources=[('UZT',)])
#Create MST Package
mst = flopy.mf6.ModflowGwtmst(
    gwt,
    sorption= 'LINEAR',
    porosity= por_dat,
    bulk_density= rb_dat,
    distcoef= kd_dat
)   
#Create IST Package
ist=flopy.mf6.ModflowGwtist(
    gwt,
    sorption= 'LINEAR',
    thetaim= sr_dat,
    zetaim = 0.001,
    bulk_density= rb_dat,
    distcoef= kd_dat
)
sim.write_simulation()
writing simulation...
  writing simulation name file...
  writing simulation tdis package...
  writing ims package baseflowmod...
  writing model baseflowmod...
    writing model name file...
    writing package dis...
    writing package npf...
    writing package ic...
    writing package oc...
    writing package obs_0...
    writing package chd_0...
    writing package rch_0...
    writing package riv_0...
    writing package riv_1...
    writing package sto...
    writing package wel_0...
    writing package vadose...
  writing model baseModelTran...
    writing model name file...
    writing package dis...
    writing package ic...
    writing package adv...
    writing package dsp...
    writing package vadose...
    writing package ssm...
    writing package mst...
    writing package ist...
emorway-usgs commented 3 months ago

@iaped I won't be able to run these commands and generate your model locally without all the additional dependencies.

However, I think I see at least one problem. When you instantiate UZT, you need to tell it the name of the corresponding flow package. That is, you need to specify FLOW_PACKAGE_NAME in the OPTIONS block of UZT. In your script, this would be FLOW_PACKAGE_NAME='vadose'. Please refer to the mf6io guide for more detail regarding this parameter.

iaped commented 3 months ago

Hi @emorway-usgs since both the UZF and UZT share the same package name the default is to look for a UZF file name with the same package name so in this case that is true. I have tried to define the flow package name but same issue occurs where modflow will not process the package data of the UZF file. All of the other packages run so the package I am struggling with is just the UZF file. In my script I reference a file name "Sources" this file just lists the row and column of the source cells. So in the UZF portion of this script it is just referencing the row and column coordinates to create the cellid to be used in the UZF package data lines. When I look at the UZF file in a text editor it appears to be the exact format required by the file and it does not flag an error, it will just stop processing the file.

emorway-usgs commented 3 months ago

@iaped, can you zip up the model and attach it? If the zip file is large, could you chop down the input written by Flopy to just 1 or two stress periods? It'll be easier to help if I can run locally.

iaped commented 3 months ago

@emorway-usgs I will be attaching 2 zip files because it is too large. Please note that the mf6 executable has been removed from both zip files and I reduced the UZF file to 1 stress period pt1.zip

iaped commented 3 months ago

pt2.zip

emorway-usgs commented 3 months ago

I was able to sort through a few issues; however, you have an issue related to your IST package that still needs to be sorted out. I'm posting a copy of the latest error message at the end of this post.

  1. You'll want to replace the SSM file with the attached, it contained errant information. May need to check up on your script so this file is written correctly.
  2. The mfsim.nam file was missing an EXCHANGE block. I'm attaching an example of one that I was able to get working
  3. MF6 no longer supports multiple models within a single IMS solution. You'll see this correction in the simulation name file and an additional IMS input file attached.
  4. The IST error message provides some further direction for you to find a solution.

fixed_files.zip

ist_err