Merging and Culling the NHDPlus Flowlines

The Mississippi Embayment region intersects six major drainage basins (vector processing units in the NHDPlus), including most of the Lower Mississippi Basin (08) and small parts of surrounding basins (fig. 1). The mapped flowlines within each basin include many more streams than are desirable to simulate in a large regional groundwater model with a coarse resolution. Within the Mississippi Alluvial Plain, most of the mapped headwater streams are ephemeral, flowing only during precipitation or periods of high water. Inspection of the hydrography overlain on satellite imagery suggested that many other mapped streams no longer exist because of land use changes; therefore, the first part of preprocessing was to merge the flowline and associated attribute data from the various drainage basins and cull the dataset to only include a reasonable density of flowlines representing the current stream network within the model area. The term “reasonable” is subjective but in this case represents a balance between including all flowlines that potentially receive groundwater discharge for at least several months of the year and maintaining a workable stream network size for the model solution and coarse grid discretization.

A Python function (cull_flowlines() in the SFRmaker Preprocessing module) was developed to perform the above tasks for a set of NHDPlus files and a polygon shapefile of the active area of the MODFLOW grid. The function outputs a single set of flowlines for the model domain and three associated attribute files (following the structure of NHDPlus version 2):

More details about these files are in McKay and others (2012).

Arbolate sum was selected as a numeric criterion for limiting stream density because it is continuous (unlike stream order) and more easily obtainable than drainage area, especially in the flat terrains of the Mississippi Alluvial Plain region. The NHDPlus flowlines were culled at various arbolate sums ranging from 4 to 30 km and overlain on satellite imagery. Flowlines in headwater areas were visually inspected at various locations around the model area for correspondence with wet channel conditions. Although this method is qualitative and depends on the time of year the satellite imagery was taken, it can provide an approximate sense for the range of appropriate arbolate sum thresholds. A threshold of 20 km was selected, meaning that streams with less than 20 km of upstream channel were excluded. From the remaining streams longer than 20 km, those identified in NHDPlus as intermittent, and with arbolate sums less than 50 km, were also culled. This resulted in a handful of isolated perennial streams that no longer had downstream connections. Inspecting these streams against satellite imagery suggested them to be minor, especially at a 1-mile or 1-km grid resolution. The isolated streams were also culled. The remaining stream network includes about 24 percent of the flowline features mapped in NHDPlus within the MERAS region (fig. 4).

Elevations and Routing at Divergences

The second part of the preprocessing focused on recasting the NHDPlus data to a format readily convertible to SFR Package input while maintaining the model grid-independent structure of the flowline features. A key part of the preprocessing was handling divergences in the stream network, where flow is routed to two or more distributaries. Many distributaries are mapped throughout the Mississippi Alluvial Plain region, but most of these features are either nonexistent or only carry flow intermittently during high water events. Although distributaries in NHDPlus are classified into “main” or “minor” paths at divergences, inspection against the satellite imagery and the most recent, lidar-based DEMs for the region (USGS, 2018c) suggested that the NHDPlus classifications are often inaccurate. The steps described below are illustrated graphically in figure 5, which shows an example location on the Quiver River near Sunflower, Mississippi. The following steps were taken to identify the main channel at each divergence:

Comparison of the sampled DEM elevations with the NHDPlus elevation attribute data revealed a high bias in many of the attribute elevations, especially near diversions. This may be a result of the upstream smoothing process described by McKay and others (2012, p. 123) when it encounters distributaries of unequal elevations such as the example shown in figure 5. To remedy this issue, the 10th percentile elevations obtained from the buffer zonal statistics were assigned to each flowline and then smoothed in the downstream direction to ensure that no flowlines had negative (uphill) streambed slopes.

Finally, routing connections to minor distributaries were removed, and arbolate sums were recomputed for the entire stream network with arbolate sums at minor distributaries starting at zero. In this way, the minor distributaries are treated like headwater streams in that they will only receive flow if the water table is greater than their assigned elevation; otherwise, they are simulated as dry and are not part of the groundwater model solution. Like headwater streams, the first 20 km of minor distributaries were trimmed from the stream network, based on inspection of satellite imagery that revealed most of them to be nonexistent or dry (fig. 5).

The above steps were implemented in the preprocess_nhdplus() function in the SFRmaker Preprocessing module, which reads the culled files produced by cull_flowlines() (or a single set of native NHDPlus version 2 files) and outputs a single shapefile with the preprocessed flowlines and attribute information including routing, channel width (described in the “Stream Reach Properties” section), and elevations. The shapefile is designed to provide a basis for an SFR package on any finite difference grid within the MERAS area.

Intersecting the Flowlines with the Model Grid

SFRmaker develops reaches by breaking the preprocessed flowlines at the model grid cell boundaries, resulting in a flowline fragment (reach) for each intersection between a flowline and a model grid cell (fig. 5a). An r-tree spatial index (Gillies, 2017) is used to speed the testing of intersections between model grid cells and flowlines. The shapely package (Gillies and others, 2017) is used to create the reaches from their parent flowlines. By default, reaches less than 5 percent of a model cell length (50 meters for 1-km cells) are discarded. The remaining reaches are grouped into segments corresponding to their parent flowlines. Stream segments are numbered sequentially with the numbering always increasing in the downstream direction (Prudic and others, 2004). Segments are routed using the attribute information developed in the preprocessing step; reaches are routed (numbered) based on geographic proximity of their end vertices (accounting for discarded reaches). Attribute information from the preprocessed flowlines is transferred to the grid-based reaches by recording the index position of model grid cells intersecting each flowline. This type of operation, where attribute information is transferred from one vector-based feature to another based on location, is often referred to as a spatial join.

Handling of Colocated Reaches

Even with the arbolate sum thresholds described in the preprocessing step, the high stream density in the Mississippi Embayment region resulted in many instances of model cells containing multiple reaches (fig. 5 shows an example). Although the SFR package is capable of computing flows between multiple colocated stream reaches and a model cell representing the groundwater flow solution, having many colocated reaches places an additional burden on the model solver and also greatly increases the size of the SFR package, which in turn slows file input and output. There is little utility in carrying these features in the model solution: all computations for a model cell are performed at the cell center, meaning the spatial component of the interactions between colocated reaches and the groundwater flow solution cannot be resolved. In other words, the level of detail that can be represented in an SFR package is limited by the model grid size. Meaningful representation of colocated streams requires refining the grid so that they are no longer colocated.

To address this issue, a “sparse” version of the SFR package was created with only one reach in each model cell. In each cell, the dominant reach (likely to carry the most flow) was identified based on the largest arbolate sum computed in the preprocessing step. The conductance term (K_vwL/m in eq. 5) was then summed for all reaches within the cell and applied to the dominant reach by increasing the hydraulic conductivity value (K_v; eq. 6). Minor reaches were then removed so that each cell contained no more than one reach, with a conductance value that represented all the original reaches. Routing connections between successive dominant reaches were established using a directed graph of the original reach connections. The one_reach_per_cell option in SFRmaker and Modflow-setup allows this “sparse reach” option to be turned on or off.

where

If values of 1 are assumed for K_v and m, equation 6 reduces to a dimensionless scaling factor (eq. 7) that can be applied to any values of leakance (K_v/m) assigned to an SFR reach. This allows leakance to be updated while maintaining the relative differences in total reach area between model cells that originally contained multiple reaches.

Example Location on the Quiver River

Figure 5 illustrates the preprocessing steps described in the “Preprocessing” section for an example location on the Quiver River west of Greenwood, Miss. Figure 5A shows the original flowlines from NHDPlus overlain on recent satellite imagery. The white numbers indicate the 10th percentile elevations sampled from the lidar-based DEM (in feet above sea level) within a 50-meter buffer of the flowlines. Inspection of the flowlines, air photo, and a shaded relief of the DEM (fig. 5B) suggests that the Quiver River historically followed a meander bend to the west (elevation of 102.27 feet) but now takes a more direct, possibly dredged cutoff (elevation of 97.39 feet). The historical meander bend is coded as the main channel in NHDPlus, but the cutoff was selected instead in the preprocessing based on a lower elevation of 97.39 feet. Inspection of the satellite imagery clearly confirms the cutoff to be the main channel. Other tributaries shown in figure 5A are mapped across what are now aquaculture ponds (green rectangles; for example, the tributary with an elevation of 105.74 feet) or row-crop fields (light brown rectangles).

Figure 5B shows the stream network after removal of flowlines with arbolate sums less than 20 km and intermittent flowlines with arbolate sums less than 50 km, overlain on a shaded relief of the DEM. The remaining (blue) lines all represent streams that are visible in the satellite imagery and DEM. The solid blue lines illustrate the dominant reach in each model cell (with the highest arbolate sum); the dashed lines illustrate minor reaches that are colocated. The red arrows illustrate routing connections between reaches as MODFLOW sees them in the “sparse” version of the SFR package. The connections extend between the centers of the finite difference cells (where all calculations are made) and only occur between successive dominant reaches. The existence of the intervening minor reaches is implied, but they are not included in the SFR package. The physical footprints of the streams in each cell (length and width; including colocated reaches) are represented implicitly via streambed hydraulic conductivity in the leakance term, as described in the “Handling of Colocated Reaches” section. A potential downside of this approach is that flux targets on minor reaches may not have simulated flows that are directly equivalent. However, the large scale and coarse discretization of the model preclude accurate simulation of these smaller streams.

Estimation of Reach Widths with Field Measurements

Stream channel widths were estimated for individual flowlines during the preprocessing step (preprocess_nhdplus() function) and applied to SFR reaches later by mapping COMIDs to reach numbers. Two methods were used to estimate channel widths for flowlines. In the first method, widths were estimated from field measurements based on their arbolate sum. In the second method, widths were sampled onto the flowlines directly from the NARWidth database (Allen and Pavelsky, 2015). The two methods were compared, and preference was given to the second method where data were available.

Field measurements of channel width recorded by USGS personnel during streamflow measurements were obtained from the National Water Information System (NWIS; USGS, 2022) for sites throughout the Mississippi Embayment region, in the period after 1970. The NWIS database includes base flow and flood measurements. Widths associated with the latter often greatly exceed the width of the stream channel during typical base flow conditions. To determine if a width measurement was made during base flow conditions, base flow periods were identified at the closest streamgage with a continuous record: an index station. Tenth percentile flows (Q90; equaled or exceeded 90 percent of the time) were computed for each index station. Width measurements were considered to be representative of base flow when the ratio between the measured flow at the index station and the index station Q90 (Q:Q90) was less than two. For locations with more than one base flow channel width measurement, the measured widths were averaged.

Each width measurement was located on the stream network developed in the preprocessing step and assigned an arbolate sum value. A power-law relation between arbolate sum and width was then developed via a curve fitting analysis, similar to Feinstein and others (2010, p. 266). Figure 6 shows the best-fit line through the log-transformed data.

Although figure 6 shows a clear relation between stream channel width and arbolate sum, considerable scatter about the best-fit line illustrates the amount of heterogeneity in channel width caused by other geomorphic factors. Additionally, streamgages are typically located in areas where the channel is straight and stable and near bridges and other artificial structures. Therefore, the measured channel widths shown in figure 6 may be biased low compared to actual widths throughout the stream network (Park, 1977; Allen and Pavelsky, 2015).

Estimation of Reach Widths from Remote Sensing Data

Widths from NARWidth were sampled onto the stream network by performing a spatial join based on the intersection of the NARWidth line features with 1-km buffer polygons around the preprocessed flowlines. To avoid erroneous overlap between mainstem NARWidth estimates and minor tributaries, flowlines with arbolate sums less than 500 km only received widths from NARWidth lines with at least 50 percent of their length inside of the 1-km buffer.

In addition to the channel widths estimated from the power law relation described above, NARWidth values were also assigned to flowlines where available. Figure 7 provides a comparison of widths from the two sources. The major rivers originating outside of the model area (Mississippi, Ohio, White, and other rivers) are visible in the clusters of points at similar arbolate sum width estimates. The drainage areas for these streams lie mostly outside of the model domain, resulting in large arbolate sums at their inflow points, and relatively small increases in arbolate sum (and therefore small increases in estimated width) as they cross the model domain (except for confluences of major rivers, where there are large step changes in arbolate sum). In other words, the clusters of points representing the largest streams in figure 7 mostly reflect the size of contributing drainage areas upstream from the model domain.

Substantial scatter in the NARWidth values for the largest streams illustrates the wide variability in width along large rivers attributable to natural geomorphic and anthropogenic factors (control structures, dredging, etc.). Arbolate sum-based width estimates, which increase monotonically downstream, cannot capture this variability. The largest scatter is at reservoirs (blue points) along the periphery of the MERAS model domain—features that vary widely in width but have low arbolate sums due to their upstream positions in the stream network.

Smaller widths falling below the 1:1 line indicate a low bias in the arbolate sum-based widths relative to NARWidth. This bias is consistent with channel widths (which form the basis for the arbolate sum estimates) measured at streamgages where the channel is narrower than average and under base flow conditions, whereas the NARWidth estimates represent mean flows. Comparing measured channel widths with NARWidth estimates for the same locations indicates the NARWidth estimates to be about 32 percent higher on average but otherwise have good correspondence (fig. 8).

Manning’s Roughness

The USGS streamflow measurement dataset used for the width estimation was also evaluated for Manning’s roughness coefficient (n). Two approaches were used: direct estimation of n from streambed bottom material type; and back calculation of n from recorded discharge, channel geometry, and an estimated slope obtained from the flowline preprocessing.

In the first approach, a representative n for each stream in the Mississippi Embayment region was estimated using channel bottom material reported by USGS personnel during streamflow measurements. Channel bottom material categories for streamflow measurements are ledge—bedrock—artificial, cobble—boulders, gravel, sand, and silt/mud. Each bottom material was assigned a representative n using the following equation (Arcement and Schneider 1989):

where

Values for n_b were estimated (table 1) for each bottom material category using the ranges listed in table 1 in Arcement and Schneider (1989) for various types of channel bottoms. Values for n₁, n₂, n₃, and n₄ were estimated (table 2) using table 2 in Arcement and Schneider (1989) and assuming:

The estimated n_b for each bottom material category is presented in table 3.

An n was calculated for each measurement location that had at least one measurement where sediment type was reported. Locations with multiple measurements where the bottom material was reported were each assigned an n, and an average of all values was used to estimate a single n for that location.

In the second approach, equation 4 was rearranged to solve for n:

Channel area recorded in the streamflow measurement dataset was divided by the recorded channel width (w) to get an average depth (y); the slope term (S₀) was estimated using the minimum and maximum smoothed elevations for the flowline associated with the streamflow measurement. Figure 9 compares the populations of roughness estimates from the two methods. The distribution of values from the first approach highlights the inherent difficulty in precisely describing the physical factors that contribute to n, an empirical parameter that is probably best estimated from actual flow observations (the second method). However, general agreement between the mean values from the two approaches, and the population densities about the means suggests that the default value of 0.037 used in SFRmaker (and the SFR packages described here) is appropriate for the MERAS region.

Length, Streambed Elevation, and Slope

SFRmaker computes reach lengths from the respective flowline fragments (linestring geometries) using the shapely package. The preprocess_nhdplus() function assigns initial streambed top elevations from the smoothed segment end elevations developed in the preprocessing by linear interpolation to the reach mid-points. Optionally, SFRmaker can also sample minimum elevations from a DEM while creating the stream reaches (typically at a higher resolution than the flowlines; Leaf and others, 2021a, b36); this can be advantageous for finer grid resolutions but was not deemed necessary here. In subsequent development of the SFR reaches, SFRmaker computes streambed slopes as the forward differences of the initial streambed tops. A minimum slope of 0.00001 was enforced to prevent excessive stream depths (eq. 4). Compared to the default value of 0.0001 in SFRmaker, a minimum slope of 0.00001 resulted in improved matches of stream stage on the lowest gradient streams in the Mississippi Delta and is also consistent with the low stream gradients present throughout the Mississippi Alluvial Plain (Leaf and others, 2023b).

Stream Depth

DEMs can be problematic for assigning streambed elevations because the DEM often represents the water surface instead of the streambed. In that case, stream depths computed by equation 4 will be added to the water surface elevation, resulting in unrealistically high simulated stages and, often, simulated leakage out of streams that is not actually occurring.

Equation 4, which does not consider continuity of energy between stream reaches, can also produce inaccurate stream depths, especially when reach properties are averaged over large areas or consolidated from multiple streams. Because streambed elevations are not considered in equation 4, computed depths can result in stages that rise in the downstream direction. These effects are most problematic for the largest streams, which can have the largest computed depths.

To address the above issues with stage, the MERAS 3 model was run once with the initial streambed top elevations to get an estimate of stream depth for each reach. Average depths from the initial run across all stress periods were then subtracted from the original streambed tops to get a revised set of streambed tops. This resulted in computed stages in subsequent runs that were generally close to the original streambed top elevations, which are more representative of actual stages. This step was performed by a third Python script (“adjust_sfr_streambed_tops_set_stages.py” in Leaf and others, 2023a) that reads SFR Package output from the model solution and writes a new SFR package with updated streambed tops.

Streambed K_v and Thickness

Although the SFR Package allows for input of both K_v and thickness (m), these parameters are often not well understood at the scale of the model, and simulation of interactions between groundwater and surface water often lump them together conceptually into a resistance (m/K_v; units of days) or leakance (K_v/m; units of days⁻¹) term that describes the total resistance to flow between the aquifer and a surface water body. Initial leakance values for the SFR packages described here were derived from the results of waterborne electrical resistivity surveys of streams in the MERAS region (Killian, 2018; Adams and others, 2019). The results of Adams and others (2019) consist of depth-integrated estimates of streambed K_v (see Killian, 2018 for methodology) and the total survey depth (the thickness of sediments contributing to the K_v values) at point locations along the surveyed streams. Leakance values were computed from these data at each survey point. The population of leakance estimates is shown in figure 10.

The point estimates of leakance within each model grid cell were identified using a spatial index and the geometric mean computed for each cell (fig. 11). K_v was then assigned to the SFR reaches by multiplying the mean leakances by the area factors from equation 7, assuming a uniform streambed thickness of 1 foot for all SFR reaches. SFR reaches with no leakance estimates were assigned the geometric mean for the whole Adams and others (2019) dataset of 1.22 days⁻¹. The initial K_v values were assigned in the parameter estimation setup for the MERAS 3 and Mississippi Delta models (“setup_parameters.py” script in Leaf and others, 2023a) and then subsequently adjusted in the parameter estimation (Leaf and others, 2023b).

Specified Inflows to the Stream Network and Streamflow Observations

The MERAS model domain receives many surface water inflows from the greater Mississippi watershed. Partial streamgage records were available for five sites near the model perimeter. Additionally, streamflow was specified at the outlets of the four flood control reservoirs that feed the Tallahatchie/Yazoo River system in the Mississippi Delta (Arkabutla, Enid, Sardis, and Grenada Lakes; Tim Rodgers, U.S. Army Corps of Engineers, written commun., March 19, 2020; fig. 1). Remaining inflows for 136 additional ungaged streams and missing flows at gaged locations dating back to July 2000 were estimated using a statistical model based on the random forest method (Dietsch and others, 2023). The random forest model compares variables including monthly precipitation, minimum and maximum temperature, evapotranspiration, catchment elevations, and geology with monthly total flow and base flow values at gaged sites throughout the Mississippi Embayment and surrounding region; and then makes predictions of monthly total flows and base flows at ungaged locations and times. A description of the random forest model is in Dietsch and others (2023). Locations of inflows to the SFR network are shown in figure 12. SFRmaker automates the process of assigning inflows to SFR reaches using NHDPlus COMIDs and proximity to the reach linestrings as grid-independent references, allowing for changes in the stream network or the application of inflows to inset models with minimal additional labor. The same process is used to locate flow or stage observations within the SFR network.

Surface Water Runoff

Comparing simulated and observed streamflow and stages, and a large runoff component in the surface water balance, indicated a need to simulate total streamflow to adequately represent surface-water availability and groundwater/surface water gradients (Leaf and others, 2023b). Total streamflow was simulated by specifying monthly mean estimates of surface water runoff from the Soil-Water-Balance code (SWB; Westenbroek and others, 2018) simulation of Nielsen and Westenbroek (2023). Gridded daily runoff output from the SWB simulation (including the “runoff” and “rejected net infiltration” variables) was aggregated to monthly mean values and NHDPlus catchment polygons (which each correspond to a COMID and a flowline) using the swb_runoff_to_csv() function in the SFRmaker Preprocessing module. swb_runoff_to_csv() produces a table of transient runoff totals for each COMID, which was then input to Modflow-setup and SFRmaker during the model build step (fig. 3). By retaining the NHDPlus COMID associated with each stream reach, SFRmaker can map runoff totals from the catchment polygons to stream reaches. Monthly runoff totals were apportioned evenly among the reaches associated with each catchment.

Source Hydrography

Figure 5 and the accompanying text highlight some of the deficiencies in the 1:100,000 scale resolution NHDPlus hydrography data for the Mississippi Embayment region and the general challenge of accurately mapping streams in a relatively flat, highly altered landscape where land use and natural geomorphic processes continue to change the stream network. NHDPlus High Resolution (Moore and others, 2019) offers similar functionality to NHDPlus version 2 at a 1:24,000 resolution, albeit at the cost of managing appreciably more files and data. The workflow described in this report could be used with NHDPlus High Resolution hydrography to produce an improved SFR network. This would most benefit inset models at resolutions that are fine enough to incorporate the additional detail.

Similarly, in other regions, improved hydrography may be available, including hydrography generated via flow accumulation on a DEM (Gardner and others, 2018). Regardless of the source, hydrography input to the SFR package creation step (fig. 3) must include information on routing. The user must supply this information if using hydrography not from NHDPlus. Future developments may relax some of these requirements. For example, routing connections could be developed from flowline start/end proximity.

Limitations of the Preprocessing Approach

The preprocessing methods described in this report are intended to reduce the source hydrography, which in this application contained many streams that no longer exist or are dry under low-flow conditions, into input that is tractable for a numerical solution covering a large area. A primary goal was to automate this process to facilitate future updates to the MERAS 3 model and the creation of inset models with minimal additional labor. Limiting stream density using an arbolate sum criterion based on ad hoc comparison to satellite imagery may result in areas within the model domain where streams are under or overrepresented. Inset models focusing on smaller streams within the network have the most potential to be affected by this shortcoming and should be checked accordingly. The arbolate sum thresholds in the preprocessing routine can be adjusted to increase or decrease the number of streams represented.

Similarly, divergences may be present where the main channel was not correctly identified. Important locations in the model should be checked by the user. A “known connections” input argument to the preprocess_nhdplus() function allows routing connections to be specified manually by the user.

Finally, the preprocessed hydrography produced by the SFRmaker preprocessing module can be edited manually in a GIS environment to correct any remaining deficiencies before building the SFR package. This can greatly facilitate rapid edits to the SFR package throughout a modeling project.

Effects of Stream Reach Consolidation

Use of the “sparse” option on a coarse model grid can result in many reaches that are only implicitly represented in connections between subsequent dominant reaches. For example, a tributary to a major stream would not have any representation (and therefore no groundwater exchange) within any cells that also include the major stream. Streamgages on the tributary within these cells therefore do not have an exact output equivalent from the model. Observed flows at such locations can only be compared to the closest upstream or downstream reach that is in the model. This problem can be eliminated by selecting a grid discretization fine enough to represent the streams of interest. Coarse regional models such as MERAS 3 are generally unsuitable for representing small streams in detail.

The stream reach consolidation described by equations 6 and 7 also violates the conceptual model for computing stream depth (eq. 4). Reaches representing multiple streams have leakance values that represent groundwater inflows or outflows for all consolidated reaches but widths and slopes that only represent the dominant reach. Some of these reaches may have simulated depths (and therefore stages) that are unrealistic. This effect was mitigated by subtracting stream depths from an initial model run from the original streambed top elevations (resulting in simulated stages that are close to the water surface elevations obtained from DEMs).

Limitations Associated with the MERAS 3 Model

The MERAS 3 application of the workflow described in this report is intended to accurately represent the effects of streams on the Mississippi River alluvial aquifer at a level of detail appropriate to the regional scale and 1-km model grid resolution. Limitations inherent to this application, but not necessarily to the automated SFR construction or MODFLOW SFR Package formulation, are discussed below.

Model Spatial and Temporal Resolution

Although more detail is available from the lidar-based DEM, waterborne resistivity surveys, and other data, the coarse grid resolution of the MERAS model only allows for a single property value to be assigned to each SFR reach representing streams within that square kilometer grid cell. At many SFR reaches, property values represent more than one physical stream. Figure 5 illustrates this reduction in detail graphically, both in terms of the consolidation of streams and streambed elevations. The outcome of this simplification is that although the model solution may adequately capture groundwater flow patterns and the overall water balance at a regional scale, the details of interactions between groundwater and surface water at scales smaller than a few grid cells cannot be resolved. This includes heterogeneity in these interactions within each model cell and the effects of pumping wells near streams. Similarly, hydrologic events occurring at scales of less than a calendar month cannot be accurately represented (Leaf and others, 2023b).

Effects of Setting Streambed Tops Using Depths Estimated from an Initial Model Run

Initial streambed top elevations derived from the lidar-based DEM were adjusted by subtracting stream depths simulated from an initial model run. Ideally, this produces simulated stream stages (streambed top plus simulated depth) that are close to the DEM elevations and that vary smoothly in space, as irregularities in simulated depths created by equation 4, are canceled in the subtraction. However, if large temporal variations in stage are simulated, the depths used to adjust streambed top need to be chosen carefully, or another method of setting streambed top might be preferable. Incorporation of stage observations can help ensure that stream stages are well represented (Leaf and others, 2023b).

Limitations Associated with the MODFLOW SFR Formulation

Although the MODFLOW SFR Package simulates continuity in flow throughout the stream network, continuity of energy is not considered. Stream stage is computed independently at each reach based on the mass balance for that reach (eq. 4). Hydraulic gradients within the stream channel and the effects of momentum are not considered. Therefore, simulated stages may rise in the downstream direction. This effect was mostly mitigated in the MERAS 3 application by subtracting stream depths from an initial model run from the original streambed top elevations. Other important limitations of the SFR formulation are discussed by Prudic and others (2004) and Langevin and others (2017).