Modelling Behaviour of the Salt Wedge in the Fraser River and Its Relationship with Climate and Man-Made Changes

Abstract

Agriculture is an important industry in the Province of British Columbia, especially in the Lower Mainland where fertile land in the Fraser River Delta combined with the enormous water resources of the Fraser River Estuary support extensive commercial agriculture, notably berry farming. However, where freshwater from inland meets saltwater from the Strait of Georgia, natural and man-made changes in conditions such as mean sea level, river discharge, and river geometry in the Fraser River Estuary could disrupt the existing balance and pose potential challenges to maintenance of the health of the farming industry. One of these challenges is the anticipated decrease in availability of sufficient freshwater from the river for irrigation purposes. The main driver for this challenge is climate change, which leads to sea level rise and to reductions in river flow at key times of the year. Dredging the navigational channel to allow bigger and deeper vessels in the river may also affect the availability of fresh water for irrigation. In this study, the salinity in the river was simulated using H3D, a proprietary three-dimensional hydrodynamic numerical model which computes the three components of velocity (u,v,w) in three dimensions (x,y,z) on a curvilinear grid developed specially for Fraser River, as well as scalar fields such as salinity and temperature. The results indicate various levels of impact to the salinity in the river and adaptive measures must be established to maintain the long-term viability of the industry. This study found that sea level rise and changes in river discharge would have a larger impact on the availability of fresh water than would channel deepening at the present sea water level. In a low river discharge regime, the impact from sea level change is more significant than in the high river discharge regime. On the other hand, the influence from changes in river discharge on withdrawal appears to increase when water level is lowered. Dredging the channel to accommodate larger vessels with deeper draft would further affect the salinity and shorten the withdrawal window; the effect of channel deepening becomes more pronounced in the lower flow period.

1. Introduction

Many of the largest economies in the world are located near or in a river estuary where saline water from the ocean meets the freshwater draining from inland: The Hudson River in New York, the Mississippi River in New Orleans, the Buffalo Bayou in Houston, the Yellow River in Shanghai, and the Pearl River in Hong Kong, just to name a few. These are not coincident, because an estuary often is of critical commercial value to nearby areas; and altogether, estuaries serve to contribute not only to human welfare and value to the world’s economy, but also provide essential ecosystem services such as food production and recycling of nutrients [1]. Among the many commercial functions it can serve, an estuary can be a conduit for transport of goods and products, and is very often an important water resource for domestic, industrial, and agricultural use. These economic benefits brought forth by an estuary are of course in addition to the environmental benefits resulting from a highly diverse ecosystem that an estuary usually sustains due to its unique hydrodynamic and geomorphologic settings.

The hydrodynamic behavior and salinity of an estuary is sensitive to climate change that could potentially affect sea level and the hydrological characteristics that dictate the amount of watershed runoff. This delicate balance is further complicated by future changes in river management such as channel deepening in order to allow navigation by larger vessels in the estuary. Understanding the complicated interaction between changes in sea level and runoff as well as bathymetry of the river and their impacts on the salinity in the river is therefore essential for developing proper and timely strategies to ensure resiliency of coastal areas against natural and anthropogenic changes in the future.

The hydrodynamics of the Fraser River estuary has been extensively studied in the past: Stronach et al. [2,3], and Halverson and Pawlowicz [4,5] studied and described the river plume dynamics of the Fraser River outflow into the Strait of Georgia. Halverson and Pawlowicz [6], Kostaschuk and Atwood [7], and Ward [8] particularly investigated the river flow and tidal forcing in the Fraser River estuary and their effects on salinity. Yin et al. [9] studied the biological significance of the salt wedge dynamics in the Fraser River estuary. Neilson-Welch and Smith [10] observed and described the effects of the upstream intrusion of saltwater.

There is limited, if any, published literature regarding specifically the effects of sea level rise (SLR) and changes in river discharge and river geometry on salinity distribution in the Fraser River, but similar studies were undertaken in other parts of the world: Krvavica et al. [11] studied the salt wedge in the Rjecina River Estuary in Croatia using a finite volume model; Funahashi et al. [12] developed a numerical model using Delft3DFlow to investigate the response of the salt wedge and salinity distribution in the river to the short-term changes in sea level and river discharge in the Yura Estuary in Japan; and Shaha et al. [13] studied the saltwater intrusion in Sumjin River Estuary in South Korea using the Finite Volume Community Ocean Model (FVCOM) and related intrusion length and river discharge by a power law.

2. Study Site

The Fraser River estuary, located south of Vancouver, is an important part of the economy in the Province of British Columbia (BC), Canada [14]. The Fraser River delta is fertile and, as such, the agricultural industry in the region is thriving, with a noticeable increase in the farmland area in the late 20th and early 21st centuries [15]. There is potential for significant growth in terms of agricultural production, provided that additional water can be withdrawn from the river for irrigation [16]. The low-salinity water withdrawn from the river for irrigation is therefore of critical importance to the continued health of the farmland and its agricultural products.

The Fraser River is connected at its downstream boundary to the Strait of Georgia (SOG), which in turn is connected to the open Pacific Ocean through the connections at Juan de Fuca Strait to the south and Johnstone Strait to the north. Therefore, water level in the SOG and the Fraser River is subject to tidal fluctuations. Figure 1 below illustrates the study location as well as the geographic and hydrographic relations between the Fraser River, the SOG, and the Pacific Ocean.

The tidal range in the Fraser River is variable, depending not only on the tidal fluctuation in the SOG, but also on distance upstream from the river mouth and on river discharge. At the river mouth, the tide is mixed, mainly semi-diurnal, and exhibits a clear alternating neap tide and spring tide pattern. The higher high water level is 2.0 m above mean sea level, which in turn is 3.1 m above chart datum, resulting in a tidal range of 5.1 m. This tidal effect, though diminishing in an upstream direction, can reach as far upstream as Mission, approximately 84 km from the river mouth.

The stratification in the SOG and the Fraser River is due to the presence of saltwater from the Pacific Ocean and freshwater from river runoff. The interaction of the freshwater and saltwater masses in the Fraser River Estuary gives rise to the formation of a salt-wedge, which intrudes from the SOG into the river in an upstream direction. While the salt wedge occupies the part of the water column from the bottom to a certain depth, the fresh river water, upon meeting the salt wedge, travels downstream on top of the salt-wedge. The extent of the intrusion varies depending on the river discharge and water level in the SOG: while the position and salinity of the salt-wedge displays an intra-day pattern in response to the semi-diurnal tidal fluctuation in the Strait, the mean daily position of the salt-wedge moves upstream and downstream in response to the changes in river discharge. The salinity in the river has a strong correlation with the salt-wedge intrusion and thus displays similar fluctuation patterns.

3. Materials and Methods

The core component of the study is the numerical simulation of the salt wedge in the Fraser River, the largest river in BC, which drains into the SOG in the Salish Sea, under a range of climate change, river flow, and river configuration scenarios. A numerical model of the Fraser River, based on the modelling platform H3D developed by Tetra Tech, has been implemented and calibrated for present conditions; however, the fundamental processes: bottom friction, vertical mixing, salt wedge migration under tidal influence, are all intrinsic processes, i.e., their formulation depends only on the specific flow and water level conditions imposed by external influences. In other words, the model is totally objective with respect to evaluating climate change. Thus, the model does not require modification for this work. However, it is important to select and properly quantify climate change scenario inputs for modelling.

3.1. Three-Dimensional Hydrodynamic Numerical Model: H3D

H3D is an implementation of the numerical model developed by Backhaus [17,18], which has had numerous applications to the European continental shelf [19,20], Arctic waters [21,22], and deep estuarine waters [23]. Locally, H3D has been used to model the temperature structure of Okanagan Lake [24], the transport of scalar contaminants in Okanagan Lake [25], sediment movement and scour/deposition in the Fraser River, circulation and wave propagation in Seymour and Capilano dams, and salinity movement in the lower Fraser River. H3D forms the basis of the model developed by Saucier and co-workers for the Gulf of St. Lawrence [26], and has been applied to the Gulf of Mexico [27].

3.1.1. Theoretical Basis

H3D is a three-dimensional time-stepping numerical model which computes the three components of velocity (u,v,w) on a structured grid in three dimensions (x,y,z), as well as scalar fields such as salinity, temperature, and contaminant concentrations. The model uses the Arakawa C-grid in space [28], and uses a two level semi-implicit scheme in the time domain. The Arakawa C-grid evaluates velocity vector components at the centres of the grid faces (i.e., u-component velocities are evaluated at the centres of the east and west grid faces, v-component velocities on the north and south grid faces, and w-component velocities on the top and bottom grid faces, while scalar quantities are calculated in the grid centres). The size of the model time step varies depending on the horizontal grid size and velocity at each of the grid cells: the model, at each time step and each grid cell, calculates the numerically-stable time step based on the u-, v-, and w-components of the velocity and the x-, y-, and z-dimensions of the grid cell; the model then updates the time step by choosing the largest time step that will satisfy the numerical stability for all grid cells. Details with respect to the derivation of the model are provided in Appendix A.

3.1.2. Grid Geometry

In the horizontal, the grid cells can be either rectilinear or orthogonal curvilinear, depending on the physical characteristics of the water body of interest. Details regarding the horizontal resolution and configuration of the models are described in Section 3.2. In the vertical, the levels near the surface are typically closely spaced to assist with resolving near-surface dynamics. Figure 2 shows a typical grid mesh utilized for the H3D model.

In addition, the model is capable of dealing with relatively large excursions in overall water level as the water level rises and falls in response to varying inflows and outflows, by allowing the number of near-surface layers to change as the water level varies. That is, as water levels rise in a particular cell, successive layers above the original layer are turned on and become part of the computational mesh. Similarly, as water levels fall, layers are turned off. This procedure has proven to be quite robust, and allows for any reasonable vertical resolution in near-surface waters. When modelling thin river plumes in areas of large tidal range, the variable number of layers approach allows for much better control over vertical resolution than does the σ-coordinate method.

In addition to tides, the model is able to capture the important response, in terms of enhanced currents and vertical mixing, to wind-driven events. This is achieved by applying wind stress to each surface grid point on each time step. Vertical mixing in the model then re-distributes this horizontal momentum throughout the water column. Similarly, heat flux through the water surface is re-distributed by turbulence and currents in temperature simulations.

3.1.3. Turbulent Closure

Turbulence modelling is important in determining the correct distribution of velocity and scalars in the model. The diffusion coefficients for momentum (A_H and A_V) and scalars (N_H and N_V) at each computational cell are dependent on the level of turbulence at that point. H3D uses a shear-dependent turbulence formulation in the horizontal presented by Smagorinsky [29]. The basic form is:

$$A_H=A_{H_0}dxdy\sqrt{{(\frac{du}{dx})}^2+{(\frac{dv}{dy})}^2+\frac{1}{2}{(\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y})}^2}$$

(1)

where the parameter A_H0 is a dimensionless tuning variable, and experience has shown it to lie in the range of 0.25 to 0.45 for most water bodies such as rivers, lakes and estuaries.

A shear and stratification dependent formulation, the Level 2 model of Mellor and Yamada [30], is used for the vertical eddy diffusivity. The basic theory for the vertical viscosity formulation is taken from an early paper by Mellor and Durbin [31]. The evaluation of length scale is based on a methodology presented by Mellor and Yamada [30].

For scalars, both horizontal and vertical eddy diffusivity are taken to be similar to their eddy viscosity counterparts, but scaled by a fixed ratio from the eddy viscosity values. Different ratios are used for the horizontal and vertical diffusivities. If data is available for calibration, these ratios can be adjusted based on comparisons between modelled and observed data. Otherwise, standard values based on experience with similar previously modelled water bodies are used. In a recent reservoir simulation with good calibration data, the ratio of vertical eddy diffusivity to vertical eddy viscosity was 0.75 and the ratio between horizontal eddy diffusivity and horizontal eddy viscosity was 1.0.

3.1.4. Scalar Transport

The scalar transport equation implements a form of the flux-corrected algorithm presented by Zalesak [32], in which all fluxes through the sides of each computational cell are first calculated using a second-order method. Although generally more accurate than a first order method, second order flux calculations can sometimes lead to unwanted high frequency oscillations in the numerical solution. To determine if such a situation is developing, the model examines each cell to see if the computed second order flux would cause a local minimum or maximum to develop. If so, then all fluxes into or out of that cell are replaced by first order fluxes, and the calculation is completed. As noted, the method is not a strict implementation of the Zalesak method, but is much faster and achieves very good performance with respect to propagation of a Gaussian distribution through a computational mesh. It does not propagate box-car distributions as well as the full Zalesak method, but achieves realistic simulations of the advection of scalars in lakes, rivers, and estuaries, which is the goal of the model. This scheme as implemented is thus a good tradeoff between precision and execution time, important since in many situations, where more than one scalar is involved, the transport-diffusion algorithm can take up more than half the execution time.

3.1.5. Heat Flux at the Air–Water Interface

The contribution of heat flux to the evolution of the water temperature (T) field is schematized in H3D as:

$$\frac{\partial T}{\partial t}=\frac{\Delta F}{\rho \ast C_P\ast h}$$

(2)

where $\Delta F$ is the net heat flux per unit area retained in a particular layer, ρ is the density of water, C_p is the heat capacity of water, and h is the layer thickness. Heat flux at the air–water interface incorporates the following terms:

• F_in: incident short wave radiation. Generally, this is not known from direct observations. Instead, it is estimated from the cloud cover and opacity observations at nearby stations, a theoretical calculation of radiation at the top of the atmosphere based on the geometry of the earth/sun system, and an empirical adjustment based on radiation measurements at Vancouver International Airport for the period 1974–1977. This procedure has worked well for many water bodies, notably Okanagan Lake and the waters of the north coast of British Columbia, in terms of allowing H3D to reproduce the observed temperature distributions in space and time. Values for albedo as a function of solar height are taken from Kondratyev [33].
• F_back: net long wave radiation. This is calculated according to Gill [34], involving the usual fourth power dependence on temperature, a factor of 0.985 to allow for the non-black body behaviour of the ocean, a factor depending on vapor pressure to allow for losses due to back radiation from moisture in the air, and a factor representing backscatter from clouds.
• F_L and F_S: Latent heat flux (F_L) is the heat carried away by the process of evaporation of water. Sensible heat flux (F_S) is driven by the air-water temperature difference and is similar to conduction, but assisted by turbulence in the air. Latent and sensible heat fluxes are described by:

  $$Q_L=1.32e^{-2}\ast L\ast U_a\ast (q_{obs}-q_{sat})\ast F_l,\mathrm{and}$$

  (3)

  $$Q_S=1.46e^{-3}\ast {\rho }_{air}\ast c_p\ast U_a\ast (T_{air}-T_{water})\ast F_s,$$

  (4)

  where q_obs and q_sat are the observed and saturated specific humidities, T_air and T_water are the air and water temperatures, L is the latent heat of evaporation of water, and C_p is the heat capacity of water, U_a is wind speed, F_l and F_s are latent and sensible factors, which are the scaling factors introduced to account for local factors, and can be adjusted, when needed, to achieve better calibration of the model. Typically, the only adjustment is that sensible factor is doubled when the air temperature is less than the water or ice surface temperature to account for increased turbulence in an unstable air column.

As light passes through the water column it is absorbed and the absorbed energy is a component of the energy balance that drives water temperature. H3D assumes that light attenuation follows an exponential decay law:

$$E(z)=E(z_0)\ast e^{-k\ast (z-z_0)},$$

(5)

where E(z) is the insolation at a distance, z, below the water surface, and E(z₀) is the incident insolation at the water surface. The model computes the energy at the top and bottom of each layer and the difference is applied to the general heat equation in that layer. The extinction coefficient (k) is related to the Secchi depth (D_s) by:

$$k=\frac{2.1}{D_s}.$$

(6)

Temperature is treated like any other scalar as far as advection and diffusion are concerned.

3.2. Model Implementation

A study of the details of the hydrodynamics of the Fraser River requires model nesting to better resolve small scale processes at and near locations of interest. The model used for this study operates in a triple-nested configuration, shown in Figure 3.

The investigation of the behaviour of the salt wedge was undertaken with a nominally 50-m resolution curvilinear model that spans the lower 41 km of the Fraser River, from Sand Heads to the Port Mann Bridge. The model uses 50-m resolution in the along-channel direction, and 20-m in the cross-channel direction. This 50-m resolution model is in turn embedded within a nominally 1-km resolution rectilinear model of the entire SOG. The vertical resolution of the models varies with depth and is 1 m for the top 10 m, 2 m for depth of 10–16 m, 3 m for 16–25 m, 10 m for 25–55 m, 25 m for 55–105 m, 50 m for 105–305 m, and 100 m for 305–505 m. The higher vertical resolution near the top is necessary to resolve the complex hydrodynamic processes and the sharp salinity gradient at the interface of the salt wedge.

Both models simulate tidal, wind-driven, and density-driven currents. Water level, velocity components, and any scalar quantities (i.e., temperature and salinity) output from the coarse grid model are passed on along the boundaries of the fine grid model and used to drive the finer-scale implementation of H3D. The fine-grid implementation provides the details of the effect of small-scale spatial variability in shorelines, depths and structures such as the tunnel cover.

The 1-km SOG model, driven by wind, temperature, and salinity as well as tidal conditions along its open boundaries bordering the northern entrance to the SOG and the western entrances to Juan de Fuca Strait, includes a coarse representation of the Fraser River, extending 41 km upstream from the river mouth at Sand Heads to the Port Mann Bridge. At that location, upstream of all salt wedge penetration, the model is dynamically coupled to a one-dimensional (1-D) model of the Fraser River, extending to Hope, which is located approximately 130 km upstream from the river mouth and is free of any tidal effects. The model coupling is two-way as the SOG model provides water level to the downstream boundary of the 1-D model, which in turn calculates the flow conditions at that same boundary given the discharge rate at the upstream boundary of the 1-D model at Hope. The initial and boundary conditions for temperature and salinity are provided by Crean and Ages [35], which documented a comprehensive temperature and salinity measurement campaign in the SOG and Juan de Fuca Strait. Tidal conditions are also specified along the open boundaries.

The 50-m lower Fraser River model is driven at its upstream end by flow boundary conditions provided by a separate, dynamically coupled 1-D Fraser River model. At the downstream end, water levels as well as temperature and salinity profiles are obtained from the 1-km SOG model, spatially interpolated from those cells of the 1-km SOG model that correspond to the boundaries of the 50-m Fraser River model.

Table 1. Modelling sequence of the study.

Seq. #	Model	Initial Conditions	Downstream Boundary Conditions	Upstream Boundary Conditions	Dynamic Model Coupling	Objective
1a	1-km SOG	Temperature and salinity profiles	Temperature, salinity and tidal level	Flow conditions provided by the 1-D Fraser River model	1-D Fraser River model for the 1-km SOG model	To provide boundary conditions to the 50-m Fraser River Model
1b	1-D Fraser River-for the 1-km SOG	Water level	Water level from the 1-km SOG model	River hydrograph at Hope	1-km SOG model	To provide upstream boundary conditions to the 1-km SOG model
2a	50-m Fraser River	Temperature and salinity profiles	Temperature, salinity and tidal level provided by the SOG model	Flow conditions provided by the 1-D Fraser River model	1-D Fraser River model for the 50-m Fraser River model	To simulate spatial and temporal salinity distributions in the Fraser River
2b	1-D Fraser River-for the 50-m Fraser River	Water level	Water level from the 1-km SOG model	River hydrograph at Hope	50-m Fraser River model	To provide upstream boundary conditions to the 50-m Fraser River model

below summarizes the modelling sequence (Seq. #), the initial conditions and boundary conditions utilized for the modelling study:

3.3. Validation of the Nested Model

The models were initialized on 1 January 2011 and run through the 6-month long harvest season period, between August and December, when the ability to withdraw freshwater from the river is crucial. The year 2011 was chosen for the model validation as this is the year for which bank-to-bank bathymetric survey data collected in the Fraser River estuary is available: the model bathymetry was constructed using this 2011 data set; the Fraser River flow rate in 2011 was used to drive the upstream boundary of the river model, and the tidal conditions in 2011 were specified along the open boundaries in the SOG model. The flow rate at the upstream boundary of the model is the combination of the flow rate at Hope and the estimated runoffs that report to the river downstream of Hope and upstream of the Port Mann Bridge. Figure 4 shows the observed river flow rate at Hope in the year 2011 [36].

The model was validated against water levels recorded at New Westminster, which is located 34 km upstream from the river mouth, and against conductivity data collected by a sensor mounted at 2 m depth from a floating platform installed at an irrigation intake near 8081 River Road, Delta, approximately 24 km upstream from the river mouth. Figure 5 shows a comparison of observed and modelled conductivity values at the water intake (location shown in Figure 6) at a depth of 2 m from 3rd to 23rd November 2011, when the leading edge of the salt wedge passed upstream and downstream through the sensor location and the conductivity signal ranged between zero and the maximum value that the sensor could record, 5500 μS/cm. Comparison of the observed and modelled results for months other than November would have been undertaken if not for the lack of reliable, continual conductivity measurement at other locations in the river (identified as data gap in Section 5.4). Also included in the figure are observed and modelled water levels for the same period. Black lines show modelled values and red lines show observed values. Since the conductivity sensor cuts off at 5500 μS/cm, the model results were similarly cut-off to facilitate comparison.

The days on which the observed conductivity is considerably lower than the modelled correspond almost perfectly with the days with rainfall in the catchment downstream of the Harrison River exceeding 15 mm/day: 7, 10–12, 17, 21–22 November. A rainfall of 15 mm/day corresponds, approximately, to an additional flow in the Fraser River of up to 1000 m³/s, depending on the hydrological processes controlling how this additional water reports to the river. This additional fresh water would undoubtedly lead to significant reduction in conductivity in the river during rainfall events.

Otherwise, the model is, in general, able to re-create the timing and the trend of the conductivity signal, as well as partially the strength of the signal (Correlation Coefficient (R) = 0.53; Root-Mean-Square-Error (RMSE) = 1571 μS/cm or 1.37 ppt), in both the daily time scale and multi-day time scale. However, there are some disagreements in the conductivity values and the exact timing of the fluctuations at the water intake. For example, the model almost always predicts an elevation of conductivity during high tides when river flow is comparatively slower and water level in the river higher (for example, on 7 November); however, the sensor at the intake did not always detect such a conductivity signal.

The complexity in the behaviour of the conductivity signal can be partly understood by considering Figure 6, showing the map of salinity at the 2-m depth, on 11 November 2011, at 7 a.m., where the modelled result appears to deviate the most from the observed value. The vectors indicate the direction and strength of the water movement at the 2-m depth. It can be seen that there is a high degree of spatial variability in the salinity field in the vicinity of the intake. Salinity can vary from 3.5 ppt (4500 μS/cm) to more than 5.0 ppt (6500 μS/cm) near the intake in a matter of metres. Further analysis of model output demonstrates that there are two mechanisms for saline water to intrude onto the relatively shallow shelf on which the intake is located: either a selective withdrawal process, whereby saltier water is drawn up onto the bench from the adjacent deeper water (on both ebb and flood), or a process whereby the toe of the salt wedge rises to the surface upstream of the bench and then falls back partially onto the bench on the ebb tide.

Besides direct comparison of observed and modelled salinity, this model validation can be considered from the perspective of water availability, which describes the onset and offset of salinity intrusion at the water intake and the time window within which river water can be safely withdrawn under the criterion salinity value of 0.35 ppt or conductivity of 400 μS/cm. Figure 7 below compares the observed and modelled number of available hours per running 24 h. The red line represents the observation and the black line represents the model results. Only a short period of record is presented to facilitate visual comparison.

The model, even though more conservative in general, was able to predict the overall trend in availability (R = 0.60; RMSE = 5.9 h). The numerical model predicts that water suitable for irrigation in November is available for 10 h per day whereas observations indicate that suitable water is available, on average, 14 h per day.

3.4. Study Methodology

The effects of changes in environmental (sea level rise and river flow changes as a result of climate change) and anthropogenic (channel deepening) conditions on the behaviour of the salt wedge and salinity distribution in the Fraser River were the main focus of this study. Climate change scenarios were chosen based on the Fraser River flow rate predicted by a Global Climate Model (GCM), selected for its demonstrated ability to closely hindcast the historical flow in the river (Section 3.1), during the irrigation period between August and October. Several sea level rise scenarios were evaluated in this study and ultimately the sea level rises appropriate for the time horizon in year 2050 and year 2100 were used (Section 3.2). Given the projected economic growth of the region, there could be a need to allow deeper draft vessels to navigate the Fraser River to provide sufficient capacity. The dredge depth which allows for most Panamax and derivative vessel classes to navigate in the Fraser River was evaluated in the model study (Section 3.3).

3.4.1. Selecting Climate Scenario and Global Climate Model for Projection of the Fraser River Hydrograph

Three climate scenarios were considered in this study. Denoted A1B, B1, and A2, these scenarios, according to the Fourth Assessment Report (AR4) published by the Intergovernmental Panel on Climate Change (IPCC) in 2007 [37], relate to how different types of societal behaviour lead to different levels of greenhouse gas emissions and subsequent climate change. The emission scenarios based on AR4 were chosen for this study because the projected Fraser River hydrograph derived by the Pacific Climate Institute Consortium, or PCIC, is based on these scenarios.

Table 2. Climate scenarios and the associated social, economic, and human theme as presented in AR4.

Climate Scenario	Description	Characteristics
A1B	This scenario is of a more integrated world	Rapid economic growth A global population that reaches 9 billion in 2050 and then quickly declines The quick spread of new and efficient technologies A convergent world—income and way of life converge between regions. Extensive social and cultural interactions worldwide Good balance on all different types of energy sources
A2	These scenarios are of a more divided world	A world of independently operating, self-reliant nations Continuously increasing population Regionally oriented economic development
B1	This scenario is of a world more integrated, and more ecologically friendly	Rapid economic growth as in A1, but with rapid changes towards a service and information economy Population rising to 9 billion in 2050 and then declining as in A1 Reductions in material intensity and the introduction of clean and resource efficient technologies. An emphasis on global solutions to economic, social, and environmental stability

summarizes the characteristics of these climate scenarios.

The selection process of a suitable projected Fraser River hydrograph involved evaluation of a set of eight GCMs, driven by the aforementioned three climate scenarios, each representing relatively low, medium, and high atmospheric greenhouse gas concentration increases in the future. These eight GCMs were selected for best matching the downscaled results of these GCMs with the historical data in the 1950–2005 period in western Canada. This method produced an ensemble of all raw projections for temperature and precipitation at a global-scale. The GCMs provided these projections in the form of average temperature and precipitation changes in grid cells, where individual grid cells are approximately 300 km × 300 km (~90,000 sq. km or 9 million hectares) in size. These grid cell values were then statistically downscaled to a watershed-scale resolution of approximately 6 km × 6 km grid cells using a technique called Bias Corrected Spatial Disaggregation before being fed into the Variable Infiltration Capacity (VIC) hydrologic model for estimating streamflow. This approach was discussed in detail in Werner [38] and Schnorbus et al. [39].

An ensemble average was calculated based on the daily flow results derived from these eight different GCMs, each appropriately downscaled to produce Fraser River flows at Hope, and averaged to produce daily flows over an indicated time period. Figure 8 illustrates the ensemble average of the daily Fraser River flow rate under a range of climate change scenarios and different time horizon (present: 2000–2014, short-to mid-term horizon: 2016–2050; long-term horizon: 2051–2098). The data results were published by PCIC [40].

The predicted average flow rate in the irrigation period, under a range of climate change scenarios in 2050 (ranging between 1142 m³/s and 2352 m³/s) and 2098 (ranging between 730 m³/s and 1540 m³/s), is less than the existing rate (ranging between 1204 m³/s and 3025 m³/s).

The climate scenario, A1B, was eventually chosen for the projected Fraser River flow in this modelling study, not only because it represents a scenario which reflects the likely evolution of human society in the next century, but also it projects a lower river discharge and likely higher salinity during the irrigation period. As the original objective of this study was to assist engineers and planners to establish appropriate adaptive measures to maximize freshwater withdrawal for irrigation, the model results based on the A1B scenario would provide the relevant parties conservative estimates of the withdrawal window such that sufficient resources can be planned and assigned to mitigate the potential impact.

With regard to selecting between the various GCMs, Figure 9 below shows the comparison of the Fraser River hydrograph at Hope between the observed and modelled values derived from the 8 GCM models for calendar years from 1998 to 2004. The level of agreement and disagreement varies both between models, and from month to month and year to year.

This figure indicates that the HADCM (HADley Centre Coupled Model), MIROC (Model for Interdisciplinary Research on Climate), and CSIRO (Commonwealth Scientific and Industrial Research Organization) models simulated the historical Fraser River flow rate closest to the observed flow rate during the irrigation period between September and December. After comparing the Root-Mean-Squared-Error (RMSE) of these models for this critical irrigation period (HADCM: RMSE = 716 m³/s; MIROC: RMSE = 694 m³/s; CSIRO: RMSE = 721 m³/s; all other models: RMSE > 730 m³/s), the MIROC model was ultimately chosen for provision of the Fraser River hydrographs as input to the hydrodynamic model for simulation in the near to mid-term future (year 2050) and long-term future (year 2100) because MIROC best matches the observed river discharge rate during the irrigation period.

The Fraser River flow that drives the upstream boundary of the river model has been derived by determining the ensemble envelope for the river hydrographs at Hope based on the MIROC model, from which ‘dry’, ‘normal’, and ‘wet’ flow years appropriate for the timeframe of interest were determined. The number of annual hydrographs chosen for formation of the ensemble for each of three time frames (present, 2050 and 2100) was limited to the ten (10) closest years to the time frame of interest, except for the ‘present’ time frame for which the ten most recent, observed annual hydrograph were used. The 10-year time window size was chosen such that the timeframe is relevant to the reference points in the timeline (e.g., 2050 and 2100). If a longer time window is chosen, the variability within the window will inevitably include the trend from the long-term climate change, not just the result from natural variability.

Table 3. Fraser River hydrograph chosen for the ensemble analysis.

Timeframe	Year of Chosen Hydrographs for Ensemble Analysis	Type of Data
Present (2014)	2005–2014	Observed
2050	2045–2055	Predicted based on MIROC
2100	2089–2098 1	Predicted based on MIROC

¹ 2098 is the end of the river flow prediction based on MIROC published by PCIC as at the time of the study.

describes the hydrographs chosen for the analysis.

Figure 10 summarizes these hydrograph ensembles and properties for present conditions, and conditions in 2050 and 2100. The darker gray areas represent daily flow values between the 90th and 10th percentile, and the lighter gray bands indicate the minimum and maximum flows. The red, green and blue lines represent respectively the dry, normal, and wet years selected for simulation. The normal, dry, and wet years are defined here respectively as the median, driest, and wettest of all the hydrographs in a certain timeframe based on total discharge volume in summer and fall. Note that daily flow in a dry year may exceed the corresponding daily flow in a wet year for short periods, which better reflects the temporal variability of the system: not every day in a dry year is ‘dry’. This characteristic illustrates the considerable variability at all time-scales in the Fraser River flow.

3.4.2. Selecting Sea Level Rise Scenarios

The change in sea level is not uniform in all parts of the ocean. Past documents and reports published by the IPCC addressed mostly sea level rise on a global scale without giving significant consideration to local effects. The relationship between global and regional mean sea level rise is not simple, owing to the complex relationship between differences in air temperature change and the resulting differential thermal expansion, feedback mechanisms involving local meteorology and wind patterns and changes in tidal hydrodynamics in coastal areas as a result of changes in water depth and circulation.

As stated by the Sea Dike Guidelines [41], the effective sea level rise of 0.5–1.0 m between 2050 and 2100 and 2.0 m between 2100 and 2200 will occur along the BC coast. Given the uncertainty in the predictions and the future evolution of the climate, the upper bound value was a prudent choice for the worst-case scenario for sea level rise. Thus, this study used a sea level rise of 1.0 m for 2050 and 2.0 m for 2100.

Along the open boundaries of the coarse resolution SOG model, the mean sea level was adjusted upward by 1.0 m for SLR of 1.0 m and by 2.0 m for SLR of 2.0 m. All other tidal constituents along the open boundaries of the 1-km SOG model were assumed the same as in the case for SLR of 0.0 m.

3.4.3. Selecting River Dredge Depth Scenarios

Vancouver serves as a gateway to Asia and the rest of North America for sea-going vessels. Previous increases in shipping demand had led to deepening of the river channel. If the shipping demand continues to increase, then there may be a need to deepen the navigation channel. Such deepening would of course be subject to stringent environmental review, because of the potential damage to the ecologically-significant wetlands along the shores of the Fraser River, particularly well-known bird and fish habitats.

Two river dredge depth scenarios were investigated for their effects on the salinity in the river: 11.5 m and 20.0 m below geodetic datum. The current dredge practice of the river bottom, nominally 11.5 m, allows vessels as large as the third-generation Panamax vessels to navigate the river channel; the hypothetical new dredge practice of 20.0-m draft would upgrade to fully allow the fourth-generation Panamax (Post-Panamax, drafting 11–13 m) vessels, the fifth-generation Panamax (Post-Panamax-Plus, drafting 13–14.5 m), and the New-Panamax (drafting 15.2 m) vessels to travel in the river safely.

4. Results

The principal issue of concern is the extent of upstream salt wedge movement and the corresponding effects on near-surface salinities. The study will mainly focus on the salinity at the three candidate sites for future additional intake installations near the existing water intake, which is located approximately 24 km upstream from the river mouth. Figure 11 shows the location of the three sites (Site 1, Site 2, and Site 3). The red lines demarcate the upstream distance along the river from the mouth.

4.1. Difference in Salinity at Different Locations Along the River

Figure 12 shows modelled time-series of salinity at the intake depth at the three sites in a normal year with no sea level rise, over a brief period in September. The green line represents Site 1, the red line represents Site 2, and the blue line represents Site 3 (see Figure 11 for site locations). Figure 12a shows a period in which the salinity progressively increases in the downstream direction from Site 1 to Site 2 and then to Site 3; whereas Figure 12b illustrates the salinity at the three sites for one of the time periods during which the surface water at Site 1 (green line) is saltier than Site 2 (red line). The salinity clearly exhibits a semi-diurnal pattern as a result of the tidal fluctuation; however, the salinity fluctuation decreases gradually in the upstream direction as the tidal fluctuation diminishes. For example, the tidal range (the elevation difference between higher high and lower low water levels) for large tide at Site 1, Site 2 and Site 3 is approximately 3.7 m, versus 4.6 m at Sand Heads, which is located at the mouth of the Fraser River. A similar trend is exhibited in other parts of the irrigation period, but only a short period in September is shown in Figure 12 to emphasize the trend.

One intriguing aspect of these model results is that salinity is occasionally lower at Site 2 than at Site 1 as indicated in the Figure 12b, even though Site 1 is approximately 4 km upstream of Site 2. This might be caused by the complex processes that govern the hydrodynamics in the shallow area at Site 2, which was demonstrated by the surface salinity contour plot in Figure 6. The figure shows the presence of high salinity waters in the midst of fresher waters, probably due to a combination of local processes of selective withdrawal and upwelling in the river.

Salinity monitoring data collected at Site 1 and Site 2 supports the occurrence of this phenomenon in that surface water at Site 1 is occasionally saltier than Site 2. This complexity renders the salinity in the area difficult to predict, thereby posing challenges to the planning process.

4.2. Effects of Sea Level Rise and River Flow

4.2.1. Time Horizon: Present

Figure 13 shows the trend in surface salinity at Site 2 during the irrigation period from August to October in a normal flow year with 0-m SLR, which represents the present case. The mean monthly flow rate during the period decreases from 3200 m³/s on 1st August to 2150 m³/s on 1st September, to 1950 m³/s on 1st October, and finally to 1600 m³/s on 31st October (Figure 10). The irrigation threshold of 0.34 ppt is shown in Figure 13 as a black line, the modelled salinity as a green line.

The salinity trend shown in Figure 13 is the combined results of (1) the semi-diurnal tidal cycles in the SOG; (2) the alternation of neap tide, when the tidal forcing from the SOG is lower, and spring tide, when the tidal forcing is stronger; as well as (3) the change in river discharge over the period. The semi-diurnal tidal cycles explain the intra-day trend in salinity; the neap and springs tide cycles have led to alternation of high and low salinity in a weekly scale; and the decrease in river discharge between August and October has led to a general rise in the surface salinity at Site 2 and reduction of water availability for irrigation.

In August, the impact of the three aforementioned factors on surface salinity and, hence, water availability is minimal as salinity remains below the criterion value because the river discharge remains sufficiently high to position the salt wedge well downstream of the intake. However, in September, as the river discharge continues to trend downward, the salinity signal appears and becomes more pronounced with time, indicating the upstream advance and near-by presence of the salt wedge; the water availability for irrigation remains near 100% at the beginning of the month, but the withdrawal window gradually shrinks with increasing surface salinity in the latter parts of the month. In October, the salinity signal continues to grow stronger and the availability window for water withdrawal expectedly becomes even narrower and very limited, entering into the drier part of the Fraser River hydrograph.

It is important to note that the river flow did not decrease at a constant, fixed rate over the time period. A short-term increase in flow rate often occurs due to passage of weather systems or occurrence of rain-on-snow events, especially in the fall season. For example, a likely rain-on-snow event, a typical occurrence in the fall season, might have caused a temporary increase in the river flow rate in early-to-mid October (the green line in Figure 10a). This spike in the flow rate directly led to a decrease in salinity in the river as shown in Figure 13c.

4.2.2. Time Horizon: Year 2050 and Year 2100

Figure 14, Figure 15 and Figure 16 respectively illustrate the time series of salinity at Site 1, Site 2, and Site 3, for the 15-day period between 20th August and 3rd September, for all three time horizons (green line for the present time with 0-m SLR; black line for the year 2050 case with 1-m SLR, and red line for the year 2100 case with 2-m SLR; the irrigation threshold of 0.34 ppt for salinity is shown with a thin black line). The river flow type (dry, normal, and wet) corresponds to the time horizon of interest (Figure 10). The difference in salinity indicated in these figures is the result of the combined effects from SLR and changes in the Fraser River hydrograph.

These figures clearly show an increase in surface salinity at all three sites as a result of sea level rise. For the normal year and wet year cases, especially, an otherwise 100% withdrawal window for the 0-m SLR case could be reduced significantly and, on some days, that window could dwindle to 0% (Site 2 on 31 August—normal flow year, for example) with 1-m SLR and 2-m SLR. This is illustrated by the surface salinity with 1-m SLR remaining continuously above the salinity threshold, whereas the surface salinity for the 0-m SLR case remains below the threshold.

4.3. Effects of Channel Deepening

Figure 17 shows the time series of salinity at Site 2 in a normal flow year and 0-m SLR and with dredging depths that accommodate vessels with 11.5-m draft (green line) and 20-m draft (blue line) for the months of August, September, and October. The irrigation threshold for salinity of 0.34 ppt is again shown in the figure as a thin black line.

Not surprisingly, salinity increases as the channel is deepened, leading to much more frequent exceedances of the salinity threshold and a shorter time window for water withdrawal. The increase in water depth would likely result in a combination of reduced bottom friction acting on the salt wedge and increase in hydrostatic pressure force that drives the salt wedge upstream, thereby allowing the salt wedge to migrate further upstream.

5. Discussion

5.1. Effect of Sea Level Rise and River Flow

Table 4. Modelled time duration (no. of hours) below salinity criterion of 0.34 ppt at site 1, site 2, and site 3 with various sea level rises and types of flow year.

	Site 1 (h)			Site 2 (h)			Site 3 (h)
Sea Level Rise	Dry	Normal	Wet	Dry	Normal	Wet	Dry	Normal	Wet
0 m	7.9	24.0	24.0	3.6	23.9	24.0	3.1	23.0	24.0
1 m	0.0	3.8	23.6	0.0	0.8	21.2	0.0	0.4	18.1
2 m	0.0	0.1	8.3	0.0	0.0	3.9	0.0	0.0	1.7

below summarizes the effects of sea level rise (0-m SLR, 1-m SLR, and 2-m SLR) and river flow (dry, normal, and wet flow years) on water availability (time duration during which salinity is 0.34 ppt) for the existing case, the year 2050 case and the year 2100 case, respectively, for the same time period in August as presented in Figure 12, Figure 13 and Figure 14. Note the hours of availability presented in Table 4 are the average time duration per day for water withdrawal during the 15-day period plotted in the above figures.

As expected, the model results show that Site 1 will provide a wider withdrawal window than Site 2, which in turn has a wider widow than Site 3. At all three sites, the window becomes narrower with an increase in sea level and with a decrease in river flow. In a lower river discharge regime, the impact of sea level rise appears to be more significant than it does in the high river discharge regime. On the other hand, the influence of river discharge on withdrawal window decreases when the sea level is higher than it does when the sea level is lower.

It is difficult to draw definitive conclusions regarding which factor (sea level rise or river flow) dominates the salt wedge dynamics from this set of model results because, as presented in Section 3.4.1 and in Figure 10, each of the sea level rise cases incorporates different flow rates associated with the projected normal year, dry year, and wet year as predicted by MIROC for its own time horizon (present, year 2050 and year 2100); simply speaking, these model runs have two varying factors, the flow rate and sea level, and it is not immediately obvious which factor is more important in governing the salinity in the river. Investigation on the dominance of these factors has not been undertaken but can be achieved by isolating these two variables in the model and keeping one factor constant and the other variable. However, it might be difficult to consider one factor without changing the other, considering that the projected sea level rise and change in the Fraser River hydrograph are likely both results of climate change.

5.2. Effect of Channel Deepening

Table 5. Time duration below salinity criterion of 0.34 ppt at site 1 and site 2 with 0-m SLR with river channel accommodating vessels with 11.5-m and 20.0-m drafts.

	Site 1		Site 2
Month	11.5-m Draft (h/Day)	20.0-m Draft (h/Day)	11.5-m Draft (h/Day)	20.0-m Draft (h/Day)
August	24.0	19.8	23.9	15.0
September	22.0	6.1	16.9	1.1
October	15.0	2.3	9.9	0.5

below summarizes the availability for water withdrawal at Site 1 and Site 2, for the existing and a deepened channel for the 0-m SLR scenarios.

Site 1 and Site 2 both display a reduction in the withdrawal window as a result of channel deepening from 11.5 m to 20.0 m, especially in the latter parts of the irrigation period in October when the flow is lower compared to the earlier parts in August. It appears that channel deepening would not affect the withdrawal availability as much as sea level rise and change in river discharge would: the withdrawal availability with deeper channel remains substantial at both sites in August with 0-m SLR. However, as the flow rate continues to decrease with time between August and October, the influence from channel deepening appears to increase as the nearly completely open withdrawal window would close significantly.

5.3. Comparison with Past Studies

Based on their model results, Krvavica et al. [11] similarly concluded that sea-level and river discharge are the dominant factors controlling salt-wedge intrusion and the relative influence from sea level rise and channeling deepening increases when the river discharge is lower. However, the study by Krvavica et al. concerned the Rječina River that drains into the Adriatic Sea with small tidal range (<1 m) and concluded that the tidal dynamics of that river is not sufficient to cause significant mixing between the freshwater and saltwater masses in the estuary. Shaha et al. [13] found that, during the wet season when river discharge in the Sumjin River is relatively high, sea level plays a less dominant role in determining the salt wedge intrusion.

5.4. Future Studies and Research

There is a significant data gap that needs to be filled in order to fully understand the salt wedge dynamics and behaviour in response to both short-term and long-term changes in climate and river geometry. There is currently a lack of continual, long-term observations of salinity and temperature profiles at multiple locations along the river. When this data gap is addressed, one can delineate the behavior of the salt wedge in the Fraser River and its relationship with changes in environmental (climate change) and man-made forcing. The observations will provide scientists and engineers the necessary data to undertake a more comprehensive modelling study. The results from such a study will assist planners and stakeholders to plot an appropriate, adaptive course of action in a timely fashion to ensure a sustainable future for industries that rely on water drawn from the river and to maintain and enhance the current ecosystem.

Future research should also pursue the relative influence of channel deepening and sea level rise in the Fraser River on salinity distributions and withdrawal windows, which was not undertaken in this study.

Additionally, a Fraser River model with a higher resolution would have resolved better the finer-scale hydrodynamic processes and features that could potentially have immense effects on the salinity distribution and therefore the withdrawal availability for irrigation and other activities.

Finally, development of a local hydrological model that provides a relationship between precipitation and runoff to the Fraser River along the section between Hope and the upstream boundary of the 50-m Fraser River model would have significantly improved the quality of the hydrodynamic model.

6. Conclusions

From this modelling study, some conclusions can be drawn regarding the evolution of salinity in the Fraser River resulting from environmental changes (sea level rise and change in river flow) and direct man-made changes (channel deepening). The following summarizes these conclusions.

• Site 1, located 4 km upstream of Site 2, has consistently shown to have a significantly wider window for water availability in all cases and for all criteria compared to Site 2, the present intake location. This indicates that salinity generally decreases in the upstream direction. However, complex hydrodynamic processes would lead to exceptions to the trend.
• The temporal projection of the withdrawal window is that, even in a wet flow year, the sea level rise of 1 m and 2 m will lead to a large reduction (minimum 85% reduction) in and, in some scenarios, complete elimination of water availability.
• This study found that sea level rise and changes in river discharge appear to have a larger impact on the withdrawal availability than does channel deepening. In a low river discharge regime, the impact from sea level change is more significant than it is in the high river discharge regime. On the other hand, the influence of changes in river discharge on withdrawal availability decreases when the sea level is higher than it does when the sea level is lower. It is difficult, however, to draw definitive conclusions regarding which factor (sea level rise or river flow) dominates the salt wedge dynamics from this set of model results because, as presented in Section 3.4.1 and in Figure 10, each of the sea level rise cases incorporates different flow rates associated with the projected normal year, dry year, and wet year as predicted by MIROC for its own time horizon (present, year 2050 and year 2100); simply speaking, these model runs have two varying factors, the flow rate and sea level, and it is not immediately obvious which factor is more important in governing the salinity in the river. However, it is difficult to consider one factor without another, considering that the sea level rise and change in the Fraser River hydrograph are both results of climate change.
• Dredging the channel to accommodate vessels with a 20 m draft will affect the salinity at the intake and will shorten the withdrawal window. The effect of channel deepening becomes more pronounced in the low flow period. However, the degree of impact from dredging on the salt wedge and on withdrawal availability under other different circumstances (i.e., different sea level rise and different dredge depths) have not been investigated.