mfsetup.discretization module

Functions related to the Discretization Package.

class mfsetup.discretization.ModflowGwfdis(*args, **kwargs)

Bases: ModflowGwfdis

property thickness

mfsetup.discretization.adjust_layers(dis, minimum_thickness=1)

Adjust bottom layer elevations to maintain a minimum thickness.

Parameters:

disflopy.modflow.ModflowDis instance

Returns:

new_layer_elevsndarray of shape (nlay, ncol, nrow)

New layer bottom elevations

mfsetup.discretization.cellids_to_kij(cellids, drop_inactive=True)

Unpack tuples of MODFLOW-6 cellids (k, i, j) to lists of k, i, j values; ignoring instances where cellid is None (unconnected cells).

Parameters:

cellidssequence of (k, i, j) tuples

drop_inactivebool

If True, drop cellids == ‘none’. If False, distribute these to k, i, j.

Returns:

k, i, j1D numpy arrays of integers

mfsetup.discretization.create_vertical_pass_through_cells(idomain)

Replaces inactive cells with vertical pass-through cells at locations that have an active cell above and below by setting these cells to -1.

Parameters:

idomainnp.ndarray with 2 or 3 dimensions. 2D arrays are returned as-is.

Returns:

revisednp.ndarray

idomain with -1s added at locations that were previous <= 0 that have an active cell (idomain=1) above and below.

mfsetup.discretization.deactivate_idomain_above(idomain, packagedata)

Sets ibound to 0 for all cells above active SFR cells.

Parameters:

packagedataMFList, recarray or DataFrame

SFR package reach data

Notes

This routine updates the ibound array of the flopy.model.ModflowBas6 instance. To produce a new BAS6 package file, model.write() or flopy.model.ModflowBas6.write() must be run.

mfsetup.discretization.fill_cells_vertically(top, botm)

In MODFLOW 6, cells where idomain < 1 are excluded from the solution. However, in the botm array, values are needed in overlying cells to compute layer thickness (cells with idomain != 1 overlying cells with idomain >= 1 need values in botm). Given a 3D numpy array with nan values indicating excluded cells, fill in the nans with the overlying values. For example, given the column of cells [10, nan, 8, nan, nan, 5, nan, nan, nan, 1], fill the nan values to make [10, 10, 8, 8, 8, 5, 5, 5, 5], so that layers 2, 5, and 9 (zero-based) all have valid thicknesses (and all other layers have zero thicknesses).

algorithm:

• given a top and botm array (top of the model and layer bottom elevations), get the layer thicknesses (accounting for any nodata values) idomain != 1 cells in thickness array must be set to np.nan

• set thickness to zero in nan cells take the cumulative sum of the thickness array along the 0th (depth) axis, from the bottom of the array to the top (going backwards in a depth-positive sense)

• add the cumulative sum to the array bottom elevations. The backward difference in bottom elevations should be zero in inactive cells, and representative of the desired thickness in the active cells.

• append the model bottom elevations (excluded in bottom-up difference)

Parameters:

top2D numpy array; model top elevations

botm3D (nlay, nrow, ncol) array; model bottom elevations

Returns:

top, botmfilled top and botm arrays

mfsetup.discretization.fill_empty_layers(array)

Fill empty layers in a 3D array by linearly interpolating between the values above and below. Layers are defined as empty if they contain all nan values. In the example of model layer elevations, this would create equal layer thicknesses between layer surfaces with values.

Parameters:

array3D numpy.ndarray

Returns:

filledndarray of same shape as array

mfsetup.discretization.find_remove_isolated_cells(array, minimum_cluster_size=10)

Identify clusters of isolated cells in a binary array. Remove clusters less than a specified minimum cluster size.

mfsetup.discretization.fix_model_layer_conflicts(top_array, botm_array, ibound_array=None, minimum_thickness=3)

Compare model layer elevations; adjust layer bottoms downward as necessary to maintain a minimum thickness.

Parameters:

top_array2D numpy array (nrow * ncol)

Model top elevations

botm_array3D numpy array (nlay * nrow * ncol)

Model bottom elevations

minimum thicknessscalar

Minimum layer thickness to enforce

Returns:

new_botm_array3D numpy array of new layer bottom elevations

mfsetup.discretization.get_highest_active_layer(idomain, null_value=-9999)

Get the highest active model layer at each i, j location, accounting for inactive and vertical pass-through cells.

mfsetup.discretization.get_layer(botm_array, i, j, elev)

Return the layers for elevations at i, j locations.

Parameters:

botm_array3D numpy array of layer bottom elevations

iscaler or sequence

row index (zero-based)

jscaler or sequence

column index

elevscaler or sequence

elevation (in same units as model)

Returns:

knp.ndarray (1-D) or scalar

zero-based layer index

mfsetup.discretization.get_layer_thicknesses(top, botm, idomain=None)

For each i, j location in the grid, get thicknesses between pairs of subsequent valid elevation values. Make a thickness array of the same shape as the model grid, assign the computed thicknesses for each pair of valid elevations to the position of the elevation representing the cell botm. For example, given the column of cells [nan nan 8. nan nan nan nan nan 2. nan], a thickness of 6 would be assigned to the second to last layer (position -2).

Parameters:

topnrow x ncol array of model top elevations

botmnlay x nrow x ncol array of model botm elevations

idomainnlay x nrow x ncol array indicating cells to be

included in the model solution. idomain=0 are converted to np.nans in the example column of cells above. (optional) If idomain is not specified, excluded cells are expected to be designated in the top and botm arrays as np.nans.

Examples

Make a fake model grid with 7 layers, but only top and two layer bottoms specified: >>> top = np.reshape([[10]]* 4, (2, 2)) >>> botm = np.reshape([[np.nan, 8., np.nan, np.nan, np.nan, 2., np.nan]]*4, (2, 2, 7)).transpose(2, 0, 1) >>> result = get_layer_thicknesses(top, botm) >>> result[:, 0, 0] array([nan 2. nan nan nan 6. nan])

example with all layer elevations specified note: this is the same result that np.diff(… axis=0) would produce; except positive in the direction of the zero axis >>> top = np.reshape([[10]] * 4, (2, 2)) >>> botm = np.reshape([[9, 8., 8, 6, 3, 2., -10]] * 4, (2, 2, 7)).transpose(2, 0, 1) >>> result = get_layer_thicknesses(top, botm) array([1., 1., 0., 2., 3., 1., 12.])

mfsetup.discretization.make_ibound(top, botm, nodata=-9999, minimum_layer_thickness=1, drop_thin_cells=True, tol=0.0001)

Make the ibound array that specifies cells that will be excluded from the simulation. Cells are excluded based on:

Parameters:

modelmfsetup.MFnwtModel model instance

Returns:

idomainnp.ndarray (int)

mfsetup.discretization.make_idomain(top, botm, nodata=-9999, minimum_layer_thickness=1, drop_thin_cells=True, tol=0.0001)

Make the idomain array for MODFLOW 6 that specifies cells that will be excluded from the simulation. Cells are excluded based on: 1) np.nans or nodata values in the botm array 2) np.nans or nodata values in the top array (applies to the highest cells with valid botm elevations; in other words, these cells have no thicknesses) 3) layer thicknesses less than the specified minimum thickness plus a tolerance (tol)

Parameters:

modelmfsetup.MF6model model instance

Returns:

idomainnp.ndarray (int)

mfsetup.discretization.make_irch(idomain)

Make an irch array for the MODFLOW 6 Recharge Package, which specifies the highest active model layer at each i, j location, accounting for inactive and vertical pass-through cells. Set all i, j locations with no active layers to 1 (MODFLOW 6 only allows valid layer numbers in the irch array).

mfsetup.discretization.make_lgr_idomain(parent_modelgrid, inset_modelgrid, ncppl)

Inactivate cells in parent_modelgrid that coincide with area of inset_modelgrid.

mfsetup.discretization.populate_values(values_dict, array_shape=None)

Given an input dictionary with non-consecutive keys, make a second dictionary with consecutive keys, with values that are linearly interpolated from the first dictionary, based on the key values. For example, given {0: 1.0, 2: 2.0}, {0: 1.0, 1: 1.5, 2: 2.0} would be returned.

Examples

>>> populate_values({0: 1.0, 2: 2.0}, array_shape=None)
{0: 1.0, 1: 1.5, 2: 2.0}
>>> populate_values({0: 1.0, 2: 2.0}, array_shape=(2, 2))
{0: array([[1., 1.],
           [1., 1.]]),
 1: array([[1.5, 1.5],
           [1.5, 1.5]]),
 2: array([[2., 2.],
           [2., 2.]])}

mfsetup.discretization.verify_minimum_layer_thickness(top, botm, isactive, minimum_layer_thickness)

Verify that model layer thickness is equal to or greater than a minimum thickness.

mfsetup.discretization.voxels_to_layers(voxel_array, z_edges, model_top=None, model_botm=None, no_data_value=0, extend_top=True, extend_botm=False, tol=0.1, minimum_frac_active_cells=0.01)

Combine a voxel array (voxel_array), with no-data values and either uniform or non-uniform top and bottom elevations, with land-surface elevations (model_top; to form the top of the grid), and additional elevation surfaces forming layering below the voxel grid (model_botm).

• In places where the model_botm elevations are above the lowest voxel elevations, the voxels are given priority, and the model_botm elevations reset to equal the lowest voxel elevations (effectively giving the underlying layer zero-thickness).

• Voxels with no_data_value(s) are also given zero-thickness. Typically these would be cells beyond a no-flow boundary, or below the depth of investigation (for example, in an airborne electromagnetic survey of aquifer electrical resisitivity). The vertical extent of the layering representing the voxel data then spans the highest and lowest valid voxels.

• In places where the model_top (typically land-surface) elevations are higher than the highest valid voxel, the voxel layer can either be extended to the model_top (extend_top=True), or an additional layer can be created between the top edge of the highest voxel and model_top (extent_top=False).

• Similarly, in places where elevations in model_botm are below the lowest valid voxel, the lowest voxel elevation can be extended to the highest underlying layer (extend_botm=True), or an additional layer can fill the gap between the lowest voxel and highest model_botm (extend_botm=False).

Parameters:

voxel_array3D numpy array

3D array of voxel data- could be zones or actually aquifer properties. Empty voxels can be marked with a no_data_value. Voxels are assumed to have the same horizontal discretization as the model_top and model_botm layers.

z_edges3D numpy array or sequence

Top and bottom edges of the voxels (length is voxel_array.shape[0] + 1). A sequence can be used to specify uniform voxel edge elevations; non-uniform top and bottom elevations can be specified with a 3D numpy array (similar to the botm array in MODFLOW).

model_top2D numpy array

Top elevations of the model at each row/column location.

model_botm2D or 3D numpy array

Model layer(s) underlying the voxel grid.

no_data_valuescalar, optional

Indicates empty voxels in voxel_array.

extend_topbool, optional

Option to extend the top voxel layer to the model_top, by default True.

extend_botmbool, optional

Option to extend the bottom voxel layer to the next layer below in model_botm, by default False.

tolfloat, optional

Depth tolerance used in comparing the voxel edges to model_top and model_botm. For example, if model_top - z_edges[0] is less than tol, the model_top and top voxel edge will be considered equal, and no additional layer will be added, regardless of extend_top. by default 0.1

minimum_frac_active_cellsfloat

Minimum fraction of cells with a thickness of > 0 for a layer to be retained, by default 0.01.

Returns:

layers3D numpy array of shape (nlay +1, nrow, ncol)

Model layer elevations (vertical edges of cells), including the model top.

Raises:

ValueError

If z_edges is not 1D or 3D

mfsetup.discretization.weighted_average_between_layers(arr0, arr1, weight0=0.5)