feat(sfr): Add picard iteration loop in sfr_fc() (#440)

Also cleanup sparse and remove as a sfrtype variables.

autotest/test_gwf_ts_sfr01.py

@@ -64,7 +64,7 @@ def build_model(ws, name, timeseries=False):
                                  nper=nper, perioddata=tdis_rc)
     # create iterative model solution and register the gwf model with it
     ims = flopy.mf6.ModflowIms(sim,
-                               print_option='NONE',
+                               print_option='SUMMARY',
                                outer_dvclose=hclose,
                                outer_maximum=nouter,
                                under_relaxation='NONE',
@@ -175,6 +175,7 @@ def build_model(ws, name, timeseries=False):
     budpth = '{}.{}.cbc'.format(name, paktest)
     cnvgpth = '{}.sfr.cnvg.csv'.format(name)
     sfr = flopy.mf6.ModflowGwfsfr(gwf,
+                                  maximum_picard_iterations=1,
                                   auxiliary=auxnames,
                                   print_input=True,
                                   budget_filerecord=budpth,

doc/ReleaseNotes/ReleaseNotes.tex

@@ -149,6 +149,7 @@ \section{History}
 		\item Add OUTER\_DVCLOSE and INNER\_DVCLOSE variables to replace OUTER\_HCLOSE and INNER\_HCLOSE variables. Warning messages will be issued if OUTER\_HCLOSE and/or INNER\_HCLOSE variables are specified. OUTER\_HCLOSE and INNER\_HCLOSE variables will eventually be deprecated. 
 		\item Add option to scale drain conductance over a user-defined range (drainage depth). Linear-conductance scaling is used with the Standard Formulation. Cubic-conductance scaling is used with the Newton-Raphson Formulation. The additional drainage depth variable is specified as an auxiliary variable and AUXDEPTHNAM is used to identify the auxiliary variable defining the drainage depth. The cubic-conductance scaling can be used as a replacement for the groundwater seepage option in the UZF Package. The scaled drainage conductance option can also be used to represent vertical seepage faces and improve model convergence in cells where simulated heads fluctuate around the elevation where the drain begins to discharge groundwater.
 		\item Add timeseries support for the reach upstream fraction variable in the SFR package.
+		\item Add Picard iterations for the SFR package to minimize differences in SFR package results between subsequent GWF Picard (non-linear) iterations as a result of non-optimal reach numbering. The number of SFR package Picard iterations can be controlled by specifying the maximum number of Picard  iteration to be used in the OPTIONS block (MAXIMUM\_PICARD\_ITERATIONS). If reaches are numbered in order, from upstream to downstream, MAXIMUM\_PICARD\_ITERATIONS can be set to 1 to reduce model run time. Specifying  MAXIMUM\_PICARD\_ITERATIONS to 1 will result in identical SFR package performance to previous versions of MODFLOW 6.
 	\end{itemize}
 
 	\textbf{\underline{BUG FIXES AND OTHER CHANGES TO EXISTING FUNCTIONALITY}} \\

doc/mf6io/mf6ivar/dfn/gwf-sfr.dfn

@@ -243,13 +243,21 @@ optional true
 longname
 description REPLACE mover {'{#1}': 'SFR'}
 
+block options
+name maximum_picard_iterations
+type integer
+reader urword
+optional true
+longname SFR picard iterations
+description value that defines the maximum number of Streamflow Routing picard iterations allowed when solving for reach stages and flows as part of the GWF formulate step. Picard iterations are used to minimize differences in SFR package results between subsequent GWF picard (non-linear) iterations as a result of non-optimal reach numbering. If reaches are numbered in order, from upstream to downstream, MAXIMUM\_PICARD\_ITERATIONS can be set to 1 to reduce model run time. By default, MAXIMUM\_PICARD\_ITERATIONS is equal to 100.
+
 block options
 name maximum_iterations
 type integer
 reader urword
 optional true
 longname SFR Newton-Raphson iterations
-description value that defines the maximum number of Streamflow Routing Newton-Raphson iterations allowed for a reach. By default, MAXSFRIT is equal to 100.
+description value that defines the maximum number of Streamflow Routing Newton-Raphson iterations allowed for a reach. By default, MAXIMUM\_ITERATIONS is equal to 100.
 
 block options
 name maximum_depth_change

doc/mf6io/mf6ivar/md/mf6ivar.md

@@ -498,7 +498,8 @@
 | GWF | SFR | OPTIONS | OBS6 | KEYWORD | keyword to specify that record corresponds to an observations file. |
 | GWF | SFR | OPTIONS | OBS6_FILENAME | STRING | name of input file to define observations for the SFR package. See the ``Observation utility'' section for instructions for preparing observation input files. Table \ref{table:obstype} lists observation type(s) supported by the SFR package. |
 | GWF | SFR | OPTIONS | MOVER | KEYWORD | keyword to indicate that this instance of the SFR Package can be used with the Water Mover (MVR) Package.  When the MOVER option is specified, additional memory is allocated within the package to store the available, provided, and received water. |
-| GWF | SFR | OPTIONS | MAXIMUM_ITERATIONS | INTEGER | value that defines the maximum number of Streamflow Routing Newton-Raphson iterations allowed for a reach. By default, MAXSFRIT is equal to 100. |
+| GWF | SFR | OPTIONS | MAXIMUM_PICARD_ITERATIONS | INTEGER | value that defines the maximum number of Streamflow Routing picard iterations allowed when solving for reach stages and flows as part of the GWF formulate step. Picard iterations are used to minimize differences in SFR package results between subsequent GWF picard (non-linear) iterations as a result of non-optimal reach numbering. If reaches are numbered in order, from upstream to downstream, MAXIMUM\_PICARD\_ITERATIONS can be set to 1 to reduce model run time. By default, MAXIMUM\_PICARD\_ITERATIONS is equal to 100. |
+| GWF | SFR | OPTIONS | MAXIMUM_ITERATIONS | INTEGER | value that defines the maximum number of Streamflow Routing Newton-Raphson iterations allowed for a reach. By default, MAXIMUM\_ITERATIONS is equal to 100. |
 | GWF | SFR | OPTIONS | MAXIMUM_DEPTH_CHANGE | DOUBLE PRECISION | value that defines the depth closure tolerance. By default, DMAXCHG is equal to $1 \times 10^{-5}$. |
 | GWF | SFR | OPTIONS | UNIT_CONVERSION | DOUBLE PRECISION | value (or conversion factor) that is used in calculating stream depth for stream reach. A constant of 1.486 is used for flow units of cubic feet per second, and a constant of 1.0 is used for units of cubic meters per second. The constant must be multiplied by 86,400 when using time units of days in the simulation. |
 | GWF | SFR | DIMENSIONS | NREACHES | INTEGER | integer value specifying the number of stream reaches.  There must be NREACHES entries in the PACKAGEDATA block. |
@@ -630,7 +631,7 @@
 | GWF | UZF | OPTIONS | SQUARE_GWET | KEYWORD | keyword specifying that groundwater ET will be simulated by assuming a constant ET rate for groundwater levels between land surface (TOP) and land surface minus the ET extinction depth (TOP-EXTDP). Groundwater ET is smoothly reduced from the PET rate to zero over a nominal interval at TOP-EXTDP. |
 | GWF | UZF | OPTIONS | SIMULATE_GWSEEP | KEYWORD | keyword specifying that groundwater discharge (GWSEEP) to land surface will be simulated. Groundwater discharge is nonzero when groundwater head is greater than land surface. |
 | GWF | UZF | OPTIONS | UNSAT_ETWC | KEYWORD | keyword specifying that ET in the unsaturated zone will be simulated as a function of the specified PET rate while the water content (THETA) is greater than the ET extinction water content (EXTWC). |
-| GWF | UZF | OPTIONS | UNSAT_ETAE | KEYWORD | keyword specifying that ET in the unsaturated zone will be simulated simulated using a capillary pressure based formulation. Capillary pressure is calculated using the Brooks-Corey retention function. |
+| GWF | UZF | OPTIONS | UNSAT_ETAE | KEYWORD | keyword specifying that ET in the unsaturated zone will be simulated using a capillary pressure based formulation. Capillary pressure is calculated using the Brooks-Corey retention function. |
 | GWF | UZF | DIMENSIONS | NUZFCELLS | INTEGER | is the number of UZF cells.  More than one UZF cell can be assigned to a GWF cell; however, only one GWF cell can be assigned to a single UZF cell. If more than one UZF cell is assigned to a GWF cell, then an auxiliary variable should be used to reduce the surface area of the UZF cell with the AUXMULTNAME option. |
 | GWF | UZF | DIMENSIONS | NTRAILWAVES | INTEGER | is the number of trailing waves.  A recommended value of 7 can be used for NTRAILWAVES.  This value can be increased to lower mass balance error in the unsaturated zone. |
 | GWF | UZF | DIMENSIONS | NWAVESETS | INTEGER | is the number of wave sets.  A recommended value of 40 can be used for NWAVESETS.  This value can be increased if more waves are required to resolve variations in water content within the unsaturated zone. |

doc/mf6io/mf6ivar/tex/gwf-sfr-desc.tex

@@ -41,7 +41,9 @@
 
 \item \texttt{MOVER}---keyword to indicate that this instance of the SFR Package can be used with the Water Mover (MVR) Package.  When the MOVER option is specified, additional memory is allocated within the package to store the available, provided, and received water.
 
-\item \texttt{maximum\_iterations}---value that defines the maximum number of Streamflow Routing Newton-Raphson iterations allowed for a reach. By default, MAXSFRIT is equal to 100.
+\item \texttt{maximum\_picard\_iterations}---value that defines the maximum number of Streamflow Routing picard iterations allowed when solving for reach stages and flows as part of the GWF formulate step. Picard iterations are used to minimize differences in SFR package results between subsequent GWF picard (non-linear) iterations as a result of non-optimal reach numbering. If reaches are numbered in order, from upstream to downstream, MAXIMUM\_PICARD\_ITERATIONS can be set to 1 to reduce model run time. By default, MAXIMUM\_PICARD\_ITERATIONS is equal to 100.
+
+\item \texttt{maximum\_iterations}---value that defines the maximum number of Streamflow Routing Newton-Raphson iterations allowed for a reach. By default, MAXIMUM\_ITERATIONS is equal to 100.
 
 \item \texttt{maximum\_depth\_change}---value that defines the depth closure tolerance. By default, DMAXCHG is equal to $1 \times 10^{-5}$.
 

doc/mf6io/mf6ivar/tex/gwf-sfr-options.dat

@@ -11,6 +11,7 @@ BEGIN OPTIONS
   [TS6 FILEIN <ts6_filename>]
   [OBS6 FILEIN <obs6_filename>]
   [MOVER]
+  [MAXIMUM_PICARD_ITERATIONS <maximum_picard_iterations>]
   [MAXIMUM_ITERATIONS <maximum_iterations>]
   [MAXIMUM_DEPTH_CHANGE <maximum_depth_change>]
   [UNIT_CONVERSION <unit_conversion>]

doc/mf6io/mf6ivar/tex/gwf-uzf-desc.tex

@@ -47,7 +47,7 @@
 
 \item \texttt{UNSAT\_ETWC}---keyword specifying that ET in the unsaturated zone will be simulated as a function of the specified PET rate while the water content (THETA) is greater than the ET extinction water content (EXTWC).
 
-\item \texttt{UNSAT\_ETAE}---keyword specifying that ET in the unsaturated zone will be simulated simulated using a capillary pressure based formulation. Capillary pressure is calculated using the Brooks-Corey retention function.
+\item \texttt{UNSAT\_ETAE}---keyword specifying that ET in the unsaturated zone will be simulated using a capillary pressure based formulation. Capillary pressure is calculated using the Brooks-Corey retention function.
 
 \end{description}
 \item \textbf{Block: DIMENSIONS}