MODFLOW stream package, lake package, stream-lake interaction, SFR vs. STR

So you want to model the interaction between lake, stream, aquifer in MODFLOW? You need the lake package and a stream package. There are two stream packages available - STR and SFR. According to my experiment  you need SFR for lake-stream interaction. STR is depreciated since MODFLOW 2000 (you can't even find STR input file documentation on USGS website. It is hidden in your MODFLOW 2000 installation folder WRDAPP\mf2k.1_19\doc\updates). Problems arise when you use software like Visual Modflow where only STR package is supported. For simple geometry, you can convert STR to SFR input file without too much trouble. This post closely examines the SFR and STR input files, compare the difference and discuss how to convert STR input file to SFR format.

The example problem is from simulation 2 of official lake package documentation. There are two lakes and three streams connection them. Flow direction is from north to south.

Some terminologies to know. GMS's wiki page explains them pretty well.
"A stream network is composed of reaches and segments. A reach is the portion of a stream that lies inside a single cell. A single cell may contain multiple reaches. A segment is a group of reaches that forms one section of the stream. The reaches within a segment are always numbered from upstream to downstream. The segments should also be numbered from upstream to downstream."

In the figure below. There are three stream segments. And Segment 1 has 8 reaches (stream covers 8 cells). Reaches are numbered sequentially from upstream to downstream.


All right, let's get to the files. Only a brief description is provided here and emphasis is on converting STR file to SFR file. Detailed descriptions of each parameter can be found in official SFR package doc page 40.
SFR package has five ways of calculating stream depth, flagged by ICALC in the input file. Manning equation with rectangle stream channel (ICALC=1) is used here. This option is also supported in STR package. Also, things get complicated when tributaries and diversions are involved. The example has no tributary or diversion involved. Stay tuned for part 2 of this post about tributary and diversion.

SFR file:
22,3,0,0,128390.4,.0001,-1,0
Explanation
22-Total # of reaches, 3-Total # of segments, these two need to be change  for each simulation
0, 0 -keep 0 if steady-steate and parameters same within a segment, 128390.4-this is the constant conversion factor if Manning equation is used (it's equal to 1.485*86400 if unit is specified in days), 0.0001-tolerance, not changed in most cases, -1, 0- output control.
1,2,5,1,1,1000.                         ;segment and reach identification and reach length
1,3,5,1,2,1000.     ; layer, row, column, segment, reach, reach length
1,4,5,1,3,1000.
1,4,6,1,4,1000.
1,4,7,1,5,500.
1,4,8,1,6,750.
1,5,8,1,7,1000.
1,6,8,1,8,1000.
1,12,9,2,1,1000.
1,13,9,2,2,1000.
1,14,10,2,3,559.
1,14,10,2,4,559.
1,15,10,2,5,1000.
1,16,10,2,6,1000.
1,22,10,3,1,1000.
1,23,10,3,2,750.
1,23,11,3,3,500.
1,23,12,3,4,1000.
1,23,13,3,5,1000.
1,23,14,3,6,1000.
1,23,15,3,7,1000.
1,23,16,3,8,1000.
3,0,0            ;3-number of segments, 0,0-read and print flags, see page 50 of doc for details
1 1 -1 0   691200.0 0.0   0.0   0.0   0.05    ;flow and environmental data for segment 1
Explanation.
1-segment number,
1-ICALC,
-1, OUTSEG, where stream is connected downstream. Negative means connection to lake, 0 means connect to nothing. In this case, it is connected to Lake 1,
0, INSEG, where stream is connected upstream. Same format as OUTSEG
691200.0, FLOW, flow into 1st reach.
0.0   0.0   0.0 - runoff, evaporation, precipitation
0.05- Manning roughness coefficient
0.50,0.50,124.5,5.0                           ;upstream stream channel information for segment 1
0.50-
0.50-stream bed thickness
124.5-elevation of top of stream bed
5.0-channel width
0.50,0.50,116.5,5.0                           ;downstream stream channel information for segment 1
2 1 -2 -1       0.0 0.0   0.0   0.0   0.05    ;flow and environmental data for segment 2
0.50,0.50,114.85,5.0                          ;upstream stream channel information for segment 2
0.50,0.50,110.65,5.0                          ;downstream stream channel information for segment
3 1  0 -2       0.0 0.0   0.0   0.0   0.05    ;flow and environmental data for segment 3
0.50,0.50,109.4286,5.0                        ;upstream stream channel information for segment 3
0.50,0.50,102.5714,5.0                        ;downstream stream channel information for segment 3  



STR file:
PARAMETER 0 0
        22         3         0         0         1 1.28383e5       154       155
Explanations
22-total number of reaches
3-total number of segments
0-total number of tributaries
0-total number of divesions
1-ICALC see explanation from above
1.28383e5-constant  = 1.486*86400 if ICALC = 1
154,155-output units
        22         0         0 ; read and printing flags, 22 is total number of reaches
    1    2    5    1    1          6.9e5        0. 4.72046e3   1.235e2    1.24e2
Explanation: lay,row,col,segment,reach,flow, stage(0=to be calculated),  bottom of streambed, top of streambed
    1    3    5    1    2             0.        0. 5.00033e3  1.2238e2  1.2288e2
    1    4    5    1    3             0.        0. 4.57631e3  1.2126e2  1.2176e2
    1    4    6    1    4             0.        0. 5.00038e3  1.2014e2  1.2064e2
    1    4    7    1    5             0.        0. 2.50019e3   1.193e2   1.198e2
    1    4    8    1    6             0.        0. 4.32848e3  1.1874e2  1.1924e2
    1    5    8    1    7             0.        0.      5.e3  1.1762e2  1.1812e2
    1    6    8    1    8             0.        0. 5.87008e3   1.165e2    1.17e2
    1   12    9    2    1           0        0. 4.74142e3    1.14e2   1.145e2
    1   13    9    2    2             0.        0. 5.00041e3 1.13222e2 1.13722e2

    1   14    9    2    3             0.        0. 3.09242e3 1.12444e2 1.12944e2
    1   14   10    2    4             0.        0. 2.47465e3 1.12056e2 1.12556e2
    1   15   10    2    5             0.        0.   5.012e3 1.11278e2 1.11778e2
    1   16   10    2    6             0.        0. 4.52491e3   1.105e2    1.11e2
    1   22   10    3    1           0.        0. 4.34095e3   1.085e2    1.09e2
    1   23   10    3    2             0.        0. 4.92948e3 1.07627e2 1.08127e2
    1   23   11    3    3             0.        0. 2.50026e3 1.07191e2 1.07691e2
    1   23   12    3    4             0.        0. 5.00051e3 1.06536e2 1.07036e2
    1   23   13    3    5             0.        0. 5.00051e3 1.05664e2 1.06164e2
    1   23   14    3    6             0.        0. 5.00051e3 1.04791e2 1.05291e2
    1   23   15    3    7             0.        0. 5.00051e3 1.03918e2 1.04418e2
    1   23   16    3    8             0.        0. 5.00051e3 1.03045e2 1.03545e2
    1   23   17    3    9             0.        0. 8.65081e1   1.025e2    1.03e2
        5.   1.12e-3     5.e-2 ; width, slope(STR only), manning roughness
        5.   1.12e-3     5.e-2
        5.   1.12e-3     5.e-2
        5.   1.12e-3     5.e-2
        5.   1.12e-3     5.e-2
        5.   1.12e-3     5.e-2
        5.   1.12e-3     5.e-2
        5.   1.12e-3     5.e-2
        5. 7.7778e-4     5.e-2
        5. 7.7778e-4     5.e-2
        5. 7.7778e-4     5.e-2
        5. 7.7778e-4     5.e-2
        5. 7.7778e-4     5.e-2
        5. 7.7778e-4     5.e-2
        5. 8.7273e-4     5.e-2
        5. 8.7273e-4     5.e-2
        5. 8.7273e-4     5.e-2
        5. 8.7273e-4     5.e-2
        5. 8.7273e-4     5.e-2
        5. 8.7273e-4     5.e-2
        5. 8.7273e-4     5.e-2
        5. 8.7273e-4     5.e-2
        5. 8.7273e-4     5.e-2


Now that we have looked at both files, let's make SFR file out of a STR
1) header line 22,3,0,0,128390.4,.0001,-1,0. Change 22, 3 according to model. Leave the printing flags as default if you don't need them.
2) add the chuck of data about location and length of each reach. The first five integers of each reach record is the same for SFR and STR, copy them over. The fifth number in the SFR file is reach length which is not specified in STR. You have to make something up. Maybe use size of the cell.

1,2,5,1,1,1000.
1,3,5,1,2,1000.   
1,4,5,1,3,1000.
1,4,6,1,4,1000.
1,4,7,1,5,500.
...

- another line about printing flags

3) 3,0,0   change 3 to appropriate total number of segments

4) now enter data about each segment. Ready to do some work... See above for parameter description.As you may have noticed already, unlike SFR there is no information about things like connection to lake, runoff, evaporation, precipitation in the STR file, so yo have to enter these information yourself  One of the important parameter for lake interaction is the 3rd one, -1 here. It means stream is connection to lake number 1. 

1 1 -1 0   691200.0 0.0   0.0   0.0   0.05
5) next line is information about upstream reach of the segment. The streambed top and bottom are included in STR file.

6) next line is information about downstream reach of the segment.

7) now repeat 4) to 6) for every other segment.

Done. Now if you comment out the STR file in the NAM file and add name of SFR file, it should just run.

Leave a comment if you have any questions.

MODFLOW convergence techniques

EDIT: if you use the new MODFLOW-NWT, you probably don't need any of these because NWT does not use wetting-drying. So go for NWT if you can. 

The biggest pain in the ass of groundwater modeling is getting your solution to converge.

 This is the due to the uniqueness of groundwater system - the aquifers can be unconfined, in which case part of the aquifer is not fully saturated with water and the rest is saturated with water. In realty, the unsaturated part still has some residual water content and is called vadose zone. However, MODFLOW does not simulate water content in the vadose zone. It treats it as completely dry. This caused a huge problem when solving the model's linear system because cells has to be constantly switched between dry and wet to satisfy the boundary conditions. And this is often accompanied with abrupt changes making it even harder to converge. The 'rewetting' option can be turned off in the input file but that means once water level in a cell drop below the bottom in a solver iteration it will never come back. You may get a solution in the end but it will not be the correct solution.

Here are a couple of tricks/techniques to better assist you with model convergence.

1. First run it with rewetting on, if the residual fluctuates but the head change is small (around 1 m) then run it again using the previous solution as initial head with rewetting of this time. You should get a solution with good mass balance and correctly reflect the groundwater system you are trying to simulate.

 2. If the above does not work, change all layer type to confined (laycon=0), then run the model. You should be able to obtain a solution with good mass balance easily. (If not, then there is something wrong with the model, check boundary conditions etc). Then using the solution as initial head, progressively change laycon of top layers to 3. For example, change layer 1 to type 3, run it, get a solution, use the new solution as initial head, change layer 1 and 2 to type 3...