Abstract
The adaptive hydraulics (AdH) numerical code was applied to study tidal propagation in the Lower Columbia River (LCR) estuary. The results demonstrate the readiness of this AdH model towards the further study of hydrodynamics in the LCR. The AdH model accurately replicated behavior of the tide as it propagated upstream into the LCR system. Results show that the MSf tidal component and the M4 overtidal component are generated in the middle LCR and contain a substantial amount of tidal energy. An analysis was performed to determine the causes of MSf tide amplification, and it was found that approximately 80% of the amplification occurs due to nonlinear interaction between the M2 and the S2 tidal components.
1. Introduction
A two-dimensional (2D) adaptive hydraulics (AdH) model was constructed to study tidal propagation through the Lower Columbia River (LCR) and the Willamette River systems. Correct application of a numerical hydrodynamic model to the LCR is extremely challenging due to several factors such as the narrow inlet, deep navigation channel, the presence of several Pile Dikes, and the hydrodynamic influence of the Bonneville Lock and Dam (L&D) and the Willamette river inflows. In addition, as the tide progresses upstream, the transfer of energy from the lower harmonic components to the higher components occurs. This conversion will be modeled incorrectly if the model is applied improperly. The AdH model was constructed and validated to several water surface gages throughout the system. A harmonic analysis of the model results and field observations was performed to ascertain the energy distribution; future analyses will include the velocity and discharge data available for this system. This harmonic analysis consisting of amplitudes as well as phases ensured that if the lower and higher harmonic constituents (overtidal components) are validated, the remaining energy in the water column must pass through as velocity head. Because the focus of this work was tidal propagation and not salinity propagation a three-dimensional (3D) analysis was not required.
2. Model Description
AdH is a US Army Corps of Engineers (USACE) developed hydrodynamics and transport code. AdH is a multiphysics numerical code and incorporates capabilities for the numerical simulation of saturated and unsaturated groundwater flow and the full Navier-Stokes equation as well as 2D/3D shallow water equations. AdH discretizes the domain on an unstructured linear triangle mesh and utilizes the implicit finite element method (FEM) to solve for primary variables. A basic description of the 2D module of the AdH model is provided below; interested readers are referred to http://adh.usace.army.mil/ for a more detailed description of model equations, solution procedures, and for downloads.
The 2D shallow-water equations are a result of the vertical integration of the equations of mass and momentum conservation for incompressible flow under the hydrostatic pressure assumption [1] Written in conservative form, the 2D shallow water equations are where where is the fluid density, is the gravitational acceleration, is the bed elevation, and is the bed shear stress drag, where the subscript () indicates the direction ( and ), is the flow depth, is the component of velocity, is the component of velocity, and the ’s are Reynolds stresses due to turbulence, where the first subscript indicates the direction, and the second indicates the face on which the stress acts.
The Reynolds stresses are determined using the Boussinesq approach to the gradient in the mean currents: where = kinematic eddy viscosity (which varies spatially).
The AdH shallow-water equations are placed in conservative form to ensure mass conservation and balance of momentum and pressure across an interface. This results in a locally mass conservative model [2, 3].
The equations are discretized using the finite element approach with the velocities and depth being represented as linear polynomials on each element. AdH utilizes a streamline-upwind Petrov-Galerkin (SUPG) scheme similar to that reported in Berger and Stockstill [4] and patterned after previous work by Hughes and Brooks [5], Moretti [6], Gabutti [7], and Steger and Warming [8]. Since the finite element scheme is not the primary focus of this paper, a more in-depth description of this method is omitted.
3. Model Application History
The AdH model has been extensively applied to various environments such as estuarine as demonstrated in Martin et al. [9], McAlpin et al. [10], and Tate et al. [11]; riverine as demonstrated in Stockstill and Vaughan [12]; dam break as demonstrated in Savant et al. [1]; and marshes as demonstrated in Sanborn et al. [13] among other.
4. Lower Columbia River Estuary
The LCR estuary consists of deep channels meandering past shallow flats and islands in a wide coastal plain type estuary. The tides in the estuary are of the mixed type with two high and two low waters per day [14, 15]. Freshwater inflows to the estuary are primarily from the Columbia River and the Willamette River; these freshwater flows range from a cumulative of 4,000 to 15,000 m3.
The hydrodynamics of the LCR estuary have been widely studied primarily to analyze the navigation channel [14, 15], Columbia River plume [16], wave effects [17], and ecology [18]. However, the tidal propagation within the estuary has not been studied extensively.
5. Mesh and Bathymetry
A numerical model mesh for the LCR and Willamette River system was generated using the surface water modeling system [19] as shown in Figure 1; this mesh extends from approximately 33 km offshore (from the Columbia River inlet) to Bonneville L&D and to Willamette Falls on the Willamette river. The mesh resolution varies significantly throughout the mesh with a higher resolution of ~10 m2 provided near the Pile Dikes (Figure 2) and the coarsest resolution of ~40,000 m2 at the ocean boundary.
The bathymetric data used in this modeling study were provided by the Portland District of the USACE and represented the latest LiDAR and survey data available. No modifications were performed on this data prior to or after inclusion in the mesh. The horizontal datum utilized for the mesh was North American Datum of 1983 (NAD 83) and State Plane Coordinate System Oregon North; the vertical datum was the North American Vertical Datum of 1988 (NAVD 88). A representation of the mesh bathymetry is provided in Figure 1.
6. Period of Simulation and Boundary Conditions
The boundary conditions used to drive the model consisted of a tidal boundary at the ocean (Figure 3), river inflows from the Columbia River at Bonneville L&D (Figure 4(a)), Willamette River at Willamette Falls (Figure 4(b)), Cowlitz River (Figure 4(c)), Sandy River (Figure 4(d)), and Louis River (Figure 4(e)). The tidal boundary consisted of the observed tidal signal obtained from the National Oceanic and Atmospheric Administration (NOAA) gage at Astoria, OR (station Id: 9439040). The Astoria gage is located inside the LCR estuary and therefore the observed signal was lowered by 0.35 m prior to application as the tidal boundary (Figure 3). The value of 0.35 m was obtained through trial and error by matching the model predicted tide with the observed tide at Astoria. The river discharges were obtained through USACE and United States Geological Survey (USGS) gages for the period of simulation.
(a)
(b)
(c)
(d)
(e)
The period of simulation was from the July 1st 2008 through December 20th 2008. This period was selected because it consisted of no discernable peaks or valleys in the inflows, allowing for the analysis of the tide without undue freshwater influences (note that the high inflows shown in Figures 4(c), 4(d), and 4(e) occur after December 20th 2008).
A Manning’s “” value of 0.025 was specified everywhere in the domain except the region designated as Pile Dikes. The bed sediment in the LCR estuary is primarily sand, comprising as much as 85% of all bed sediment [14]; the Manning’s “” value of 0.025 is close to the value attributed to sand. The energy loss due to Pile Dikes was represented using the unsubmerged vegetation formulation described in [20]. The pile stem density was specified as 3.3 piles/m2, with a pile diameter of 0.3 m and an equivalent roughness height of 0.07 m. Pile density and diameter was determined through drawings provided by the Portland District of the USACE.
The kinematic eddy viscosity was specified using a Smagorinski formulation [21] with a coefficient of 0.05 in the following form: where kinematic eddy viscosity, = Smagorinski coefficient, = element size dependent length scale, direction velocity component and direction velocity component, and “” indicates the multiplicative operator.
7. Model Validation and Results
AdH results for the period of simulation were compared to the observed water surface elevations at several NOAA and USGS gages in the system (Figure 1). A total of 9 gage stations were available (Figure 1). Detailed results will be provided for 4 gages in the system (covering the system ocean to the upstream boundary). Results for other gages will only be in the form of harmonic constituent analysis and relative errors for the sake of compactness.
Figures 5(a), 5(b), 5(c), and 5(d) present the water surface elevation comparisons for the Astoria, Long View, Saint Helens, and Vancouver gages. Figures 6(a), 6(b), 6(c), and 6(d) present the box form of comparisons, in these figures the closer the data points are to the diagonal line, the better the representation of the field by the model.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
Figures 7(a), 7(b), 7(c), and 7(d) present the amplitudes of the dominant tidal components for the Astoria, Long View, Saint Helens, and Vancouver gages, respectively. Figures 8(a), 8(b), 8(c), and 8(d) present the phases for the dominant tidal signals at the same 4 gage locations. The phases computed from ADH results are very close to the dominant M2, K1, S2, and O1 tidal components throughout the system with the largest error of ~10° in the M2 component occurring at Vancouver.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
Table 1 lists the amplitudes for the dominant M2, K1, S2, O1, M4, and Msf components for all gages in the numerical domain. The largest absolute error of 0.05 m in the amplitude is computed at Vancouver for the M2 component and represents approximately 6% of the tidal amplitude at Vancouver.
Table 2 presents the results of a statistical analysis done on the AdH model results and the field observations. The statistical analysis consisted of determining the root mean square error (RMSE), the correlation coefficient, the Nash-Sutcliffe coefficient, and the Willmott coefficient. The RMSE is an oft-used indicator of the differences between values predicted by a model and the values actually observed in the field [22]. The value of RMSE varies from 0 to 1, indicating excellent replication and no replication, respectively. The correlation coefficient provides a single number to represent the closeness of relation between two related variables [23]. The value of the correlation coefficient varies between +1 and −1, with a value of +1 representing direct correlation and −1 representing no correlation. The Nash-Sutcliffe coefficient is a function of the ratio of the variation in the data not explained by the model to the variation in the observed data [24]. The value of the Nash-Sutcliffe coefficient varies between and +1, with +1 indicating perfect model replication of the field observations. The Willmott coefficient is very similar to the Nash-Sutcliffe coefficient except that it includes the variation in the model about the observed mean value whereas the Nash-Sutcliffe only involves the variation of the observed values about their mean [25].
All of the gages show good statistical parameter comparisons, with the results at Skamokawa showing the worst performance. This gage is on an off-channel and not highly resolved in the model; this could explain the slightly worse performance when compared to the other gages. There was one USGS water surface elevation gage available on the Willamette River at Portland. The statistical analysis on this gage shows good statistical comparisons with the ADH model results. This gage is close to the Willamette and Columbia River confluence and indicates that the tidal exchange is being correctly represented between the Columbia and Willamette Rivers.
8. Analysis of Tidal Transformation
This section will provide an analysis of the tidal transformation and distortion as it propagates into the LCR estuary system.
At Astoria the water surface exhibits a relatively pure symmetric tidal signal; however, by Longview the tide starts to become more complex due to shallow water effects as well as the influence of the river inflows. The transfer of energy from the lower order tidal components to the higher order constituents results in this tidal complexity; this effect is more obvious in the troughs as the asymmetric effects of friction are more significant at shallower depths that occur at ebb tide. A result of this is the transfer of energy from the M2 component to the M4 component; this should be evident in a harmonic decomposition of the tide. Figures 7(a), 7(b), 7(c), and 7(d) show an amplitude analysis of the tide at the four aforementioned stations, and it shows that the amplitude of the M4 component increases as the tide propagates up through the system compared to the amplitude of the M2 component. In addition, the amplitude plots also show the magnitudes of the components when compared to the observed; all the dominant signals show a close match between the primary components as well as the overtides. The close match of the overtides is a further indication that the model is representing the field adequately and that the friction values specified in the model are correct and the energy is being divided between the water head and the velocity head correctly.
The celerity of the tidal wave is a function of the flow depth implying that at ebb the wave travels at a slower speed when compared to celerity of the flood tide wave. According to [26] the energy loss is proportional to the square of the total current speed; this indicates that energy loss during flood will be greater than the energy loss during the ebb tide when the wave celerity is slower. This is seen in Figures 5(b), 5(c), and 5(d) where the trough has a higher water surface elevation when compared to the water surface elevation at Astoria (Figure 5(a)).
At Astoria the model computes a higher amplitude M2 tide when compared to the observed; this is an artifact of the 2D treatment of the system. For locations where river systems interact with the tide, the behavior of the river and the saline tidal water is governed by 3D processes. During ebb the freshwater and the tidal propagation are in the same direction. Because the energy loss is directly proportional to square of the speed of the tidal current, during ebb there is a higher energy loss compared to the flood where the direction of the freshwater and tidal propagation is opposite resulting in a lower flood velocity. The 2D model cannot account for this; hence, the energy that would have been transferred to the higher tidal components from M2 in a 3D representation stays within the M2 component resulting in higher M2 amplitude. At other stations the dominant amplitudes differ by less than 0.025 m between the model and the observed indicating that the 3D effects are primarily important only in the lower LCR system.
Figures 9(a), 9(b), 9(c), 9(d), 10(a), and 10(b) show the behavior of the dominant tidal constituents as the tide progresses upstream into the LCR estuary. The four predominant constituents at the mouth of LCR estuary are M2, K1, S2, and O1, respectively. As the tide progresses inland the tide is quickly distorted and energy transferred to the MSf (or the shallow water equivalent MS) tide as well as the M4 overtidal component. However, the MSf component is clearly more amplified when compared to the M4 overtides.
(a)
(b)
(c)
(d)
(a)
(b)
The expected amplitude of the MSf tidal component is described in [27] as where = amplitude of MSf tide, = angular speed of MSf tide (1.0159°/solar hour [26], and = geographical latitude (45.5° for LCR)).
Utilizing the above expression, the amplitude for MSf is computed to be very close to 0.0 m compared to 0.1 m observed and modeled in LCR estuary.
This amplification of the MSf tide can be attributed to one or a combination of three forcings (a) resonance, (b) nonlinear interaction between M2 and S2 tides, and (c) bottom friction.
Amplification of certain tidal constituents is often observed in estuaries due to resonance caused by similar tidal constituent wavelength and basin length. A resonance in the MSf constituent is in most circumstances through quarter wavelength resonance [27] given as where time period of oscillation (14.76 days for MSf), basin length, average basin depth, and acceleration due to gravity.
The average depth of the LCR estuary is ~17 m and it is ~180 km long. For quarter wavelength resonance to occur the estuary length would have to be ~4200 km compared to the 180 km length of the LCR estuary. Therefore, quarter wavelength resonance is not the reason for MSf amplification.
In estuaries nonlinear interactions between the M2 and S2 tides can generate MS tides (referred to as MSf in this manuscript), whose amplitudes are considerably larger than the theoretical amplitude of ~0 computed previously. Pugh [28] provides four harmonic components generated through nonlinear interaction of M2 and S2 as where location specific constant, tidal component amplitude (subscript indicates tidal component name), and lunar and solar angular speeds (radians), respectively, and = time.
The component represented by the speed that is equal to the difference in speed between the lunar and solar tides is the MSf tide and is . The speed of the MSf component is , with or a period of where time period corresponding to , time period corresponding to , and = time period corresponding to . Representing in terms of mean solar day results in and days. This provides a days. This period is the representative of the spring-neap modulations on the tidal amplitudes [26]. Thus, the MSf amplification is a consequence of the nonlinear interaction of M2 with S2 or vice versa or tidal pumping. Utilizing the above relationship, amplitude of 0.08 m is obtained for MSf at St. Helens; this is close to the observed and predicted value of 0.1 m.
Another explanation for this amplification that has held true in other estuaries is the shallowness of the estuary; however, the LCR estuary is relatively deep compared to others [29], where this has been true. So, excessive shallowness is not a factor and the MSf amplification is attributed to tidal pumping.
9. Conclusions
A numerical model was validated to the tidal hydrodynamics in the Lower Columbia River estuary.
The AdH model-simulated water surface elevations at all observation locations throughout the system were comparable to the observed. The average error metrics values were 0.17 (RMSE), 0.90 (Nash-Sutcliff), 0.97 (Willmott), and 0.95 (Correlation coefficient). The average error in harmonic constituents was 0.003 m (M2), 0.016 m (K1), 0.004 m (S2), 0.009 m (O1), 0.011 m (M4), and 0.008 m (MSf). The largest phase error for all the dominant harmonic constituents (M2, K1, S2, and O1) was ~10°.
The numerical model accurately predicted the amplification of the M4 overtidal component due to the transfer of energy from lower tidal components. The model also accurately predicted the amplification of the MSf tidal component, and analysis of the generation of MSf indicated that the amplification occurs due to the nonlinear interaction between the M2 and S2 components.
The correct representation of the main tidal components, M4 overtide (dependent upon transfer of energy due to friction) and the MSf tide (dependent upon representation of nonlinear processes) indicates that the model is ready for utilization towards the study of tidal propagation in the LCR estuary.
A limitation of the 2D representation of the LCR estuary in the lower LCR was also realized. The 2D representation, probably, misrepresents the energy distribution between the M2 tide and the higher components because of the homogeneous representation of the water column resulting in a lesser transfer of energy from M2. In the absence of conclusive evidence the hypothesis that the 2D representation misrepresents the energy distribution/transfer requires further investigation.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The work described and results presented in this paper were obtained through research sponsored by the US Army Corps of Engineers Portland District and the US Army Corps of Engineers System Wide Water Resources Program (SWWRP). Permission to publish this paper was granted by the Chief of Engineers, US Army Corps of Engineers.