MODFLOW 2005 Models and Packages

class Modflow:
    """MODFLOW 2005 model container and execution manager"""
    def __init__(
        self,
        modelname: str = 'modflowtest',
        namefile_ext: str = 'nam',
        version: str = 'mf2005',
        exe_name: str = 'mf2005',
        structured: bool = True,
        listunit: int = 2,
        model_ws: str = '.',
        external_path: str = None,
        verbose: bool = False,
        **kwargs
    ): ...
    
    def write_input(
        self,
        SelPackList: list = None,
        check: bool = True
    ) -> None:
        """Write MODFLOW 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 the MODFLOW model and return success status and output"""
        ...
    
    def load(
        cls,
        f: str,
        version: str = 'mf2005',
        exe_name: str = None,
        verbose: bool = False,
        model_ws: str = '.',
        load_only: list = None,
        forgive: bool = False,
        check: bool = True
    ) -> 'Modflow':
        """Load existing MODFLOW model from files"""
        ...
    
    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"""
        ...
    
    def get_package_list(self, ftype: str = None) -> list[str]:
        """Get list of packages"""
        ...
    
    def set_model_units(self, iunit: dict = None) -> None:
        """Set model unit numbers"""
        ...
    
    @property
    def modelgrid(self) -> object:
        """Model grid object"""
        ...
    
    @property
    def packages(self) -> list[object]:
        """List of model packages"""
        ...
    
    @property
    def packagelist(self) -> list[str]:
        """List of package names"""
        ...

class ModflowBas:
    """MODFLOW basic package"""
    def __init__(
        self,
        model: Modflow,
        ibound: ArrayData = 1,
        strt: ArrayData = 1.0,
        ifrefm: bool = True,
        ixsec: bool = False,
        ichflg: bool = False,
        stoper: float = None,
        hnoflo: float = 999.0,
        extension: str = 'bas',
        unitnumber: int = None,
        filenames: str = None,
        **kwargs
    ): ...

class ModflowDis:
    """MODFLOW discretization package"""
    def __init__(
        self,
        model: Modflow,
        nlay: int = 1,
        nrow: int = 2,
        ncol: int = 2,
        nper: int = 1,
        delr: ArrayData = 1.0,
        delc: ArrayData = 1.0,
        laycbd: ArrayData = 0,
        top: ArrayData = 1.0,
        botm: ArrayData = 0.0,
        perlen: ArrayData = 1.0,
        nstp: ArrayData = 1,
        tsmult: ArrayData = 1.0,
        steady: ArrayData = True,
        itmuni: int = 4,
        lenuni: int = 2,
        extension: str = 'dis',
        unitnumber: int = None,
        filenames: str = None,
        **kwargs
    ): ...

class ModflowOc:
    """MODFLOW output control package"""
    def __init__(
        self,
        model: Modflow,
        ihedfm: int = 0,
        iddnfm: int = 0,
        chedfm: str = None,
        cddnfm: str = None,
        cboufm: str = None,
        compact: bool = True,
        stress_period_data: dict = None,
        extension: str = 'oc',
        unitnumber: int = None,
        filenames: str = None,
        **kwargs
    ): ...

class ModflowLpf:
    """MODFLOW layer property flow package"""
    def __init__(
        self,
        model: Modflow,
        laytyp: ArrayData = 0,
        layavg: ArrayData = 0,
        chani: ArrayData = 1.0,
        layvka: ArrayData = 0,
        laywet: ArrayData = 0,
        ipakcb: int = None,
        hdry: float = -1e30,
        iwdflg: int = 0,
        wetfct: float = 0.1,
        iwetit: int = 1,
        ihdwet: int = 0,
        hk: ArrayData = 1.0,
        hani: ArrayData = 1.0,
        vka: ArrayData = 1e-6,
        ss: ArrayData = 1e-5,
        sy: ArrayData = 0.15,
        vkcb: ArrayData = 0.0,
        wetdry: ArrayData = -0.01,
        storagecoefficient: bool = False,
        constantcv: bool = False,
        thickstrt: bool = False,
        nocvcorrection: bool = False,
        novfc: bool = False,
        extension: str = 'lpf',
        unitnumber: int = None,
        filenames: str = None,
        **kwargs
    ): ...

class ModflowBcf:
    """MODFLOW block-centered flow package"""
    def __init__(
        self,
        model: Modflow,
        ipakcb: int = None,
        intercellt: int = 0,
        laycon: ArrayData = 3,
        trpy: ArrayData = 1.0,
        hdry: float = -1e30,
        iwdflg: int = 0,
        wetfct: float = 0.1,
        iwetit: int = 1,
        ihdwet: int = 0,
        tran: ArrayData = 1.0,
        hy: ArrayData = 1.0,
        vcont: ArrayData = 1.0,
        sf1: ArrayData = 1e-5,
        sf2: ArrayData = 0.15,
        wetdry: ArrayData = -0.01,
        extension: str = 'bcf',
        unitnumber: int = None,
        filenames: str = None,
        **kwargs
    ): ...

class ModflowUpw:
    """MODFLOW upstream weighting package"""
    def __init__(
        self,
        model: Modflow,
        laytyp: ArrayData = 0,
        layavg: ArrayData = 0,
        chani: ArrayData = 1.0,
        layvka: ArrayData = 0,
        laywet: ArrayData = 0,
        ipakcb: int = None,
        hdry: float = -1e30,
        npupw: int = 1,
        iphdry: int = 0,
        hk: ArrayData = 1.0,
        hani: ArrayData = 1.0,
        vka: ArrayData = 1e-6,
        ss: ArrayData = 1e-5,
        sy: ArrayData = 0.15,
        vkcb: ArrayData = 0.0,
        wetdry: ArrayData = -0.01,
        extension: str = 'upw',
        unitnumber: int = None,
        filenames: str = None,
        **kwargs
    ): ...

class ModflowHfb:
    """MODFLOW horizontal flow barriers package"""
    def __init__(
        self,
        model: Modflow,
        nhfbnp: int = 0,
        mxfb: int = 0,
        nhfbx1: int = 0,
        nhfbx2: int = 0,
        hfb_data: list = None,
        nactp: int = 0,
        nactcbhfb: int = 0,
        no_print: bool = False,
        options: list = None,
        extension: str = 'hfb',
        unitnumber: int = None,
        filenames: str = None,
        **kwargs
    ): ...

class ModflowWel:
    """MODFLOW well package"""
    def __init__(
        self,
        model: Modflow,
        ipakcb: int = None,
        stress_period_data: StressPeriodData = None,
        dtype: np.dtype = None,
        extension: str = 'wel',
        unitnumber: int = None,
        filenames: str = None,
        options: list = None,
        binary: bool = False,
        **kwargs
    ): ...

class ModflowRch:
    """MODFLOW recharge package"""
    def __init__(
        self,
        model: Modflow,
        nrchop: int = 3,
        ipakcb: int = None,
        rech: ArrayData = 0.001,
        irch: ArrayData = 0,
        extension: str = 'rch',
        unitnumber: int = None,
        filenames: str = None,
        **kwargs
    ): ...

class ModflowEvt:
    """MODFLOW evapotranspiration package"""
    def __init__(
        self,
        model: Modflow,
        nevtop: int = 3,
        ipakcb: int = None,
        surf: ArrayData = 0.0,
        evtr: ArrayData = 0.0004,
        exdp: ArrayData = 1.0,
        ievt: ArrayData = 0,
        extension: str = 'evt',
        unitnumber: int = None,
        filenames: str = None,
        **kwargs
    ): ...

class ModflowGhb:
    """MODFLOW general head boundary package"""
    def __init__(
        self,
        model: Modflow,
        ipakcb: int = None,
        stress_period_data: StressPeriodData = None,
        dtype: np.dtype = None,
        no_print: bool = False,
        options: list = None,
        extension: str = 'ghb',
        unitnumber: int = None,
        filenames: str = None,
        binary: bool = False,
        **kwargs
    ): ...

class ModflowRiv:
    """MODFLOW river package"""
    def __init__(
        self,
        model: Modflow,
        ipakcb: int = None,
        stress_period_data: StressPeriodData = None,
        dtype: np.dtype = None,
        extension: str = 'riv',
        unitnumber: int = None,
        filenames: str = None,
        options: list = None,
        binary: bool = False,
        **kwargs
    ): ...

class ModflowDrn:
    """MODFLOW drain package"""
    def __init__(
        self,
        model: Modflow,
        ipakcb: int = None,
        stress_period_data: StressPeriodData = None,
        dtype: np.dtype = None,
        extension: str = 'drn',
        unitnumber: int = None,
        filenames: str = None,
        options: list = None,
        binary: bool = False,
        **kwargs
    ): ...

class ModflowChd:
    """MODFLOW constant head package"""
    def __init__(
        self,
        model: Modflow,
        stress_period_data: StressPeriodData = None,
        dtype: np.dtype = None,
        extension: str = 'chd',
        unitnumber: int = None,
        filenames: str = None,
        options: list = None,
        **kwargs
    ): ...

class ModflowLak:
    """MODFLOW lake package"""
    def __init__(
        self,
        model: Modflow,
        nlakes: int = 1,
        ipakcb: int = None,
        theta: float = -1.0,
        nssitr: int = 0,
        sscncr: float = 0.0,
        surfdepth: float = 0.0,
        stages: ArrayData = 1.0,
        stage_range: ArrayData = None,
        lakarr: ArrayData = None,
        bdlknc: ArrayData = None,
        sill_data: dict = None,
        flux_data: dict = None,
        extension: str = 'lak',
        unitnumber: int = None,
        filenames: str = None,
        **kwargs
    ): ...

class ModflowSfr2:
    """MODFLOW streamflow routing package"""
    def __init__(
        self,
        model: Modflow,
        nstrm: int = 0,
        nss: int = 0,
        const: float = 86400.0,
        dleak: float = 0.0001,
        ipakcb: int = None,
        istcb2: int = None,
        isfropt: int = 1,
        nstrail: int = 10,
        isuzn: int = 1,
        nsfrsets: int = 30,
        irtflg: int = 1,
        numtim: int = 2,
        weight: float = 0.75,
        flwtol: float = 0.0001,
        reach_data: np.recarray = None,
        segment_data: dict = None,
        channel_geometry_data: dict = None,
        channel_flow_data: dict = None,
        dataset_5: dict = None,
        reachinput: bool = False,
        transroute: bool = False,
        tabfiles: bool = False,
        tabfiles_dict: dict = None,
        extension: str = 'sfr',
        unitnumber: int = None,
        filenames: str = None,
        **kwargs
    ): ...

class ModflowUzf1:
    """MODFLOW unsaturated zone flow package"""
    def __init__(
        self,
        model: Modflow,
        nuztop: int = 1,
        iuzfopt: int = 0,
        irunflg: int = 0,
        ietflg: int = 0,
        ipakcb: int = None,
        iuzfcb2: int = None,
        ntrail2: int = 10,
        nsets2: int = 20,
        nuzgag: int = 0,
        surfdep: float = 1.0,
        iuzfbnd: ArrayData = 1,
        irunbnd: ArrayData = 0,
        vks: ArrayData = 1e-6,
        eps: ArrayData = 3.5,
        thts: ArrayData = 0.35,
        thtr: ArrayData = 0.15,
        thti: ArrayData = 0.20,
        row_col_iftunit_iuzopt: list = None,
        extension: str = 'uzf',
        unitnumber: int = None,
        filenames: str = None,
        **kwargs
    ): ...

class ModflowSub:
    """MODFLOW subsidence package"""
    def __init__(
        self,
        model: Modflow,
        ipakcb: int = None,
        isuboc: int = 0,
        acc: float = 1e-5,
        itmin: int = 5,
        idsave: int = None,
        idrest: int = None,
        nndb: int = 1,
        ndb: int = 1,
        nmz: int = 1,
        nn: int = 20,
        ac1: float = 0.0,
        ac2: float = 1.0,
        itmin_sub: int = 5,
        ln: ArrayData = 0,
        rnb: ArrayData = 1.0,
        hc: ArrayData = 1.0,
        sfe: ArrayData = 1e-4,
        sfv: ArrayData = 1e-3,
        com: ArrayData = 1e-3,
        dp: list = None,
        dz: list = None,
        dhc: list = None,
        dcom: list = None,
        dz_mult: list = None,
        ids15: list = None,
        ids16: list = None,
        extension: str = 'sub',
        unitnumber: int = None,
        filenames: str = None,
        **kwargs
    ): ...

class ModflowSwi2:
    """MODFLOW seawater intrusion package"""
    def __init__(
        self,
        model: Modflow,
        nsrf: int = 1,
        istrat: int = 1,
        nobs: int = 0,
        iswizt: int = None,
        ipakcb: int = None,
        iswibd: int = 0,
        ihcofadd: int = 0,
        izcfl: int = 0,
        izcfm: int = 0,
        iglfl: int = 0,
        iglfm: int = 0,
        toeslope: float = 0.05,
        tipslope: float = 0.05,
        alpha: float = None,
        beta: float = 0.1,
        nadptmx: int = 1,
        nadptmn: int = 1,
        adptfct: float = 1.0,
        nu: ArrayData = 0.025,
        zeta: ArrayData = None,
        ssz: ArrayData = 0.25,
        isource: ArrayData = 0,
        obsnam: list = None,
        obslrc: list = None,
        extension: str = 'swi2',
        unitnumber: int = None,
        filenames: str = None,
        **kwargs
    ): ...

class ModflowPcg:
    """MODFLOW preconditioned conjugate gradient solver"""
    def __init__(
        self,
        model: Modflow,
        mxiter: int = 50,
        iter1: int = 30,
        npcond: int = 1,
        hclose: float = 1e-5,
        rclose: float = 1e-5,
        relax: float = 1.0,
        nbpol: int = 0,
        iprpcg: int = 0,
        mutpcg: int = 3,
        damp: float = 1.0,
        dampt: float = 1.0,
        ihcofadd: int = 0,
        extension: str = 'pcg',
        unitnumber: int = None,
        filenames: str = None,
        **kwargs
    ): ...

class ModflowSip:
    """MODFLOW strongly implicit procedure solver"""
    def __init__(
        self,
        model: Modflow,
        mxiter: int = 200,
        nparm: int = 5,
        accl: float = 1.0,
        hclose: float = 1e-5,
        ipcalc: int = 1,
        wseed: float = 0.0,
        iprsip: int = 0,
        extension: str = 'sip',
        unitnumber: int = None,
        filenames: str = None,
        **kwargs
    ): ...

class ModflowDe4:
    """MODFLOW direct solver"""
    def __init__(
        self,
        model: Modflow,
        itmx: int = 50,
        mxup: int = 0,
        mxlow: int = 0,
        mxbw: int = 0,
        ifreq: int = 3,
        mutd4: int = 0,
        accl: float = 1.0,
        hclose: float = 1e-5,
        iprd4: int = 1,
        extension: str = 'de4',
        unitnumber: int = None,
        filenames: str = None,
        **kwargs
    ): ...

class ModflowNwt:
    """MODFLOW Newton solver"""
    def __init__(
        self,
        model: Modflow,
        headtol: float = 0.01,
        fluxtol: float = 500.0,
        maxiterout: int = 100,
        thickfact: float = 1e-5,
        linmeth: int = 1,
        iprnwt: int = 0,
        ibotav: int = 0,
        options: str = 'COMPLEX',
        Continue: bool = False,
        dbdtheta: float = 0.4,
        dbdkappa: float = 1e-5,
        dbdgamma: float = 0.0,
        momfact: float = 0.1,
        backflag: int = 1,
        maxbackiter: int = 50,
        backtol: float = 1.1,
        backreduce: float = 0.7,
        maxitinner: int = 50,
        ilumethod: int = 2,
        levfill: int = 5,
        stoptol: float = 1e-10,
        msdr: int = 15,
        iacl: int = 2,
        norder: int = 1,
        level: int = 5,
        north: int = 7,
        iredsys: int = 0,
        rrctols: float = 0.0,
        idroptol: int = 1,
        epsrn: float = 1e-4,
        hclosexmd: float = 1e-4,
        mxiterxmd: int = 50,
        extension: str = 'nwt',
        unitnumber: int = None,
        filenames: str = None,
        **kwargs
    ): ...

class ModflowGmg:
    """MODFLOW geometric multigrid solver"""
    def __init__(
        self,
        model: Modflow,
        mxiter: int = 50,
        iiter: int = 20,
        iadamp: int = 3,
        hclose: float = 1e-5,
        rclose: float = 1e-5,
        relax: float = 0.8,
        ioutgmg: int = 0,
        iunitmhc: int = 0,
        ism: int = 0,
        isc: int = 0,
        damp: float = 1.0,
        dup: float = 0.75,
        dlow: float = 0.01,
        chglimit: float = 0.5,
        extension: str = 'gmg',
        unitnumber: int = None,
        filenames: str = None,
        **kwargs
    ): ...

class ModflowHob:
    """MODFLOW head observation package"""
    def __init__(
        self,
        model: Modflow = None,
        iuhobsv: int = 1051,
        hobdry: float = 0,
        tomulth: float = 1.0,
        obs_data: list = None,
        hobname: list = None,
        extension: str = 'hob',
        no_print: bool = False,
        options: list = None,
        unitnumber: int = None,
        filenames: str = None,
        **kwargs
    ): ...

class ModflowFlwob:
    """MODFLOW flow observation package"""
    def __init__(
        self,
        model: Modflow,
        nqfb: int = 0,
        nqcfb: int = 0,
        nqtfb: int = 0,
        iufbobsv: int = 1052,
        tomultfb: float = 1.0,
        nqobfb: list = None,
        nqclfb: list = None,
        obsnam: list = None,
        irefsp: list = None,
        toffset: list = None,
        flwsim: list = None,
        flwobs: list = None,
        extension: str = 'flwob',
        unitnumber: int = None,
        filenames: str = None,
        **kwargs
    ): ...

import flopy
import numpy as np

# Create model
mf = flopy.modflow.Modflow(modelname='basic_model', exe_name='mf2005')

# Add discretization
nlay, nrow, ncol = 3, 50, 50
dis = flopy.modflow.ModflowDis(
    mf,
    nlay=nlay, nrow=nrow, ncol=ncol,
    delr=100.0, delc=100.0,
    top=100.0,
    botm=[50.0, 0.0, -50.0],
    nper=2,
    perlen=[1.0, 365.0],
    nstp=[1, 12],
    steady=[True, False]
)

# Add basic package
ibound = np.ones((nlay, nrow, ncol), dtype=int)
ibound[:, 0, :] = -1  # Constant head boundary
ibound[:, -1, :] = -1
bas = flopy.modflow.ModflowBas(
    mf,
    ibound=ibound,
    strt=75.0
)

# Add layer property flow
hk = [10.0, 5.0, 1.0]  # Different K for each layer
lpf = flopy.modflow.ModflowLpf(
    mf,
    hk=hk,
    vka=0.1,
    sy=0.2,
    ss=1e-5,
    laytyp=[1, 1, 0]  # Top two layers convertible
)

# Add recharge
rch = flopy.modflow.ModflowRch(mf, rech=0.001)

# Add wells
wel_sp_data = {
    0: [[0, 25, 25, -1000.0]],  # Steady state pumping well
    1: [[0, 25, 25, -2000.0]]   # Increased pumping in transient
}
wel = flopy.modflow.ModflowWel(mf, stress_period_data=wel_sp_data)

# Add output control
oc = flopy.modflow.ModflowOc(mf)

# Add solver
pcg = flopy.modflow.ModflowPcg(mf)

# Write input files
mf.write_input()

# Run model
success, buff = mf.run_model()

import flopy
import numpy as np

# Create regional model
mf = flopy.modflow.Modflow(
    modelname='regional_model',
    exe_name='mf2005',
    version='mf2005'
)

# Large grid for regional analysis
nlay, nrow, ncol = 5, 200, 300
dis = flopy.modflow.ModflowDis(
    mf,
    nlay=nlay, nrow=nrow, ncol=ncol,
    delr=500.0, delc=500.0,
    top=200.0,
    botm=[150.0, 100.0, 50.0, 0.0, -100.0],
    nper=12,  # Monthly stress periods
    perlen=30.25,  # Average month length
    nstp=5,
    tsmult=1.2,
    steady=False
)

# Basic package with complex boundary conditions
ibound = np.ones((nlay, nrow, ncol), dtype=int)
# Set no-flow boundaries
ibound[:, 0, :] = 0
ibound[:, -1, :] = 0
ibound[:, :, 0] = 0
# Set constant head boundaries along one edge
ibound[:, :, -1] = -1

bas = flopy.modflow.ModflowBas(
    mf,
    ibound=ibound,
    strt=150.0
)

# Layer property flow with heterogeneous properties
hk_layers = [50.0, 20.0, 10.0, 5.0, 1.0]
lpf = flopy.modflow.ModflowLpf(
    mf,
    hk=hk_layers,
    vka=[5.0, 2.0, 1.0, 0.5, 0.1],
    sy=0.25,
    ss=1e-5,
    laytyp=1,
    hdry=-999.0
)

# Seasonal recharge pattern
rech_data = {}
for kper in range(12):
    # Higher recharge in winter months (Nov-Mar)
    if kper in [10, 11, 0, 1, 2]:
        rech_rate = 0.005
    else:
        rech_rate = 0.001
    rech_data[kper] = rech_rate

rch = flopy.modflow.ModflowRch(mf, rech=rech_data)

# Multiple rivers
riv_data = {}
for kper in range(12):
    # River stages vary seasonally
    stage_mult = 1.0 + 0.2 * np.sin(2 * np.pi * kper / 12)
    rivers = []
    # Main river along row 50
    for col in range(50, 250):
        rivers.append([0, 50, col, 145.0 * stage_mult, 1000.0, 140.0])
    # Tributary along column 100
    for row in range(60, 120):
        rivers.append([0, row, 100, 140.0 * stage_mult, 500.0, 135.0])
    riv_data[kper] = rivers

riv = flopy.modflow.ModflowRiv(mf, stress_period_data=riv_data)

# General head boundaries for regional flow
ghb_data = {}
for kper in range(12):
    ghbs = []
    # Eastern boundary with regional gradient
    for row in range(nrow):
        head = 120.0 + (row / nrow) * 30.0  # North-south gradient
        ghbs.append([0, row, ncol-1, head, 100.0])
    ghb_data[kper] = ghbs

ghb = flopy.modflow.ModflowGhb(mf, stress_period_data=ghb_data)

# Municipal wells with seasonal pumping
wel_data = {}
well_locations = [(0, 75, 150), (1, 125, 200), (0, 160, 100)]
for kper in range(12):
    wells = []
    # Higher pumping in summer months
    if kper in [5, 6, 7, 8]:
        pump_mult = 1.5
    else:
        pump_mult = 1.0
    
    base_rates = [-5000.0, -3000.0, -2000.0]
    for i, (lay, row, col) in enumerate(well_locations):
        rate = base_rates[i] * pump_mult
        wells.append([lay, row, col, rate])
    wel_data[kper] = wells

wel = flopy.modflow.ModflowWel(mf, stress_period_data=wel_data)

# Output control for detailed output
stress_period_data = {}
for kper in range(12):
    stress_period_data[kper] = ['save head', 'save budget']

oc = flopy.modflow.ModflowOc(
    mf,
    stress_period_data=stress_period_data,
    compact=True
)

# Newton solver for better convergence
nwt = flopy.modflow.ModflowNwt(
    mf,
    headtol=0.01,
    fluxtol=500.0,
    maxiterout=100,
    linmeth=2  # GMRES
)

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

# MODFLOW 2005 specific types
ArrayData = Union[int, float, np.ndarray, list, str]
StressPeriodData = dict[int, list[list]]

# Boundary condition data formats
WelData = list[list[int, int, int, float]]  # [layer, row, col, flux]
RivData = list[list[int, int, int, float, float, float]]  # [layer, row, col, stage, cond, rbot]
GhbData = list[list[int, int, int, float, float]]  # [layer, row, col, head, cond]
DrnData = list[list[int, int, int, float, float]]  # [layer, row, col, elev, cond]
ChdData = list[list[int, int, int, float, float]]  # [layer, row, col, shead, ehead]

# Discretization arrays
DelrDelc = Union[float, list[float], np.ndarray]
Botm = Union[float, list[float], np.ndarray]
Perlen = Union[float, list[float], np.ndarray]

# Package unit numbers
UnitNumber = Union[int, None]