sutra2.ModPath_Well module

class sutra2.ModPath_Well.ModPathWell(schematisation: HydroChemicalSchematisation, workspace: Optional[str] = None, modelname: Optional[str] = None, mf_exe='mf2005.exe', mp_exe='mpath7.exe', bound_left: str = 'xmin', bound_right: str = 'xmax', bound_top: str = 'top', bound_bot: str = 'bot', bound_north: str = 'ymin', bound_south: str = 'ymax', trackingdirection='forward')

Bases: object

Compute travel time distribution using MODFLOW and MODPATH.

MP7modelrun(mf_namfile=None, mp_exe=None)

Modpath model run.

assign_material(schematisation: dict, dict_keys: Optional[dict] = None)

Assign grid material using subkeys in schematisation_dict.

assign_wellloc(schematisation: dict, dict_key: str = 'well_parameters', well_names: list = 'None', discharge_parameter='well_discharge')

Determine the location of the pumping wells and the relative discharge per cell.

Boundaries of the well locations are to be included in the dictionaries with keys “dict_keys” [list], based on presence of “discharge_parameter” [str (default = “Q”)]:

bound_left: str = “xmin” bound_right: str = “xmax” bound_top: str = “top” bound_bot: str = “bot” bound_north: str = “ymin” bound_south: str = “ymax”

# Output: # Creates a dict with list of tuples representing the well locations well_loc = {“well1”: [(5,0,0),(6,0,0)], “well2”:}

# Calculate cumulative KD per well [m2/day] (to calculate relative discharge per cell) KD_well = {“well1”: 500., “well2”: 100.,…} # Create stress period data [list of lists per stress period]: # wel_spd = {0: [[lay,row,col,discharge1],[lay,row,col,discharge2]]

# Total discharge per day per well Qwell_day = = {“well1”: -1000.., “well2”: -1.,…}

axisym_correction(grid: array, dtype: Optional[str] = None, theta=6.283185307179586)

Correct modflow parameters to correctly calculate axisymmetric flow. Adjust the array ‘grid’ if it already exists, using multiplier theta and xmid (nrow == 1 or nrow == 2) or ymid (ncol == 1 or ncol == 2).

calc_flux_cell(frf, flf, fff, loc, flux_direction='total')

Calculate the total volume flux along the cell boundary, using frf, flf and fff.

Parameters

  • frf (np.array) – flux right face (= positive to the right) [m^3 d-1]

  • flf (np.array) – flux lower face (= positive in downward direction) [m^3 d-1]

  • fff (np.array) – flux front face (= positive in ‘southward’ direction) [m^3 d-1]

  • loc (tuple) – cell location (lay,row,col)

  • flux_direction (str) – cell flux along boundary [‘total’,’west’,’east’,’north’,’south’,’top’,’bottom’, ‘vertical’, ‘horizontal’] (NB ‘total’: ‘magnitude’ of flow in cell) NB2: ‘vertical –> corrected for additional inflow horizontally NB3: ‘vertical –> corrected for additional inflow vertically

calc_traveltime_vadose_analytical(vadose_parameters: dict = {}, dict_key='vadose_zone', distance=None, gw_level=None)

Inherit functionality from Analytical_Well module

cell_bounds(schematisation: dict, dict_key: str = 'None', dict_subkey: str = 'None', model_type='axisymmetric')

Create cell boundaries of box created using following boundaries: self.bound_left: str = “xmin” self.bound_right: str = “xmax” self.bound_top: str = “top” self.bound_bot: str = “bot” self.bound_north: str = “ymin” # not required for axisymmetric or 2D model self.bound_south: str = “ymax” # not required for axisymmetric or 2D model # Return boundary indices of row and columns plus parameter value return layidx_min,layidx_max,rowidx_min,rowidx_max,colidx_min,colidx_max

create_MNW(schematisation: dict, dict_key: str = 'well_parameters', well_names: Optional[list] = None, discharge_parameter='well_discharge')

Add multi nodal well Package to the MODFLOW model. # Function allows for a single row, col combination in current version per identified well. Documented in flopy docs - MNW2.

Parameters

  • wellid (str or int) – name of the well

  • nnodes (int) – nr of filter screens per vertical well (if nnodes < 0). Allows for input of ztop and zbot per filter screen.

  • ztop (float) – top elevation of open intervals of vertical well. [m]

  • zbotm (float) – bottom elevation of open intervals of vertical well. [m]

  • rw (float) – well radius of the vertical well(s) [m]

create_MP7object(mp_exe='mpath7', mf_model=None, headfilename=None, budgetfilename=None)

Load modpath7 model, using input from modflow model ‘mf_model’ (DIS), results from head calculations, and cell-by-cell data.

create_bas()

Add basic Package to the MODFLOW model

create_dis(per_nr=0)

Add dis Package to the MODFLOW model

create_lpf()

Add lpf Package to the MODFLOW model

create_mfobject(mf_exe='mf2005.exe', fname_nam=None)

create mf-object (if fname_nam equals None). Otherwise: load mf-model using nam-file ‘fname_nam’.

create_modflow_input(**kwargs)

Generate flopy files.

create_modflow_packages(**kwargs)

Create modflow packages used in the model run.

create_modpath_packages(mp_exe='mpath7', **kwargs)

Create modpath packages used in the model run.

create_mpbas()

Add BAS Package to the ModPath model

create_oc(spd_oc={(0, 0): ['save head', 'save budget']}, per_nr=0)

Add output control Package to the MODFLOW model.

create_pcg(hclose=0.0001, rclose=0.001)

Add PCG Package to the MODFLOW model

create_rch(rech=None)

Add recharge Package to the MODFLOW model.

create_wel(spd_wel=None)

Add well Package to the MODFLOW model.

extract_vadose_zone_parameters(schematisation: dict, schematisation_type: str, dict_key='geo_parameters', parameter: str = 'vadose', remove_layer=None)

Extract vadose_zone interval from geo_parameters dictionary ‘dict_key’, if ‘parameter’ value in the ‘dict_key’ (e.g. “vadose”) is True. remove_layer: –> remove the key from schematisation[dict_key] if it is a vadose_zone (T/F).

fill_df_flowline(df_particle, model_cbc)

Fill df_flowline dataframe

Parameters

df_particle (pd.DataFrame) –

Returns

df_flowline – Column ‘flowline_id’: Integer Column ‘flowline_type’: string Described the type of contamination associated with the flowline, either ‘diffuse_source’ or ‘point_source’. Column ‘flowline_discharge’: float. Discharge associated with the flowline, [m3/d]. Column ‘particle_release_day’: float Column ‘input_concentration’ Column ‘endpoint_id’: Integer ID of Well (or drain) where the flowline ends. Column ‘well_discharge’: float Column ‘substance’: string Column ‘removal_function’: string

Return type

pandas.DataFrame

fill_df_particle(particle_group, pg_nodes, parm_list: list = [], mppth=None)

Fill the df_particle, ordered by flowline_id and total_travel_time

Parameters

  • particle_group (str) – particle group name of released particles

  • pg_nodes (list of tuple or list of list) – locations at which particles were started [(lay,row,col),(…)]

  • parm_list – list of columns to include in df_particle (to be returned)

  • mppth (str or Path) – ‘modpath_name.pth’ file location

Returns

df_particle – Holder for df_particle, filled by the flowline_id

flowline_id: int

ID of the flowline

xcoord: float

x coordinate of particle

ycoord: float

y coordinate of particle

zcoord: float

z coordinate of particle

total_travel_time:

travel time since start of model [days]

material: str

identifier for geological layer or other material/object wherein particle resides

redox: str

redox condition [‘suboxic’,’anoxic’,’deeply_anoxic’]

temp_water: float

Water temperature [degrees celcius]

porosity: float

effective porosity [-]

dissolved_organic_carbon:

dissolved organic carbon (DOC) fraction [-]

pH: float

pH of the water [-]

fraction_organic_carbon

organic carbon fraction (foc) [-]

solid_density: float

Bulk density of the material [kg / L]

Return type

pandas.dataframe

fill_grid(schematisation: dict, dict_keys: Optional[list] = None, parameter: str = 'None', grid: Optional[array] = None, dtype: str = 'float')

Assign values to ‘grid’ [np.array] for parameter name ‘parameter’ [str], using the keys dict_keys [list] in schematisation dictionary ‘self.schematisation’. ‘dtype’ [grid dtype] –> grid dtype [str] is obtained from grid if initial array is given. The schematisation type refers to the model_type[“axisymmetric”,”2D” or “3D”]

Boundaries of the parameter values are to be included in the dictionaries:

self.bound_left: str = “xmin” self.bound_right: str = “xmax” self.bound_top: str = “top” self.bound_bot: str = “bot” self.bound_north: str = “ymin” self.bound_south: str = “ymax”

The function returns: - grid # grid [np.array] filled with (numeric) values for parameter ‘parameter’.

generate_modflow_files()

Write package data MODFLOW model.

get_node_ID(locs)

Obtain/return model node index (int) belonging to layer ‘iLay’, row ‘iRow’ and column ‘iCol’. for point locations ‘locs’.

get_node_indices(xyz_nodes, particle_list: Optional[list] = None)

Obtain/return layer,row,column idx (“node_indices”) as dict of np.arrays corresponding to xyz-coördinates of type dict per tracked particle (as key) with nodes “xyz_nodes” (np.array).

Method requires center points xmid, ymid, zmid to be predefined.

make_discretisation(schematisation: dict, dict_keys=None, model_type='axisymmetric')

Generate spatial grid for model_type choices: ‘axisymmetric’, ‘2D’ or ‘3D’.

Parameter ‘schematisation’ is of type dict with (sub)dictionaries with keys ‘dict_keys’. The subdictionaries should contain specific keyword arguments to generate the grids. The function indirectly uses “_assign_cellboundaries” to obtain grid data outputs.

# Layer data assignment (model_type: axi-symmetric | 2D | 3D) Required keys: - “bot” # bottom of local grid refinement - “top” # top of local grid refinement Optional keyword argument(s): - “nlayers” # number of layers within local grid refinement

  • self.delv: layer depths of the model layers [np.array]

    rounded to two decimals [cm scale].

  • self.zmid: z-coordinates (middle) of the model layers [np.array]

  • self.nlay: number of model layers

  • self.top: model top [float or np.array]

  • self.bot: model bottoms per layer [1D | 2D | 3D np.array of floats]

# Column data assignment (model_type: axi-symmetric | 2D | 3D) Required keys: - “xmin” # left side of local grid refinement - “xmax” # right side of local grid refinement Optional keyword argument(s): - “ncols” # number of columns within local grid refinement

# Column data outputs: - self.delr: column widths of the model columns [np.array] rounded to three decimals [mm scale]. - self.xmid: x-coordinates (middle) of the model columns [np.array] - self.ncol: number of model columns

# Row data assignment (model_type: 3D) # N.B. for modeltype: (axisymmetric | 2D) the rows have a predefined width [2 rows, 1 m width]

Required keys: - “ymin” # ‘northern’ side of local grid refinement - “ymax” # ‘southern’ side of local grid refinement Optional keyword argument(s): - “nrows” # number of rows within local grid refinement

# Row data outputs: - self.delc: row widths of the model rows [np.array] - self.ymid: y-coordinates (middle) of the model rows [np.array] - self.nrow: number of model rows

Returns

  • empty_grid (np.array) – Empty numpy array with size (self.nlay,nrow,ncol)

  • lay_bounds (np.array) – Upper and lower boundaries of the model grid cells [1D-array]

  • row_bounds (np.array) – Left and right boundaries of the model grid cells [1D-array]

  • col_bounds (np.array) – North-south boundaries of the grid-cells [1D-array]

mfmodelrun()

Main model run code.

model_cbc

Initialize the AnalyticalWell object by making the dictionaries and adding the schematisation (and therefore attributes of the schematisation) to the object.

Return type

Dictionaries (Modflow) with all parameters as an attribute of the function.

modpath_simulation(mp_model=None, trackingdirection='backward', simulationtype='combined', stoptimeoption='specified', particlegroups=None, stoptime=100000.0, zones=None)

input MODPATH Simulation File Package Class, see flopy docs.

mpbas_input(prsity=0.3, defaultiface={'RECHARGE': 6})

Read model prsity and iface values into object.

particle_data(partlocs=None, structured=True, particleids=[0], localx=None, localy=0.5, localz=None, timeoffset=0.0, drape=None, pgname='Recharge', pg_filename='Recharge.sloc', trackingdirection='forward', releasedata=0.0)

Class to create the most basic particle data type (starting location input style 1). Input style 1 is the most general input style and provides the highest flexibility in customizing starting locations, see flopy docs.

plot_age_distribution(df_particle: DataFrame, vmin=0.0, vmax=1.0, orientation={'row': 0}, fpathfig=None, figtext=None, x_text=0, y_text=0, lognorm=True, xmin=0.0, xmax=None, ymin=None, ymax=None, line_dist=1, dpi=192, trackingdirection='forward', cmap='viridis_r', show_vadose=True)

Create pathline plots with residence times using colours as indicator. with: - df_particle: dataframe containing xyzt-points of the particle paths. - fpathfig = output location of plots figtext: figure text to show within plot starting at x-location ‘x_text’ and y-location ‘y_text’. lognorm: True/False (if vmin <= 0. –> vmin = 1.e-5) xmin: left boundary figure xmax: right boundary figure ymin: lower boundary figure ymax: upper boundary figure line_dist: min. distance between well pathlines at source (m) line_freq: show every ‘X’ pathlines dpi: pixel density trackingdirection: direction of calculating flow along pathlines” cmap: Uses colormap ‘viridis_r’ (viridis reversed as default) show_vadose: T/F (if True, include vadose zone flow lines; default = True)

read_binarycbc_flow(fname=None)

Read binary cell budget file (fname). This is modflow output.

return frf, flf, fff.

read_hdsobj(fname=None, time=None)

Return head data from file. If time = -1 –> return final time and head grid, elif time = ‘all’ –> return all time values and head grids, elif time = [1.,2.,time_n]–> Return head grids for prespecified times.

read_pathlinedata(fpth, nodes)

Read pathlinedata from file fpth (extension: ‘.mppth’), given the particle release node index ‘nodes’ obtained from a tuple or list of tuples (iLay,iRow,iCol).

Parameters

  • fpth (str) – filepath of pathline data

  • nodes – Cell node indices of starting locations. Can be obtained using the method get_node_ID((iLay,iRow,iCol))

Returns

  • xyz_nodes – The xyz coördinates of each node per tracked particle

  • dist – The distance between each node per tracked particle [L]

  • tdiff – The travel time (difference) between each node per particle [T]

  • dist_tot – Total distance covered by each tracked particle [L]

  • time_tot – Total duration between release and ending of each particle [T]

  • pth_data (np.recarray) – Complete recarray is returned to see what is in the file fpth.

run_ModPathmod()

Run ModPath model

run_model(xll=0.0, yll=0.0, perlen: dict = 18250.0, nstp: dict = 1, nper: int = 1, steady: dict = True, run_mfmodel=True, run_mpmodel=True)

Run the combined modflow and modpath model using one of four possible schematisation types: - “Phreatic” - “Semi-confined” - “Recharge basin (BAR)” - “River bank filtration (RBF)” Currently (13-7-2021) only the Phreatic schematisation is supported.

run_modflowmod()

Run modflow model

write_input_mp()

Write package data ModPath model.

xyz_to_layrowcol(xyz_point, decimals=3)

obtain lay, row, col index of an xyz_point list [x,y,z].