0
# File I/O and Post-processing
1

2
This module provides comprehensive file readers for MODFLOW output files, budget analysis, observation processing, and result visualization. These utilities enable extraction, analysis, and manipulation of model results from various MODFLOW-family codes including binary and text output formats.
3

4
## Binary Output File Readers
5

6
### HeadFile
7

8
Reader for MODFLOW head output files (binary format).
9

10
```python { .api }
11
class HeadFile:
12
    """MODFLOW head file reader"""
13
    def __init__(
14
        self,
15
        filename: str,
16
        text: str = 'head',
17
        precision: str = 'single',
18
        verbose: bool = False,
19
        **kwargs
20
    ): ...
21
    
22
    def get_data(
23
        self,
24
        kstpkper: tuple[int, int] = None,
25
        idx: int = None,
26
        totim: float = None
27
    ) -> np.ndarray:
28
        """Get head data for specific time step and stress period"""
29
        ...
30
    
31
    def get_alldata(
32
        self,
33
        mflay: int = None,
34
        nodata: float = -9999.0
35
    ) -> np.ndarray:
36
        """Get all head data for all time steps"""
37
        ...
38
    
39
    def get_ts(
40
        self,
41
        idx: tuple = None,
42
        kstpkper: list = None,
43
        totim: list = None
44
    ) -> np.ndarray:
45
        """Get time series of head at specific locations"""
46
        ...
47
    
48
    def get_kstpkper(self) -> list[tuple[int, int]]:
49
        """Get all available time step/stress period pairs"""
50
        ...
51
    
52
    def get_times(self) -> list[float]:
53
        """Get all simulation times"""
54
        ...
55
    
56
    @property
57
    def nlay(self) -> int:
58
        """Number of layers"""
59
        ...
60
    
61
    @property
62
    def nrow(self) -> int:
63
        """Number of rows"""
64
        ...
65
    
66
    @property
67
    def ncol(self) -> int:
68
        """Number of columns"""
69
        ...
70
    
71
    @property
72
    def shape(self) -> tuple[int, int, int]:
73
        """Array shape (nlay, nrow, ncol)"""
74
        ...
75
    
76
    def close(self) -> None:
77
        """Close the file"""
78
        ...
79
```
80

81
### CellBudgetFile
82

83
Reader for MODFLOW cell budget files containing flow data.
84

85
```python { .api }
86
class CellBudgetFile:
87
    """MODFLOW cell budget file reader"""
88
    def __init__(
89
        self,
90
        filename: str,
91
        precision: str = 'single',
92
        verbose: bool = False,
93
        **kwargs
94
    ): ...
95
    
96
    def get_data(
97
        self,
98
        kstpkper: tuple[int, int] = None,
99
        idx: int = None,
100
        totim: float = None,
101
        text: str = None,
102
        full3D: bool = False
103
    ) -> np.ndarray:
104
        """Get budget data for specific time and budget component"""
105
        ...
106
    
107
    def get_alldata(
108
        self,
109
        text: str = None,
110
        mflay: int = None
111
    ) -> list[np.ndarray]:
112
        """Get all data for specified budget text"""
113
        ...
114
    
115
    def get_ts(
116
        self,
117
        idx: tuple = None,
118
        text: str = None,
119
        kstpkper: list = None,
120
        totim: list = None
121
    ) -> np.ndarray:
122
        """Get time series for specific cells and budget component"""
123
        ...
124
    
125
    def get_budget(
126
        self,
127
        names: list[str] = None,
128
        zones: np.ndarray = None,
129
        net: bool = False,
130
        masked_values: list = None
131
    ) -> np.recarray:
132
        """Get budget summary by zones"""
133
        ...
134
    
135
    def list_unique_records(self) -> list[str]:
136
        """List all unique budget record names"""
137
        ...
138
    
139
    def get_unique_record_names(
140
        self,
141
        decode: bool = True
142
    ) -> list[str]:
143
        """Get unique budget record names"""
144
        ...
145
    
146
    @property
147
    def recordarray(self) -> np.ndarray:
148
        """Record array with budget metadata"""
149
        ...
150
```
151

152
### UcnFile
153

154
Reader for MT3DMS concentration files.
155

156
```python { .api }
157
class UcnFile:
158
    """MT3DMS concentration file reader"""
159
    def __init__(
160
        self,
161
        filename: str,
162
        text: str = 'concentration',
163
        precision: str = 'single',
164
        verbose: bool = False,
165
        **kwargs
166
    ): ...
167
    
168
    def get_data(
169
        self,
170
        kstpkper: tuple[int, int] = None,
171
        idx: int = None,
172
        totim: float = None
173
    ) -> np.ndarray:
174
        """Get concentration data"""
175
        ...
176
    
177
    def get_alldata(
178
        self,
179
        mflay: int = None,
180
        nodata: float = -999.0
181
    ) -> np.ndarray:
182
        """Get all concentration data"""
183
        ...
184
    
185
    def get_ts(
186
        self,
187
        idx: tuple = None,
188
        kstpkper: list = None
189
    ) -> np.ndarray:
190
        """Get concentration time series"""
191
        ...
192
```
193

194
### HeadUFile
195

196
Reader for unformatted head files (MODFLOW-USG and others).
197

198
```python { .api }
199
class HeadUFile:
200
    """Unformatted head file reader"""
201
    def __init__(
202
        self,
203
        filename: str,
204
        text: str = 'headu',
205
        precision: str = 'single',
206
        verbose: bool = False,
207
        **kwargs
208
    ): ...
209
    
210
    def get_data(
211
        self,
212
        kstpkper: tuple[int, int] = None,
213
        idx: int = None,
214
        totim: float = None
215
    ) -> np.ndarray:
216
        """Get unstructured head data"""
217
        ...
218
```
219

220
### FormattedHeadFile
221

222
Reader for formatted (ASCII) head files.
223

224
```python { .api }
225
class FormattedHeadFile:
226
    """Formatted head file reader"""
227
    def __init__(
228
        self,
229
        filename: str,
230
        **kwargs
231
    ): ...
232
    
233
    def get_data(
234
        self,
235
        kstpkper: tuple[int, int] = None
236
    ) -> np.ndarray:
237
        """Get formatted head data"""
238
        ...
239
```
240

241
## List File Readers
242

243
### MfListBudget
244

245
Reader for MODFLOW list file water budgets.
246

247
```python { .api }
248
class MfListBudget:
249
    """MODFLOW list file budget reader"""
250
    def __init__(
251
        self,
252
        filename: str,
253
        timeunit: str = 'days',
254
        start_datetime: str = None,
255
        **kwargs
256
    ): ...
257
    
258
    def get_budget(
259
        self,
260
        names: list[str] = None,
261
        times: list[float] = None
262
    ) -> pd.DataFrame:
263
        """Get water budget data as DataFrame"""
264
        ...
265
    
266
    def get_model_runtime(self) -> dict:
267
        """Get model run time statistics"""
268
        ...
269
    
270
    def get_solver_info(self) -> dict:
271
        """Get solver convergence information"""
272
        ...
273
    
274
    @property
275
    def data(self) -> pd.DataFrame:
276
        """Budget data DataFrame"""
277
        ...
278
    
279
    @property
280
    def flux_names(self) -> list[str]:
281
        """Available flux component names"""
282
        ...
283
```
284

285
### Mf6ListBudget
286

287
Reader for MODFLOW 6 list file budgets.
288

289
```python { .api }
290
class Mf6ListBudget:
291
    """MODFLOW 6 list file budget reader"""
292
    def __init__(
293
        self,
294
        filename: str,
295
        timeunit: str = 'days',
296
        **kwargs
297
    ): ...
298
    
299
    def get_budget(
300
        self,
301
        names: list[str] = None,
302
        times: list[float] = None
303
    ) -> pd.DataFrame:
304
        """Get MODFLOW 6 budget data"""
305
        ...
306
    
307
    def get_model_runtime(self) -> dict:
308
        """Get model execution statistics"""
309
        ...
310
```
311

312
### MfusgListBudget
313

314
Reader for MODFLOW-USG list file budgets.
315

316
```python { .api }
317
class MfusgListBudget:
318
    """MODFLOW-USG list file budget reader"""
319
    def __init__(
320
        self,
321
        filename: str,
322
        timeunit: str = 'days',
323
        **kwargs
324
    ): ...
325
```
326

327
### MtListBudget
328

329
Reader for MT3DMS list file mass budgets.
330

331
```python { .api }
332
class MtListBudget:
333
    """MT3DMS list file budget reader"""
334
    def __init__(
335
        self,
336
        filename: str,
337
        timeunit: str = 'days',
338
        **kwargs
339
    ): ...
340
    
341
    def get_budget(
342
        self,
343
        names: list[str] = None,
344
        times: list[float] = None
345
    ) -> pd.DataFrame:
346
        """Get mass budget data"""
347
        ...
348
```
349

350
### SwrListBudget
351

352
Reader for SWR (Surface Water Routing) list file budgets.
353

354
```python { .api }
355
class SwrListBudget:
356
    """SWR list file budget reader"""
357
    def __init__(
358
        self,
359
        filename: str,
360
        timeunit: str = 'days',
361
        **kwargs
362
    ): ...
363
```
364

365
### SwtListBudget
366

367
Reader for SEAWAT list file budgets.
368

369
```python { .api }
370
class SwtListBudget:
371
    """SEAWAT list file budget reader"""
372
    def __init__(
373
        self,
374
        filename: str,
375
        timeunit: str = 'days',
376
        **kwargs
377
    ): ...
378
```
379

380
## Observation File Readers
381

382
### Mf6Obs
383

384
Reader for MODFLOW 6 observation files.
385

386
```python { .api }
387
class Mf6Obs:
388
    """MODFLOW 6 observation file reader"""
389
    def __init__(
390
        self,
391
        filename: str,
392
        itype: str = 'head',
393
        **kwargs
394
    ): ...
395
    
396
    def get_data(
397
        self,
398
        obsname: str = None,
399
        totim: float = None,
400
        start_datetime: str = None
401
    ) -> pd.DataFrame:
402
        """Get observation data as DataFrame"""
403
        ...
404
    
405
    def get_ts(
406
        self,
407
        obsname: list[str] = None
408
    ) -> pd.DataFrame:
409
        """Get time series for specified observations"""
410
        ...
411
    
412
    @property
413
    def obs_names(self) -> list[str]:
414
        """Available observation names"""
415
        ...
416
    
417
    @property
418
    def data(self) -> pd.DataFrame:
419
        """All observation data"""
420
        ...
421
```
422

423
### HydmodObs
424

425
Reader for HYDMOD observation files.
426

427
```python { .api }
428
class HydmodObs:
429
    """HYDMOD observation file reader"""
430
    def __init__(
431
        self,
432
        filename: str,
433
        **kwargs
434
    ): ...
435
    
436
    def get_data(self) -> pd.DataFrame:
437
        """Get HYDMOD observation data"""
438
        ...
439
```
440

441
### SwrObs
442

443
Reader for SWR observation files.
444

445
```python { .api }
446
class SwrObs:
447
    """SWR observation file reader"""
448
    def __init__(
449
        self,
450
        filename: str,
451
        **kwargs
452
    ): ...
453
```
454

455
## Specialized Output File Readers
456

457
### SfrFile
458

459
Reader for SFR (Streamflow Routing) output files.
460

461
```python { .api }
462
class SfrFile:
463
    """SFR output file reader"""
464
    def __init__(
465
        self,
466
        filename: str,
467
        **kwargs
468
    ): ...
469
    
470
    def get_data(
471
        self,
472
        kstpkper: tuple[int, int] = None,
473
        idx: int = None
474
    ) -> np.ndarray:
475
        """Get SFR data for specific time"""
476
        ...
477
    
478
    def get_ts(
479
        self,
480
        iseg: int = None,
481
        ireach: int = None
482
    ) -> pd.DataFrame:
483
        """Get time series for specific reach/segment"""
484
        ...
485
```
486

487
### SwrBudget
488

489
Reader for SWR budget output files.
490

491
```python { .api }
492
class SwrBudget:
493
    """SWR budget output reader"""
494
    def __init__(
495
        self,
496
        filename: str,
497
        **kwargs
498
    ): ...
499
```
500

501
### SwrExchange
502

503
Reader for SWR exchange output files.
504

505
```python { .api }
506
class SwrExchange:
507
    """SWR exchange output reader"""
508
    def __init__(
509
        self,
510
        filename: str,
511
        **kwargs
512
    ): ...
513
```
514

515
### SwrFlow
516

517
Reader for SWR flow output files.
518

519
```python { .api }
520
class SwrFlow:
521
    """SWR flow output reader"""
522
    def __init__(
523
        self,
524
        filename: str,
525
        **kwargs
526
    ): ...
527
```
528

529
### SwrStage
530

531
Reader for SWR stage output files.
532

533
```python { .api }
534
class SwrStage:
535
    """SWR stage output reader"""
536
    def __init__(
537
        self,
538
        filename: str,
539
        **kwargs
540
    ): ...
541
```
542

543
### SwrStructure
544

545
Reader for SWR structure output files.
546

547
```python { .api }
548
class SwrStructure:
549
    """SWR structure output reader"""
550
    def __init__(
551
        self,
552
        filename: str,
553
        **kwargs
554
    ): ...
555
```
556

557
## MODPATH Output File Readers
558

559
### EndpointFile
560

561
Reader for MODPATH endpoint files.
562

563
```python { .api }
564
class EndpointFile:
565
    """MODPATH endpoint file reader"""
566
    def __init__(
567
        self,
568
        filename: str,
569
        **kwargs
570
    ): ...
571
    
572
    def get_data(
573
        self,
574
        partid: int = None,
575
        groupname: str = None
576
    ) -> np.ndarray:
577
        """Get endpoint data"""
578
        ...
579
    
580
    def get_alldata(self) -> np.ndarray:
581
        """Get all endpoint data"""
582
        ...
583
    
584
    def get_destination_data(
585
        self,
586
        dest_cells: list = None
587
    ) -> np.ndarray:
588
        """Get endpoints reaching specific destinations"""
589
        ...
590
```
591

592
### PathlineFile
593

594
Reader for MODPATH pathline files.
595

596
```python { .api }
597
class PathlineFile:
598
    """MODPATH pathline file reader"""
599
    def __init__(
600
        self,
601
        filename: str,
602
        **kwargs
603
    ): ...
604
    
605
    def get_data(
606
        self,
607
        partid: int = None,
608
        groupname: str = None
609
    ) -> list[np.ndarray]:
610
        """Get pathline data for specific particles"""
611
        ...
612
    
613
    def get_alldata(self) -> list[np.ndarray]:
614
        """Get all pathline data"""
615
        ...
616
    
617
    def intersect_polygon(
618
        self,
619
        polygon: list[tuple[float, float]]
620
    ) -> list[np.ndarray]:
621
        """Get pathlines intersecting polygon"""
622
        ...
623
```
624

625
### TimeseriesFile
626

627
Reader for MODPATH timeseries files.
628

629
```python { .api }
630
class TimeseriesFile:
631
    """MODPATH timeseries file reader"""
632
    def __init__(
633
        self,
634
        filename: str,
635
        **kwargs
636
    ): ...
637
    
638
    def get_data(
639
        self,
640
        partid: int = None,
641
        totim: float = None
642
    ) -> np.ndarray:
643
        """Get timeseries data"""
644
        ...
645
```
646

647
## Zone Budget Analysis
648

649
### ZoneBudget
650

651
Zone budget analysis for MODFLOW models.
652

653
```python { .api }
654
class ZoneBudget:
655
    """Zone budget analysis for MODFLOW"""
656
    def __init__(
657
        self,
658
        cbc_file: str,
659
        z: np.ndarray,
660
        kstpkper: tuple[int, int] = None,
661
        totim: float = None,
662
        aliases: dict = None,
663
        **kwargs
664
    ): ...
665
    
666
    def get_budget(
667
        self,
668
        names: list[str] = None,
669
        zones: list[int] = None,
670
        net: bool = False,
671
        pivot: bool = False
672
    ) -> pd.DataFrame:
673
        """Calculate zone budget"""
674
        ...
675
    
676
    def get_model_budget(
677
        self,
678
        f_out: str = None,
679
        **kwargs
680
    ) -> pd.DataFrame:
681
        """Get overall model budget"""
682
        ...
683
    
684
    def export(
685
        self,
686
        f_out: str,
687
        df: pd.DataFrame = None,
688
        **kwargs
689
    ) -> None:
690
        """Export budget results"""
691
        ...
692
    
693
    @property
694
    def budget(self) -> pd.DataFrame:
695
        """Zone budget DataFrame"""
696
        ...
697
```
698

699
### ZoneBudget6
700

701
Zone budget analysis for MODFLOW 6.
702

703
```python { .api }
704
class ZoneBudget6:
705
    """Zone budget analysis for MODFLOW 6"""
706
    def __init__(
707
        self,
708
        model: object = None,
709
        flowja: np.ndarray = None,
710
        budget_file: str = None,
711
        zone_file: str = None,
712
        kstpkper: tuple[int, int] = None,
713
        **kwargs
714
    ): ...
715
    
716
    def get_budget(
717
        self,
718
        zones: list[int] = None,
719
        net: bool = False
720
    ) -> pd.DataFrame:
721
        """Calculate MODFLOW 6 zone budget"""
722
        ...
723
```
724

725
### ZoneFile6
726

727
Reader for MODFLOW 6 zone files.
728

729
```python { .api }
730
class ZoneFile6:
731
    """Zone file handling for MODFLOW 6"""
732
    def __init__(
733
        self,
734
        zonefile: str,
735
        **kwargs
736
    ): ...
737
    
738
    def get_zones(self) -> np.ndarray:
739
        """Get zone array"""
740
        ...
741
```
742

743
### ZBNetOutput
744

745
Zone budget network output processor.
746

747
```python { .api }
748
class ZBNetOutput:
749
    """Zone budget network output"""
750
    def __init__(
751
        self,
752
        zbobj: ZoneBudget,
753
        **kwargs
754
    ): ...
755
    
756
    def write_network(
757
        self,
758
        filename: str,
759
        **kwargs
760
    ) -> None:
761
        """Write network file"""
762
        ...
763
```
764

765
## Binary Data Handling
766

767
### BinaryHeader
768

769
Binary file header information.
770

771
```python { .api }
772
class BinaryHeader:
773
    """Binary file header handling"""
774
    def __init__(
775
        self,
776
        bintype: str = None,
777
        precision: str = 'single',
778
        **kwargs
779
    ): ...
780
    
781
    @classmethod
782
    def set_dtype(
783
        cls,
784
        bintype: str = None,
785
        precision: str = 'single'
786
    ) -> np.dtype:
787
        """Set numpy dtype for binary data"""
788
        ...
789
```
790

791
### FlopyBinaryData
792

793
Container for binary data with metadata.
794

795
```python { .api }
796
class FlopyBinaryData:
797
    """Binary data container"""
798
    def __init__(
799
        self,
800
        array: np.ndarray,
801
        bintype: str = None,
802
        **kwargs
803
    ): ...
804
    
805
    def plot(self, **kwargs) -> object:
806
        """Plot binary data"""
807
        ...
808
    
809
    def to_shapefile(
810
        self,
811
        filename: str,
812
        **kwargs
813
    ) -> None:
814
        """Export to shapefile"""
815
        ...
816
```
817

818
## Usage Examples
819

820
### Basic Head and Budget Analysis
821

822
```python
823
import flopy.utils as fpu
824
import numpy as np
825
import matplotlib.pyplot as plt
826

827
# Read head file
828
hds = fpu.HeadFile('model.hds')
829

830
# Get head data for specific times
831
head_final = hds.get_data(kstpkper=(0, 0))  # First time step, first stress period
832
head_all = hds.get_alldata()  # All time steps
833

834
print(f"Head array shape: {head_final.shape}")
835
print(f"All heads shape: {head_all.shape}")
836

837
# Time series at specific locations
838
# Wells at (layer, row, col) locations
839
well_locs = [(0, 25, 25), (1, 40, 60), (0, 15, 35)]
840
head_ts = hds.get_ts(idx=well_locs)
841

842
# Plot head time series
843
times = hds.get_times()
844
plt.figure(figsize=(10, 6))
845
for i, (lay, row, col) in enumerate(well_locs):
846
    plt.plot(times, head_ts[:, i], label=f'Well {i+1} (L{lay},R{row},C{col})')
847
plt.xlabel('Time (days)')
848
plt.ylabel('Head (m)')
849
plt.legend()
850
plt.title('Head Time Series at Monitoring Wells')
851
plt.show()
852

853
# Read budget file
854
cbb = fpu.CellBudgetFile('model.cbb')
855

856
# List available budget records
857
records = cbb.list_unique_records()
858
print("Available budget records:", records)
859

860
# Get specific budget components
861
storage = cbb.get_data(text='STORAGE')
862
wells = cbb.get_data(text='WELLS')
863
recharge = cbb.get_data(text='RECHARGE')
864

865
print(f"Storage shape: {storage.shape if storage is not None else 'Not found'}")
866
print(f"Wells shape: {wells.shape if wells is not None else 'Not found'}")
867

868
# Budget time series at specific cells
869
if wells is not None:
870
    well_flow_ts = cbb.get_ts(idx=well_locs, text='WELLS')
871
    
872
    plt.figure(figsize=(10, 6))
873
    for i, (lay, row, col) in enumerate(well_locs):
874
        plt.plot(times, well_flow_ts[:, i], label=f'Well {i+1}')
875
    plt.xlabel('Time (days)')
876
    plt.ylabel('Well Flow Rate (m³/d)')
877
    plt.legend()
878
    plt.title('Well Flow Time Series')
879
    plt.show()
880

881
# Close files
882
hds.close()
883
cbb.close()
884
```
885

886
### Zone Budget Analysis
887

888
```python
889
import flopy.utils as fpu
890
import numpy as np
891
import pandas as pd
892

893
# Define zone array
894
nlay, nrow, ncol = 3, 50, 100
895
zones = np.ones((nlay, nrow, ncol), dtype=int)
896

897
# Create different zones based on model features
898
# Zone 1: Natural recharge area (rows 0-20)
899
zones[:, 0:20, :] = 1
900

901
# Zone 2: Agricultural area (rows 20-35)  
902
zones[:, 20:35, :] = 2
903

904
# Zone 3: Urban area (rows 35-50)
905
zones[:, 35:50, :] = 3
906

907
# Different zones by layer (geological units)
908
zones[0, :, :] *= 10  # Multiply by 10 for top layer
909
zones[1, :, :] *= 20  # Multiply by 20 for middle layer  
910
zones[2, :, :] *= 30  # Multiply by 30 for bottom layer
911

912
# Zone budget analysis
913
zb = fpu.ZoneBudget('model.cbb', zones, kstpkper=(0, 1))
914

915
# Get budget for all zones
916
budget_df = zb.get_budget(
917
    names=['WELLS', 'RECHARGE', 'ET', 'RIVERS', 'STORAGE'],
918
    zones=[10, 20, 30, 40, 50, 60],  # Top layer zones
919
    net=False  # Show inflows and outflows separately
920
)
921

922
print("Zone Budget Summary:")
923
print(budget_df)
924

925
# Calculate net budget (IN - OUT)
926
net_budget = zb.get_budget(
927
    names=['WELLS', 'RECHARGE', 'ET', 'RIVERS'],
928
    zones=[10, 20, 30],  # Selected zones
929
    net=True
930
)
931

932
print("\nNet Zone Budget:")
933
print(net_budget)
934

935
# Export results
936
budget_df.to_csv('zone_budget_detailed.csv')
937
net_budget.to_csv('zone_budget_net.csv')
938

939
# Model-wide budget
940
model_budget = zb.get_model_budget()
941
print("\nModel-wide Budget:")
942
print(model_budget)
943
```
944

945
### MODFLOW 6 Advanced File Reading
946

947
```python
948
import flopy
949
import flopy.utils as fpu
950
import pandas as pd
951

952
# Read MODFLOW 6 simulation results
953
sim = flopy.mf6.MFSimulation.load(sim_ws='mf6_model')
954
gwf = sim.get_model('gwf')
955

956
# Head file
957
hds = fpu.HeadFile('gwf.hds')
958
head_data = hds.get_alldata()
959

960
# Budget file with detailed analysis
961
cbc = fpu.CellBudgetFile('gwf.cbc')
962

963
# Get specific discharge vectors for flow visualization
964
spdis = cbc.get_data(text='DATA-SPDIS')  # Specific discharge
965
if spdis is not None:
966
    qx, qy, qz = spdis[0], spdis[1], spdis[2]
967
    print(f"Specific discharge shapes: qx={qx.shape}, qy={qy.shape}, qz={qz.shape}")
968

969
# MODFLOW 6 observations
970
try:
971
    obs_head = fpu.Mf6Obs('gwf.obs.head.csv', itype='head')
972
    head_obs_df = obs_head.get_data()
973
    
974
    print("Available head observations:", obs_head.obs_names)
975
    print("\nHead observation data:")
976
    print(head_obs_df.head())
977
    
978
    # Plot observations vs time
979
    import matplotlib.pyplot as plt
980
    
981
    plt.figure(figsize=(12, 6))
982
    for obs_name in obs_head.obs_names[:5]:  # First 5 observations
983
        obs_ts = obs_head.get_ts([obs_name])
984
        plt.plot(obs_ts['totim'], obs_ts[obs_name], label=obs_name)
985
    
986
    plt.xlabel('Time (days)')
987
    plt.ylabel('Head (m)')
988
    plt.legend()
989
    plt.title('Head Observations Time Series')
990
    plt.show()
991
    
992
except FileNotFoundError:
993
    print("Observation files not found")
994

995
# MODFLOW 6 list budget
996
try:
997
    mf6_lst = fpu.Mf6ListBudget('gwf.lst')
998
    budget_df = mf6_lst.get_budget()
999
    
1000
    print("\nMODFLOW 6 Budget Summary:")
1001
    print(budget_df.head())
1002
    
1003
    # Runtime information
1004
    runtime_info = mf6_lst.get_model_runtime()
1005
    print(f"\nModel runtime: {runtime_info}")
1006
    
1007
except FileNotFoundError:
1008
    print("List file not found")
1009

1010
# MODFLOW 6 zone budget
1011
try:
1012
    # Create zone array
1013
    grid = gwf.modelgrid
1014
    zones = np.ones(grid.shape, dtype=int)
1015
    
1016
    # Different zones based on K values or other criteria
1017
    k_array = gwf.npf.k.array
1018
    zones[k_array > 10.0] = 2  # High K zone
1019
    zones[k_array < 1.0] = 3   # Low K zone
1020
    
1021
    zb6 = fpu.ZoneBudget6(
1022
        model=gwf,
1023
        budget_file='gwf.cbc',
1024
        zone_file=None  # Will use zones array directly
1025
    )
1026
    
1027
    mf6_budget = zb6.get_budget(zones=[1, 2, 3])
1028
    print("\nMODFLOW 6 Zone Budget:")
1029
    print(mf6_budget)
1030
    
1031
except Exception as e:
1032
    print(f"Zone budget error: {e}")
1033
```
1034

1035
### Transport Results Analysis
1036

1037
```python
1038
import flopy.utils as fpu
1039
import numpy as np
1040
import matplotlib.pyplot as plt
1041

1042
# Read MT3DMS concentration results
1043
ucn = fpu.UcnFile('MT3D001.UCN')
1044

1045
# Get concentration data
1046
conc_final = ucn.get_data()  # Final time step
1047
conc_all = ucn.get_alldata()  # All time steps
1048

1049
print(f"Final concentration shape: {conc_final.shape}")
1050
print(f"All concentrations shape: {conc_all.shape}")
1051

1052
# Concentration time series at monitoring points
1053
monitor_locs = [(0, 10, 25), (0, 20, 50), (0, 30, 75)]
1054
conc_ts = ucn.get_ts(idx=monitor_locs)
1055

1056
# Plot concentration breakthrough curves
1057
times = ucn.get_times()
1058
plt.figure(figsize=(10, 6))
1059
for i, (lay, row, col) in enumerate(monitor_locs):
1060
    plt.semilogy(times, conc_ts[:, i], 'o-', label=f'MW-{i+1} (R{row},C{col})')
1061
plt.xlabel('Time (days)')
1062
plt.ylabel('Concentration (mg/L)')
1063
plt.legend()
1064
plt.title('Concentration Breakthrough Curves')
1065
plt.grid(True)
1066
plt.show()
1067

1068
# Calculate mass distribution over time
1069
total_mass = []
1070
for i in range(len(times)):
1071
    conc_3d = ucn.get_data(idx=i)
1072
    # Assume porosity, cell volumes for mass calculation
1073
    porosity = 0.25
1074
    cell_volume = 100.0 * 100.0 * 10.0  # delr * delc * thickness
1075
    mass = np.sum(conc_3d * porosity * cell_volume)
1076
    total_mass.append(mass)
1077

1078
plt.figure(figsize=(10, 6))
1079
plt.plot(times, total_mass, 'b-', linewidth=2)
1080
plt.xlabel('Time (days)')
1081
plt.ylabel('Total Mass in Domain (kg)')
1082
plt.title('Mass Conservation Over Time')
1083
plt.grid(True)
1084
plt.show()
1085

1086
# MT3DMS list budget for mass balance
1087
try:
1088
    mt_lst = fpu.MtListBudget('MT3D.LST')
1089
    mass_budget = mt_lst.get_budget()
1090
    
1091
    print("\nMT3DMS Mass Budget:")
1092
    print(mass_budget)
1093
except FileNotFoundError:
1094
    print("MT3DMS list file not found")
1095

1096
ucn.close()
1097
```
1098

1099
### MODPATH Results Analysis
1100

1101
```python
1102
import flopy.utils as fpu
1103
import numpy as np
1104
import matplotlib.pyplot as plt
1105

1106
# Read MODPATH endpoint file
1107
ep = fpu.EndpointFile('model.mpend')
1108
endpoint_data = ep.get_alldata()
1109

1110
print(f"Number of particles: {len(endpoint_data)}")
1111
print(f"Endpoint data columns: {endpoint_data.dtype.names}")
1112

1113
# Analyze particle destinations
1114
# Status codes: 1=active, 2=terminated at boundary, 3=terminated at sink, etc.
1115
status_counts = {}
1116
for status in np.unique(endpoint_data['status']):
1117
    count = np.sum(endpoint_data['status'] == status)
1118
    status_counts[status] = count
1119
    print(f"Status {status}: {count} particles")
1120

1121
# Plot particle starting vs ending locations
1122
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))
1123

1124
# Starting locations
1125
ax1.scatter(endpoint_data['x0'], endpoint_data['y0'], 
1126
           c=endpoint_data['particlegroup'], alpha=0.6)
1127
ax1.set_xlabel('X Coordinate')
1128
ax1.set_ylabel('Y Coordinate')
1129
ax1.set_title('Particle Starting Locations')
1130
ax1.grid(True)
1131

1132
# Ending locations colored by travel time
1133
travel_time = endpoint_data['time'] - endpoint_data['time0']
1134
scatter = ax2.scatter(endpoint_data['x'], endpoint_data['y'], 
1135
                     c=travel_time, cmap='viridis', alpha=0.6)
1136
ax2.set_xlabel('X Coordinate')
1137
ax2.set_ylabel('Y Coordinate')
1138
ax2.set_title('Particle Ending Locations (colored by travel time)')
1139
ax2.grid(True)
1140

1141
cbar = plt.colorbar(scatter, ax=ax2)
1142
cbar.set_label('Travel Time (days)')
1143

1144
plt.tight_layout()
1145
plt.show()
1146

1147
# Read pathline file for detailed tracking
1148
try:
1149
    pth = fpu.PathlineFile('model.mppth')
1150
    pathline_data = pth.get_alldata()
1151
    
1152
    print(f"Number of pathlines: {len(pathline_data)}")
1153
    
1154
    # Plot selected pathlines
1155
    plt.figure(figsize=(12, 8))
1156
    
1157
    # Plot first 20 pathlines
1158
    for i, pathline in enumerate(pathline_data[:20]):
1159
        plt.plot(pathline['x'], pathline['y'], 'b-', alpha=0.5, linewidth=1)
1160
        # Mark starting point
1161
        plt.plot(pathline['x'][0], pathline['y'][0], 'go', markersize=4)
1162
        # Mark ending point
1163
        plt.plot(pathline['x'][-1], pathline['y'][-1], 'ro', markersize=4)
1164
    
1165
    plt.xlabel('X Coordinate')
1166
    plt.ylabel('Y Coordinate')
1167
    plt.title('Particle Pathlines (first 20)')
1168
    plt.legend(['Pathlines', 'Start', 'End'])
1169
    plt.grid(True)
1170
    plt.show()
1171
    
1172
    # Calculate pathline statistics
1173
    path_lengths = []
1174
    travel_times = []
1175
    
1176
    for pathline in pathline_data:
1177
        # Calculate path length
1178
        dx = np.diff(pathline['x'])
1179
        dy = np.diff(pathline['y'])
1180
        dz = np.diff(pathline['z'])
1181
        segments = np.sqrt(dx**2 + dy**2 + dz**2)
1182
        path_length = np.sum(segments)
1183
        path_lengths.append(path_length)
1184
        
1185
        # Travel time
1186
        travel_time = pathline['time'][-1] - pathline['time'][0]
1187
        travel_times.append(travel_time)
1188
    
1189
    # Plot statistics
1190
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))
1191
    
1192
    ax1.hist(path_lengths, bins=20, alpha=0.7, edgecolor='black')
1193
    ax1.set_xlabel('Path Length (m)')
1194
    ax1.set_ylabel('Frequency')
1195
    ax1.set_title('Distribution of Path Lengths')
1196
    ax1.grid(True)
1197
    
1198
    ax2.hist(travel_times, bins=20, alpha=0.7, edgecolor='black')
1199
    ax2.set_xlabel('Travel Time (days)')
1200
    ax2.set_ylabel('Frequency')
1201
    ax2.set_title('Distribution of Travel Times')
1202
    ax2.grid(True)
1203
    
1204
    plt.tight_layout()
1205
    plt.show()
1206
    
1207
    print(f"Average path length: {np.mean(path_lengths):.1f} m")
1208
    print(f"Average travel time: {np.mean(travel_times):.1f} days")
1209
    
1210
except FileNotFoundError:
1211
    print("Pathline file not found")
1212
```
1213

1214
## Common Types
1215

1216
```python { .api }
1217
# File I/O types
1218
FilePath = Union[str, os.PathLike]
1219
Precision = Literal['single', 'double']
1220
TextRecord = str
1221

1222
# Time and indexing types
1223
TimeStepStressPeriod = tuple[int, int]  # (kstp, kper)
1224
CellIndex = tuple[int, int, int]  # (layer, row, col) or node number
1225
TimeValue = float
1226

1227
# Data array types
1228
HeadArray = np.ndarray  # Shape: (nlay, nrow, ncol) or (nodes,)
1229
BudgetArray = np.ndarray  # Shape varies by budget type
1230
ConcentrationArray = np.ndarray  # Shape: (nlay, nrow, ncol) or (nodes,)
1231
ZoneArray = np.ndarray  # Integer array for zone definitions
1232

1233
# Budget and observation types
1234
BudgetRecord = str  # Budget text identifier
1235
ObservationName = str
1236
FluxComponentNames = list[str]
1237

1238
# Output data types
1239
TimeSeriesData = np.ndarray  # Shape: (ntimes, nlocations)
1240
BudgetData = pd.DataFrame
1241
ObservationData = pd.DataFrame
1242
PathlineData = list[np.ndarray]
1243
EndpointData = np.ndarray
1244

1245
# File format specifications  
1246
BinaryFileType = Literal['head', 'budget', 'concentration', 'drawdown']
1247
ListFileType = Literal['modflow', 'modflow6', 'mt3dms', 'seawat', 'swr']
1248
```
1249

1250
This comprehensive documentation covers the complete file I/O and post-processing API for FloPy including all major file readers, budget analysis tools, and usage patterns. The examples demonstrate basic to advanced result analysis scenarios including time series extraction, zone budget analysis, and particle tracking visualization.