MODPATH Particle Tracking

class Modpath7:
    """MODPATH 7 particle tracking model"""
    def __init__(
        self,
        modelname: str = 'modpathtest',
        flowmodel: object = None,
        exe_name: str = 'mp7',
        modflowmodel: object = None,
        dis_file: str = None,
        dis_unit: int = 87,
        head_file: str = None,
        budget_file: str = None,
        model_ws: str = '.',
        external_path: str = None,
        verbose: bool = False,
        load: bool = True,
        listunit: int = 7,
        **kwargs
    ): ...
    
    def write_input(self) -> None:
        """Write MODPATH 7 input files"""
        ...
    
    def run_model(
        self,
        silent: bool = False,
        pause: bool = False,
        report: bool = False,
        normal_msg: str = 'normal termination'
    ) -> tuple[bool, list[str]]:
        """Run MODPATH 7 model and return success status and output"""
        ...
    
    @classmethod
    def create_mp7(
        cls,
        modelname: str,
        trackdir: str = 'forward',
        flowmodel: object = None,
        exe_name: str = 'mp7',
        model_ws: str = '.',
        rowcelldivisions: int = 1,
        columncelldivisions: int = 1,
        layercelldivisions: int = 1,
        **kwargs
    ) -> 'Modpath7':
        """Create basic MODPATH 7 model"""
        ...
    
    def add_package(self, p: object) -> None:
        """Add package to model"""
        ...
    
    def remove_package(self, pname: str) -> None:
        """Remove package from model"""
        ...
    
    def get_package(self, name: str) -> object:
        """Get package by name"""
        ...
    
    @property
    def modelgrid(self) -> object:
        """Model grid object"""
        ...

class Modpath6:
    """MODPATH 6 particle tracking model"""
    def __init__(
        self,
        modelname: str = 'modpathtest',
        exe_name: str = 'mp6',
        modflowmodel: object = None,
        dis_file: str = None,
        head_file: str = None,
        budget_file: str = None,
        model_ws: str = '.',
        external_path: str = None,
        verbose: bool = False,
        load: bool = True,
        **kwargs
    ): ...
    
    def write_input(self) -> None:
        """Write MODPATH 6 input files"""
        ...
    
    def run_model(
        self,
        silent: bool = False,
        pause: bool = False,
        report: bool = False
    ) -> tuple[bool, list[str]]:
        """Run MODPATH 6 model"""
        ...

class Modpath7Bas:
    """MODPATH 7 basic package"""
    def __init__(
        self,
        model: Modpath7,
        hnoflo: float = 1e30,
        hdry: float = -1e30,
        def_face_ct: int = 6,
        bud_loc: int = 0,
        bud_label: str = None,
        def_iface: ArrayData = None,
        laytyp: ArrayData = None,
        ibound: ArrayData = None,
        porosity: ArrayData = 0.1,
        retfac: ArrayData = 1.0,
        izone: ArrayData = 1,
        extension: str = 'mpbas',
        **kwargs
    ): ...

class Modpath7Sim:
    """MODPATH 7 simulation package"""
    def __init__(
        self,
        model: Modpath7,
        simulationtype: str = 'combined',
        trackingdirection: str = 'forward',
        weaksinkoption: str = 'pass_through',
        weaksourceoption: str = 'pass_through',
        budgetoutputoption: str = 'no',
        traceparticledata: list = None,
        budgetcellnumbers: list = None,
        referencetime: float = 0.0,
        stoptimeoption: str = 'extend',
        stoptime: float = None,
        timepointdata: list = None,
        zonedataoption: str = 'no',
        zones: list = None,
        retardationfactoroption: str = 'no',
        retardation: list = None,
        particlegroups: list = None,
        extension: str = 'mpsim',
        **kwargs
    ): ...

class Modpath6Bas:
    """MODPATH 6 basic package"""
    def __init__(
        self,
        model: Modpath6,
        hnoflo: float = 1e30,
        hdry: float = -1e30,
        def_face_ct: int = 6,
        bud_loc: int = 0,
        bud_label: str = None,
        def_iface: ArrayData = None,
        laytyp: ArrayData = None,
        ibound: ArrayData = None,
        porosity: ArrayData = 0.1,
        extension: str = 'mpbas',
        **kwargs
    ): ...

class Modpath6Sim:
    """MODPATH 6 simulation package"""
    def __init__(
        self,
        model: Modpath6,
        option_flags: list = None,
        ref_time: float = 0.0,
        ref_time_per_stp: list = None,
        stop_time: float = None,
        group_name: list = None,
        group_placement: list = None,
        release_times: list = None,
        group_region: list = None,
        mask_nlay: list = None,
        mask_layer: list = None,
        mask_1lay: list = None,
        face_ct: list = None,
        ifaces: list = None,
        part_ct: list = None,
        time_ct: list = None,
        release_time_incr: list = None,
        time_pts: list = None,
        particle_cell_cnt: list = None,
        cell_bd_ct: list = None,
        bud_grp: list = None,
        trace_id: list = None,
        stop_zone: list = None,
        zone: list = None,
        retard_fac: list = None,
        retard_fcCB: list = None,
        sf_water_cont: list = None,
        extension: str = 'mpsim',
        **kwargs
    ): ...

class ParticleData:
    """Base particle data class"""
    def __init__(
        self,
        partlocs: ArrayData,
        structured: bool = True,
        particleids: ArrayData = None,
        **kwargs
    ): ...
    
    @property
    def particlecount(self) -> int:
        """Number of particles"""
        ...
    
    def to_coords(self, grid: object) -> list[tuple[float, float, float]]:
        """Convert particle locations to coordinates"""
        ...
    
    def write(self, f: object) -> None:
        """Write particle data to file"""
        ...

class LRCParticleData(ParticleData):
    """Layer-row-column particle data for structured grids"""
    def __init__(
        self,
        subdivisiondata: list = None,
        lrcregions: list = None,
        **kwargs
    ): ...
    
    @classmethod
    def create_subdivisiondata(
        cls,
        grid: object,
        subdivisions: tuple[int, int, int] = (1, 1, 1),
        localx: ArrayData = 0.5,
        localy: ArrayData = 0.5,
        localz: ArrayData = 0.5
    ) -> list:
        """Create subdivision data for uniform particle placement"""
        ...
    
    @classmethod  
    def create_lrcregions(
        cls,
        lrcregions: list[tuple[int, int, int, int, int, int]]
    ) -> list:
        """Create LRC regions for particle placement"""
        ...

class NodeParticleData(ParticleData):
    """Node-based particle data for unstructured grids"""
    def __init__(
        self,
        subdivisiondata: list = None,
        nodes: ArrayData = None,
        **kwargs
    ): ...

class CellDataType:
    """Cell-based particle starting locations"""
    def __init__(
        self,
        drape: int = 0,
        columncelldivisions: int = 1,
        rowcelldivisions: int = 1,
        layercelldivisions: int = 1,
        **kwargs
    ): ...

class FaceDataType:
    """Face-based particle starting locations"""
    def __init__(
        self,
        drape: int = 0,
        verticaldivisions: int = 1,
        horizontaldivisions: int = 1,
        **kwargs
    ): ...

class ParticleGroup:
    """Particle group container"""
    def __init__(
        self,
        particlegroupname: str,
        particledata: ParticleData,
        filename: str = None,
        releasedata: object = None,
        **kwargs
    ): ...
    
    @property
    def particlecount(self) -> int:
        """Total number of particles in group"""
        ...
    
    def set_releaseoption(
        self,
        option: str = 'specified',
        releasedata: list = None
    ) -> None:
        """Set particle release timing"""
        ...
    
    def write(self, f: object) -> None:
        """Write particle group to file"""
        ...

class ParticleGroupLRCTemplate:
    """LRC-based particle group template"""
    def __init__(
        self,
        particlegroupname: str,
        particledata: LRCParticleData,
        filename: str = None,
        **kwargs
    ): ...

class ParticleGroupNodeTemplate:
    """Node-based particle group template"""
    def __init__(
        self,
        particlegroupname: str,
        particledata: NodeParticleData,
        filename: str = None,
        **kwargs
    ): ...

import flopy
import numpy as np

# Create or load MODFLOW model first
mf = flopy.modflow.Modflow(modelname='flow_model')
# ... MODFLOW setup and execution ...

# Create MODPATH 7 model
mp = flopy.modpath.Modpath7(
    modelname='basic_tracking',
    flowmodel=mf,  # Link to MODFLOW model
    exe_name='mp7',
    model_ws='.'
)

# Basic package - define porosity and zones
bas = flopy.modpath.Modpath7Bas(
    mp,
    porosity=0.25,  # Effective porosity
    retfac=1.0,     # No retardation
    izone=1         # Single zone
)

# Create particle starting locations
# Place particles at well locations for capture zone analysis
nlay, nrow, ncol = mf.dis.nlay, mf.dis.nrow, mf.dis.ncol

# Define particle locations around pumping wells
particledata = []
well_locations = [(0, 25, 25), (1, 40, 60)]  # Layer, row, col

for lay, row, col in well_locations:
    # Create 3x3x1 grid of particles around each well
    for dr in [-1, 0, 1]:
        for dc in [-1, 0, 1]:
            r = max(0, min(nrow-1, row + dr))
            c = max(0, min(ncol-1, col + dc))
            particledata.append([lay, r, c, 0.5, 0.5, 0.5])

# Create particle data object
pdata = flopy.modpath.LRCParticleData(
    subdivisiondata=particledata
)

# Create particle group
pg = flopy.modpath.ParticleGroup(
    particlegroupname='wells',
    particledata=pdata,
    filename='wells.sloc'
)

# Simulation package for forward tracking
sim = flopy.modpath.Modpath7Sim(
    mp,
    simulationtype='combined',  # Both pathlines and endpoints
    trackingdirection='forward',
    weaksinkoption='pass_through',
    weaksourceoption='pass_through',
    referencetime=0.0,
    stoptimeoption='extend',
    particlegroups=[pg]
)

# Write and run MODPATH
mp.write_input()
success, buff = mp.run_model()

# Results will be in endpoint and pathline files

import flopy
import numpy as np

# MODFLOW model (abbreviated)
mf = flopy.modflow.Modflow(modelname='source_tracking')
# ... model setup ...

# MODPATH 7 for backward tracking
mp = flopy.modpath.Modpath7(
    modelname='backward_track',
    flowmodel=mf
)

# Basic package
bas = flopy.modpath.Modpath7Bas(
    mp,
    porosity=0.30,
    retfac=1.0
)

# Backward tracking from stream discharge points
# Place particles at all stream cells
stream_cells = []
# Assume stream along row 30
for col in range(20, 80):  # Stream extent
    stream_cells.append([0, 30, col, 0.5, 0.5, 0.9])  # Near top of cell

pdata = flopy.modpath.LRCParticleData(
    subdivisiondata=stream_cells
)

pg = flopy.modpath.ParticleGroup(
    particlegroupname='stream_discharge',
    particledata=pdata
)

# Backward simulation with time limit
sim = flopy.modpath.Modpath7Sim(
    mp,
    simulationtype='pathline',
    trackingdirection='backward',
    referencetime=0.0,
    stoptimeoption='specified',
    stoptime=3650.0,  # 10 years maximum travel time
    particlegroups=[pg]
)

mp.write_input()
success, buff = mp.run_model()

import flopy
import numpy as np

# Complex particle tracking scenario
mf = flopy.modflow.Modflow(modelname='complex_tracking')
# ... transient MODFLOW model setup ...

mp = flopy.modpath.Modpath7(
    modelname='complex_mp',
    flowmodel=mf
)

# Enhanced basic package with multiple zones
# Different porosity and retardation by layer
porosity_3d = np.ones((nlay, nrow, ncol)) * 0.25
porosity_3d[0, :, :] = 0.35  # Higher porosity in top layer
porosity_3d[2, :, :] = 0.15  # Lower porosity in bottom layer

retardation_3d = np.ones((nlay, nrow, ncol))
retardation_3d[1, :, :] = 2.0  # Retardation in middle layer

# Zone array for different geological units
zones = np.ones((nlay, nrow, ncol), dtype=int)
zones[0, :, :] = 1  # Alluvium
zones[1, :, :] = 2  # Clay layer
zones[2, :, :] = 3  # Bedrock

bas = flopy.modpath.Modpath7Bas(
    mp,
    porosity=porosity_3d,
    retfac=retardation_3d,
    izone=zones
)

# Multiple particle groups with different release strategies

# Group 1: Continuous release at contamination source
source_particles = []
source_row, source_col = 15, 40
# Dense particle grid at source
for i in range(5):
    for j in range(5):
        localx = 0.1 + i * 0.2
        localy = 0.1 + j * 0.2
        source_particles.append([0, source_row, source_col, localx, localy, 0.5])

pdata1 = flopy.modpath.LRCParticleData(subdivisiondata=source_particles)
pg1 = flopy.modpath.ParticleGroup(
    particlegroupname='continuous_source',
    particledata=pdata1
)

# Group 2: Periodic releases at injection wells
injection_particles = []
injection_locations = [(0, 20, 25), (0, 30, 45), (1, 25, 35)]
for lay, row, col in injection_locations:
    # Single particle per injection point
    injection_particles.append([lay, row, col, 0.5, 0.5, 0.5])

pdata2 = flopy.modpath.LRCParticleData(subdivisiondata=injection_particles)
pg2 = flopy.modpath.ParticleGroup(
    particlegroupname='periodic_injection',
    particledata=pdata2
)

# Group 3: Area source with uniform distribution
area_particles = []
# Distributed source over agricultural area
for row in range(10, 20):
    for col in range(50, 70):
        if (row - 10) % 2 == 0 and (col - 50) % 3 == 0:  # Sparse distribution
            area_particles.append([0, row, col, 0.5, 0.5, 0.8])

pdata3 = flopy.modpath.LRCParticleData(subdivisiondata=area_particles)
pg3 = flopy.modpath.ParticleGroup(
    particlegroupname='area_source',
    particledata=pdata3
)

# Complex simulation with time points and budget output
time_points = [30.0, 90.0, 365.0, 1825.0, 3650.0]  # Monthly to 10 years

sim = flopy.modpath.Modpath7Sim(
    mp,
    simulationtype='combined',
    trackingdirection='forward',
    weaksinkoption='stop',  # Stop at weak sinks
    weaksourceoption='pass_through',
    budgetoutputoption='summary',
    referencetime=0.0,
    stoptimeoption='specified',
    stoptime=7300.0,  # 20 years maximum
    timepointdata=time_points,
    zonedataoption='on',
    zones=[1, 2, 3],  # Track zone transitions
    particlegroups=[pg1, pg2, pg3]
)

mp.write_input()
success, buff = mp.run_model()

# Post-process results
try:
    # Read pathline data
    pthobj = flopy.utils.PathlineFile(mp.model_ws + '/complex_mp.mppth')
    pathlines = pthobj.get_alldata()
    
    # Read endpoint data
    epobj = flopy.utils.EndpointFile(mp.model_ws + '/complex_mp.mpend')
    endpoints = epobj.get_alldata()
    
    print(f"Number of pathlines: {len(pathlines)}")
    print(f"Number of endpoints: {len(endpoints)}")
    
except:
    print("Could not read MODPATH output files")

import flopy
import numpy as np

# High-resolution particle tracking in specific zones
mf = flopy.modflow.Modflow(modelname='detailed_tracking')
# ... MODFLOW setup ...

mp = flopy.modpath.Modpath7(
    modelname='detailed_mp',
    flowmodel=mf
)

bas = flopy.modpath.Modpath7Bas(mp, porosity=0.25)

# Create detailed subdivision data for high-density particle placement
# Method 1: Uniform subdivision across selected cells
nlay, nrow, ncol = mf.dis.nlay, mf.dis.nrow, mf.dis.ncol

# Define region of interest (contamination plume area)
roi_layers = [0, 1]  # Top two layers
roi_rows = range(20, 40)  # Central portion
roi_cols = range(30, 70)  # Extended width

subdivisiondata = []
for lay in roi_layers:
    for row in roi_rows:
        for col in roi_cols:
            # 4x4x2 subdivision in each cell (32 particles per cell)
            for k_div in range(2):
                for i_div in range(4):
                    for j_div in range(4):
                        localx = (j_div + 0.5) / 4.0
                        localy = (i_div + 0.5) / 4.0
                        localz = (k_div + 0.5) / 2.0
                        subdivisiondata.append([lay, row, col, localx, localy, localz])

# Method 2: Using built-in subdivision functionality
subdivision_specs = []
for lay in roi_layers:
    for row in roi_rows:
        for col in roi_cols:
            # Each cell subdivided into 3x3x2 = 18 particles
            subdivision_specs.append([lay, row, col, 3, 3, 2])

# Create particle data with manual subdivision
pdata_manual = flopy.modpath.LRCParticleData(
    subdivisiondata=subdivisiondata
)

# Create using built-in subdivision
grid = mf.modelgrid
pdata_auto = flopy.modpath.LRCParticleData.create_subdivisiondata(
    grid=grid,
    subdivisions=(2, 3, 3),  # z, y, x subdivisions
    localx=np.array([0.25, 0.75]),  # Two columns
    localy=np.array([0.2, 0.5, 0.8]),  # Three rows
    localz=np.array([0.25, 0.75])  # Two layers
)

# Create particle groups
pg_manual = flopy.modpath.ParticleGroup(
    particlegroupname='detailed_manual',
    particledata=pdata_manual
)

# Simulation for detailed analysis
sim = flopy.modpath.Modpath7Sim(
    mp,
    simulationtype='pathline',
    trackingdirection='forward',
    referencetime=0.0,
    stoptimeoption='specified',
    stoptime=1825.0,  # 5 years
    particlegroups=[pg_manual]
)

mp.write_input()
success, buff = mp.run_model()

print(f"Total particles tracked: {pdata_manual.particlecount}")

import flopy
import numpy as np

# MODFLOW-USG or MODFLOW 6 unstructured model
# Example assumes MODFLOW 6 with DISV
sim = flopy.mf6.MFSimulation(sim_name='unstructured_sim')
# ... simulation setup ...

gwf = flopy.mf6.ModflowGwf(sim, modelname='gwf_unstructured')
# ... unstructured grid model setup with DISV ...

# MODPATH 7 for unstructured grid
mp = flopy.modpath.Modpath7(
    modelname='unstructured_mp',
    flowmodel=gwf
)

# Basic package for unstructured grid
ncpl = gwf.disv.ncpl.array  # Cells per layer
nlay = gwf.disv.nlay.array
porosity = np.ones(ncpl * nlay) * 0.25

bas = flopy.modpath.Modpath7Bas(
    mp,
    porosity=porosity
)

# Node-based particle data for unstructured grid
# Select nodes around specific features
source_nodes = []
target_nodes = []

# Assume contamination source at specific nodes
contamination_nodes = [150, 151, 152, 180, 181, 182]  # Node numbers
for node in contamination_nodes:
    # Multiple particles per node with different local positions
    for i in range(9):  # 3x3 pattern
        localx = 0.2 + (i % 3) * 0.3
        localy = 0.2 + (i // 3) * 0.3
        localz = 0.5
        source_nodes.append([node, localx, localy, localz])

# Discharge area nodes for backward tracking
discharge_nodes = [800, 801, 802, 820, 821, 822]
for node in discharge_nodes:
    target_nodes.append([node, 0.5, 0.5, 0.9])  # Near cell top

# Create node particle data
pdata_source = flopy.modpath.NodeParticleData(
    subdivisiondata=source_nodes
)

pdata_target = flopy.modpath.NodeParticleData(
    subdivisiondata=target_nodes
)

# Create particle groups
pg_forward = flopy.modpath.ParticleGroup(
    particlegroupname='source_forward',
    particledata=pdata_source
)

pg_backward = flopy.modpath.ParticleGroup(
    particlegroupname='target_backward',
    particledata=pdata_target
)

# Simulation with both forward and backward tracking
sim_mp = flopy.modpath.Modpath7Sim(
    mp,
    simulationtype='combined',
    trackingdirection='forward',  # Will be overridden per group
    particlegroups=[pg_forward, pg_backward]
)

mp.write_input()
success, buff = mp.run_model()

# Particle tracking types
TrackingDirection = Literal['forward', 'backward']
SimulationType = Literal['endpoint', 'pathline', 'timeseries', 'combined']
SinkSourceOption = Literal['pass_through', 'stop']
StopTimeOption = Literal['extend', 'specified']

# Particle location data
ParticleLocation = list[Union[int, float]]  # [lay/node, row, col, localx, localy, localz]
SubdivisionData = list[ParticleLocation]
RegionData = list[tuple[int, int, int, int, int, int]]  # (lay1, row1, col1, lay2, row2, col2)

# Grid and array types
ArrayData = Union[float, list[float], np.ndarray]
GridType = Literal['structured', 'unstructured', 'vertex']
NodeNumbers = list[int]

# Time and release data
TimeData = list[float]
ReleaseOption = Literal['specified', 'periodic', 'continuous']
ReleaseData = Union[float, list[float]]

# Output options
BudgetOutput = Literal['no', 'summary', 'record_summary']
ZoneOption = Literal['no', 'on']
RetardationOption = Literal['no', 'on']

# File types
ParticleFile = str
PathlineFile = str
EndpointFile = str
TimeseriesFile = str