Automatic option for moving the water based on DEM in MODFLOW 6

Question

Is there any automatic option in MODFLOW6-FloPy to route the rejected infiltration from the UZF cells to the SFR reaches based on the DEM and grid cells elevations?

Answer 1

Everything that follows assumes you are familiar with python and flopy.

I think the python below will help you out, I wrote it to essentially accomplish what you are asking for. You'll likely need to make a few edits to get inputs and outputs into the form your particular application requires.

It requires a couple of inputs:

• elev_arr: this is a 2D numpy array that contains the land surface elevations. For the model that I used this in, I had a saved file that had all of the land surface elevations called 'top1.txt'. I used this file directly in MODFLOW array. If, in your model, your land surface elevations occur in different layers, you'll need to account for that.
• ibnd: gives the extent of the active cells in the elevation array
• sfr_dat: is a list of tuples, where each tuple contains (lay, row, col). The order of this list matters! This list should be in the same order as the data appearing in the packagedata block of the SFR input file.
• sfrlayout_conn_candidate_elev: a value that is greater than all of your active cell's maximum elevation

The function determine_runoff_conns_4mvr generates and returns a 2D numpy array that is analogous to the IRUNBND array used in UZF1. I put this function into its own python script called mvr_conns and then imported it. The value at each position of the numpy array returned by determine_runoff_conns_4mvr corresponds to the SFR reach (identified by its index in the SFR PackageData block) that rejected infiltration will be routed to, which I think is what you want.

A note about lakes: if you have any, the script below will require some modification. I only had 1 lake and I dealt with it in a very non-robust manner (because I could get away with that). I have to leave this fix to you as an exercise if you have more than one lake.

So, for the MVR package, you'll need to process the array that is returned by determine_runoff_conns_4mvr. Here is how I did it, but you may need to alter for your purposes. Basically, the information contained in static_mvrperioddata can be passed to the perioddata argument when instantiating the MVR package with flopy.

# This function uses the top elevation array and SFR locations to calculate the irunbnd array from the UZF1 package.
import mvr_conns
# Of course in MF6, MVR now handles these runoff connections
irunbnd = mvr_conns.determine_runoff_conns_4mvr(pth)  # at this point, irunbnd is 1-based, compensate below

iuzno = 0
k     = 0             # Hard-wire the layer no.
for i in range(0, iuzfbnd.shape[0]):
    for j in range(0,iuzfbnd.shape[1]):
        if irunbnd[i,j] > 0:           # This is a uzf -> sfr connection
            static_mvrperioddata.append(['UZF-1', iuzno, 'SFR-1', irunbnd[i, j] - 1, 'FACTOR', 1.])
            iuzno += 1
        elif irunbnd[i, j] < 0:        # This is a uzf -> lak connection
            static_mvrperioddata.append(['UZF-1', iuzno, 'RES-1', abs(irunbnd[i, j]) - 1, 'FACTOR', 1.])
            iuzno += 1

Here is the base function: def determine_runoff_conns_4mvr(pth):

    elev_arr = np.loadtxt(os.path.join(pth,'dis_support','top1.txt'), dtype=np.float)
    ibnd     = np.loadtxt(os.path.join(pth,'bas_support','ibnd1.txt'), dtype=np.int)
    
    # Get the sfr information stored in a companion script
    sfr_dat  = sfrinfo.get_sfrlist()
    sfrlayout = np.zeros_like(ibnd)
    for i, rchx in enumerate(sfr_dat):
        row = rchx[1]
        col = rchx[2]
        sfrlayout[row - 1, col - 1] = i + 1
    
    sfrlayout_new = sfrlayout.copy()
    
    nrow = 64
    ncol = 133
    stop_candidate = False
    
    for i in np.arange(0, nrow):
        for j in np.arange(0, ncol):
        
            # Check to ensure current cell is active
            if ibnd[i, j] == 0:
                continue
            
            # Check to make sure it is not a stream cell
            if not sfrlayout[i, j] == 0:
                continue
            
            # Recursively trace path by steepest decent back to a stream
            curr_i = i
            curr_j = j
            
            sfrlayout_conn_candidate_elev = 10000.
            while True:
                direc = 0
                min_elev = elev_arr[curr_i, curr_j]
                
                # Look straight left
                if curr_j > 0:
                    if not sfrlayout[curr_i, curr_j - 1] == 0 and not ibnd[curr_i, curr_j - 1] == 0:   # Step in if neighbor is a stream cell
                        if elev_arr[curr_i, curr_j - 1] > 0 and (elev_arr[curr_i, curr_j - 1] < elev_arr[curr_i, curr_j] and
                                                                 elev_arr[curr_i, curr_j - 1] < sfrlayout_conn_candidate_elev):
                            sfrlayout_conn_candidate = sfrlayout[curr_i, curr_j - 1]
                            sfrlayout_conn_candidate_elev = elev_arr[curr_i, curr_j - 1]
                            stop_candidate = True
                    
                    elif not elev_arr[curr_i, curr_j - 1] == 0 and not ibnd[curr_i, curr_j - 1] == 0:  # Step here if neighbor is not an sfr cell
                        if elev_arr[curr_i, curr_j - 1] < elev_arr[curr_i, curr_j] and elev_arr[curr_i, curr_j - 1] < min_elev:
                            elevcm1 = elev_arr[curr_i, curr_j - 1]
                            min_elev = elevcm1
                            direc = 2
                
                # Look up and left
                if curr_j > 0 and curr_i > 0:
                    if not sfrlayout[curr_i - 1, curr_j - 1] == 0 and not ibnd[curr_i - 1, curr_j - 1] == 0:   # Step in if neighbor is a stream cell
                        if elev_arr[curr_i - 1, curr_j - 1] > 0 and (elev_arr[curr_i - 1, curr_j - 1] < elev_arr[curr_i, curr_j] and
                                                                     elev_arr[curr_i - 1, curr_j - 1] < sfrlayout_conn_candidate_elev):
                            sfrlayout_conn_candidate = sfrlayout[curr_i - 1, curr_j - 1]
                            sfrlayout_conn_candidate_elev = elev_arr[curr_i - 1, curr_j - 1]
                            stop_candidate = True

                    elif not elev_arr[curr_i - 1, curr_j - 1] == 0 and not ibnd[curr_i - 1, curr_j - 1] == 0:   # Step here if neighbor is not an sfr cell
                        if elev_arr[curr_i - 1, curr_j - 1] < elev_arr[curr_i, curr_j] and elev_arr[curr_i - 1, curr_j - 1] < min_elev:
                            elevrm1cm1 = elev_arr[curr_i - 1, curr_j - 1]
                            min_elev = elevrm1cm1
                            direc = 5

                
                # Look straight right
                if curr_j < ncol - 1:
                    if not sfrlayout[curr_i, curr_j + 1] == 0 and not ibnd[curr_i, curr_j + 1] == 0:   # Step in if neighbor is a stream cell
                        if elev_arr[curr_i, curr_j + 1] > 0 and (elev_arr[curr_i, curr_j + 1] < elev_arr[curr_i, curr_j] and
                                                                 elev_arr[curr_i, curr_j + 1] < sfrlayout_conn_candidate_elev):
                            sfrlayout_conn_candidate = sfrlayout[curr_i, curr_j + 1]
                            sfrlayout_conn_candidate_elev = elev_arr[curr_i, curr_j + 1]
                            stop_candidate = True
                    
                    elif not elev_arr[curr_i, curr_j + 1] == 0 and not ibnd[curr_i, curr_j + 1] == 0:  # Step here if neighbor is not an sfr cell
                        if elev_arr[curr_i, curr_j + 1] < elev_arr[curr_i, curr_j] and elev_arr[curr_i, curr_j + 1] < min_elev:
                            elevcm1 = elev_arr[curr_i, curr_j + 1]
                            min_elev = elevcm1
                            direc = 4
                
                # Look straight right and down
                if curr_i < nrow - 1 and curr_j < ncol - 1:
                    if not sfrlayout[curr_i + 1, curr_j + 1] == 0 and not ibnd[curr_i + 1, curr_j + 1] == 0:   # Step in if neighbor is a stream cell
                        if elev_arr[curr_i + 1, curr_j + 1] > 0 and (elev_arr[curr_i + 1, curr_j + 1] < elev_arr[curr_i, curr_j] and
                                                                     elev_arr[curr_i + 1, curr_j + 1] < sfrlayout_conn_candidate_elev):
                            sfrlayout_conn_candidate = sfrlayout[curr_i + 1, curr_j + 1]
                            sfrlayout_conn_candidate_elev = elev_arr[curr_i + 1, curr_j + 1]
                            stop_candidate = True
                    
                    elif not elev_arr[curr_i + 1, curr_j + 1] == 0 and not ibnd[curr_i + 1, curr_j + 1] == 0:   # Step here if neighbor is not an sfr cell
                        if elev_arr[curr_i + 1, curr_j + 1] < elev_arr[curr_i, curr_j] and elev_arr[curr_i + 1, curr_j + 1] < min_elev:
                            elevrp1cp1 = elev_arr[curr_i + 1, curr_j + 1]
                            min_elev = elevrp1cp1
                            direc = 7
                
                
                # Look straight up
                if curr_i > 0:
                    if not sfrlayout[curr_i - 1, curr_j] == 0 and not ibnd[curr_i - 1, curr_j] == 0:   # Step in if neighbor is a stream cell
                        if elev_arr[curr_i - 1, curr_j] > 0 and (elev_arr[curr_i - 1, curr_j] < elev_arr[curr_i, curr_j] and
                                                                 elev_arr[curr_i - 1, curr_j] < sfrlayout_conn_candidate_elev):
                            sfrlayout_conn_candidate = sfrlayout[curr_i - 1, curr_j]
                            sfrlayout_conn_candidate_elev = elev_arr[curr_i - 1, curr_j]
                            stop_candidate = True
                    
                    elif not elev_arr[curr_i - 1, curr_j] == 0 and not ibnd[curr_i - 1, curr_j] == 0:   # Step here if neighbor is not an sfr cell
                        if elev_arr[curr_i - 1, curr_j] < elev_arr[curr_i, curr_j] and elev_arr[curr_i - 1, curr_j] < min_elev:
                            elevcm1 = elev_arr[curr_i - 1, curr_j]
                            min_elev = elevcm1
                            direc = 3
                    
                
                # Look up and right
                if curr_i > 0 and curr_j < ncol - 1:
                    if not sfrlayout[curr_i - 1, curr_j + 1] == 0 and not ibnd[curr_i - 1, curr_j + 1] == 0:   # Step in if neighbor is a stream cell
                        if elev_arr[curr_i - 1, curr_j + 1] > 0 and (elev_arr[curr_i - 1, curr_j + 1] < elev_arr[curr_i, curr_j] and
                                                                 elev_arr[curr_i - 1, curr_j + 1] < sfrlayout_conn_candidate_elev):
                            sfrlayout_conn_candidate = sfrlayout[curr_i - 1, curr_j + 1]
                            sfrlayout_conn_candidate_elev = elev_arr[curr_i - 1, curr_j + 1]
                            stop_candidate = True
                    
                    elif not elev_arr[curr_i - 1, curr_j + 1] == 0 and not ibnd[curr_i - 1, curr_j + 1] == 0:   # Step here if neighbor is not an sfr cell
                        if elev_arr[curr_i - 1, curr_j + 1] < elev_arr[curr_i, curr_j] and elev_arr[curr_i - 1, curr_j + 1] < min_elev:
                            elevrm1cp1 = elev_arr[curr_i - 1, curr_j + 1]
                            min_elev = elevrm1cp1
                            direc = 6
                
                # Look straight down
                if curr_i < nrow - 1:
                    if not sfrlayout[curr_i + 1, curr_j] == 0 and not ibnd[curr_i + 1, curr_j] == 0:   # Step in if neighbor is a stream cell
                        if elev_arr[curr_i + 1, curr_j] > 0 and (elev_arr[curr_i + 1, curr_j] < elev_arr[curr_i, curr_j] and
                                                                 elev_arr[curr_i + 1, curr_j] < sfrlayout_conn_candidate_elev):
                            sfrlayout_conn_candidate = sfrlayout[curr_i + 1, curr_j]
                            sfrlayout_conn_candidate_elev = elev_arr[curr_i + 1, curr_j]
                            stop_candidate = True
                    
                    elif not elev_arr[curr_i + 1, curr_j] == 0 and not ibnd[curr_i + 1, curr_j] == 0:   # Step here if neighbor is not an sfr cell
                        if elev_arr[curr_i + 1, curr_j] < elev_arr[curr_i, curr_j] and elev_arr[curr_i + 1, curr_j] < min_elev:
                            elevrp1 = elev_arr[curr_i + 1, curr_j]
                            min_elev = elevrp1
                            direc = 1
                    
                # Look down and left
                if curr_i < nrow - 1 and curr_j > 0:
                    if not sfrlayout[curr_i + 1, curr_j - 1] == 0 and not ibnd[curr_i + 1, curr_j - 1] == 0:   # Step in if neighbor is a stream cell
                        if elev_arr[curr_i + 1, curr_j - 1] > 0 and (elev_arr[curr_i + 1, curr_j - 1] < elev_arr[curr_i, curr_j] and
                                                                 elev_arr[curr_i + 1, curr_j - 1] < sfrlayout_conn_candidate_elev):
                            sfrlayout_conn_candidate = sfrlayout[curr_i + 1, curr_j - 1]
                            sfrlayout_conn_candidate_elev = elev_arr[curr_i + 1, curr_j - 1]
                            stop_candidate = True
                    
                    elif not elev_arr[curr_i + 1, curr_j - 1] == 0 and not ibnd[curr_i + 1, curr_j - 1] == 0:   # Step here if neighbor is not an sfr cell
                        if elev_arr[curr_i + 1, curr_j - 1] < elev_arr[curr_i, curr_j] and elev_arr[curr_i + 1, curr_j - 1] < min_elev:
                            elevrp1cm1 = elev_arr[curr_i + 1, curr_j - 1]
                            min_elev = elevrp1cm1
                            direc = 8
                
                # if stop candidate found, don't move the cell indices
                if not stop_candidate:
                    # Direc corresponds to:
                    #  |----------------------
                    #  |  5  |    3    |  6  |
                    #  |----------------------
                    #  |  2  | cur_cel |  4  |
                    #  |----------------------
                    #  |  8  |    1    |  7  |
                    #  |----------------------
                    if direc == 0:
                        break
                    elif direc == 1:
                        curr_i += 1
                    elif direc == 2:
                        curr_j -= 1
                    elif direc == 3:
                        curr_i -= 1
                    elif direc == 4:
                        curr_j += 1
                    elif direc == 5:
                        curr_i -= 1
                        curr_j -= 1
                    elif direc == 6:
                        curr_i -= 1
                        curr_j += 1
                    elif direc == 7:
                        curr_i += 1
                        curr_j += 1
                    elif direc == 8:
                        curr_i += 1
                        curr_j -= 1
                
                if stop_candidate:
                    sfrlayout_new[i, j] = sfrlayout_conn_candidate
                    stop_candidate = False
                    break  # Bust out of while loop
                elif not stop_candidate:
                    # Check if encountered ibnd == 0, which may be a lake or boundary that drains out of model
                    if ibnd[curr_i, curr_j] == 0:
                        # This condition is dealt with after looping through all cells,
                        # see comment that starts, "Last step is set..."
                        break
                    pass  # Commence next downstream cell search

    # Last step is set the 0's in the vicinity of the lake equal to the negative of the lake connection
    for i in np.arange(0, nrow):
        for j in np.arange(0, ncol):
            if sfrlayout_new[i, j] == 0 and ibnd[i, j] > 0:
                sfrlayout_new[i, j] = -1
    
    # Once all cells are filled, save to array
    np.savetxt(os.path.join(pth, 'uzf_support','irunbnd_mf6.txt'), sfrlayout_new, fmt='%5d')
    
    return sfrlayout_new