FloPy - Lake Stage Observation starting at top of model

[patjhayes]

Hi,

I'm running a simple lake model where I'm trying to simulate vertical infiltration only through a storage device into a shallow groundwater model with a fixed groundwater level.

I'm running this in flopy as I'm interested in running many simulations with many variables. I've managed to run these types of simulations ok with ModelMuse.

However, with my FloPy code, I get the stage relationship in the lake starting at the top of my model despite setting initial lake stage at the bottom.

#!/usr/bin/env python
# coding: utf-8

# In[1]:


import os
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import flopy as fp


# In[2]:


# Define lake properties

Stormtrap_Width = [2,3]
Stormtrap_Length = [7,8]
path = os.getcwd()


# In[3]:


for Width in Stormtrap_Width:
    for Length in Stormtrap_Length:
        # Update the model name and workspace based on the current Stormtrap_Length
        modelname = f'ST_len_{Length}_wid_{Width}'
        modelws = f'./{modelname}'

            #grid properties
        Lx = 50
        Ly = 50
        ncol = 100
        nrow = 100
        delc = Lx / ncol
        delr = Ly / nrow

        #Layers
        nlay = 3
        top = 3
        botm_stormtrap = 2
        watertable = 1.0
        base_aquifer = -10


        botm = [botm_stormtrap,watertable,base_aquifer]

        k = 1.16E-4


        #Aquifer properties
        k_hor = [k,k,k]
        k_ver = [k,k,k]
        porosity = 0.2

        #Boundary conditions

        recharge = 0.0
        recharge_lake = 1e-5
        recharge_lake1 = 1e-99
        lak_etrate = 0


        L = Length*2
        W = Width*2

        lake_col_start = (ncol//2) -L//2
        lake_col_end = lake_col_start + L
        lake_row_start = (nrow//2) - W//2
        lake_row_end = lake_row_start + W
        lake_top_elev = top  # Example
        lake_bot_elev = botm_stormtrap  # Example

        # Define lake properties
        lake_id = 0
        lake_stage = lake_bot_elev
        lake_row, lake_col = nrow // 2, ncol // 2  # Center of the model

        lake_cells_amount = (lake_col_end - lake_col_start) * (lake_row_end - lake_row_start)

        # Lake package data
        packagedata = [(lake_id, lake_stage, lake_cells_amount, "LAKE")]
        leakance = 1

        runoff_rate =  0.000001
        runoff_rate1 = 0

        perioddata = [[0, "runoff", runoff_rate]]

        connectiondata = []

        connection_idx = 0
        for r in range(lake_row_start, lake_row_end):
            for c in range(lake_col_start, lake_col_end):
                conn = (lake_id, connection_idx, (0, r, c), "VERTICAL", leakance, lake_bot_elev, lake_top_elev, 0.5, 0.5)
                connectiondata.append(conn)
                connection_idx += 1


        sim = fp.mf6.MFSimulation(sim_name=modelname,
                             version='mf6.exe',
                             exe_name=path + '/mf6.exe',
                             sim_ws=modelws,
                             )


        ims = fp.mf6.ModflowIms(sim,
                            complexity = 'MODERATE',
                           linear_acceleration='bicgstab',
                           relaxation_factor = 0.97,
                            outer_maximum = 100, 
                            outer_dvclose = 0.01,
                            rcloserecord=0.1, 
                            inner_maximum=30,
                           inner_dvclose=1.e-6,
                           )

        gwf = fp.mf6.ModflowGwf(sim,
                           modelname=modelname,
                           model_nam_file=f"{modelname}.nam",
                           save_flows=True)

        dis = fp.mf6.ModflowGwfdis(gwf,
                              length_units='METERS',
                              nlay=nlay,
                              nrow=nrow,
                              ncol=ncol,
                              delr=delr,
                              delc=delc,
                              top=top,
                              botm=botm,
                              )


        npf = fp.mf6.ModflowGwfnpf(gwf,
                              save_specific_discharge=True,
                              icelltype=[1,0,0],
                              k=k_hor,
                              k33=k_ver,
                              )

        ic = fp.mf6.ModflowGwfic(gwf,
                            strt=watertable
                            )

        tdis_rc = [(1.0,1,1.0)] #length, number of steps and timestep
        tdis = fp.mf6.ModflowTdis(sim,
                                     time_units='SECONDS',
                                     perioddata=tdis_rc,
                                     )

        chd_spd = []
        for r in range(nrow):
                chd_spd.append([(2, r, 0), watertable])

        # Right boundary
        for r in range(nrow):
                chd_spd.append([(2, r, ncol - 1), watertable])

        # Top boundary
        for c in range(1, ncol - 1):
                chd_spd.append([(2, 0, c), watertable])

        # Bottom boundary
        for c in range(1, ncol - 1):
                chd_spd.append([(2, nrow - 1, c), watertable])

        chd = fp.mf6.ModflowGwfchd(gwf,stress_period_data=chd_spd,)

        rch = fp.mf6.ModflowGwfrcha(gwf,
                               recharge = recharge, 
                               pname = 'rch'
                               )

            # Create LAK package
        lak = fp.mf6.ModflowGwflak(gwf, pname='lake', boundnames=True, time_conversion=86400.0,
                                     nlakes=1, noutlets=0, 
                                      packagedata=packagedata, 
                                      connectiondata=connectiondata,
                                 perioddata=perioddata)

        sto = fp.mf6.ModflowGwfsto(
            gwf, 
            ss=1e-5,
            transient=True,)



        tdis_rc = [(120, 120, 1)]
        tids = fp.mf6.ModflowTdis(sim,
                                     time_units='SECONDS',
                                     nper=1,
                                     perioddata=tdis_rc,
                                     )

        obs_file = f"{modelname}.lak.obs"

        oc = fp.mf6.ModflowGwfoc(gwf,
                            budget_filerecord=f"{modelname}.cbc",
                            head_filerecord=f"{modelname}.hds",
                            saverecord=[("HEAD","ALL"), ("BUDGET","ALL")])



        csv_file = obs_file + ".csv"

        # Define additional observation points outside the lake
        additional_obs = [
            ["head_obs1", "stage", (0, 60, 60)],  # Replace with actual cell location
            # Add as many observations as needed
        ]

        obslist = [["conductance", "conductance", (0,50,50)], 
              ["stage", "stage", (0,50,50)],
              ["Runoff", "Runoff", (0,50,50)],
              ["volume", "volume", (0,50,50)], 
              ["WITHDRAWAL", "WITHDRAWAL", (0,50,50)], 
              ["LAK", "LAK", (0,50,50)], 
              ["storage", "storage", (0,50,50)],
              ["surface-area", "surface-area", (0,50,50)],
              ["inflow", "inflow", (0,50,50)]]

        
        obslist = obslist + additional_obs
        
        obsdict = {csv_file: obslist}

        # Define the observation

        lakobs = fp.mf6.ModflowUtlobs(lak, filename=obs_file, print_input=False, continuous=obsdict)

        sim.write_simulation()
        sim.run_simulation();

        # Ensure the directory exists, if not, create it
        if not os.path.exists(modelws):
            os.mkdir(modelws)


plt.rcParams['figure.figsize']=(5,5)
plt.rcParams["figure.autolayout"] = True


# 
# Check input

# In[4]:


import flopy.utils.binaryfile as bf

# Read in the head data
head_file = os.path.join(modelws, f"{modelname}.hds")
headobj = bf.HeadFile(head_file)
head = headobj.get_data(totim=headobj.get_times()[-1])  # Get the head for the last time step

# Check the range of your head data
head_min = head.min()
head_max = head.max()



# In[8]:


print(head_max)
print(head_min)


# In[5]:


# Create the plot
pmv = fp.plot.PlotMapView(model=gwf, layer=1)
pmv.plot_bc('CHD', color='b')  # Plot CHD boundaries
pmv.plot_bc('LAK', plotAll=True, color='y')  # Plot lake boundaries
pmv.plot_grid(colors='silver', lw=1)  # Plot the grid

# Plot contours of hydraulic head
levels = np.linspace(head.min(), head.max(), num=21)  # Adjust the number of levels and range as needed
contours = pmv.contour_array(head, levels=levels, colors='blue', linewidths=1)  # Add contours

# Optionally, add a color bar
plt.colorbar(contours, shrink=0.0001, aspect=5)

plt.show()  # Show the plot


# In[ ]:


plt.figure(figsize=(5,3))
pxs = fp.plot.PlotCrossSection(model=gwf, line={"Row": nrow//2})
pxs.plot_bc('CHD')
#pxs.plot_bc('WEL')
pxs.plot_bc('LAK', color = 'y')
pxs.plot_grid(colors='silver',lw=1)


# Solve Model

# READ AND VISUALISE RESULTS

# In[ ]:


print(f'FINISHED {modelname}')

`

Stormtrap2.zip

[emorway-usgs]

@patjhayes, just a couple of quick things to try. When you instantiate Sto, you didn't specify specific yield but did specify that the model is transient. Also, icelltype looks to be zero or the bottom two layers which directs the model to hold them as saturated. Also, it doesn't look as though you specified idomain which means all cells will specified as active. However, it looks like you incised the infiltration pond (LAK) into the model, so you have both a lake and aquifer cells active in the same location.

[patjhayes]

@emorway-usgs Thanks Eric,

I've added specific yield, made icelltype 1 in all layers and yet the stage observation sill starts at 3.

Is it incorrect to incise the infiltration pond? I was using the flopy examples as a guide and it seems that they are done this way also?

[emorway-usgs]

No, it's not incorrect to incise an infiltration basin. As best I can tell from reading the script (I haven't run the script locally), a lake spans multiple cells, but it doesn't look like idomain has been deactivated where the lake cells exist? Just something to check. Maybe post the .lak input file to this discussion and maybe that'll help me understand how the lakes are being set up.

[patjhayes]

I don't have a seperate .lak input file. I've now tested setting the idomain for the lake cells to 0 but it doesn't like it giving an error:

MFDataException: An error occurred in data element "connectiondata" model "ST_len_7_wid_2" package "lake". The error occurred while setting data in the "__init__" method.
Additional Information:
(1) Cellid (0, 48, 43) is outside of the active model grid.
(2) Error occurred while adding dataset "connectiondata" to block "connectiondata"