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]) self.K.append(float(lineInput[i])) 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) # grid origin self.thickMin = float(lineInput) self.zMax = float(lineInput) self.dz = float(lineInput) 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.
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).
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.
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.