Hydrodynamic Modeling of the St. Lawrence Fluvial Estuary. I: Model Setup, Calibration, and Validation

Matte, Pascal; Secretan, Yves; Morin, Jean

doi:10.1061/(ASCE)WW.1943-5460.0000397

Open access

Technical Papers

Mar 27, 2017

Hydrodynamic Modeling of the St. Lawrence Fluvial Estuary. I: Model Setup, Calibration, and Validation

This article has been corrected.

VIEW CORRECTION

Authors: Pascal Matte https://orcid.org/0000-0003-0968-507X [email protected];[email protected], Yves Secretan [email protected], and Jean Morin [email protected]Author Affiliations

Publication: Journal of Waterway, Port, Coastal, and Ocean Engineering

Volume 143, Issue 5

https://doi.org/10.1061/(ASCE)WW.1943-5460.0000397

PDF

Abstract

In this study, a high-resolution, two-dimensional (2D), time-dependent hydrodynamic model of the St. Lawrence fluvial estuary was developed with the objective of documenting the tidal hydrodynamics of this complex yet poorly understood region. The hydrodynamic model solves the shallow-water equations over a finite-element–discretized domain, with an average spatial resolution of 50 m, and includes a drying–wetting component for the treatment of shallow intertidal areas. The numerical terrain model is composed of high-density topographic data and detailed friction fields associated with bottom substrate and macrophytes. Calibration and validation were carried out using recently acquired data for water level and velocity. Results show very good accuracy in water levels, with prediction skills higher than 0.99 at all stations (where a skill of 1 means perfect agreement between model and observations in terms of their relative average error) and root-mean-square errors (RMSEs) less than 5% of local tidal ranges downstream; at upstream stations where tidal ranges are significantly reduced, RMSEs were lower than 6 cm. Discharges were reproduced with similarly good accuracy, with errors lower than 6% of the maximum observed discharges at 11 of the 13 surveyed transects; the two remaining sections are subject to larger interpolation and bathymetric uncertainties. In this paper, critical aspects of model development are discussed, including the 2D approximation, temporal and spatial resolution, bathymetric uncertainty, error in the boundary conditions, and calibration under nonstationary conditions. This work is the first part of a two-part investigation serving as a methodological framework for model setup, calibration, and validation in large tidal rivers.

Introduction

The St. Lawrence River is the third largest river in North America, with a drainage basin of ∼1.6 × 10⁶ km² and an average freshwater discharge of 12,200 m³/s at Québec. It connects the Great Lakes with the Atlantic Ocean and is the primary drainage of the Great Lakes basin, one of the most industrialized regions of the world. This river enables diverse economic activity for both Canada and the United States, involving commercial navigation, numerous industries, recreational activities, and tourism. It is simultaneously the drinking-water source and effluent receptor of major cities, and it encompasses great aquatic habitat diversity.

The St. Lawrence fluvial estuary (SLFE) spans 180 river kilometers (rkm) from the exit of Lake Saint-Pierre to the eastern tip of Orleans Island, located at the upper limit of the saline intrusion (Fig. 1). The circulation of the SLFE is characterized by vertically well-mixed freshwater (Simons et al. 2010), and it is driven by strong tidal and river flows. Ocean tides are amplified as they enter the St. Lawrence until they reach their highest level (∼7 m in range) in the upper estuary at Saint-Joseph-de-la-Rive, hereafter referred to as rkm 0. Approximately 66 rkm upstream, tidal ranges still exceed 6 m during the largest spring tides at Saint-François (i.e., the downstream limit of the SLFE). Increases in water levels of more than 1 m/h can be observed at these locations during the rising tide, leading to rapid changes in flow conditions and in the wetted areas. This generates strong current reversals with peak tidal discharges up to five times larger than the daily average in both upstream and downstream directions. The tidal signal is increasingly distorted and damped as it propagates upstream as a result of frictional effects (Godin 1999; Matte et al. 2014a); the ebb tides are lengthened, and the flood tides are steepened and shortened. The limit where the flow becomes unidirectional (i.e., only one slack water) moves between Grondines (rkm 179.5) and Bécancour (rkm 217) as a function of tidal range and river flow. At Trois-Rivières (rkm 231), the fortnightly modulation of mean water levels (MWLs) induced by the neap-spring cycle exceeds in amplitude the semidiurnal tide (LeBlond 1979), whose range is 0.2 m for a mean tide. Most of the short-period tide (i.e., diurnal, semidiurnal, etc.) is damped in Lake Saint-Pierre (rkm 264), but long-period oscillations are still noticeable as far as Montréal (rkm 360). These flow properties exhibit both lateral and longitudinal variations that were confirmed by field measurements (Matte et al. 2014b).

Fig. 1. (Color) Map of the SLFE (Québec) showing locations of the DFO’s tide gauges (dark-blue circles), pressure sensors (red stars), and measurement transects (light-blue lines) [Note: rkm are indicated for each station (see Table 3 for a list of station names); major tributaries are the St-Maurice (SM), Batiscan (BA), Sainte-Anne (SA), Jacques-Cartier (JC), and Chaudière (CH) rivers; the global simulation domain is partitioned into two overlapping domains, the upstream (dotted line) and downstream (dashed line) models]

River forcing in the SLFE comes from the freshwater outflows of Lake Ontario, Ottawa River, and other tributaries along its course. Although the average discharge is 12,200 m³/s at Québec, observed minimum and maximum daily net discharges in the St. Lawrence amounted to 7,000 and 32,700 m³/s over the 1960–2010 period, taking into account the contribution of all tributaries and drainage areas upstream of Québec (Bouchard and Morin 2000). The effects of such variations on MWL and tidal range are severe, particularly in the upper portion of the SLFE (Matte et al. 2014a).

Combined with the strong tidal and fluvial forcing are the effects of variable bathymetry, characterized by deep channels (> 60 m); a section width that varies from less than 1 km to more than 15 km; extensive intertidal areas, river bends, shoals, and islands; and high spatial heterogeneity in substrate composition, distribution of macrophytes, and ice and wind conditions. Together, these characteristics provide a unique environment for the study of tidal hydrodynamics.

Previous numerical studies conducted in the SLFE only provided limited knowledge on the circulation complexities of the river. El-Sabh and Murty (1990) reviewed the early modeling efforts made for the St. Lawrence, from the gulf entrance up to Montréal in some instances. Numerous analytical or tidal models have been developed [e.g., Chassé et al. (1993); Godin (1999); Marche and Partenscky (1974); Matte et al. (2014a); Partenscky and Warmoes (1970); Vincent (1965)], providing valuable insights into the propagation and distortion of tides. Development and application of one-dimensional (1D) numerical models of the St. Lawrence have also been conducted [e.g., Bourgault and Koutitonsky (1999); Cheylus and Ouellet (1971); Godin (1971); Kamphuis (1968); Morse (1990); Prandle (1970, 1971); Prandle and Crookshank (1972, 1974); Robert et al. (1992)]. These models were able to reproduce the main tidal and fluvial characteristics encountered in the system, although only qualitatively in the upstream portion of the SLFE for most of them, partly because of imprecise discharge conditions at the boundaries. Among the models, the ONE-D model (Dailey and Harleman 1972; Morse 1990) is currently run in operational mode, fed by the 30-day outflow forecast from Lake Ontario and the Ottawa River and by the 48-h wind forecast of Environment Canada (Meteorological Service of Canada) at the downstream boundary (Lefaivre et al. 2009, 2016). This operational model meets the need for water-level prediction in the SLFE, but it cannot account for lateral exchanges [e.g., Matte et al. (2014b)], which are fundamentally two-dimensional (2D).

A number of 2D (depth- or laterally averaged) numerical models have been developed for the St. Lawrence [e.g., De Borne de Grandpré et al. (1981); Leclerc et al. (1990); Lévesque (1977); Lévesque et al. (1979); Ouellet and Cerceau (1975); Prandle and Crookshank (1972); Prandle and Crookshank (1974); Tee and Lim (1987)]. Three-dimensional (3D) models have also been presented [e.g., Gagnon (1994); Saucier and Chassé (2000); Saucier et al. (2003); Simons et al. (2010)], notably leading to the production of an atlas of tidal currents in the St. Lawrence Estuary (Saucier et al. 1997, 1999). Most of these models do not include the SLFE upstream of Orleans Island. When they do, they usually suffer from a lack of validation, and their spatial resolution is generally too coarse (≥ 200 m) to account for local variations (both lateral and longitudinal) in topography, friction, and hydrodynamic properties. In parallel, model development for specific applications has been conducted by private or government agencies [e.g., Doyon (2011)], but these models are usually restricted to smaller areas of interest and thereby do not provide a complete description of the system hydrodynamics. The SLFE thus remains incompletely documented.

This paper is Part I of a two-part investigation presenting the development of a high-resolution, 2D, time-dependent hydrodynamic model of the SLFE, with the objective of providing detailed spatial and temporal descriptions of the hydrodynamics in response to tidal and fluvial forcings. The model components, underlying data, and setup properties are described in detail in Part I. The finite-element program H2D2, used to run the simulations, is presented, along with its drying–wetting component, which allows water in intertidal areas to be cyclically stored and evacuated. The hydrodynamic model was thoroughly calibrated and validated using an extensive data set composed of water-level data from 29 tide gauges and cross-sectional water-level and velocity data collected along 13 transects, repeatedly surveyed during semidiurnal tidal cycles (Matte et al. 2014a). Tide gauge and transect data are used in Part I to quantitatively evaluate the agreement between modeled and observed water levels and discharges during the calibration phase, whereas transect data are further exploited in Matte et al. (2017) (hereafter referred to as Part II) to validate the spatial and temporal processes observed in the system.

The 2D (vertically integrated) shallow-water equations remain a valid approximation of the tidal hydrodynamics throughout the SLFE because of its small depth-to-width ratio, vertically well-mixed freshwater, and generally small vertical velocities, as confirmed by field measurements (Matte et al. 2014a). However, the model not only encompasses the SLFE but extends further downstream into the estuarine transition zone (Simons et al. 2010), down to Île aux Coudres (Fig. 1), where salinity and vertical stratification start to appear. The inclusion of this segment is meant to allow the placement of a water-level boundary both aligned with a permanent tide gauge and far enough from Orleans Island, where important flow exchanges occur as the river splits into two branches. Although the 2D barotropic model cannot capture the effects of density currents and stratification in this transition zone, model predictability in the SLFE per se (i.e., upstream of the salinity intrusion limit) remains very accurate. For a given computational cost, making the transition to a 3D model would otherwise result in a loss of horizontal resolution, thus compromising the ability of the model to represent complex topographic features.

This work is the first step toward the development of a comprehensive model of the SLFE ecosystem, thus filling the gap between the 2D upstream model (Montréal–Trois-Rivières reach) (Morin and Champoux 2006) and the 3D ocean model of the estuary and gulf (Smith et al. 2013), which are both operational.

Field Campaigns

Hydrodynamic Measurement Campaign

A field campaign, the most extensive to date, was conducted in the SLFE during the summer of 2009 for the measurement of water levels and velocities (Matte et al. 2014b). The objective was twofold: (1) document the hydrodynamics and improve current knowledge of the system, and (2) obtain recent and detailed data for the calibration and validation of the hydrodynamic model.

A set of 15 pressure sensors was installed from June to October 2009 for water-level measurements, adding to the 14 tide gauges operated by Canada’s Department of Fisheries and Oceans (DFO) already available (Fig. 1). Boat surveys were also conducted over semidiurnal tidal periods along 13 cross sections of 1 to 4 km in width of the SLFE, representative of the longitudinal variability in geomorphological and tidal-fluvial properties. At each transect, the boat repeatedly moved back and forth across the channel following the same transect line, perpendicular to the mean flow direction. Water levels and velocities were measured during each crossing with the use of a mounted real-time kinematic global positioning system (RTK GPS) and an acoustic Doppler current profiler (ADCP), respectively. The measurements were carried out between June 15, 2009, and August 25, 2009, so that most transects were surveyed at different phases on the neap-spring tidal cycle. Daily-averaged net discharges at Québec varied between 11,100 and 14,600 m³/s during that period. Details of the data acquisition and analyses have been provided by Matte et al. (2014b, c).

Data from both the extended tide gauge network and boat measurements were used to calibrate and validate the model presented herein, whereas they are used in Part II to further assess model performance based on its capability to reproduce tidal-fluvial processes in the SLFE.

LIDAR Campaign

A light detection and ranging (LIDAR) campaign was conducted during the summer of 2012 to collect topographic data in shallow areas. The topography of all intertidal zones within a 200-km reach of the St. Lawrence, extending from Islet-sur-Mer (rkm 30) to Trois-Rivières (rkm 231), was measured around low tide. The survey was carried out from June 16, 2012, to July 7, 2012. This time period was chosen after the spring freshet and at a time where aquatic plants are not at their peak in terms of density and coverage. A total of 8 days of measurements was necessary to cover the entire domain. The LIDAR system used during the campaign (Optech ALTM Gemini No. 07sen209; Teledyne Optech, Vaughan, ON, Canada) allowed a scanning frequency of 37 Hz, with an angle of 20°. Flight height was 1,350 m aboveground, yielding a scanning width of 877 m. Line spacing was 614 m, allowing a 30% overlap.

The low tide propagates landward at a speed of approximately 15–30 km/h, depending on channel depth, whereas the speed of the plane is approximately 10 times faster (i.e., 256 km/h). The propagation times and heights of low water also vary longitudinally, both over the neap-spring cycle and as a function of river discharge (Godin 1984; LeBlond 1979). Reliable information on the times of arrival of low water was thus essential to the planning of the LIDAR campaign for the measurement to be synchronized with the low tide. Tidal predictions from the ONE-D model of the St. Lawrence (Lefaivre et al. 2009) were used to determine the times of arrival of low water at different stations in the SLFE, under the discharge and tidal conditions prevailing during the period of the survey. The flight lines and time of the survey were determined according to these predictions. Therefore, the extent of dry areas covered was highly dependent on the prediction accuracy of the model and the concordance between the planned and executed flight lines.

Because of the very energetic tidal environment, measurements taken slightly before or after low tide represent a significant loss in the dry areas surveyed. A vertical tolerance of 20 cm around the low tide was considered in the delimitation of zones of measurement. At stations of moderate to strong tidal range (i.e., between Grondines and Saint-François), this tolerance was respected for time intervals from 20 min before low tide to 10 min after low tide (the rising tide being more abrupt than the falling tide). This gave the operator some latitude in the determination of flight lines and times. The flight lines were delineated in such a way as to cover the entire area from the chart datum elevation in the river to the bank limit inland, based on the approximate limit of a 1,000-yr recurrence flood. The limits were further extended by 500 m inside the major tributaries. During measurements, the flight lines were executed in a predetermined order, from the channel to the bank, to ensure that the most critical areas (intertidal zone, shoals, etc.) were surveyed as close as possible to the low tide.

Model Setup

Simulation Domain and Boundaries

The simulation domain of the SLFE model extends beyond the actual limits of the freshwater estuary, typically defined from the exit of Lake Saint-Pierre to the eastern tip of Orleans Island (Fig. 1). The downstream boundary is, in fact, located in the estuarine transition zone of the St. Lawrence Estuary (Simons et al. 2010). It is positioned along a 15-km-wide cotidal line of constant phase, a few kilometers downstream from the nearest tide gauge, located on the north shore at Saint-Joseph-dela-Rive (rkm 0). This station is the closest permanent tide gauge to the east of Orleans Island, past Saint-François (rkm 66). It was chosen as the downstream boundary, far enough from Orleans Island to allow exchanges between the north and south arms to occur freely. Because of the section width and the absence of a tide station on the south shore of the boundary, water-level distribution across the section was determined by assimilation, using differences between observed and modeled water levels at nearby stations to orient changes in the imposed water levels (discussed further under Assimilation).

The upstream boundary is positioned in alignment with the Port Saint-François tide station (rkm 241), at the exit of Lake Saint-Pierre (i.e., the upstream limit of the SLFE). Measured water levels, rather than estimated tidal discharges, were specified at the boundary to ensure accurate results because a tide of 20–30 cm in range still subsists at the exit of Lake Saint-Pierre. To actually impose the freshwater discharge would require extending the model beyond the head of the tide, to the closest station where continuous discharge measurements are taken (Lasalle station, more than 100 rkm upstream). Although work is currently under way to connect the present model with the upstream 2D operational model (Morin and Champoux 2006), the current focus is on the SLFE reach.

Major tributaries were included in the model to allow water to be cyclically stored and evacuated as a function of the tide: the Saint-Maurice (SM), Batiscan (BA), Sainte-Anne (SA), Jacques-Cartier (JC), and Chaudière (CH) rivers (Fig. 1). The boundaries were positioned at upstream locations removed from tidal influence for the imposition of daily-averaged discharges uninfluenced by the tide. The latter were reconstructed by adding the discharge measured at an upstream station to the estimated lateral inflow, consisting of surface water runoff and groundwater inflow. Because virtually no data are available for groundwater inflow, only surface water runoff was considered and was approximated based on gauged areas. For ungauged areas, the inflow was estimated from the runoff coefficient of an adjoining gauged area. Relations for each tributary to the St. Lawrence were developed by Morse (1990) and adapted by Bouchard and Morin (2000).

The global simulation domain was partitioned into two smaller overlapping domains, each sharing the same grid configuration and data as the global model (Fig. 1). This separation was necessary to dissociate the calibration process from the assimilation process. Hence, calibration in the upstream model was performed on a reduced domain extending from Port Saint-François (rkm 241) to Québec (rkm 106.5), where water-level conditions are defined at the boundaries from tide gauge data. The downstream model was calibrated on a reduced domain extending from Neuville (rkm 138) to Saint-Joseph-dela-Rive (rkm 0), where assimilation of the downstream boundary conditions was carried out to determine the optimal water-level distribution. Once the two submodels were calibrated and assimilated, simulations were performed using the global domain for validation.

Numerical Terrain Model

Development of the numerical terrain model (NTM) was made using MODELEUR software, a geographic information system (GIS) adapted to fluvial hydrodynamics (Secretan and Leclerc 1998; Secretan et al. 2001). The NTM is horizontally positioned in a Universal Transverse Mercator (UTM) projection in the NAD83 geodetic reference system. The vertical datum is CGVD28. The different components composing the NTM are described in the following sections.

Topography

The main explanatory factor of the hydrodynamics is the topography, which includes channel bathymetry, floodplains, and engineering structures. Vertical accuracy of a few centimeters can be achieved in shallow-water coastal environments with today’s high-resolution multibeam echosoundings, coupled with high-accuracy kinematic GPS positioning (Ernstsen et al. 2006). In the St. Lawrence, new bathymetric soundings are routinely made in the navigational channel, but several regions outside the channel were only surveyed decades ago, and shallow regions are generally not covered. Not only measurement precision may vary in time as the technology improves, but morphological changes may have occurred since the surveys were conducted. Significant offsets may also appear as a result of inaccuracies in the local chart datum, typically derived from tide measurements at fixed stations not necessarily representative of the spatial variations in tidal amplitude, particularly as the estuary width increases. Therefore, quantifying these errors is difficult; they can be as high as 1 m, especially if century-old data are included (Burningham and French 2011). Data density and associated errors are therefore a function of the year of acquisition and location. Approximately 42 million bathymetric data points were obtained from the Canadian Hydrographic Service (CHS). They were reduced from chart datum to CGVD28 using a kriging grid based on known conversions at the tide gauges.

To complete the topography in shallow areas, data from the LIDAR campaign were used. Water lines were identified from the data and used to separate points on the water from those on the ground. Data were validated with supplementary LIDAR and bathymetric data sets in overlapping regions. The accuracy of LIDAR data is of the order of 15 cm on average but can vary as a function of land cover [e.g., Schmid et al. (2011)]. In all, 420 million LIDAR data points were integrated to the NTM.

As for the tributaries, the geometry of the river shores was extracted using geospatial data from Natural Resources Canada’s GeoBase website (http://www.geobase.ca/). Because of the lack of data, their bathymetry was represented by a regularly shaped trapezoidal channel of constant slope and variable depth and width. Channel depth was determined by calculating the depth needed to discharge the average river flow at an approximate mean velocity of 1 m/s, given the local width of the river.

Engineering structures, such as bridge pillars, piers, marinas, and ports, were also included in the NTM. A free-slip condition was used along vertical walls, but most lateral boundaries are mobile and controlled by the drying–wetting component. The topography as described by the SLFE model is presented in Fig. 2.

Fig. 2. (Color) Topography of the SLFE as projected on the finite-element mesh, with zoom-ins at Grondines and at the junction of Orleans Island

Friction

Information on bottom substrate and macrophytes was added to the NTM for friction description. Substratum composition was defined for homogeneous regions based on 6,400 substratum observation points obtained from the CHS. They were converted into Manning coefficients following Morin et al. (2000a). Friction resulting from macrophytes was included as an additional layer in regions where aquatic plants are observed during summer. Friction coefficients were adjusted during calibration, within the range of values provided by Morin et al. (2000b). Finally, a constant friction coefficient was used inside the tributaries. The distribution map of Manning coefficients used in the SLFE during the calibration period is shown in Fig. 3. They range from 0.0125 to 0.0625, the highest values accounting for very rough floodplain elements (e.g., trees). Surface friction by wind or ice was set to 0.

Fig. 3. (Color) Manning coefficients in the SLFE, based on substrate composition and macrophyte distribution, with a zoom-in on Gentilly shoal (Bécancour)

Finite-Element Mesh

Both topography and friction data were assembled onto a computing grid, in this case, a 2D finite-element mesh composed of triangular P₁–P₁isoP₂ elements of continuity C⁰ (Heniche et al. 2000). The mesh for the global model forms a triangular irregular network composed of 1,347,515 nodes and 662,934 elements. The calibration and assimilation domains include 585,798 and 849,480 nodes, and 286,300 and 419,318 elements, respectively. The mesh was built following the river morphology in such a way as to reduce errors in regions of strong variability and to represent the terrain adequately. Hence, average grid resolution is 50 m, with refinements down to ∼1 m near engineering structures (e.g., bridge pillars, piers) and over regions of complex topography (e.g., steep bathymetry, narrow channels), as evidenced by the bathymetric data. These refinements are localized (i.e., in the near surroundings of strong topographic variations) and are meant to capture the effects of these topographic elements on large-scale processes, such as tidal propagation and lateral exchanges between the shallow intertidal areas and the main channel. This also allows for calibrating and validating the model with a similar level of detail as provided by the available hydrodynamic data (Matte et al. 2014b, 2017). An even higher spatial resolution would be necessary, however, with an accordingly increasing computational cost, to capture finer structures of the flow propagating away from their originating point; this was not the objective here. The finite-element mesh is shown in Fig. 4.

Hydrodynamic Model

Hydrodynamic simulations were performed using H2D2, which allows for robust, distributed, and shared-memory computing for large systems and time-dependent problems. It solves the 2D depth-integrated Navier–Stokes equations over a finite-element discretized domain, with special treatment of drying–wetting areas. The shallow-water model, as implemented in H2D2, is based on the assumptions of incompressibility, hydrostatic pressure, and stable riverbed [for a derivation of the shallow-water equations, see, e.g., Dronkers (1964) and Bois (2000)]. The weak variational formulation and finite-element discretization have been detailed by Heniche et al. (2000) and Dhatt et al. (2005). A brief summary of the model equations is provided in the Appendix. Further details on the various modules of H2D2 can be found online (http://www.gre-ehn.ete.inrs.ca/H2D2). Descriptions of the drying–wetting model, temporal scheme, and application to the SLFE are presented in the following sections.

Drying–Wetting Model

A Eulerian approach is used for the treatment of drying–wetting areas (Heniche et al. 2000), in which the water level can plunge under the bed level and generate both positive (i.e., wet) and negative (i.e., dry) water depths. The effective depth in the dry area is, however, limited to a minimum depth H_min (Fig. 5), allowing only a thin layer of water to subsist that mimics a groundwater flow. With the goal being to freeze the flow in the dry area, for mass to be conserved, the Manning coefficient n is drastically increased (n ∼ 0.5), whereas in the wet area, it is set in accordance with local flow properties. Increased viscosity is also imposed in the dry area to reduce to a minimum the contribution of velocities to momentum conservation. Moreover, a Darcy viscosity (via the hydraulic conductivity δ) is added to further smooth the free surface in dry areas (Appendix). In general, to ensure a smoother transition, changes in the parameters between the wet and dry areas are made over a certain distance, delimited by the position of H_min and H_threshold (Fig. 5). Within this transitory depth zone, all the parameters (e.g., Manning’s n) vary linearly from their wet value to their dry value. Furthermore, this transition follows a hysteresis loop, identified by the dashed lines in Fig. 5, so that the passage of a node from the wet to the dry state is made differently than inversely. This feature is meant to reduce the rigidity of the system in allowing faster convergences.

Resolution

The discretized shallow-water equations are integrated in time using an implicit Euler time scheme. The resulting nonlinear algebraic system is solved by Newton–Raphson iterations (Dhatt et al. 2005). A spin-up of approximately one tidal cycle (12h25) is needed to dissipate the effects of initial conditions and to propagate the tidal wave from one end of the domain to the other [average propagation time from Saint-Joseph-dela-Rive (rkm 0) to Port Saint-François (rkm 241) is ∼8 h]. This rather short spin-up time is consistent with the type of boundary conditions used. Because the model is forced by water levels at both upstream and downstream boundaries, the low-frequency (fortnightly) tide, whose amplitude exceeds that of the semidiurnal tide upstream, is directly prescribed at the boundaries so that the effects of initialization are rapidly dissipated. Extending the model beyond the head of the tide and specifying discharge at the upstream boundary would, however, require a much longer spin-up time (∼1–2 weeks) for the fortnightly tide to reach equilibrium.

Implicit temporal schemes do not constrain the time-step size. Simulations are thus run at a 5-min nominal time step, with a time-step halving algorithm in case of nonconvergence. This time-step size is fine enough to accurately represent the tidal dynamics of the SLFE but also acts as a temporal filter for small-scale processes by smoothing any motion whose frequency is higher than that prescribed by the time step. Numerical experiments using explicit, semi-implicit, and implicit time integration schemes, with time steps varying between 1 and 300 s, were made (not presented) to assess the impact of the chosen spatial and temporal resolutions. Sensitivity analyses were also conducted to adjust turbulent and numerical viscosities through the mixing length and Peclet number, respectively (Appendix), within their limit of stability and/or applicability. Values between 0.01 and 1 for the mixing length coefficient and between 1 and 5 for the Peclet number were tested. Overall, velocity fields were minimally sensitive to these changes when compared with the reference simulation (5-min time step, mixing length coefficient of 1, Peclet number equal to 1). This is consistent with the fact that refinements in the mesh are very localized and do not allow small-scale structures to develop and propagate, even at low viscosity. Only larger-scale structures are observed, which remain unaffected by a reduction of the time step.

Using a relatively large time step is also consistent with the hydrostatic assumption of the 2D shallow-water model, by low-pass filtering of small-scale, and hence nonhydrostatic, flow features. Smaller intratidal transverse flow structures caused by channel curvature (e.g., at Portneuf; rkm 163.5), rapidly varying bottom slopes (e.g., at Deschambault; rkm 168), and channel bifurcation (e.g., at the junction of Orleans Island; see Fig. 1) can lead to stronger nonhydrostatic effects. However, the effect of nonhydrostatic pressure on the large-scale hydrodynamics (i.e., tidal variations of water levels and depth-averaged velocities) is negligible, as demonstrated by Wang et al. (2009).

Calibration

Calibration was conducted for two 15-day simulation periods starting on 14 June 2009 00:00:00 EDT and 19 August 2009 00:00:00 EDT, respectively, corresponding to the periods of transect surveys. The run time for the two simulations was approximately 40 and 31 wall-clock days, respectively, on two AMD (Sunnyvale, California) Opteron 2384 (four cores, 2.7 GHz) processors with 16 GB of memory.

Calibration consisted of model adjustments to reduce errors associated with (1) topography, (2) substrate, and (3) friction. In rivers, topography is the primary factor controlling the hydrodynamics. Bathymetric uncertainties are often neglected in the calibration and validation of hydrodynamic models, but their effect on the computed velocity field can be significant (Cea and French 2012), especially over intertidal flats (Wang et al. 2009). Large uncertainty may arise from surveying errors, georeferencing, digitization, datum conversions, morphological changes, data density and coverage, model discretization, and resolution, among other factors. Failure to adequately capture topographic features and gradients will most likely result in erroneous velocity patterns. Furthermore, the tidal phase is also sensitive to bathymetric errors because wave celerity is a function of water depth [Eq. (7)]. Ensuring that modeled tidal signals are in phase with observations is an indirect way to validate the bathymetry. In the present case, the tidal phase was constrained by having two water-level boundary conditions at the upstream and downstream limit of the model. Thus, topographic errors were mainly detected and corrected by looking for the presence of incoherent velocity patterns or large errors in water levels and by making an educated guess on where the error was likely to be maximal, based on the error sources mentioned previously and on intercomparisons between data sets of different sources.

In a second step, substrate composition was locally controlled to ensure that mean surface slopes and tidal ranges were adequately reproduced by the model. Finally, further adjustments in Manning friction coefficients were made, especially in shallow regions where macrophytes are present, until errors in water levels were minimized. Increases in bottom friction act on the tides by reducing their amplitude and increasing the MWL upstream of the modification. Because the model is controlled in water levels at both upstream and downstream boundaries, tidal amplitudes are also slightly increased downstream of the modification when friction is increased. Experimentation showed that a performance criterion based on these two metrics (MWL and tidal range) can be more informative in tidal rivers than a calibration solely based on tidal amplitudes and phases of major astronomical and shallow-water constituents. In fact, river tides are nonstationary, meaning that the tidal constituent amplitudes and phases vary as a function of time, mostly with river discharge. Harmonic analysis of nonstationary records can only extract the average frequency content of the tide and provides no information about its temporal variability. As a result, a large part of the original signal variance gets lost in the residual even though significant (unresolved tidal or nontidal) energy remains. With short records, even nonstationary harmonic analysis methods [e.g., Matte et al. (2013)] or wavelet analysis techniques [e.g., Buschman et al. (2009); Jay and Kukulka (2003); Sassi et al. (2011)] present limited descriptive capabilities.

Model calibration was therefore performed based on raw water levels and filtered MWLs and tidal ranges, allowing for easier isolation of the separate effects of friction on the low-frequency fluctuations of mean surface slope and tidal range. Such an approach is also believed to be more versatile than harmonic-based validation alone because it can be applied to a variety of nonstationary contexts other than tidal.

MWLs were calculated using a 3-day moving average, similar to the tidal eliminator introduced by Godin (1972). Tidal ranges were extracted from the high-pass-filtered data by calculating the difference between higher high water and lower low water using a 27-h moving window, similar to the tidal range filter used by Kukulka and Jay (2003).

To assess model performance, two criteria were used, namely, root-mean-square error (RMSE) and skill (Willmott et al. 1985), based on quantitative comparisons between simulated and observed values; they are respectively given by

RMSE = \sqrt{\frac{1}{n} \sum_{n} {(X_{sim} - X_{obs})}^{2}}

(1)

and

Skill = 1 - \frac{\sum_{n} {(X_{sim} - X_{obs})}^{2}}{\sum_{n} {(| X_{sim} - \bar{X_{obs}} | + | X_{obs} - \bar{X_{obs}} |)}^{2}}

(2)

where X = variable being compared; and

\bar{X}

= time average of n values. RMSE is a measure of the average error between the model and observations in the units of the quantities assessed. Skill is a unitless measure of the relative average error between the model and observations, where a skill of 1 means perfect agreement, and a skill of 0 means that the model is equivalent to the mean of the observations.

Assimilation

Imposed water levels along the downstream boundary were derived from observations at Saint-Joseph-dela-Rive (rkm 0) rather than tidal harmonic constituents to allow nontidal (residual) variations to be propagated in the system. Because a permanent station exists only on the north side of the boundary, water levels across the 15-km-wide section were optimized to determine the distribution that would best reproduce observations at upstream tide gauges, located on both sides of the river. During assimilation, time shifts were iteratively applied on the imposed water levels until simulated and observed signals were in phase at the stations. Lateral gradients (both positive and negative) in MWL and tidal range were also tested along the boundary to assess model sensitivity to the imposed conditions. The simulated variations were generally not sufficient to explain the differences between observed and simulated water levels at the stations, so the decision was made to specify the same (time-shifted) water-level time series on every node of the downstream boundary, with no lateral gradient. Assimilation of the boundary condition was thus limited to (1) determining the optimal time shift needed for the tide to propagate in phase with observations at the stations and (2) defining the corrections in MWL (i.e., vertical shift) and tidal range (i.e., amplitude factor) that would best reproduce the tidal signal at Saint-Joseph-dela-Rive. More precisely, the imposed water levels h were lagged temporally and updated as follows: h_i(t) = a + b*h_i–1(t – τ), where i = iteration number; t = time; τ = computed time lag; and a and b = constants iteratively determined from model ouputs at Saint-Joseph-dela-Rive.

A slightly different assimilation strategy was used by Lefaivre et al. (2016) to optimize the downstream boundary conditions of a 1D model, partitioned into three branches in the St. Lawrence estuarine transition zone. Differences between the two strategies mostly arise from the fact that the 2D model comprises the complete bathymetry, which is particularly complex in the area, whereas the 1D model needs to compensate for these bathymetric variations, partly through its boundary conditions.

Validation

Several authors have described the process of model evaluation, calibration, and validation [e.g., Dias and Lopes (2006); Espino et al. (2007); Hsu et al. (1999); Jung et al. (2012); Liu et al. (2009); Umgiesser et al. (2004); Warner et al. (2005); Willmott (1982); Willmott et al. (1985)]. Classic calibration and validation may be considered as two independent exercises. In practice, however, iterations are generally made between the two processes so that validation most often becomes an integral part of the calibration process. Validation is defined generically as the process of determining the degree to which a model is an accurate representation of the real world. In general, there can be four types of validation (Arhonditsis and Brett 2004; Power 1993): (1) replicative validation, where the agreement between predicted values and observational data from the real system is quantified during the calibration phase, based on the notions of goodness of fit or distributional similarity; (2) predictive validation, where model performance is evaluated against independent data sets, or once the data are acquired from the real system after a forecast has been issued; (3) structural validation, where the ability of the model to reproduce the real system behavior is assessed, with regard to its operational characteristics, spatial and temporal patterns, and relative magnitudes of derived quantities [the expression scientific validation is sometimes used to refer to this type of validation (Biondi et al. 2012); herein and in Part II, the term process-based validation is adopted]; and (4) model transferability, which is a measure of the performance of a specific model structure in different systems.

In this first part, a replicative validation is performed, focusing on a quantification of model accuracy in water levels and discharges (detailed in the following discussion). In Part II, a process-based validation is conducted, aiming at the reproduction of the dynamic patterns observed in the field data.

Two simulation periods were tested here (June and August 2009), for which detailed measurements are available. Calibration was mainly based on a comparison of modeled and observed water levels at the tide gauges, using tidal ranges and MWLs as metrics for these changes (as previously detailed under Calibration). Most of the calibration was performed based on the simulation in June. The result accuracy was then verified for the simulation in August, and minor adjustments were made.

During the calibration phase, the domain was also divided into two subdomains to separate the calibration from the assimilation process; the calibration and assimilation separately performed on the subdomains were then validated on the global model for the two simulation periods, based on a comparison of observed and modeled water levels, MWLs, and tidal ranges at the tide gauges.

To validate the frequency content of the tide, harmonic properties (i.e., tidal amplitudes and phases of major tidal and shallow-water constituents) were computed from observed and modeled tidal signals by performing classical harmonic analyses (Pawlowicz et al. 2002) at all stations of the global model. For each tidal constituent, RMSE and skill were calculated by summing results at all 29 stations, using Eqs. (1) and (2).

Finally, the simulated discharges were controlled at each of the 13 measurement transects (Fig. 1) to ensure that the right amount of flow was in the system and that the overall dynamics were accurately reproduced.

Results and Discussion

Statistics on Water Levels

Results from the calibrated model are presented separately for the upstream and downstream segments (Tables 1 and 2, respectively) and for the global model (Table 3) for the two simulation periods of June 14–29, 2009, and August 19–September 3, 2009. Because friction acts on both MWLs and tidal ranges, statistics on these two quantities were also calculated. They were used during calibration to orient the local modifications toward an increase or a decrease in friction. In the upstream model (141-km-long reach), RMSEs between observed and simulated water levels were lower than 10 cm at all 16 stations (Table 1). Modeled water levels were found to be in almost perfect agreement with observations at the stations, with prediction skills varying between 0.993 and 1.000. The highest RMSE (0.094 m) was obtained at Grondines and was mostly dominated by errors in MWL. These are likely related to topographic errors in the intertidal flats between Deschambault (rkm 168) and Grondines (rkm 179.5), combined with an overestimated friction.

Table 1. Calibration Results for the Upstream Model: RMSE and Skill for Modeled Water Levels (Raw), MWLs, and Tidal Range at the Stations (Fig. 1) for Each 15-Day Simulation Period

rkm	Station	June 14–29, 2009				August 19–September 3, 2009
		RMSE (m)			Skill	RMSE (m)			Skill
		Raw	MWL	Range	Skill	Raw	MWL	Range	Skill
115	Québec Bridge	0.053	0.028	0.057	1.000	0.051	0.022	0.085	1.000
124	Saint-Nicolas	0.027	0.005	0.020	1.000	0.030	0.005	0.015	1.000
138	Neuville	0.043	0.015	0.083	1.000	0.042	0.013	0.061	1.000
146	Sainte-Croix-Est	0.046	0.013	0.052	1.000	0.057	0.018	0.082	0.999
157	Cap-Santé	0.046	0.019	0.077	1.000	0.047	0.021	0.051	1.000
161	Pointe-Platon	0.042	0.014	0.048	1.000	0.043	0.020	0.022	1.000
163.5	Portneuf	0.061	0.050	0.033	0.999	0.076	0.066	0.027	0.999
168	Deschambault	0.062	0.025	0.124	0.999	0.067	0.030	0.193	0.999
179.5	Leclercville	—	—	0.021	—	—	—	0.020	—
179.5	Grondines	0.094	0.090	0.021	0.995	0.087	0.082	0.017	0.996
186	Cap-à-la-Roche	0.028	0.006	0.032	0.999	0.037	0.022	0.027	0.999
199	Batiscan	0.045	0.039	0.053	0.996	0.041	0.033	0.051	0.997
213	Champlain	0.054	0.047	0.068	0.990	0.045	0.032	0.063	0.995
217	Bécancour	0.016	0.006	0.014	0.999	0.050	0.047	0.022	0.993
231	Trois-Rivières	0.028	0.026	0.007	0.989	0.014	0.006	0.008	0.999
241	Port Saint-François	0.008	0.007	0.007	0.999	0.007	0.006	0.006	1.000

Table 2. Calibration and Assimilation Results for the Downstream Model: RMSE and Skill for Modeled Water Levels (Raw), MWLs, and Tidal Range at the Stations (Fig. 1) for Each 15-Day Simulation Period

rkm	Station	June 14–29, 2009				August 19–September 3, 2009
		RMSE (m)			Skill	RMSE (m)			Skill
		Raw	MWL	Range	Skill	Raw	MWL	Range	Skill
0	Saint-Joseph-de-la-Rive	0.095	0.014	0.066	0.999	0.091	0.020	0.069	0.999
30	Islet-sur-Mer	0.130	0.026	0.292	0.998	0.143	0.047	0.304	0.998
38	Rocher Neptune	0.106	0.069	0.198	0.999	—	—	—	—
45	Ile-aux-Grues South	0.114	0.059	0.213	0.999	0.143	0.102	0.196	0.998
46	Ile-aux-Grues North	0.160	0.119	0.297	0.997	0.133	0.092	0.274	0.998
54	Banc du Cap Brûlé	0.102	0.043	0.192	0.999	0.113	0.042	0.189	0.999
66	Saint-François	0.120	0.066	0.257	0.999	0.120	0.052	0.257	0.998
78	Saint-Jean	0.155	0.131	0.233	0.998	0.106	0.048	0.224	0.999
97	Beauport	0.197	0.173	0.047	0.996	0.157	0.122	0.068	0.997
100	Lauzon	0.163	0.128	0.039	0.997	0.124	0.071	0.047	0.998
104	Saint-Charles River	0.157	0.120	0.029	0.997	0.126	0.070	0.039	0.998
106.5	Lévis	0.148	0.106	0.040	0.997	0.119	0.061	0.043	0.998
106.5	Québec	0.151	0.111	0.037	0.997	0.119	0.065	0.036	0.998
115	Québec Bridge	0.151	0.103	0.127	0.997	0.117	0.061	0.168	0.998
124	Saint-Nicolas	0.082	0.056	0.086	0.999	0.067	0.041	0.098	0.999
138	Neuville	0.005	0.002	0.006	1.000	0.006	0.003	0.005	1.000

Table 3. Validation of the Global Model: RMSE and Skill for Modeled Water Levels, MWLs, and Tidal Range at All Stations (Fig. 1) for Each 15-Day Simulation Period

rkm	Station	June 14–29, 2009				August 19–September 3, 2009
		RMSE (m)			Skill	RMSE (m)			Skill
		Raw	MWL	Range	Skill	Raw	MWL	Range	Skill
0	Saint-Joseph-de-la-Rive	0.095	0.013	0.064	0.999	0.091	0.021	0.066	0.999
30	Islet-sur-Mer	0.130	0.031	0.303	0.998	0.143	0.051	0.319	0.998
38	Rocher Neptune	0.103	0.067	0.198	0.999	—	—	—	—
45	Ile-aux-Grues South	0.110	0.057	0.212	0.999	0.145	0.107	0.193	0.998
46	Ile-aux-Grues North	0.154	0.115	0.293	0.998	0.135	0.098	0.271	0.998
54	Banc du Cap Brûlé	0.100	0.047	0.195	0.999	0.113	0.049	0.192	0.999
66	Saint-François	0.124	0.078	0.243	0.999	0.121	0.061	0.249	0.998
78	Saint-Jean	0.166	0.147	0.194	0.997	0.108	0.058	0.202	0.999
97	Beauport	0.221	0.198	0.064	0.995	0.160	0.120	0.058	0.997
100	Lauzon	0.185	0.155	0.074	0.996	0.133	0.076	0.084	0.998
104	Saint-Charles River	0.182	0.147	0.102	0.996	0.136	0.075	0.099	0.998
106.5	Lévis	0.175	0.137	0.102	0.997	0.133	0.069	0.086	0.998
106.5	Québec	0.177	0.141	0.086	0.997	0.132	0.072	0.088	0.998
115	Québec Bridge	0.177	0.140	0.040	0.996	0.131	0.069	0.079	0.998
124	Saint-Nicolas	0.138	0.112	0.087	0.998	0.104	0.062	0.078	0.999
138	Neuville	0.135	0.107	0.135	0.997	0.107	0.057	0.116	0.998
146	Sainte-Croix-Est	0.123	0.096	0.062	0.997	0.111	0.061	0.103	0.998
157	Cap-Santé	0.125	0.102	0.109	0.997	0.107	0.066	0.095	0.998
161	Pointe-Platon	0.118	0.094	0.083	0.997	0.109	0.066	0.085	0.998
163.5	Portneuf	0.089	0.062	0.056	0.999	0.121	0.086	0.076	0.997
168	Deschambault	0.106	0.069	0.139	0.997	0.126	0.074	0.222	0.996
179.5	Leclercville	—	—	0.029	—	—	—	0.042	—
179.5	Grondines	0.070	0.054	0.029	0.997	0.111	0.087	0.044	0.993
186	Cap-à-la-Roche	0.070	0.055	0.036	0.996	0.080	0.050	0.061	0.996
199	Batiscan	0.040	0.027	0.037	0.996	0.060	0.044	0.063	0.994
213	Champlain	0.046	0.035	0.074	0.992	0.058	0.037	0.056	0.992
217	Bécancour	0.027	0.010	0.042	0.997	0.056	0.048	0.036	0.992
231	Trois-Rivières	0.025	0.022	0.019	0.992	0.015	0.005	0.008	0.999
241	Port Saint-François	0.008	0.007	0.007	0.999	0.007	0.006	0.005	1.000

Results for the downstream model are presented in Table 2. Errors in water levels were larger than those in the upstream segment, with the maximum RMSE reaching nearly 20 cm at Beauport. In an attempt to reduce these errors by assimilation, lateral variations in MWLs and tidal ranges were imposed along the downstream boundary, with little success; variations in modeled water levels were not sufficient to explain the observed differences at the stations. The error is likely attributable to an underestimated friction in the intertidal zones of Orleans Island, particularly in the north arm, as a result of the presence of aquatic plants during summer. This error has a repercussion on upstream water levels, and efforts should be made in the future to reduce it below 10 cm, or even 5 cm, to meet similar quality standards as the upstream operational model of the St. Lawrence (Morin and Champoux 2006). Nonetheless, prediction skills were found to be very high (> 0.997) at all stations, indicating an overall good performance of the model.

Table 3 presents statistics on water levels obtained with the global model. The errors made in the downstream segment propagated upstream, increasing the gaps between observed and modeled water levels by several centimeters at upstream stations compared with the results shown in Table 1. Errors in water levels were attributed mostly to underestimated tidal ranges in the first 78 rkm. They were dominated by errors in MWL at Beauport (rkm 97) and relatively balanced between errors in MWL and tidal ranges upstream, depending on the station and simulation period considered. Despite these discrepancies, prediction skills were still found to be higher than 0.992 at all stations. These very high skills cn be explained by the large tidal ranges characterizing the region, which often exceed 5 m. In fact, for stations downstream of Batiscan (rkm 199), the RMSEs correspond to less than 5% of the local tidal ranges. At upstream stations, the ratio of the error to tidal range increases because of the rapidly decreasing tidal ranges, but RMSEs were lower than 6 cm. This confirms that longitudinal variations in friction were overall well captured by the model.

Harmonic Properties

Statistics from the observed and modeled amplitudes and phases of major tidal and shallow-water constituents were computed for each 15-day simulation period. Results for the 9 dominant constituents, among the 17 resolved ones, are presented in Table 4. Classical harmonic analysis implies that tides are stationary. In reality, daily discharges at Québec varied by 1,000 m³/s (12,400 m³/s average) and 1,200 m³/s (11,300 m³/s average) during the two simulation periods, respectively. Therefore, only the average frequency content of the tides could be retrieved from the analyses. For longer signals, however, nonstationary harmonic analyses could be performed to extract the time-varying amplitudes and phases of the tides [e.g., Matte et al. (2013, 2014a)].

Table 4. Statistics from Observed and Modeled Amplitudes and Phases of Principal Tidal Constituents at All Stations of the Global Model for Each 15-Day Simulation Period

Constituent	June 14–29, 2009			August 19–September 3, 2009
	Amplitude		Phase	Amplitude		Phase
	RMSE (m)	Skill	RMSE (degrees)	RMSE (m)	Skill	RMSE (degrees)
MSf	0.020	0.981	5.5	0.034	0.966	4.1
O₁	0.009	0.997	2.9	0.011	0.992	2.6
K₁	0.019	0.988	2.7	0.008	0.997	4.4
M₂	0.066	0.998	1.5	0.059	0.998	1.9
S₂	0.010	0.997	7.4	0.035	0.993	4.3
M₃	0.010	0.975	2.4	0.005	0.980	9.2
M₄	0.015	0.996	7.2	0.019	0.996	6.7
M₆	0.009	0.986	17.8	0.016	0.982	15.5
M₈	0.006	0.991	15.8	0.011	0.961	23.4

The results shown in Table 4 highlight the capacity of the model to represent the tidal frequency content of the observed signals with very high accuracy, with skills greater than 0.961 for all tidal constituents. MWLs were well accounted for by the model, based on statistics from the low-frequency MSf component. The larger RMSE was obtained for the M₂ (dominant) constituent, whose amplitude exceeded 2 m at the most downstream stations. The relative error for M₂ was therefore the lowest, presenting the highest skills among analyzed constituents. Phases were also well reproduced by the model, indicating a good synchronism of the tides with observations. Furthermore, tidal asymmetry, which can be assessed through the relative importance of M₂ and M₄ components, is expected to be well accounted for by the model. In fact, the skills associated with these two components were among the highest (≥0.996). Phase errors degraded with constituents of higher frequency, namely, M₆ and M₈. However, the lowest errors for these components were observed at upstream locations, where they were the most significant, because they are generated by energy transfers from M₂ through nonlinear frictional interactions.

Statistics on Discharges

Validation of the tidal discharges was performed at each measurement transect of Fig. 1. Statistics are presented in Table 5. The RMSEs generally increased with the strength of the tidal flow, but they remained below 6% of the maximum discharge (i.e., peak ebb discharge) on average. Only at the sections of Beauport (rkm 97) and Lauzon (rkm 100) was the relative RMSE larger than 10%. Because the discharges downstream, in the north and south arms of Orleans Island, were quite accurate (relative RMSEs of 6.3% and 5.3% at Château-Richer and Saint-Jean, respectively), these larger RMSEs at the island junction can most likely be attributed to the larger interpolation errors of the transect data between each crossing (Matte et al. 2014b) in combination with possible local bathymetric errors responsible for larger discrepancies in modeled water levels at these stations (cf. Table 3).

Table 5. Statistics at the Validation Transects (Fig. 1) for the Global Model: RMSE and Skill for Modeled Discharges

rkm	Transect	Date	Length of record (h)	RMSE [m³/s (%)]	Skill
78	Saint-Jean	24 Aug 09	7.6	3,221 (5.3)	0.997
79	Château-Richer	25 Aug 09	8.7	639 (6.3)	0.996
97	Beauport	24 Jun 09	10.7	946 (12.0)	0.985
100	Lauzon	24 Jun 09	11.4	6,332 (11.1)	0.993
106.5	Québec	15 Jun 09	11.3	2,934 (6.1)	0.997
124	Saint-Nicolas	21 Aug 09	8.8	1,217 (2.7)	1.000
138	Neuville	25 Jun 09	9.1	1,166 (3.1)	0.999
163.5	Portneuf	26 Jun 09	6.7	1,164 (5.2)	0.998
168	Deschambault	20 Aug 09	10.1	719 (3.6)	0.997
179.5	Grondines	19 Jun 09	8.9	901 (5.5)	0.988
199	Batiscan	23 Jun 09	7.6	856 (5.6)	0.994
217	Bécancour	18 Jun 09	5.0	650 (4.6)	0.826
231	Trois-Rivières	18 Jun 09	0.8	680 (5.4)	0.124

Note: Values in parentheses are errors relative to the maximum observed ebb discharge.

Prediction skills, overall, were found to be very high considering all transects, with the exception of Bécancour (rkm 217) and Trois-Rivières (rkm 231), where the length of the discharge records was less than half the tidal period. At Trois-Rivières, the low skill can also be explained by the fact that tides are considerably weaker there than downstream, leading to only a very slight increase in the observed discharge during the measurement period. The model, however, reproduced the average flow with an accuracy of approximately 5% near the upstream boundary, meaning that the right amount of water entered into the system through the upstream (water-level) boundary. Elsewhere, the very high skills indicate that the variability in tidal discharges and the overall dynamics were accurately reproduced, with errors evenly distributed throughout the system.

Validity of the 2D Approximation

The geometry and size of the SLFE (notably, the small depth relative to horizontal width) combined with the effects of strong tides in a freshwater environment make the assumptions leading to a 2D approximation valid. Furthermore, recent velocity measurements made in the SLFE (Matte et al. 2014b) have shown that virtually no vertical structures are present in the flow, which is dominantly 2D. In fact, current reversals generally occur in a very short period of time and are nearly simultaneous within the water column. Secondary circulation, conversely, was observed only at Portneuf, a region of high channel curvature. Although the vertical variability in velocity is not accounted for by the 2D shallow-water model, the results for this region were very accurate (cf. Table 3). In the estuarine transition zone, between Orleans Island and Île aux Coudres (Fig. 1), the presence of density currents and stratification was observed by Simons et al. (2010), which varies at both semidiurnal and fortnightly timescales. Water-column stratification may cause a neap-spring periodicity in the bottom drag coefficient and affect tidal-height prediction (Spitz and Klinck 1998; Zhong and Li 2006). Although the 2D vertically averaged model developed here cannot resolve the vertical flow structure or density currents in this transition zone, the impact on water levels is believed to be small, as evidenced by high model accuracy in tidal levels, both in the estuarine transition zone and upstream (cf. Table 3).

Error Sources and Resolution

High-resolution meshes allow for a decrease in the overall discretization error and improvement in, for example, the prediction of subtidal velocities, compared with coarser grids, particularly over important topographic features (Rayson et al. 2015). However, accurate predictions are strongly dependent on the precision of tidal and river forcing at the boundaries, bathymetry, and bottom roughness. In fact, in regions of complex bathymetry, an increase in the grid resolution cannot be obtained without an increase in the bathymetric resolution, and greater care must be taken in parameterizing bottom friction. Moreover, nonhydrostatic effects may become important for small-scale flow structures, but they have a small influence on the large-scale dynamics (Wang et al. 2009). Because our main focus here was to model processes at the tidal scale reported in recently acquired field data, the model was not configured to reproduce phenomena at scales that are influenced by the nonhydrostatic pressure. Rather, local refinements of the mesh were intended for a better inclusion of topographic elements influencing the mean flow, including the navigational channel and shallow intertidal areas.

Model results in the downstream segment presented larger errors than in the upstream portion of the SLFE. Past Orleans Island, the most significant error sources were related to the downstream boundary conditions and the particularly complex bathymetry, characterized by a large number of islands. The strong tides typical of this area are therefore much more likely to be influenced by these boundary conditions and bathymetry than by spatial variations in bottom friction. The actual lateral distribution of water levels along the downstream boundary could not be optimized entirely. The assimilation performed here only allowed for the definition of laterally invariant time lags, vertical shifts in water levels, and amplitude corrections. Coupling with a 3D ocean model of the St. Lawrence Estuary and gulf (Saucier et al. 2003; Saucier et al. 2009; Smith et al. 2013) may allow for an improvement in water-level forcing at the downstream boundary. Numerical experimentation showed this exercise to be nontrivial, stressing the need for numerical investigations to be complemented with field data. Moreover, despite the detailed LIDAR topographic data that were acquired for shallow areas, bathymetry in the Montmagny archipelago, downstream of Orleans Island, is very complex, and the only available data in some areas were measured several decades ago. Also, data reduction from chart datum to mean sea level is not trivial in this region because of the high spatial variability in tidal properties. This therefore introduced further errors in the model, which are likely to be more pronounced during low spring tides (Falcão et al. 2013). One possible way to improve model performance would be to treat the bathymetric error as a calibration parameter (Cea and French 2012).

Broad Implications

The SLFE is among the largest tidal rivers in North America, based on its size, the volume of freshwater it discharges to the ocean, and the amplitude of its tides. It therefore constitutes an ideal natural laboratory for studying the complex interaction between tides, river flow, and system geometry, the implications of which can be translated to many similar energetic systems worldwide. The methodological framework provided here should apply to the setup, calibration, and validation of numerical models in freshwater tidal environments in general and forms the basis for the analysis presented in Part II. The following elements of the modeling process should be considered to build confidence in model results:

1.

Local refinements of the mesh should be intended for a better inclusion of topographic elements influencing the mean flow, such as the navigational channel, engineering structures, shoals, and shallow intertidal areas. Further refinements may allow for resolution of smaller-scale circulation patterns if combined with shorter time steps and appropriate viscosity parameterization.

2.

The accuracy and density of topographic and bathymetric data may present significant spatial variations; combined with grid resolution, these can directly affect the computed velocity field, especially the flooding–drying processes occurring over intertidal regions (cf. Part II).

3.

Bed friction associated with substrate composition, bedforms, and/or macrophytes directly acts on tidal ranges and MWLs (used herein as metrics to assess model performance) and on tidal propagation and modulation properties (cf. Part II). However, adjustment of the friction coefficients alone may not suffice to calibrate velocities effectively. In this case, bathymetric errors can be treated as a calibration parameter in areas where their effects on the computed fields are likely (or known) to be significant.

4.

Errors in boundary conditions directly translate into errors in modeled hydrodynamic fields. Imposed conditions should be assimilated when necessary, within the range prescribed by the estimated error and variability in the data.

5.

Calibration and validation should be an iterative process during which model accuracy is assessed under a variety of conditions corresponding to distinct periods of analysis and locations throughout the system. A versatile approach (cf. Part II) consists of model evaluation based on a series of physical variables conjointly used for comparisons between events occurring at different frequencies and scales in time (e.g., intratidal, fortnightly, seasonal) and in space (e.g., cross-sectional, longitudinal).

Conclusion

A 2D hydrodynamic model of the SLFE was developed based on a finite-element grid with an average spatial resolution of 50 m, far denser than previous/existing models. The model was calibrated and validated using water-level and velocity data collected in the summer of 2009, constituting the most detailed data set used to date in this section of the St. Lawrence. Results showed good agreement between modeled and observed water levels, with prediction skills higher than 0.99 at all stations and RMSEs corresponding to less than 5% of the local tidal ranges in the first 186 rkm; at upstream stations where tidal ranges are significantly reduced, RMSEs were lower than 6 cm. Similarly, errors in discharge remained within 6% of the maximum observed discharges at 11 of the 13 surveyed transects; larger relative errors at the two remaining sections can mainly be attributed to interpolation errors in the transect data and local bathymetric uncertainties. In Part II, the ability of the model to reproduce tidal and flow features observed in the field data is demonstrated. Together, these results confirm that the terrain description, boundary conditions, and parameterizations used in the model were well defined.

This research provides, for the first time, a detailed 2D description of the tidal hydrodynamics of this complex region, thoroughly validated with recent field data. This is the first step toward the development of a comprehensive model of the SLFE ecosystem that will be run in operational mode and include variables for the assessment of habitat and water quality [e.g., Morin et al. (2003)].

Appendix: Shallow-Water Equations

The shallow-water equations as implemented in H2D2 can be written as follows (Heniche et al. 2000) for mass conservation [Eq. (3)] and momentum conservation [Eqs. (4a) and (4b)]:

\frac{\partial h}{\partial t} + \frac{\partial q_{x}}{\partial x} + \frac{\partial q_{y}}{\partial y} - (γ + δ) (\frac{\partial^{2} h}{\partial x^{2}} + \frac{\partial^{2} h}{\partial y^{2}}) = 0

(3)

\frac{\partial q_{x}}{\partial t} + \frac{\partial}{\partial x} (\frac{q_{x} q_{x}}{H}) + \frac{\partial}{\partial y} (\frac{q_{x} q_{y}}{H}) + c^{2} \frac{\partial h}{\partial x} - \frac{1}{ρ} \times [\frac{\partial}{\partial x} (H τ_{x x}) + \frac{\partial}{\partial y} (H τ_{x y}) + τ_{x}^{s} - τ_{x}^{b}] - f_{c} q_{y} = 0

(4a)

\frac{\partial q_{y}}{\partial t} + \frac{\partial}{\partial x} (\frac{q_{y} q_{x}}{H}) + \frac{\partial}{\partial y} (\frac{q_{y} q_{y}}{H}) + c^{2} \frac{\partial h}{\partial x} - \frac{1}{ρ} \times [\frac{\partial}{\partial x} (H τ_{y x}) + \frac{\partial}{\partial y} (H τ_{y y}) + τ_{y}^{s} - τ_{y}^{b}] - f_{c} q_{y} = 0

(4b)

where x(x, y) = east and north Cartesian coordinates (m); t = time (s); h = water level (m); and q(q_x, q_y) = specific discharge (m²/s), defined as

q = u H

(5)

where u(u_x, u_y) = water velocity (m/s); and H = water depth (m), defined as

H = h - z_{b}

(6)

where z_b = bed level with respect to the mean sea level; and c = celerity of waves (m/s), defined as

c = \sqrt{g H}

(7)

where g = gravitational acceleration (9.81 m/s); ρ = density of water (10³ kg/m³); and γ = Lapidus coefficient, defined as (Heniche et al. 2000; Lapidus 1967)

γ = γ_{0} Δ^{2} \sqrt{{(\frac{\partial h}{\partial x})}^{2} + {(\frac{\partial h}{\partial y})}^{2}}

(8)

where γ₀ = a constant value set to 2.5 × 10⁻⁵; Δ = local element size (this added viscosity damps the solution only in regions with a high water-level gradient, thus preventing oscillations of the free surface); δ = hydraulic conductivity (set to 0 in wet areas and to 1 in dry areas); and

τ_{ij}

= Reynolds stress (kg/s²/m), with each stress component being defined as

\frac{1}{ρ} [\begin{matrix} τ_{x x} & τ_{x y} \\ τ_{y x} & τ_{y y} \end{matrix}] = (ν_{l} + ν_{t} + ν_{n}) [\begin{matrix} 2 \frac{\partial u_{x}}{\partial x} & (\frac{\partial u_{x}}{\partial y} + \frac{\partial u_{y}}{\partial x}) \\ (\frac{\partial u_{x}}{\partial y} + \frac{\partial u_{y}}{\partial x}) & 2 \frac{\partial u_{y}}{\partial y} \end{matrix}]

(9)

where ν_l = laminar viscosity; ν_t = turbulent viscosity; and ν_n = numerical viscosity; and ν_t is expressed either as a constant viscosity or as a function of the flow gradient, derived from the mixing length theory as (Rodi 1984)

ν_{t} = l_{m}^{2} \sqrt{2 {(\frac{\partial u_{x}}{\partial x})}^{2} + 2 {(\frac{\partial u_{y}}{\partial y})}^{2} + {(\frac{\partial u_{x}}{\partial y} + \frac{\partial u_{y}}{\partial x})}^{2}}

(10)

where l_m = mixing length and is defined as (Ouellet et al. 1986; Soulaïmani 1983)

l_{m} = λ H

(11)

with a calibration coefficient λ set to 1 [this representation of the mixing length differs from Smagorinsky’s (1963) subgrid model in that the length scale at which turbulent processes can be represented is limited by a fraction of the local depth instead of the element size Δ; in the present application, the mean depth is overall smaller than the average element size, although spatial variations of the depth do not necessarily follow the variations in mesh resolution]; ν_n is controlled by the Peclet number P (Zienkiewicz et al. 2014), as follows:

ν_{n} = | \frac{q}{H} | \frac{Δ}{P}

(12)

τ_{i}^{s}

= surface friction (N/m), defined as

τ_{i}^{s} = ρ_{a} C_{w} | w | w_{i}

(13)

where ρ_a = air density; C_w = wind drag coefficient (cf. Heniche et al. 2000); and w(w_x, w_y) = wind velocity;

τ_{i}^{b}

= bottom friction (N/m²), defined as

τ_{i}^{b} = (α + \frac{ρ g n^{2} | q |}{H^{7 / 3}}) q_{i}

(14)

where α = a damping coefficient acting as a linear friction parameter (set to 0 in the present application); n = the Manning coefficient, which defines quadratic resistance to flow by substrate, macrophytes, macrorugosity, and so forth, with each being quadratically added as follows (Boudreau et al. 1994):

n^{2} = \sum_{m} n_{m}^{2}

(15)

f_c = Coriolis factor (s⁻¹), defined as

f_{c} = 2 ω \sin ϕ

(16)

where

ω

= Earth’s rotational rate; and

ϕ

= latitude.

The weak form of the shallow-water equations is derived via the Galerkin weighted residuals method (Dhatt et al. 2005; Heniche et al. 2000). The higher-order terms and the continuity equation are integrated by parts, leading to a natural condition of impermeability on solid boundaries where no explicit condition is specified (i.e., q_n = 0, where q_n is the specific discharge normal to the boundary).

Acknowledgments

Work by Pascal Matte was supported by scholarships from the Natural Sciences and Engineering Research Council of Canada and Fonds de recherche du Québec—Nature et technologies. The authors thank Environment Canada (Meteorological Service of Canada) for financial support; the Ministère du Développement Durable, de l’Environnement et de la Lutte contre les Changements Climatiques (MDDELCC) for financing the LIDAR campaign; and the Canadian Hydrographic Service (CHS) and Ministère des Transports du Québec (MTQ) for providing bathymetric and topographic data. Special thanks go to Olivier Champoux, Patrice Fortin, Jimmy Poulin, Charles Gignac, and Alain Soucy for their contribution to this work. The authors also thank Daniel Bourgault for his valued comments on a previous version of the manuscript and two anonymous reviewers for their constructive comments. Tide gauge data are available for download from the DFO’s Canadian Tides and Water Levels Data Archive (http://www.meds-sdmm.dfo-mpo.gc.ca/isdm-gdsi/twl-mne/index-eng.htm). The H2D2 software and source code are freely available for download, upon request to Yves Secretan, at http://www.gre-ehn.ete.inrs.ca/H2D2/contenu_download.

References

Arhonditsis, G. B., and Brett, M. T. (2004). “Evaluation of the current state of mechanistic aquatic biogeochemical modeling.” Mar. Ecol. Prog. Ser., 271, 13–26.

Abstract

Introduction

Field Campaigns

Hydrodynamic Measurement Campaign

LIDAR Campaign

Model Setup

Simulation Domain and Boundaries

Numerical Terrain Model

Topography

Friction

Finite-Element Mesh

Hydrodynamic Model

Drying–Wetting Model

Resolution

Calibration

Assimilation

Validation

Results and Discussion

Statistics on Water Levels

Harmonic Properties

Statistics on Discharges

Validity of the 2D Approximation

Error Sources and Resolution

Broad Implications

Conclusion

Appendix: Shallow-Water Equations

Acknowledgments

References

Information

Published In

Copyright

History

Authors

Affiliations

Metrics

Citations

Download citation

Cited by

Figures

Other

Share

Copy the content Link

Share with email

Share

Request Username

Create a new account

Change Password

Password Changed Successfully

Verify Phone

Congrats!