fix(SfrCrossSectionUtils): fix for n-point x-section area calculation (…

…MODFLOW-USGS#1172)

* fix(SfrCrossSectionUtils): fix for n-point x-section area calculation

* Changes that accounted for the vertically wetted sides in the wide rectangular channel approach should not have been made.

* updates I originally forgot to transfer when redoing the PR (originally MODFLOW-USGS#1134)

* now that the areas are being calculated correctly, really small widths are being multiplied by really small depths on the intial guess at channel depth. As a result, the iterative stage calculation blows up leading to a temporary fix.

* one more fix I should've included earlier

* goes with e579d10

* More work to do on this PR after conversation with @jdhughes-usgs

* clean-up python to go with reworked fortran changes.

* Bolster discussion in 'Streamflow Routing Package Cross-SectionTable Input File' section as discussed.

* Update doc/mf6io/gwf/sfr.tex

Co-authored-by: langevin-usgs <[email protected]>

* remove errant "("

* more changes after checking answers for a more complex x-section; reran black

---------

Co-authored-by: langevin-usgs <[email protected]>

autotest/cross_section_functions.py

@@ -12,7 +12,7 @@ def calculate_rectchan_mannings_discharge(
 
  """
  area = width * depth
- return conversion_factor * area * depth**mpow * slope**0.5 / roughness
+ return conversion_factor * area * depth ** mpow * slope ** 0.5 / roughness
 
 
 # n-point cross-section functions
@@ -69,17 +69,20 @@ def get_wetted_perimeter(
 
  # -- calculate the wetted perimeter for the segment
  xlen = x1 - x0
+ dlen = 0.0
  if xlen > 0.0:
- if depth > hmax:
+ if depth >= hmax:
  dlen = hmax - hmin
  else:
  dlen = depth - hmin
+ xlen = x1 - x0
  else:
  if depth > hmin:
  dlen = min(depth, hmax) - hmin
  else:
  dlen = 0.0
- return np.sqrt(xlen**2.0 + dlen**2.0)
+
+ return np.sqrt(xlen ** 2.0 + dlen ** 2.0)
 
 
 def get_wetted_area(x0, x1, h0, h1, depth):
@@ -95,8 +98,12 @@ def get_wetted_area(x0, x1, h0, h1, depth):
  if depth > hmax:
  area = xlen * (depth - hmax)
  # -- add the area below zmax
- if hmax != hmin and depth > hmin:
- area += 0.5 * (depth - hmin)
+ if hmax != hmin:
+ if depth >= hmax:
+ area += 0.5 * (hmax - hmin) * xlen
+ elif depth > hmin:
+ x0, x1 = get_wetted_station(x0, x1, h0, h1, depth)
+ area += 0.5 * (depth - hmin) * (x1 - x0)
  return area
 
 
@@ -114,9 +121,6 @@ def wetted_area(
  x0, x1 = x[idx], x[idx + 1]
  h0, h1 = h[idx], h[idx + 1]
 
- # get station data
- x0, x1 = get_wetted_station(x0, x1, h0, h1, depth)
-
  # get wetted area
  a = get_wetted_area(x0, x1, h0, h1, depth)
  area += a
@@ -172,20 +176,24 @@ def manningsq(
  q = 0.0
  for i0 in range(x.shape[0] - 1):
  i1 = i0 + 1
- perimeter = get_wetted_perimeter(x[i0], x[i1], h[i0], h[i1], depth)
+
+ # get station data
+ x0, x1 = get_wetted_station(x[i0], x[i1], h[i0], h[i1], depth)
+
+ perimeter = get_wetted_perimeter(x0, x1, h[i0], h[i1], depth)
  area = get_wetted_area(x[i0], x[i1], h[i0], h[i1], depth)
  if perimeter > 0.0:
  radius = area / perimeter
  q += (
- conv * area * radius**mpow * slope**0.5 / roughness[i0]
+ conv * area * radius ** mpow * slope ** 0.5 / roughness[i0]
  )
  else:
  perimeter = wetted_perimeter(x, h, depth)
  area = wetted_area(x, h, depth)
  radius = 0.0
  if perimeter > 0.0:
  radius = area / perimeter
- q = conv * area * radius**mpow * slope**0.5 / roughness[0]
+ q = conv * area * radius ** mpow * slope ** 0.5 / roughness[0]
  return q
 
 

autotest/test_gwf_sfr_evap.py

@@ -401,11 +401,11 @@ def eval_results(sim):
  # Establish known answer:
  stored_strm_evap = np.array(
  [
- -66.01388942706727,
- -65.99953139524196,
- -65.98517173673817,
- -65.97081044940869,
- -65.95644753110243,
+ -62.17272623,
+ -62.15731943,
+ -62.14191043,
+ -62.12649925,
+ -62.11108587,
  ]
  )
 

autotest/test_gwf_sfr_npoint02.py

@@ -5,6 +5,7 @@
 import pytest
 from framework import TestFramework
 from simulation import TestSimulation
+from cross_section_functions import get_depths
 
 paktest = "sfr"
 
@@ -46,14 +47,14 @@
  "h": np.array([0.0, 0.0], dtype=float),
  }
 
+
 # depth as a function of flow for a wide cross-section
 def flow_to_depth_wide(rwid, q):
  return ((q * roughness) / (conversion_fact * rwid * np.sqrt(slope))) ** 0.6
 
 
 #
 def build_model(idx, ws):
-
  # build MODFLOW 6 files
  name = ex[idx]
  sim = flopy.mf6.MFSimulation(
@@ -219,10 +220,22 @@ def eval_npointdepth(sim):
  obs["INFLOW"], np.abs(obs["OUTFLOW"])
  ), "inflow not equal to outflow"
 
- d = flow_to_depth_wide(
- obs["WIDTH"],
- inflow,
- )
+ d = []
+ for n in range(nper):
+ x0 = 0.0
+ x1 = rwid * (n + 1) # generates absolute widths generated above
+ x = np.array([x0, x1])
+ cdepth = get_depths(
+ inflow,
+ x=x,
+ h=np_data[n]["h"],
+ roughness=roughness,
+ slope=slope,
+ conv=1.0,
+ dd=1e-4,
+ verbose=False,
+ )
+ d.append(cdepth[0])
 
  assert np.allclose(
  obs["DEPTH"], d

autotest/test_gwf_sfr_npoint03.py

@@ -59,8 +59,16 @@
  "r": np.array([10.0, 1.0, 10.0, 10.0], dtype=float),
 }
 
-
+# Cross section depiction
+#
+# | | | |
+# | (left) | (channel) | (right) |
+# | 10.0 | 1.0 | 10.0 | <- "MANFRACTION"
+# | | | |
+# +-----------+-----------+-----------+ y: 0, 0, 0, 0
+# x: 0 1/3 2/3 1
 #
+
 def build_model(idx, ws, base=False):
 
  if base:

doc/ReleaseNotes/v6.5.0.tex

@@ -26,6 +26,7 @@
  \item The input for some stress packages is read in a list format consisting of a cellid, the form of which depends on the type of discretization package, and stress information on each line. The cellid is checked upon reading to ensure that the cell is within the model grid. If the cell is outside the model grid, the program issues an error message and terminates. This cellid check was not implemented when the list was provided from an OPEN/CLOSE binary input file. The program was modified to include this check for both text and binary input.
  \item In some cases, unrecognized keywords and invalid auxiliary input did not terminate with a useful error message. The program was corrected to provide error handling for these cases.
  \item Based on the LAK package input-output instructions (mf6io.pdf), the variable ``connlen must be greater than zero for a HORIZONTAL, EMBEDDEDH, or EMBEDDEDV lake-GWF connection.'' However, a value of zero could be specified and the model would run with no LAK-groundwater exchange. A minor fix was made to enforce connlen to be strictly greater than zero per the input instructions. The error message thrown when connlen is specified as zero was augmented with additional information for assisting the user.
+ \item An SFR channel defined with the n-point cross-section option was calculating the wetted cross-sectional area incorrectly. The cross sectional area for the area of a triangle was being calculated as one-half multiplied by the depth of the channel, as opposed to one-half multiplied by the base width multiplied by the height. As a result, the units in the mannings equation were not correct owing to the missing dimension in the area calculation. The change in the area calculation will slightly alter the solution found using Manning's equation since the cross-sectional area term appears in it. As a result, existing models may reflect slightly different answers in groundwater\/surface-water exchange amounts owing to slight differences in the calculated stream stage. In addition to the fix, some clarifying text, including a new figure, was added to mf6io.pdf.
  % \item xxx
  \end{itemize}
 

doc/mf6io/gwf/sfr.tex

@@ -94,6 +94,16 @@ \subsection{Streamflow Routing Package Cross-Section Table Input File} \label{se
  \label{fig:sfr-n-point}
 \end{figure}
 
+Where irregular cross sections are used to define cross-sectional stream geometries, the wetted perimeter used in Manning’s equation [Equation 7-7, \cite{modflow6gwf}] depends on the number of points defining the cross section and the simulated stage. Using only the minimum number of points (i.e., 2-point cross section), \mf does not include perimeter lengths above the uppermost defined points in the wetted perimeter calculations. For example, the 2-point cross sections shown in fig.~\ref{fig:sfr-n-point-wp}A-C depict the cross-sectional areas (light blue) and wetted perimeters (orange) calculated by the SFR package and used in \cite{modflow6gwf} (Equation 7-7). In applications where the intent is for the wetted perimeter to include the entire lengths of wetted sides, additional points above the maximum anticipated stage should be defined (fig.~\ref{fig:sfr-n-point-wp}D). Note that when the simulated stream stage rises above the points representing the top of the channel, the additional cross-sectional flow area above the defined points will be accounted for but the corresponding wetted perimeter will not extend above the defined points (fig.~\ref{fig:sfr-n-point-wp}E,F).
+
+
+\begin{figure}[ht]
+ \centering
+ \includegraphics[scale=1.0]{../Figures/n-point-cross-section-wetted-perimeter}
+ \caption[Illustrations of variously defined n-point cross-sections that show how wetted perimeter will vary depending on the stage and the number of points used to define the cross-section]{Example irregular cross-section geometries showing the corresponding wetted perimeter based on the number of points that define a cross-section and the simulated stage. (A-C) Wetted perimeters (orange lines) for variously configured 2-point cross-sections. (D-F) Wetted perimeters for variously configurated 4-point cross-sections}
+ \label{fig:sfr-n-point-wp}
+\end{figure}
+
 Cross-Section tables are specified by including file names in the CROSSSECTIONS or PERIOD blocks of the SFR Package for specific reaches. These file names correspond to a Streamflow Routing cross-section table input file. The format of the Streamflow Routing cross-section table input file is described here.
 
 \vspace{5mm}

src/Model/ModelUtilities/SfrCrossSectionUtils.f90

@@ -382,7 +382,11 @@ subroutine get_cross_section_areas(npts, stations, heights, d, a)
  !
  ! -- add the area below dmax
  if (dmax /= dmin .and. d > dmin) then
- a(n) = a(n) + DHALF * (d - dmin)
+ if (d < dmax) then
+ a(n) = a(n) + DHALF * (d - dmin) * xlen
+ else
+ a(n) = a(n) + DHALF * (dmax - dmin) * xlen
+ end if
  end if
  end if
  end do