0
# MT3DMS Transport Modeling
1

2
This module provides comprehensive support for MT3DMS and MT3D-USGS transport models with advection, dispersion, reactions, and specialized packages for complex transport scenarios. MT3DMS is the multi-species transport model for simulating contaminant fate and transport in groundwater systems.
3

4
## Core Model Class
5

6
### Mt3dms
7

8
Main MT3DMS model class that serves as a container for transport packages.
9

10
```python { .api }
11
class Mt3dms:
12
    """MT3DMS transport model container and execution manager"""
13
    def __init__(
14
        self,
15
        modelname: str = 'mt3dtest',
16
        namefile_ext: str = 'nam',
17
        modflowmodel: object = None,
18
        ftlfilename: str = 'mt3d_link.ftl',
19
        ftlfree: bool = False,
20
        version: str = 'mt3dms',
21
        exe_name: str = 'mt3dms5',
22
        structured: bool = True,
23
        listunit: int = 16,
24
        ftlunit: int = 10,
25
        model_ws: str = '.',
26
        external_path: str = None,
27
        verbose: bool = False,
28
        load: bool = True,
29
        silent: int = 0,
30
        **kwargs
31
    ): ...
32
    
33
    def write_input(self, SelPackList: list = None, check: bool = True) -> None:
34
        """Write MT3DMS input files"""
35
        ...
36
    
37
    def run_model(
38
        self,
39
        silent: bool = False,
40
        pause: bool = False,
41
        report: bool = False,
42
        normal_msg: str = 'program completed'
43
    ) -> tuple[bool, list[str]]:
44
        """Run the MT3DMS model and return success status and output"""
45
        ...
46
    
47
    @classmethod
48
    def load(
49
        cls,
50
        f: str,
51
        version: str = 'mt3dms',
52
        exe_name: str = 'mt3dms5',
53
        verbose: bool = False,
54
        model_ws: str = '.',
55
        load_only: list = None,
56
        forgive: bool = False,
57
        check: bool = True
58
    ) -> 'Mt3dms':
59
        """Load existing MT3DMS model from files"""
60
        ...
61
    
62
    def add_package(self, p: object) -> None:
63
        """Add package to model"""
64
        ...
65
    
66
    def remove_package(self, pname: str) -> None:
67
        """Remove package from model"""
68
        ...
69
    
70
    def get_package(self, name: str) -> object:
71
        """Get package by name"""
72
        ...
73
    
74
    @property
75
    def modelgrid(self) -> object:
76
        """Model grid object"""
77
        ...
78
    
79
    @property
80
    def packages(self) -> list[object]:
81
        """List of model packages"""
82
        ...
83
    
84
    @property
85
    def ncomp(self) -> int:
86
        """Number of chemical components"""
87
        ...
88
    
89
    @property
90
    def nlay(self) -> int:
91
        """Number of layers"""
92
        ...
93
    
94
    @property
95
    def nrow(self) -> int:
96
        """Number of rows"""
97
        ...
98
    
99
    @property
100
    def ncol(self) -> int:
101
        """Number of columns"""
102
        ...
103
```
104

105
## Core Transport Packages
106

107
### Mt3dBtn
108

109
Basic transport package that defines fundamental transport parameters and component setup.
110

111
```python { .api }
112
class Mt3dBtn:
113
    """MT3DMS basic transport package"""
114
    def __init__(
115
        self,
116
        model: Mt3dms,
117
        nlay: int = 1,
118
        nrow: int = 2,
119
        ncol: int = 2,
120
        nper: int = 1,
121
        ncomp: int = 1,
122
        mcomp: int = 1,
123
        tunit: str = 'D',
124
        lunit: str = 'M',
125
        munit: str = 'KG',
126
        laycon: ArrayData = 0,
127
        delr: ArrayData = 1.0,
128
        delc: ArrayData = 1.0,
129
        htop: ArrayData = 1.0,
130
        dz: ArrayData = 1.0,
131
        prsity: ArrayData = 0.30,
132
        icbund: ArrayData = 1,
133
        sconc: ArrayData = 0.0,
134
        cinact: float = -1e30,
135
        thkmin: float = 0.01,
136
        ifmtcn: int = 0,
137
        ifmtnp: int = 0,
138
        ifmtrf: int = 0,
139
        ifmtdp: int = 0,
140
        savucn: bool = True,
141
        nprs: int = 0,
142
        timprs: ArrayData = None,
143
        obs: ArrayData = None,
144
        nprobs: int = 1,
145
        chkmas: bool = True,
146
        nprmas: int = 1,
147
        dt0: float = 0,
148
        mxstrn: int = 50000,
149
        ttsmult: float = 1.0,
150
        ttsmax: float = 0,
151
        species_names: list[str] = None,
152
        extension: str = 'btn',
153
        unitnumber: int = None,
154
        filenames: str = None,
155
        **kwargs
156
    ): ...
157
```
158

159
**Parameters:**
160
- `ncomp` (int): Number of mobile species
161
- `mcomp` (int): Total number of species (mobile + immobile)
162
- `tunit` (str): Time unit ('S', 'M', 'H', 'D', 'Y')
163
- `lunit` (str): Length unit ('CM', 'M', 'KM', 'FT', 'MI')
164
- `munit` (str): Mass unit ('G', 'KG', 'LB', 'TON')
165
- `prsity` (ArrayData): Effective porosity
166
- `icbund` (ArrayData): Boundary condition array
167
- `sconc` (ArrayData): Starting concentration
168

169
## Process Packages
170

171
### Mt3dAdv
172

173
Advection package for simulating solute transport by groundwater flow.
174

175
```python { .api }
176
class Mt3dAdv:
177
    """MT3DMS advection package"""
178
    def __init__(
179
        self,
180
        model: Mt3dms,
181
        mixelm: int = 3,
182
        percel: float = 0.75,
183
        mxpart: int = 800000,
184
        nadvfd: int = 1,
185
        itrack: int = 3,
186
        wd: float = 0.5,
187
        dceps: float = 1e-5,
188
        nplane: int = 2,
189
        npl: int = 10,
190
        nph: int = 40,
191
        npmin: int = 5,
192
        npmax: int = 80,
193
        nlsink: int = None,
194
        npsink: int = None,
195
        dchmoc: float = 1e-3,
196
        extension: str = 'adv',
197
        unitnumber: int = None,
198
        filenames: str = None,
199
        **kwargs
200
    ): ...
201
```
202

203
**Parameters:**
204
- `mixelm` (int): Advection solution method
205
  - -1: Implicit finite difference
206
  - 0: Method of characteristics (MOC)
207
  - 1: Modified method of characteristics (MMOC)
208
  - 2: Hybrid method of characteristics (HMOC)
209
  - 3: Third-order TVD scheme (default)
210
- `percel` (float): Courant number for advection
211
- `mxpart` (int): Maximum number of particles
212

213
### Mt3dDsp
214

215
Dispersion package for modeling mechanical dispersion and molecular diffusion.
216

217
```python { .api }
218
class Mt3dDsp:
219
    """MT3DMS dispersion package"""
220
    def __init__(
221
        self,
222
        model: Mt3dms,
223
        al: ArrayData = 0.01,
224
        trpt: ArrayData = 0.1,
225
        trpv: ArrayData = 0.01,
226
        dmcoef: ArrayData = 0.0,
227
        multiDiff: bool = False,
228
        extension: str = 'dsp',
229
        unitnumber: int = None,
230
        filenames: str = None,
231
        **kwargs
232
    ): ...
233
```
234

235
**Parameters:**
236
- `al` (ArrayData): Longitudinal dispersivity
237
- `trpt` (ArrayData): Ratio of horizontal transverse to longitudinal dispersivity
238
- `trpv` (ArrayData): Ratio of vertical transverse to longitudinal dispersivity
239
- `dmcoef` (ArrayData): Molecular diffusion coefficient
240

241
### Mt3dSsm
242

243
Source/sink mixing package for handling mass loading from boundary conditions.
244

245
```python { .api }
246
class Mt3dSsm:
247
    """MT3DMS source/sink mixing package"""
248
    def __init__(
249
        self,
250
        model: Mt3dms,
251
        crch: ArrayData = None,
252
        cevt: ArrayData = None,
253
        mxss: int = None,
254
        stress_period_data: dict = None,
255
        dtype: np.dtype = None,
256
        extension: str = 'ssm',
257
        unitnumber: int = None,
258
        filenames: str = None,
259
        **kwargs
260
    ): ...
261
```
262

263
**Parameters:**
264
- `crch` (ArrayData): Concentration in recharge
265
- `cevt` (ArrayData): Concentration in evapotranspiration
266
- `stress_period_data` (dict): Source/sink data by stress period
267
  - Format: {period: [[k, i, j, css, itype, (cssm)], ...]}
268
  - `css`: Source/sink concentration
269
  - `itype`: Source/sink type
270

271
### Mt3dRct
272

273
Chemical reactions package for simulating decay, sorption, and other reactions.
274

275
```python { .api }
276
class Mt3dRct:
277
    """MT3DMS chemical reactions package"""
278
    def __init__(
279
        self,
280
        model: Mt3dms,
281
        isothm: int = 0,
282
        ireact: int = 0,
283
        igetsc: int = 0,
284
        rhob: ArrayData = 1.8,
285
        prsity2: ArrayData = None,
286
        srconc: ArrayData = None,
287
        sp1: ArrayData = None,
288
        sp2: ArrayData = None,
289
        rc1: ArrayData = None,
290
        rc2: ArrayData = None,
291
        extension: str = 'rct',
292
        unitnumber: int = None,
293
        filenames: str = None,
294
        **kwargs
295
    ): ...
296
```
297

298
**Parameters:**
299
- `isothm` (int): Sorption isotherm type
300
  - 0: No sorption
301
  - 1: Linear isotherm
302
  - 2: Freundlich isotherm  
303
  - 3: Langmuir isotherm
304
  - 4: First-order kinetic sorption
305
- `ireact` (int): Kinetic rate reaction type
306
- `rhob` (ArrayData): Bulk density
307
- `sp1` (ArrayData): First sorption parameter
308
- `sp2` (ArrayData): Second sorption parameter
309
- `rc1` (ArrayData): First-order reaction rate
310

311
## Solver Package
312

313
### Mt3dGcg
314

315
Generalized conjugate gradient solver for the matrix equations.
316

317
```python { .api }
318
class Mt3dGcg:
319
    """MT3DMS generalized conjugate gradient solver"""
320
    def __init__(
321
        self,
322
        model: Mt3dms,
323
        mxiter: int = 1,
324
        iter1: int = 50,
325
        isolve: int = 3,
326
        ncrs: int = 0,
327
        accl: float = 1.0,
328
        cclose: float = 1e-5,
329
        iprgcg: int = 0,
330
        extension: str = 'gcg',
331
        unitnumber: int = None,
332
        filenames: str = None,
333
        **kwargs
334
    ): ...
335
```
336

337
**Parameters:**
338
- `mxiter` (int): Maximum outer iterations
339
- `iter1` (int): Maximum inner iterations
340
- `isolve` (int): Solver type
341
  - 1: Jacobi
342
  - 2: SSOR
343
  - 3: Modified incomplete Cholesky (MIC)
344
- `cclose` (float): Closure criterion
345
- `accl` (float): Acceleration parameter
346

347
## Observation Package
348

349
### Mt3dTob
350

351
Transport observation package for model calibration and analysis.
352

353
```python { .api }
354
class Mt3dTob:
355
    """MT3DMS transport observations package"""
356
    def __init__(
357
        self,
358
        model: Mt3dms,
359
        outnam: str = 'mt3d_obs.out',
360
        CScale: float = 1.0,
361
        FluxGroups: int = 0,
362
        FScale: float = 1.0,
363
        iOutFlux: int = 0,
364
        obs_data: list = None,
365
        extension: str = 'tob',
366
        no_print: bool = False,
367
        unitnumber: int = None,
368
        filenames: str = None,
369
        **kwargs
370
    ): ...
371
```
372

373
## Advanced Transport Packages
374

375
### Mt3dSft
376

377
Streamflow transport package for stream-aquifer solute exchange.
378

379
```python { .api }
380
class Mt3dSft:
381
    """MT3DMS streamflow transport package"""
382
    def __init__(
383
        self,
384
        model: Mt3dms,
385
        nsfinit: int = 0,
386
        mxsfbc: int = 0,
387
        icbcsf: int = 0,
388
        ioutobs: int = 0,
389
        ietsfr: int = 0,
390
        isfsolv: int = 1,
391
        cclosesf: float = 1e-6,
392
        mxitersf: int = 10,
393
        crntsf: float = 1.0,
394
        iprtxmd: int = 0,
395
        coldsf: ArrayData = 0.0,
396
        dispsf: ArrayData = 0.0,
397
        nobssf: int = 0,
398
        obs_sf: list = None,
399
        sf_stress_period_data: dict = None,
400
        extension: str = 'sft',
401
        unitnumber: int = None,
402
        filenames: str = None,
403
        **kwargs
404
    ): ...
405
```
406

407
### Mt3dLkt
408

409
Lake transport package for lake-aquifer solute exchange.
410

411
```python { .api }
412
class Mt3dLkt:
413
    """MT3DMS lake transport package"""
414
    def __init__(
415
        self,
416
        model: Mt3dms,
417
        nlkinit: int = 0,
418
        mxlkbc: int = 0,
419
        icbclk: int = 0,
420
        ietlak: int = 0,
421
        coldlk: ArrayData = 0.0,
422
        lk_stress_period_data: dict = None,
423
        extension: str = 'lkt',
424
        unitnumber: int = None,
425
        filenames: str = None,
426
        **kwargs
427
    ): ...
428
```
429

430
### Mt3dUzt
431

432
Unsaturated zone transport package.
433

434
```python { .api }
435
class Mt3dUzt:
436
    """MT3DMS unsaturated zone transport package"""
437
    def __init__(
438
        self,
439
        model: Mt3dms,
440
        iuzfbnd: ArrayData = None,
441
        cuzinf: ArrayData = 0.0,
442
        cuzet: ArrayData = 0.0,
443
        cgwet: ArrayData = 0.0,
444
        extension: str = 'uzt',
445
        unitnumber: int = None,
446
        filenames: str = None,
447
        **kwargs
448
    ): ...
449
```
450

451
### Mt3dPhc
452

453
PHC3D reactions package for advanced geochemical modeling.
454

455
```python { .api }
456
class Mt3dPhc:
457
    """MT3DMS PHC3D reactions package"""
458
    def __init__(
459
        self,
460
        model: Mt3dms,
461
        os: int = 2,
462
        temp: ArrayData = 25.0,
463
        asbin: ArrayData = None,
464
        Prime: ArrayData = None,
465
        extension: str = 'phc',
466
        unitnumber: int = None,
467
        filenames: str = None,
468
        **kwargs
469
    ): ...
470
```
471

472
## Usage Examples
473

474
### Basic Single-Species Transport Model
475

476
```python
477
import flopy
478
import numpy as np
479

480
# Create MODFLOW model first (required for transport)
481
mf = flopy.modflow.Modflow(modelname='transport_base')
482

483
# Basic MODFLOW setup
484
nlay, nrow, ncol = 1, 50, 100
485
dis = flopy.modflow.ModflowDis(
486
    mf, nlay=nlay, nrow=nrow, ncol=ncol,
487
    delr=10.0, delc=10.0, top=10.0, botm=0.0
488
)
489
bas = flopy.modflow.ModflowBas(mf, ibound=1, strt=5.0)
490
lpf = flopy.modflow.ModflowLpf(mf, hk=1.0, sy=0.1, ss=1e-5)
491
pcg = flopy.modflow.ModflowPcg(mf)
492
oc = flopy.modflow.ModflowOc(mf)
493

494
# Write and run MODFLOW
495
mf.write_input()
496
success, buff = mf.run_model()
497

498
# Create MT3DMS model
499
mt = flopy.mt3d.Mt3dms(
500
    modelname='transport',
501
    modflowmodel=mf,
502
    ftlfilename='mt3d_link.ftl'
503
)
504

505
# Basic transport package
506
btn = flopy.mt3d.Mt3dBtn(
507
    mt,
508
    ncomp=1,  # Single species
509
    prsity=0.25,
510
    icbund=1,  # Active transport cells
511
    sconc=0.0,  # Starting concentration
512
    timprs=np.array([100.0, 200.0, 500.0, 1000.0]),  # Output times
513
    nprs=4,
514
    chkmas=True,
515
    nprmas=1
516
)
517

518
# Advection package
519
adv = flopy.mt3d.Mt3dAdv(
520
    mt,
521
    mixelm=3,  # TVD scheme
522
    percel=0.75,
523
    nadvfd=1
524
)
525

526
# Dispersion package  
527
dsp = flopy.mt3d.Mt3dDsp(
528
    mt,
529
    al=0.5,     # Longitudinal dispersivity
530
    trpt=0.1,   # Transverse/longitudinal ratio
531
    dmcoef=1e-9  # Molecular diffusion
532
)
533

534
# Source/sink mixing - contaminant source
535
ssm_data = {}
536
# Injection well at center with contamination
537
ssm_data[0] = [
538
    [0, 25, 50, 1000.0, 2]  # Layer, row, col, concentration, type (2=well)
539
]
540
ssm = flopy.mt3d.Mt3dSsm(mt, stress_period_data=ssm_data)
541

542
# Solver
543
gcg = flopy.mt3d.Mt3dGcg(mt, mxiter=1, iter1=50, cclose=1e-5)
544

545
# Write and run MT3DMS
546
mt.write_input()
547
success, buff = mt.run_model()
548
```
549

550
### Multi-Species Reactive Transport Model
551

552
```python
553
import flopy
554
import numpy as np
555

556
# Create base MODFLOW model (abbreviated)
557
mf = flopy.modflow.Modflow(modelname='reactive_transport')
558
# ... MODFLOW setup code ...
559

560
# Multi-species MT3DMS model
561
mt = flopy.mt3d.Mt3dms(
562
    modelname='reactive_mt',
563
    modflowmodel=mf
564
)
565

566
# Multiple species setup
567
ncomp = 3  # Parent compound + 2 decay products
568
species_names = ['PCE', 'TCE', 'DCE']
569

570
btn = flopy.mt3d.Mt3dBtn(
571
    mt,
572
    ncomp=ncomp,
573
    mcomp=ncomp,
574
    species_names=species_names,
575
    prsity=0.30,
576
    icbund=1,
577
    sconc=[[0.0] * ncol * nrow] * ncomp,  # Zero initial concentration
578
    timprs=np.logspace(0, 3, 20),  # Logarithmic time spacing
579
    savucn=True
580
)
581

582
# Advection
583
adv = flopy.mt3d.Mt3dAdv(
584
    mt,
585
    mixelm=0,  # MOC for accuracy with reactions
586
    percel=0.5,
587
    mxpart=100000
588
)
589

590
# Dispersion - different for each species
591
al = [1.0, 0.8, 0.6]  # Decreasing dispersivity with molecular size
592
dsp = flopy.mt3d.Mt3dDsp(
593
    mt,
594
    al=al,
595
    trpt=0.1,
596
    dmcoef=[1e-9, 1.2e-9, 1.5e-9]  # Increasing diffusion
597
)
598

599
# Reactions - sequential decay chain
600
# PCE -> TCE -> DCE
601
rct = flopy.mt3d.Mt3dRct(
602
    mt,
603
    isothm=1,  # Linear sorption
604
    ireact=1,  # First-order decay
605
    rhob=1.8,  # Bulk density
606
    sp1=[[0.5, 0.3, 0.1]],  # Kd values (decreasing sorption)
607
    rc1=[
608
        [0.01, 0.0, 0.0],   # PCE decay rate
609
        [0.0, 0.02, 0.0],   # TCE decay rate  
610
        [0.0, 0.0, 0.05]    # DCE decay rate
611
    ],
612
    rc2=[
613
        [0.0, 0.01, 0.0],   # PCE -> TCE yield
614
        [0.0, 0.0, 0.02],   # TCE -> DCE yield
615
        [0.0, 0.0, 0.0]     # DCE -> end products
616
    ]
617
)
618

619
# Source with declining concentration
620
ssm_data = {}
621
for kper in range(10):  # 10 stress periods
622
    # Exponentially declining source
623
    conc = 1000.0 * np.exp(-0.1 * kper)
624
    ssm_data[kper] = [
625
        [0, 10, 20, conc, 2, 0.0, 0.0]  # Only PCE injected
626
    ]
627

628
ssm = flopy.mt3d.Mt3dSsm(mt, stress_period_data=ssm_data)
629

630
# Solver
631
gcg = flopy.mt3d.Mt3dGcg(
632
    mt,
633
    mxiter=1,
634
    iter1=100,
635
    isolve=3,
636
    cclose=1e-6
637
)
638

639
# Observations at multiple locations
640
obs_data = []
641
# Monitoring wells at different distances
642
for i, dist in enumerate([50, 100, 200]):
643
    row = 25
644
    col = 20 + dist // 10
645
    obs_data.append(['MW_{}'.format(i+1), 0, row, col])
646

647
tob = flopy.mt3d.Mt3dTob(
648
    mt,
649
    obs_data=obs_data,
650
    outnam='reactive_obs.out'
651
)
652

653
# Write and run
654
mt.write_input()
655
success, buff = mt.run_model()
656
```
657

658
### Stream-Aquifer Transport Interaction
659

660
```python
661
import flopy
662
import numpy as np
663

664
# MODFLOW model with SFR package (abbreviated setup)
665
mf = flopy.modflow.Modflow(modelname='stream_transport')
666
# ... discretization, basic, lpf packages ...
667

668
# Add SFR package for stream routing
669
# Stream reaches along model diagonal
670
reach_data = []
671
segment_data = {}
672
for i in range(50):  # 50 reaches
673
    row = i
674
    col = i
675
    reach_data.append([
676
        1, i+1, row, col, 5.0, 1.0, 0.001, 2.0, 1.0, 0.03, 0.5, 0.5
677
    ])
678

679
segment_data[0] = {
680
    1: [1, 1, 0, 100.0, 0.0, 0.0, 0.0, 'stream_seg']
681
}
682

683
sfr = flopy.modflow.ModflowSfr2(
684
    mf,
685
    nstrm=50,
686
    nss=1,
687
    reach_data=np.array(reach_data, dtype=flopy.modflow.ModflowSfr2.get_default_reach_dtype()),
688
    segment_data=segment_data
689
)
690

691
# Run MODFLOW
692
mf.write_input()
693
success, buff = mf.run_model()
694

695
# MT3DMS with stream transport
696
mt = flopy.mt3d.Mt3dms(
697
    modelname='stream_mt',
698
    modflowmodel=mf
699
)
700

701
# Basic transport
702
btn = flopy.mt3d.Mt3dBtn(
703
    mt,
704
    ncomp=1,
705
    prsity=0.25,
706
    sconc=0.0,
707
    timprs=np.array([10.0, 30.0, 100.0, 365.0])
708
)
709

710
# Advection and dispersion
711
adv = flopy.mt3d.Mt3dAdv(mt, mixelm=3)
712
dsp = flopy.mt3d.Mt3dDsp(mt, al=2.0, trpt=0.1)
713

714
# Stream transport package
715
sf_data = {}
716
# Initial stream concentration and dispersion
717
nsfinit = 50  # Number of stream reaches
718
coldsf = np.ones(nsfinit) * 500.0  # Initial stream concentration
719
dispsf = np.ones(nsfinit) * 10.0   # Stream dispersion coefficient
720

721
# Stream boundary conditions by stress period
722
sf_stress_data = {}
723
for kper in range(4):
724
    # Upstream inflow concentration varies seasonally
725
    if kper == 0:  # Spring - high concentration
726
        inflow_conc = 800.0
727
    elif kper == 1:  # Summer - medium concentration  
728
        inflow_conc = 400.0
729
    elif kper == 2:  # Fall - low concentration
730
        inflow_conc = 200.0
731
    else:  # Winter - medium concentration
732
        inflow_conc = 500.0
733
    
734
    sf_stress_data[kper] = [
735
        [1, 1, inflow_conc, 2]  # Segment, reach, concentration, type
736
    ]
737

738
sft = flopy.mt3d.Mt3dSft(
739
    mt,
740
    nsfinit=nsfinit,
741
    mxsfbc=20,
742
    coldsf=coldsf,
743
    dispsf=dispsf,
744
    sf_stress_period_data=sf_stress_data
745
)
746

747
# Source/sink mixing for groundwater-stream interaction
748
ssm_data = {}
749
# Point source in aquifer
750
ssm_data[0] = [
751
    [0, 10, 10, 2000.0, 2]  # High concentration injection
752
]
753
ssm = flopy.mt3d.Mt3dSsm(mt, stress_period_data=ssm_data)
754

755
# Solver
756
gcg = flopy.mt3d.Mt3dGcg(mt)
757

758
# Run transport model
759
mt.write_input()
760
success, buff = mt.run_model()
761
```
762

763
### Advanced Geochemical Modeling with PHC3D
764

765
```python
766
import flopy
767
import numpy as np
768

769
# Complex reactive transport with PHC3D
770
mf = flopy.modflow.Modflow(modelname='geochemical_model')
771
# ... MODFLOW setup ...
772

773
# Multi-component reactive system
774
mt = flopy.mt3d.Mt3dms(modelname='geochemical_mt', modflowmodel=mf)
775

776
# Multiple chemical species
777
ncomp = 6  # Multiple interacting species
778
species_names = ['Ca', 'CO3', 'H', 'CaCO3', 'CaHCO3', 'HCO3']
779

780
btn = flopy.mt3d.Mt3dBtn(
781
    mt,
782
    ncomp=ncomp,
783
    mcomp=ncomp,
784
    species_names=species_names,
785
    prsity=0.25,
786
    # Initial concentrations from equilibrium calculations
787
    sconc=[
788
        [1.0e-3] * (nrow * ncol),  # Ca concentration
789
        [1.0e-4] * (nrow * ncol),  # CO3 concentration  
790
        [1.0e-7] * (nrow * ncol),  # H concentration (pH 7)
791
        [0.0] * (nrow * ncol),     # CaCO3 (mineral)
792
        [1.0e-4] * (nrow * ncol),  # CaHCO3 complex
793
        [1.0e-3] * (nrow * ncol)   # HCO3 concentration
794
    ]
795
)
796

797
# Transport processes
798
adv = flopy.mt3d.Mt3dAdv(mt, mixelm=3)
799
dsp = flopy.mt3d.Mt3dDsp(mt, al=1.0, trpt=0.1)
800

801
# PHC3D for geochemical reactions
802
# Temperature field
803
temp = np.ones((nlay, nrow, ncol)) * 25.0  # 25°C
804

805
# Specify which species are in equilibrium vs kinetic
806
os_flag = 2  # Mixed equilibrium-kinetic system
807

808
phc = flopy.mt3d.Mt3dPhc(
809
    mt,
810
    os=os_flag,
811
    temp=temp
812
)
813

814
# Source with varying chemistry
815
ssm_data = {}
816
for kper in range(5):
817
    # Injection with different pH values
818
    ca_conc = 2.0e-3
819
    co3_conc = 5.0e-4
820
    h_conc = 1.0e-6 * (10 ** (kper - 2))  # pH 4-8 range
821
    
822
    ssm_data[kper] = [
823
        [0, 25, 50, ca_conc, 2, co3_conc, h_conc, 0.0, 0.0, 0.0]
824
    ]
825

826
ssm = flopy.mt3d.Mt3dSsm(mt, stress_period_data=ssm_data)
827

828
# Specialized solver for reactive transport
829
gcg = flopy.mt3d.Mt3dGcg(
830
    mt,
831
    mxiter=1,
832
    iter1=200,  # More iterations for convergence
833
    isolve=3,
834
    cclose=1e-7  # Tighter convergence
835
)
836

837
# Write and run
838
mt.write_input()
839
success, buff = mt.run_model()
840
```
841

842
## Common Types
843

844
```python { .api }
845
# Transport-specific types
846
ConcentrationArray = Union[float, list[float], np.ndarray]
847
ReactionParameter = Union[float, list[float], np.ndarray]
848
SpeciesData = list[ConcentrationArray]
849

850
# Boundary condition types
851
SourceSinkData = list[list[Union[int, float]]]
852
StreamTransportData = dict[int, list[list[Union[int, float]]]]
853
ObservationData = list[list[Union[str, int, float]]]
854

855
# Reaction types
856
IsothermType = Literal[0, 1, 2, 3, 4]  # None, Linear, Freundlich, Langmuir, Kinetic
857
ReactionType = int
858
SolverType = Literal[1, 2, 3]  # Jacobi, SSOR, MIC
859

860
# Advection methods
861
AdvectionMethod = Literal[-1, 0, 1, 2, 3]  # FD, MOC, MMOC, HMOC, TVD
862

863
# Time and mass units
864
TimeUnit = Literal['S', 'M', 'H', 'D', 'Y']
865
LengthUnit = Literal['CM', 'M', 'KM', 'FT', 'MI']  
866
MassUnit = Literal['G', 'KG', 'LB', 'TON']
867

868
# File format types
869
FileFormat = Literal['formatted', 'unformatted']
870
OutputFormat = int
871
```
872

873
This comprehensive documentation covers the complete MT3DMS transport modeling API including all packages, reaction types, and usage patterns. The examples demonstrate basic to advanced transport scenarios including multi-species systems, reactive transport, and coupled surface water-groundwater transport.