Introduction

This is a brief follow-on to my post on December 28, 2017 which concerned simulating the formation of an alluvial fan with a modified sediment transport modeling approach. In that model, as streams with different discharge intensities meander across the sloping terrain, varying sediment size class distributions are deposited at different locations, augmented by occasional debris flows. The result of these processes is the slow accumulation of sediments with an implied three-dimensional hydraulic conductivity structure.

As a subsequent step to the fan formation model, I wrote short python script to parse the simulated alluvial fan into layers, and to assign a lumped hydraulic conductivity value to each cell in each layer, based on the local sediment composition. The objective was to enable modeling the flow of groundwater through the fan, based on simple imposed boundary conditions, using MODFLOW. This relatively short script consists of the following code:

from __future__ import print_function
from numpy import *
from pandas import *
import matplotlib.pyplot as plt

options.mode.chained_assignment = None

class Sediment:

    def __init__(self):
        # read sediment properties
        self.name = []
        self.K = []
        lineInput = []
        inputFile = open('sed_props.txt','r')
        for line in inputFile: lineInput.append(line.split())
        inputFile.close()
        for i in xrange(1, len(lineInput)):
            self.name.append(lineInput[i][0])
            self.K.append(float(lineInput[i][1]))
        self.K = array(self.K)
        print('Read sediment properties.')

    def CalcK(self, fractions):
        # sediment size class weighted hydraulic conductivity
        K = dot(self.K, fractions)
        return 10**K

class Params:

    def __init__(self):
         # grid parameters
        lineInput = []
        inputFile = open('model_params.txt','r')
        for line in inputFile: lineInput.append(line.split())
        inputFile.close()
        self.numLayers = int(lineInput[0][1])                # grid origin
        self.thickMin = float(lineInput[1][1])
        self.zMax = float(lineInput[2][1])
        self.dz = float(lineInput[3][1])
        print('Read model grid constraints.')

class Fan:

    def __init__(self, grid, sed):
        # import output of alluvial fan generator and convert to data frame
        self.fan_df = read_csv('sed_distribution.csv', sep=',')
        self.fan_df = self.fan_df[self.fan_df['z']<=grid.zMax]
        self.xSet = list(set(self.fan_df['x']))
        self.ySet = list(set(self.fan_df['y']))
        self.xSet.sort()
        self.ySet.sort()
        self.numRows = len(self.ySet)
        logK = dot(self.fan_df[sed.name], log10(sed.K))
        self.fan_df['K'] = 10.**logK
        self.epsilon = 0.01                      # factor used to correct layer numbering scheme
        print('Read fan model output.')               

    def Layers(self, grid):
        # return dataframes based on fan_df, but vertically aggregated into grid.numLayers
        print('Processing ...')
        aggStatus = False
        for i, xCol in enumerate(self.xSet):
            for j, yRow in enumerate(self.ySet):
                column_df = self.fan_df[(self.fan_df['x']==xCol) & (self.fan_df['y']==yRow)]
                if len(column_df) < grid.numLayers+1:                     # add extra sediment slabs to this location                     for k in xrange(len(column_df), grid.numLayers+1):                         column_df = column_df.append(column_df.tail(1), ignore_index=True)                         column_df['z'][k] = column_df['z'][k-1] + grid.dz                 z0 = column_df['z'].min()                 zExtent = column_df['z'].max() - z0 + self.epsilon                 dz = zExtent/grid.numLayers                  column_df['layer'] = (column_df['z']-z0)/dz                 column_df['layer'] = column_df['layer'].astype('int64')                        column_df = column_df.groupby('layer').mean()                 column_df.reset_index(inplace=True)                          column_df['row'] = self.numRows-j 		# renumber to match MODFLOW row numbering scheme                 column_df['col'] = i+1                 column_df['zBase'] = z0 + column_df['layer']*dz                                          column_df['zTop'] = column_df['zBase'] + dz                                          column_df['ibound'] = bool(zExtent>=grid.thickMin) 	# fan areas that are too thin are marked as inactive
                column_df = column_df[['col', 'row', 'layer', 'zBase', 'zTop', 'K', 'ibound']]
                if not aggStatus:
                    layers_df = column_df.copy()
                    aggStatus = True
                else: layers_df = concat([layers_df, column_df], axis=0)
        layers_df.reset_index(inplace=True)
        layers_df['layer'] = grid.numLayers - layers_df['layer']
        hydro_df = layers_df[['col', 'row', 'layer', 'K', 'ibound']]        # hydrology props in 3-D
        print('Generated hydrology data set.')
        layerStruct_df = layers_df[layers_df['layer']==1]
        zb = array(layerStruct_df['zBase'])
        topModel = array(layerStruct_df['zTop'])
        layerStruct_df = layerStruct_df[['col', 'row']]
        layerStruct_df['top_model'] = topModel
        layerStruct_df['base_1'] = zb
        for i in xrange(2, grid.numLayers+1):
            nextLayer_df = layers_df[layers_df['layer']==i]
            zb = array(nextLayer_df['zBase'])
            layerStruct_df['base_' + str(i)] = zb
        print('Generated layer structure data set.')
        return hydro_df, layerStruct_df

### main code ###

def FanGrid():

    # set up sediment class
    sed = Sediment()    

    # read grid parameter constraints
    grid = Params()

    # import and process prior fan model output
    fan = Fan(grid, sed)

    # collapse vertical sediment stacks into layers
    hydro_df, layerStruct_df = fan.Layers(grid)
    active_df = hydro_df[hydro_df['ibound']==1]
    logK = log10(active_df['K'])		# plot log K histogram (over all active cells)
    plt.figure()
    ax = logK.plot.hist(bins=40)
    ax.set_xlabel('Log K')
    plt.show()

    print('Writing output.')
    hydro_df.to_csv('hydro.csv', index=False, sep='\t')
    layerStruct_df.to_csv('layers_struct.csv', index=False, sep='\t')

    print('Done.')

The script employ the pandas library to streamline coding and to speed up execution (in place of for-loops, etc.). The input files for the script consist of (1) the raw comma-delimited text files generated by the alluvial fan simulator, (2) a file containing a short list of constraints (e.g., number of layers, maximum elevation clipping value) for gridding, and (3) a file containing the names of textural classes of the simulated alluvial sediments, along with their corresponding user-provided hydraulic conductivities. A logarithmic weighting scheme is used to posit a resulting hydraulic conductivity for each cell. These files, along with the python script itself, can be viewed and downloaded from my GitHub repository.

histogram
Distribution of hydraulic conductivity (log m/day) across all active cells in the parsed model grid to be used as MODFLOW input. Values reflect the distributions of sediment size classes in each cell and the posited end-member conductivities per textural class (e.g., 0.00001 m/day for clay, 0.01 m/day for silt, 10 m/day for fine sand, 100 m/day for coarse sand, and 500 m/day for gravel).

The script’s output consists of two files that can be used directly to set up a MODFLOW simulation (in this case, via the U.S. Geological Survey’s ModelMuse pre-processor): (1) a layer structure file indicating the model top surface and individual layer bottom elevations for all column and row numbers comprising the grid, and (2) the distribution of hydraulic conductivity and active cell labels (i.e., boolean flag) for every cell in the system (by column, row, and layer number).

Example Application

The example synthetic alluvial fan presented in my earlier post consists of a 30 x 50 grid of cells with respect to the horizontal plane, with each cell a square, 30 meters on a side. This small-scale alluvial fan is fed by two sediment sources. The parsing script shown above was used to divide the modeled fan into eight layers, varying in thickness from the sediment influx locations to the more distal portions of the resulting fan. Using the layer structure and property distribution output files, I created a steady-state MODFLOW model for the fan, based on the following set of assumptions:

Specified head boundary conditions are employed at each of the two sediment influx locations (matching the original alluvial fan generation model) and along a portion of the distal edge of the fan where the total sediment elevation is approximately 4 meters or less.

  • All layers are convertible (i.e., confined or unconfined conditions, as applicable).
  • Vertical recharge from stream flow or precipitation is not considered.
  • The system is not stressed by any other fluxes, such as pumping wells.
  • The vertical-to-horizontal hydraulic conductivity anisotropy ratio is 1:10.

Once the steady-state head distribution was calculated, particle tracking with MODPATH  was used to delineate groundwater flow lines away from the two sediment influx locations. The relatively coarse grid implied by the original fan generation model was maintained, although this could have been refined (1) within the above-listed script itself by introducing additional columns, rows, and or layers and interpolating the property values, or (2) by adjusting the grid directly using within ModelMuse using its various grid refinement tools.

grid_sections
MODFLOW grid representing the simulated alluvial fan (active cells). The extent of the fan in the x-direction is approximately 1 km. A relatively coarse grid is used for this demonstration.

layers_1_and_2

layers_3_and_4

layers_5_and_6

layers_7_and_8
Synthetic hydraulic conductivity distributions per layer (left) and corresponding MODPATH particle tracks representing steady-state flow from the two sources along the left boundary. Warmer colors among the particle tracks correspond to the longest travel times and hence represent slow groundwater flow pathways. The slices, of variable thickness and sequenced from the model surface to match MODFLOW’s layer numbering scheme, correspond to the same synthetic alluvial fan shown in my earlier post.

In summary, the model results indicate a complex flow environment, with the manifestation of preferential flow pathways corresponding to areas of higher hydraulic conductivity as well as flow between layers (as implied, for example, in the top two layers where portions of the aquifer are dry). The coarse grid, and the selection of the sediment transport momentum term employed in the fan generation model, results in some pathlines maintaining uni-directional flow over long sections of the model, a somewhat unnatural-looking result that could be remedied by a finer grid. On the other hand, a subsequent step I have in mind is to explore the possible impacts of sediment heterogeneity (e.g., clay-rich versus sand-rich) on groundwater chemistry. This will entail a reactive transport model to address ion exchange and other processes. A coarse mesh will offer the advantage of reduced computation time when exploring this issue.

More later.