ArchPy2Modflow module¶

This module propose several functions a class to interface ArchPy with MODFLOW 6

class ap_mf.archpy2modflow(T1, sim_name='sim_test', model_dir='workspace', model_name='test', exe_name='mf6', vb=1)¶

Bases: object

Class to convert an ArchPy table to a MODFLOW 6 model

Parameters:

T1 (base.Arch_table object) – ArchPy table object to convert
sim_name (str) – name of the simulation
model_dir (str) – directory where the model will be saved
model_name (str) – name of the model
exe_name (str) – path to the mf6 executable

create_reduced_child_object(iu=0, ifa=0, ip=0, properties=None)¶: This function creates a child object of the current archpy2modflow object. where the contained archpy model has only a limited number of units, facies and properties simulations It is a copy of the current object but with a reduced number of units, facies and properties.

create_sim(grid_mode='archpy', iu=0, lay_sep=1, modflowgrid_props=None, xorigin=0, yorigin=0, angrot=0, factor_x=None, factor_y=None, factor_z=None, unit_limit=None, surface_layer=False, surface_thickness=0.05, top=None, botm=None, idomain=None, load_existing_model=False)¶

Create a modflow simulation from an ArchPy table

Parameters:

grid_mode (str) – “archpy” : use the grid defined in the ArchPy table “layers” : use the surfaces of each unit to define the grid “new_resolution” : use factors to change the resolution of the grid In this case, factor_x, factor_y and factor_z must be provided “disv” : use a disv grid, modflowgrid_props must contain the properties of the disv grid “disu” : use a disu grid, modflowgrid_props must contain the properties of the disu grid “given_layers” : use given top, botm and idomain arrays to define the grid
iu (int) – index of the unit to use when grid_mode is “layers”
lay_sep (int or list of int of size nlay) – if grid_mode is layers, lay_sep indicates the number of layers to separate each unit
factor_x (float) – factor to change the resolution of the grid in the x direction. e.g. 2 means that the resolution will be divided by 2
factor_y (float) – factor to change the resolution of the grid in the y direction.
factor_z (float) – factor to change the resolution of the grid in the z direction.
unit_limit (unit name) – unit under which cells are considered as inactive
surface_layer (bool) – if grid_mode is layers, surface_layer indicates whether to add a thin layer (of thickness surface_thickness) on top of all others within the model domain
surface_thickness (float) – if grid_mode is layers and surface_layer is True, indicates the thickness of the additional surface layer

create_sim_energy(strt_temp=10, ktw=0.56, kts=2.5, al=1, ath1=1, prsity=0.2, cpw=4186, cps=840, rhow=1000, rhos=2650, lhv=2260000.0)¶: Create a simulation object for an energy model with some default values

create_sim_transport(strt_conc=0.0, porosity=0.2, diff=1e-09, adv_scheme='upstream', alh=None, ath1=None, decay=0.0, decay_0=False, decay_1=False)¶: Create a simulation object for a transport model with some default values

create_ssm_e(sourcerecarray)¶

Create the ssm package for the energy model

Parameters:: sourcerecarray (list of tuples) – list of tuples where each indicate a gwf package which has auxiliar variables (package name, AUX, auxname)

create_ssm_t(sourcerecarray)¶

Create the ssm package for the transport model

Parameters:: sourcerecarray (list of tuples) – list of tuples where each indicate a gwf package which has auxiliar variables (package name, AUX, auxname)

get_gw_energy()¶: Get the MODFLOW 6 groundwater energy object

get_gw_transport()¶: Get the MODFLOW 6 groundwater energy object

get_gwf()¶: Get the MODFLOW 6 groundwater flow object

get_heads(kstpkper=(0, 0))¶: Get the heads of the simulation

get_list_active_cells()¶: Get the list of active cells in the modflow model

get_mp()¶: Get the modpath object

get_sim()¶: Get the MODFLOW 6 simulation object

get_sim_energy()¶: Get the MODFLOW 6 energy simulation object

get_sim_prt()¶: Get the MODFLOW 6 particle tracking simulation object

get_sim_transport()¶: Get the MODFLOW 6 energy simulation object

gmean()¶

mp_create(mpexe, trackdir='forward', locs=None, rowcelldivisions=1, columncelldivisions=1, layercelldivisions=1, one_group_per_cells=False, list_p_coords=None, **kwargs)¶

Create a modpath simulation from an ArchPy table

Parameters:

mpexe (str) – path to the modpath executable
trackdir (str) – direction of tracking
locs (list of tuples) – list of cells to put particles. List of int where values correspond to the indices of the cells. Each cell in locs will be then spatially divided and will output rowcelldivisions * columncelldivisions * layercelldivisions particles
rowcelldivisions (int) – number of row divisions. If 1, no division, if 2, divide each cell in 2 parts, etc.
columncelldivisions (int) – number of column divisions
layercelldivisions (int) – number of layer divisions
one_group_per_cells (bool) – if True, create one particle group per cell
list_p_coords (list of tuples) – list of particles coordinates. Each tuple must have 3 values (xp, yp, zp) corresponding to the coordinates of the particle

mp_get_endpoints_object()¶: Get the endpoints object from the modpath simulation

mp_get_facies_path_particle(i_particle, fac_time=1.1574074074074073e-05, iu=0, ifa=0)¶: Function to retrieve the facies sequence along a pathline

mp_get_pathlines_object()¶: Get the pathlines object from the modpath simulation

mp_run(silent=False)¶: Write MODPATH files and run the particle tracking simulation using MODPATH 7

plot_3D_heads(kstpkper=(0, 0), ax=None, **kwargs)¶: Plot the heads of the simulation

prt_create(prt_name='test', workspace='./', trackdir='forward', list_p_coords=None, cellids=None, particle_names=None, **kwargs)¶

Create a particle tracking simulation from a list of coordinates

Parameters:

prt_name (str) – name of the particle tracking simulation
workspace (str) – path to the workspace where the simulation will be saved
trackdir (str) – tracking direction. Can be “forward” or “backward”
list_p_coords (list of tuples) – list of particles coordinates. Each tuple must have 3 values (xp, yp, zp) corresponding to the coordinates of the particle
particle_names (list of str) – list of particle names. If None, particles will be named “P1”, “P2”, etc.
kwargs (dict) – additional arguments to pass to the ModflowPrt simulation

prt_get_facies_path_particle(i_particle=1, fac_time=1.1574074074074073e-05, iu=0, ifa=0, facies_archpy=False)¶

Function to retrieve the facies sequences along a pathline

Parameters:

i_particle (int) – index of the particle to retrieve the pathline for.
fac_time (float) – factor to convert time from seconds to days (default is 1/86400)
iu (int) – index of the unit realization to retrieve the facies from (default is 0)
ifa (int) – index of the facies realization to retrieve the facies from (default is 0)
facies_archpy (bool) – if True, retrieve the facies from the the fine ArchPy model (not upscaled)

prt_get_pathlines(i_particle=None)¶: Retrieve the raw pathlines from the particle tracking simulation

prt_run(silent=False)¶: Write PRT files and run the particle tracking simulation using MODFLOW 6

set_cnd(xt3d_off=None, xt3d_rhs=None, alh=1, alv=None, ath1=1, ath2=None, atv=None, ktw=0.56, kts=2.5, vb=1, **kwargs)¶

Set the conduction and dispersion package

Parameters:

xt3d_off (bool) – flag indicating whether the XT3D package is active
xt3d_rhs (bool) – flag indicating whether the XT3D package is using the right-hand side formulation
alh (float or array of model dimension or string) – longitudinal dispersivity (if flow is horizontal) If a string is provided, it is assumed to be the name of the property in the ArchPy table
alv (float or array of model dimension or string) – longitudinal dispersivity (if flow is vertical)
ath1 (float or array of model dimension or string) – horizontal (y direction) transverse dispersivity (if flow is directed in the x-direction)
ath2 (float or array of model dimension or string) – horizontal (z direction) transverse dispersivity (if flow is directed in the x-direction)
atv (float or array of model dimension or string) – transverse (x-y plane) dispersivity (if flow is vertical)
ktw (float or array of model dimension or string) – thermal conductivity of water
kts (float or array of model dimension or string) – thermal conductivity of solid
vb (int (0 or 1)) – flag to print messages

set_dsp(alh=None, ath1=None, alv=None, ath2=None, atv=None, diffc=1e-09, xt3d_off=False, xt3d_rhs=False, vb=1)¶

Set the dispersion package for the transport model Note that this package reinitalize pre assigned values

Parameters:

alh (float) – longitudinal dispersivity if flow is horizontal
ath1 (float) – transverse dispersivity (in y direction) if flow is horizontal
ath2 (float) – transverse dispersivity (in z direction) if flow is horizontal
alv (float) – longitudinal dispersivity if flow is vertical
atv (float) – transverse dispersivity (in xy direction) if flow is vertical
diffc (float) – molecular diffusion coefficient
xt3d_off (bool) – xt3d formulation

set_est(cpw=4184.0, cps=840, rhow=1000, rhos=2500, lhv=2453500.0, porosity=None, save_flows=True, vb=1, **kwargs)¶

Set the heat transport mass storage package

Parameters:

cpw (float) – specific heat capacity of the fluid
cps (float or array of model dimension or string) – specific heat capacity of the solid if a string is provided, it is assumed to be the name of the property object (see base.Prop)
rhow (float) – density of the fluid
rhos (float or array of model dimension or string) – density of the solid if a string is provided, it is assumed to be the name of the property in the ArchPy table
lhv (float) – latent heat of vaporization
porosity (float or array of model dimension or string) – porosity. Note that it is preferable to update the porosity using the set_porosity() method. As this method only update porosity in the est package but not in the other packages (mp, prt, …) if a string is provided, it is assumed to be the name of the property in the ArchPy table
save_flows (bool) – flag to save flows in the cell budget file
vb (int (0 or 1)) – flag to print messages

set_grid_layers_mode(iu=0, unit_limit=None, lay_sep=1)¶: Change the ArchPy unit simulation used to defined the layers with the Layers mode. Note that this preserves the modflow model and do not recreate a new simulation object. For this use create_sim() method.

set_imsgwe(**kwargs)¶: Create the ims package for the energy model

set_imsgwt(**kwargs)¶: Create the ims package for the transport model

set_k(k_key='K', iu=0, ifa=0, ip=0, log=False, k=None, k22=None, k33=None, k_average_method='anisotropic', average_facies=False, upscaling_method='simplified_renormalization', xt3doptions=None, surface_k=1e-06)¶

Set the hydraulic conductivity for a specific facies

Parameters:

k_key (str) – name of the property to use for the hydraulic conductivity, default is “K”
iu (int) – index of the unit simulation to use when grid_mode is “layers”
ifa (int) – index of the facies simulation to use when grid_mode is “layers”
ip (int) – index of the property simulation to use when grid_mode is “layers”
log (bool) – if True, the hydraulic conductivity property is assumed to be in log10 scale, default is False
k (np.ndarray, optional) – hydraulic conductivity in x direction if provided, the hydraulic conductivity will be set to this value, default is None if k is None, the hydraulic conductivity will be computed from the property defined in the ArchPy table
k22 (np.ndarray, optional) – hydraulic conductivity in y direction
k33 (np.ndarray, optional) – hydraulic conductivity in z direction
k_average_method (str) – method to use for averaging the hydraulic conductivity when grid_mode is “layers” can be “arithmetic”, “harmonic” or “anisotropic”, default is “arithmetic”
average_facies (bool) – if True, the facies will be upscaled and saved in the model if False, the facies will not be upscaled, default is False
surface_k (float) – hydraulic conductivity to set for the surface layer, default is -6 #FIXME only used when grid_mode is ‘layers’

set_model_dir(model_dir)¶: Set the model directory for the modflow simulation

set_mst(porosity=None, decay=None, decay_0=False, decay_1=False, decay_sorbed=None, sorption=None, bulk_density=None, distcoef=None, sp2=None, vb=1)¶

Set the mass storage package for the transport model To assign a property from an Archpy model, use the name of the property as a string

Parameters:

porosity (float, array or string) – porosity of the medium
decay (float, array or string) – decay coefficient. A negative value indicates production
decay_0 (bool) – if True, decay is a zero order decay
decay_1 (bool) – if True, decay is a first order decay
decay_sorbed (float, array or string) – decay coefficient for the sorbed phase
sorption (string) – sorption keyword to indicate that soprtion is used. Options are: linear, langmuir, freundlich
bulk_density (float, array or string) – bulk density of the medium in mass/length^3 Bulk density corresponds to the amount of solid per unit volume of the medium
distcoef (float, array or string) – distribution coefficient for the equilibrium sorption. unit is length^3/mass
sp2 (float, array or string) – exponent for the freundlich isotherm

set_oc_t(saverecord, printrecord, concentrationprintrecord)¶

Set the output control for the transport model

Parameters:

budget_filerecord (str) – name of the budget file
concentration_filerecord (str) – name of the concentration file

set_porosity(iu=0, ifa=0, ip=0, porosity=None, prop_key='porosity', packages=None)¶

Set the porosity of the model

Parameters:

iu (int) – unit simulation index
ifa (int) – facies simulation index
ip (int) – property simulation index
porosity (float) – porosity value, if None, the porosity is taken from the table according to the prop_key
prop_key (str) – key of the property in the table
packages (list) – list of packages to update the porosity. Can be “mp”, “prt” and “energy”. If None, all existing packages are updated

set_strt(heads=None)¶: Set the starting heads

set_strt_conc(strt_conc)¶

Set the initial concentration for the transport model

Parameters:: strt_conc (float or array of model dimension) – initial concentration

set_strt_temp(strt_temp)¶

Set the initial temperature of the energy model

Parameters:: strt_temp (float or array of model dimension) – initial temperature

set_tdisgwe(perioddata)¶

Create the tdis package for the energy model

Parameters:: perioddata (list of tuples) – list of tuples where each tuple is a period data (perlen, nstp, tsmult)

set_tdisgwt(perioddata)¶

Create the tdis package for the transport model

Parameters:: perioddata (list of tuples) – list of tuples where each tuple is a period data (perlen, nstp, tsmult)

upscale_prop(prop_key, iu=0, ifa=0, ip=0, method='arithmetic', log=False)¶

ap_mf.array2cellids(array, idomain)¶: Convert an array to a list of cellids array : 3D array, array to convert idomain : 3D array, idomain array which determine if a cell is active or not (1 active, 0 inactive)

ap_mf.cellidBD(idomain, layer=0)¶: extract the cellids at the boundary of the domain at a given layer idomain : 3D array, idomain array which determine if a cell is active or not (1 active, 0 inactive) layer : int, layer on which the boundary cells are extract

ap_mf.check_thk(top, botm)¶: check if all the cells in a modflow model (given top and botm array) which have a thickness <= 0 for and this for each layer input : top (the top surface) and botm (botom of each layer) output : lst of bool (false mean everything’s okay in that specific layer !)

ap_mf.get_locs(nodes, nx, ny)¶

ap_mf.get_nodes(locs, nx, ny)¶

ap_mf.hmean_speed(a)¶

ap_mf.mask_below_unit(T1, unit, iu=0)¶

ap_mf.mean_speed(a)¶

ap_mf.plot_particle_facies_sequence(arch_table, df, plot_time=False, plot_distance=False, proportions=False, fac_time=1.1574074074074073e-05)¶

Plot the facies sequence of a particle in time and distance

Parameters:

arch_table (base.Arch_table object) – ArchPy table object
df (pandas dataframe) – dataframe with the particle data. Must contain the columns “facies” or “facies_prop_i” where i is the facies id and “time” or “distance” or “cum_distance”.
plot_time (bool, optional) – if True, plot the facies sequence in time. The default is False.
plot_distance (bool, optional) – if True, plot the facies sequence in distance. The default is False.
proportions (bool, optional) – if True, plot the proportions of each facies in time or distance. The default is False. df must contain the columns “facies_prop_i” where i is the facies id.
time_unit (str, optional) – unit of time in the column “time”. The default is “s”.
resampling (str, optional) – resampling frequency for the time series. The default is “d”. Check pandas documentation for more information.
fac_time (float, optional) – factor to apply to time column in case number are too large and cause an overflow error in pandas

ap_mf.points2grid_index(points, grid)¶

Convert a list of points to a list of grid indices

Parameters:

points (list of tuples) – list of points to convert. Each tuple must have 3 values (xp, yp, zp) corresponding to the coordinates of the point
grid (modflow grid object) – modflow grid object

Returns:

list of tuples with the indices of the points in the grid. Depends on the type of grid if grid is disv, return (ilay, icell) if grid is disu, return (icell)

Return type:

list of tuples