Proof-of-Concept of a Quasi-2D Water-Quality Modelling Approach to Simulate Transverse Mixing in Rivers

: A quasi-two-dimensional (quasi-2D) modelling approach is introduced to mimic transverse mixing of an inﬂow into a river from one of its banks, either an industrial outfall or a tributary. The concentrations of determinands in the inﬂow vary greatly from those in the river, leading to very long mixing lengths in the river downstream of the inﬂow location. Ideally, a two-dimensional (2D) model would be used on a small scale to capture the mixing of the two ﬂow streams. However, for large-scale applications of several hundreds of kilometres of river length, such an approach demands too many computational resources and too much computational time, especially if the application will at some point require ensemble input from climate-change scenario data. However, a one-dimensional (1D) model with variables varying in the longitudinal ﬂow direction but averaged across the cross-sections is too simple of an approach to capture the lateral mixing between different ﬂow streams within the river. Hence, a quasi-2D method is proposed in which a simpliﬁed 1D solver is still applied but the discretisation of the model setup can be carried out in such a way as to enable a 2D representation of the model domain. The quasi-2D model setup also allows secondary channels and side lakes in ﬂoodplains to be incorporated into the discretisation. To show proof-of-concept, the approach has been tested on a stretch of the lower Athabasca River in Canada ﬂowing through the oil sands region between Fort McMurray and Fort MacKay. A dye tracer and suspended sediments are the constituents modelled in this test case.


Background
Much of the process water from oil sands surface mining operations is recycled and managed in tailing ponds. However, the capacity for storage is approaching unmanageable and unsustainable levels; hence, some release of treated process water into the Athabasca River is anticipated as early as 2025 [1]. The release of treated process water may pose a risk to aquatic species and to humans who harvest and consume these species, in particular fish. Therefore, effective models to describe the transport and fate of oil sands related substances are required [2]. These substances can be transported downstream and deposited in Lake Athabasca and the Peace-Athabasca-Delta (PAD); in addition, secondary channels and lakes within the floodplain along the lower river reach may also trap released sediment and associated constituents. An important objective of this research is to determine the fate of such effluent within these features using computer modelling.

Quasi-Two-Dimensional Modelling
This paper describes a unique quasi-two-dimensional representation of river hydraulics that is particularly suited to the application of the WASP model to accurately represent the configuration of the lower Athabasca River with a high level of computational efficiency. In the current study, we introduce a novel approach to modelling transverse mixing in a river with secondary channels and side lakes to study the water quality along the area of the Athabasca River with extensive oil sands development.
In order to maintain short computational times, a one-dimensional (1D) modelling approach is necessary. However, the transverse mixing along the river requires modelling with at least a two-dimensional (2D) representation, especially since the lengths of complete mixing along this river are relatively long (>100 km). Hence, the use of a quasi-twodimensional (quasi-2D) approach is proposed, in which flow is simulated in 1D, but in such a way to allow a 2D discretisation of the domain.
Quasi-2D water-quality modelling has been carried out in the past for other applications. For off-channel storage facilities (polders) along the Elbe River in Germany, Lindenschmidt et al. [5] modelled dissolved oxygen and nutrient dynamics for various flow regimes (low, medium, high (flood) flows). Deposition of sediments and heavy metals in the off-channel storage basins was captured using the quasi-2D method [6][7][8]. The quasi-2D approach has also been used to capture flows between main river channels and their floodplains, in particular through dike breaches [9] and capping flood peaks using side-channel storage [10,11]. Flow between the Mekong River and its delta [12] and between the Po River and a portion of its floodplain [13] were simulated using quasi-2D models. Sediment transport was included in a quasi-2D model of the Rhine River's main channel and floodplain [14]. In this study, we extend the quasi-2D approach to model transverse mixing in rivers.

Objectives
The objectives of this study are to:

1.
Develop a quasi-2D modelling approach to characterise transverse mixing of two water streams of different sediment concentrations.

2.
Assess the role of drawing on the output of additional models to complement the implementation of models with different complexity and spatial scale.
As proof of concept of the approach, the sediment transport along the lower Athabasca River was quantified. The model domain includes a secondary channel and side lake connected hydraulically to the river's main stem.

Site Description
The Athabasca River (see Figure 1) is a northern river (ice-covered in winter) in western Canada that originates on the eastern slopes of the Rocky Mountains and flows approximately 1500 km in a north-easterly direction to empty into Lake Athabasca. It is the longest unregulated river in Alberta. For the last 80 km, the river flows through the Athabasca Delta. The Athabasca Delta is part of the larger Peace-Athabasca Delta (PAD), which is the largest inland delta in North America. The PAD is an important ecosystem for many aquatic and terrestrial plant and animal species, and the area has been named both a Ramsar and United Nations Educational, Scientific and Cultural Organization (UNESCO) World Heritage site. From a human perspective, the PAD is an important social, cultural and economic entity for many of the Aboriginal communities in the area. Many conservation efforts have been carried out to maintain the ecological integrity of the aquatic and terrestrial systems.
Referring to Figure 1, the average discharge at the Water Survey of Canada (WSC) gauge at Embarras (gauge 07DD001-Athabasca River at Embarras Airport), which is approximately 80 km upstream of the PAD, is around 850 m 3 /s. The total catchment area draining at the same gauge is approximately 156,000 km 2 . Much of the basin's land use consists of forests (89%) with some agricultural lands (8%). The basin is sparsely populated (≈200,000 inhabitants), with the population concentrated in a few urban centres, which release effluent from municipal wastewater treatment plants into the river.
Water 2021, 13, x FOR PEER REVIEW 3 of 20 3. As proof of concept of the approach, the sediment transport along the lower Athabasca River was quantified. The model domain includes a secondary channel and side lake connected hydraulically to the river's main stem.

Site Description
The Athabasca River (see Figure 1) is a northern river (ice-covered in winter) in western Canada that originates on the eastern slopes of the Rocky Mountains and flows approximately 1500 km in a north-easterly direction to empty into Lake Athabasca. It is the longest unregulated river in Alberta. For the last 80 km, the river flows through the Athabasca Delta. The Athabasca Delta is part of the larger Peace-Athabasca Delta (PAD), which is the largest inland delta in North America. The PAD is an important ecosystem for many aquatic and terrestrial plant and animal species, and the area has been named both a Ramsar and United Nations Educational, Scientific and Cultural Organization (UNESCO) World Heritage site. From a human perspective, the PAD is an important social, cultural and economic entity for many of the Aboriginal communities in the area. Many conservation efforts have been carried out to maintain the ecological integrity of the aquatic and terrestrial systems.
Referring to Figure 1, the average discharge at the Water Survey of Canada (WSC) gauge at Embarras (gauge 07DD001-Athabasca River at Embarras Airport), which is approximately 80 km upstream of the PAD, is around 850 m 3 /s. The total catchment area draining at the same gauge is approximately 156,000 km 2 . Much of the basin's land use consists of forests (89%) with some agricultural lands (8%). The basin is sparsely populated (≈200,000 inhabitants), with the population concentrated in a few urban centres, which release effluent from municipal wastewater treatment plants into the river.   The oil sand region covers a large area of Northern Alberta and has one of the richest deposits of petroleum in the world. For a 65 km stretch of the lower Athabasca River, between Fort McMurray and downstream of Fort MacKay (see Figure 1), the river flows through the active oil sand surface mines. Much of the oil is extracted through open-pit mining and processed using water from the Athabasca River. Almost 5% of the river's average flow has been allocated for anthropogenic usage. Although approximately half of that amount is allocated for surface and in situ mining activities, less than 1% is used [15].
The lower Athabasca River contains many islands, secondary channels, wetlands, and floodplain lakes. Many of these secondary channels freeze to the bottom or have low oxygen levels during winter [16]. Sediment deposition areas occur downstream of tributaries (confluence bars), in mid-channel bars, in secondary channels and in side lake and wetland areas along the river. During high flow events, substantial amounts of oil sands material can be transported from tributaries into the river to be deposited in these depositional areas. Over time, these areas of high oil sands material are mixed and diluted with water and sediment coming from upstream.
To show proof-of-concept, the domain modelled in this study is only a short stretch of 15 km. It is part of a larger system (>200 km) that is to be modelled in the future, with additional outfalls and sensitive depositional areas (e.g., secondary channels and side lakes). The model domain is indicated in Figure 2 and was chosen based on the location of proposed outfalls from oilsands operations. One particular side lake of interest in this study is Shipyard Lake, which is located near the Suncor oil facilities on the east side of the Athabasca River and is within the model domain. The lake, shown in Figure 3, has a surface area of 21.3 ha and a maximum depth of 2.2 m [17]. It is essentially a large wetland area that is flooded by higher flows from the Athabasca River [17]. Numerous small creeks feed the lake and water from the Athabasca River flows into Shipyard Lake over a low levee between the watercourse and water body. Flow begins to overflow the levee at an Athabasca River flow of 2800 m 3 /s [18], which corresponds to a water surface elevation of 237.7 m a.s.l. just upstream of the lake [18]. This corresponds to a water level elevation of approximately 240.1 m a.s.l. at the Athabasca River gauge downstream of Fort McMurray [19], the location of which is shown in Figure 2. The oil sand region covers a large area of Northern Alberta and has one of the richest deposits of petroleum in the world. For a 65 km stretch of the lower Athabasca River, between Fort McMurray and downstream of Fort MacKay (see Figure 1), the river flows through the active oil sand surface mines. Much of the oil is extracted through open-pit mining and processed using water from the Athabasca River. Almost 5% of the river's average flow has been allocated for anthropogenic usage. Although approximately half of that amount is allocated for surface and in situ mining activities, less than 1% is used [15].
The lower Athabasca River contains many islands, secondary channels, wetlands, and floodplain lakes. Many of these secondary channels freeze to the bottom or have low oxygen levels during winter [16]. Sediment deposition areas occur downstream of tributaries (confluence bars), in mid-channel bars, in secondary channels and in side lake and wetland areas along the river. During high flow events, substantial amounts of oil sands material can be transported from tributaries into the river to be deposited in these depositional areas. Over time, these areas of high oil sands material are mixed and diluted with water and sediment coming from upstream.
To show proof-of-concept, the domain modelled in this study is only a short stretch of 15 km. It is part of a larger system (>200 km) that is to be modelled in the future, with additional outfalls and sensitive depositional areas (e.g., secondary channels and side lakes). The model domain is indicated in Figure 2 and was chosen based on the location of proposed outfalls from oilsands operations. One particular side lake of interest in this study is Shipyard Lake, which is located near the Suncor oil facilities on the east side of the Athabasca River and is within the model domain. The lake, shown in Figure 3, has a surface area of 21.3 ha and a maximum depth of 2.2 m [17]. It is essentially a large wetland area that is flooded by higher flows from the Athabasca River [17]. Numerous small creeks feed the lake and water from the Athabasca River flows into Shipyard Lake over a low levee between the watercourse and water body. Flow begins to overflow the levee at an Athabasca River flow of 2800 m 3 /s [18], which corresponds to a water surface elevation of 237.7 m a.s.l. just upstream of the lake [18]. This corresponds to a water level elevation of approximately 240.1 m a.s.l. at the Athabascsa River gauge downstream of Fort McMurray [19], the location of which is shown in Figure 2.    A secondary channel is also included in the model domain, also shown in Figure 3. Indicated in the figure is the location of a cross-section of the Athabasca River main stem and the secondary channel, which is shown in Figure 4. Some water does flow into the secondary channel from the Athabasca River at a flow of 600 m 3 /s.  A secondary channel is also included in the model domain, also shown in Figure 3. Indicated in the figure is the location of a cross-section of the Athabasca River main stem and the secondary channel, which is shown in Figure 4. Some water does flow into the secondary channel from the Athabasca River at a flow of 600 m 3 /s. A secondary channel is also included in the model domain, also shown in Figure 3. Indicated in the figure is the location of a cross-section of the Athabasca River main stem and the secondary channel, which is shown in Figure 4. Some water does flow into the secondary channel from the Athabasca River at a flow of 600 m 3 /s.  Total suspended solids (TSS) data from four stations along the Athabasca River were used in the analysis: Empirical relationships between TSS concentrations and river discharge are provided in Appendix A, from which equations were established and are shown in Figure 6 for a lower flow regime. Sediment transport was simulated only for the 2800 m 3 /s flow since, in this scenario, Athabasca River water flows into both the secondary channel and Shipyard Lake. A total flow of 2800 m 3 /s at the Athabasca River gauge below Fort McMurray corresponds to a flow of 250 m 3 /s for the Clearwater River (from Figure   The hydraulic model HEC-RAS [20], developed by the U.S Army Corps of Engineers, was used to generate the segmentation for the WASP water quality model. The water level elevations of the HEC-RAS model at flows 600 and 2800 m 3 /s, which formed the geometrical basis for the WASP segmentation, were verified with flow simulations using the hydraulic model ONE-D, a hydrodynamic model that uses a finite difference implicit scheme for the solution ( [21,22] both referenced in [23]). Both mass and momentum of sub-critical fluid motion are conserved for open-channel flow. The model is robust and well-tested with many consultancy applications (e.g., [24,25]). Bathymetry data were obtained from Alberta Environment and Parks (AEP) [26].  The hydraulic model HEC-RAS [20], developed by the U.S Army Corps of Engineers, was used to generate the segmentation for the WASP water quality model. The water level elevations of the HEC-RAS model at flows 600 and 2800 m 3 /s, which formed the geometrical basis for the WASP segmentation, were verified with flow simulations using the hydraulic model ONE-D, a hydrodynamic model that uses a finite difference implicit scheme for the solution ( [21,22] both referenced in [23]). Both mass and momentum of sub-critical fluid motion are conserved for open-channel flow. The model is robust and welltested with many consultancy applications (e.g., [24,25]). Bathymetry data were obtained from Alberta Environment and Parks (AEP) [26].

Water-Quality Modelling Setup
The Water Quality Analysis Simulation Program (WASP 8.32) was used to model the water quality of the river-lake system. WASP is a general dynamic model that applies a segmentation network to solve the conservation of momentum, energy, and mass equations and simulate contaminant and sediment transport. WASP was initially developed in the 1980s and has undergone many upgrades since then. WASP is a widely used model, especially in North America, for addressing various environmental and water quality concerns. The WASP segmentation can be structured in 1D for streams, or 2D for rivers with branching channels or 3D for lakes.
The WASP stream transport module, TOXI, is able to calculate the flow of water, sediment and dissolved constituents through branched and ponded segments and is coupled with flow routing for free flow streams, ponded segments, and backwater reaches. The kinematic wave was used for the 1D flow routing, which is based on solutions of onedimensional continuity equations and a simplified form of the momentum equation that considers effects of gravity and friction, and calculates variations in velocities, widths and depths throughout the network. Flow through the ponded segments is calculated based on a sharp-crested weir equation, which calculates outflow based on water elevation above the weir crest, using kinematic wave flow for a river and ponded weir overflow for a lake. For the current study, boundary and initial conditions were defined based on water quality data of the closest station to the study area by spatial linear interpolation.
Input data for the model includes channel geometry, flow routing, boundary conditions, environmental time functions, loads and initial conditions for segments. Geometry for the WASP segments (average depths, widths and slopes) was extracted from the fiftyeight HEC-RAS cross-sections, which were spaced at a distance of 493 m apart. Segments are a series of discretized components used to estimate the model state variables at every

Water-Quality Modelling Setup
The Water Quality Analysis Simulation Program (WASP 8.32) was used to model the water quality of the river-lake system. WASP is a general dynamic model that applies a segmentation network to solve the conservation of momentum, energy, and mass equations and simulate contaminant and sediment transport. WASP was initially developed in the 1980s and has undergone many upgrades since then. WASP is a widely used model, especially in North America, for addressing various environmental and water quality concerns. The WASP segmentation can be structured in 1D for streams, or 2D for rivers with branching channels or 3D for lakes.
The WASP stream transport module, TOXI, is able to calculate the flow of water, sediment and dissolved constituents through branched and ponded segments and is coupled with flow routing for free flow streams, ponded segments, and backwater reaches. The kinematic wave was used for the 1D flow routing, which is based on solutions of one-dimensional continuity equations and a simplified form of the momentum equation that considers effects of gravity and friction, and calculates variations in velocities, widths and depths throughout the network. Flow through the ponded segments is calculated based on a sharp-crested weir equation, which calculates outflow based on water elevation above the weir crest, using kinematic wave flow for a river and ponded weir overflow for a lake. For the current study, boundary and initial conditions were defined based on water quality data of the closest station to the study area by spatial linear interpolation.
Input data for the model includes channel geometry, flow routing, boundary conditions, environmental time functions, loads and initial conditions for segments. Geometry for the WASP segments (average depths, widths and slopes) was extracted from the fiftyeight HEC-RAS cross-sections, which were spaced at a distance of 493 m apart. Segments are a series of discretized components used to estimate the model state variables at every time step. The sediment transport model is an independent set of routines with associated parameters for the model.
The domain was modelled with a network of river segments, a secondary channel and a side lake. Water quality modelling of such a complex system requires special consideration of the varying hydrodynamics within the system. The primary segments of the system are 493 m long; as mentioned before, right stream segments accommodate 20% of the width of cross-sections produced by HEC-RAS for the 600 m 3 /s scenario and 10% of the width for the 2800 m 3 /s flow scenario.
Bathymetric information for Shipyard Lake was obtained from a bathymetry survey carried out by Golder [17]. The lake was divided into two segments: • a deep part with a surface area of 14 ha and an average depth of 1.5 m and • a shallow part with a surface area of 7.3 ha and an average depth of 0.75 m.
As mentioned above, the mixing length for the stream is long, and the geographical complexity of the study area cannot be represented with a 1D modelling approach alone. On the other hand, a 2D approach would require higher data computational resources for such a long stream. Quasi-2D models can be a reasonable alternative by using the simplicity of a 1D solution but still capturing 2D discretisation of the substance transport. The approach is purely a balance of water and mass between each adjoining segments. This includes the exchange of water and mass between each adjacent-lying, left and right stream segments. The exchange is a small percentage of the water flowing into each of those segments from their upstream segments. The exchange alternates from right-to-left and left-to-right stream segments along the model domain. The flow exchange was calibrated in such a way that the constituents within the two streams mixed together more and more as we moved further downstream along the domain. The flow transfers between adjacent segments were calibrated so that the concentrations matched the observed dispersion pattern calculated with the Athabasca River Model (ARM). Concentration gradients were generally higher between the left and right stream segments at the upstream boundary leading to more rapid mixing in this area; as we move further downstream, the mixing between the adjacent segments of each stream led to a decrease in their concentration gradients and a decrease in the mixing rate.
In this study, results from one model, ARM, were used to calibrate the mixing simulated with WASP. However, field data sampled not only longitudinally along the flow direction of the river, but also transversely across the river can be used to calibrate the model, an example of which is provided in a follow-up study to simulate the transport of vanadium [27].
The stretches between every two cross-sections of the HEC-RAS model were divided into two stream tubes. The flow from the Clearwater River controlled the segment on the right side (right stream), and the main Athabasca River system controlled the segments on the left side (left stream), as shown in Figure 8. Flow for the secondary channel was supplied through the right side stream. It was estimated that 2% of flow into Shipyard Lake was required to maintain water volume continuity. The flow through Shipyard Lake was anticipated to flow back into the right side stream just downstream of the lake. A schematic of the segmentation is shown in Figure 9. In order to characterize erosion and sedimentation, benthic segments were added under each surface water segment. Solids can be separated into clay, silt and sand, which can refine simulation results of solids. The user has four choices for simulating solids in WASP: solid flow fields, descriptive, the mechanistic Van Rijn equation, and the mechanistic Robert equation. Currently, WASP has two kinetic modules for modelling sediment transport. For the current study, the advanced toxicant module was used because it is considered to be the best module to simulate oil-sands associated substances.

WASP Mixing Calibration (Q = 600 m 3 /s)
Calibration of the mixing of the WASP model was done using simulation results from the Athabasca River Model (ARM) [16,17]. ARM is a vertically averaged, two-dimensional mixing model used to predict substance concentrations varying across the width and length of the lower Athabasca River for various flow conditions. The model's calculations In order to characterize erosion and sedimentation, benthic segments were added under each surface water segment. Solids can be separated into clay, silt and sand, which can refine simulation results of solids. The user has four choices for simulating solids in WASP: solid flow fields, descriptive, the mechanistic Van Rijn equation, and the mechanistic Robert equation. Currently, WASP has two kinetic modules for modelling sediment transport. For the current study, the advanced toxicant module was used because it is considered to be the best module to simulate oil-sands associated substances.

WASP Mixing Calibration (Q = 600 m 3 /s)
Calibration of the mixing of the WASP model was done using simulation results from the Athabasca River Model (ARM) [16,17]. ARM is a vertically averaged, two-dimensional mixing model used to predict substance concentrations varying across the width and length of the lower Athabasca River for various flow conditions. The model's calculations are analytical/explicit in nature, so there is no discretization other than differentiating between different reaches with different mixing and hydraulic parameters. ARM applies analytical solutions for river dispersion equations under steady-state conditions. Mid-field mixing is represented by dispersion equations described by Fischer et al. [28]. The same representation is used for passive ambient mixing in the Cornel Mixing Zone Expert System (CORMIX) [29] except that, for ARM, transverse mixing coefficients have been calibrated based on a tracer-dye study completed in 2003 [30,31] and validated using a range of data sources [32]. The calibration of ARM is described in [30,31].
A hypothetical conservative tracer was chosen to calibrate the mixing path of the WASP model to the ARM predictions. A load of 5670 kg/day, which led to the chemical concentration of 1 mg/l in the most-upstream segment (Segment 1) was added. The flow system included two upstream boundary conditions: the Athabasca River water for the left stream boundary and the Clearwater River water for the right stream boundary.
As expected, the conservative tracer was transported to the left stream when a flow was designated between them. Segments divided by islands (segment combinations 16-17, 18-19 and 20-21) had no flow between them. WASP simulated the concentrations of the tracer in each segment. These values were compared with values predicted by the ARM model and shown in Figure 10.
Water 2021, 13, x FOR PEER REVIEW 12 of 20 system included two upstream boundary conditions: the Athabasca River water for the left stream boundary and the Clearwater River water for the right stream boundary. As expected, the conservative tracer was transported to the left stream when a flow was designated between them. Segments divided by islands (segment combinations 16-17, 18-19 and 20-21) had no flow between them. WASP simulated the concentrations of the tracer in each segment. These values were compared with values predicted by the ARM model and shown in Figure 10. Good agreement for both left and right streams was observed. As can be seen, both models predicted the same trend for both streams: a reduction in chemical concentrations in the right stream and an increase in chemical concentrations in the left stream. Contour plots of the tracer concentrations modelled in WASP and ARM are shown in Figure 11. Good agreement for both left and right streams was observed. As can be seen, both models predicted the same trend for both streams: a reduction in chemical concentrations in the right stream and an increase in chemical concentrations in the left stream. Contour plots of the tracer concentrations modelled in WASP and ARM are shown in Figure 11. Good agreement for both left and right streams was observed. As can be seen, both models predicted the same trend for both streams: a reduction in chemical concentrations in the right stream and an increase in chemical concentrations in the left stream. Contour plots of the tracer concentrations modelled in WASP and ARM are shown in Figure 11.

WASP Mixing Validation (Q = 2800 m 3 /s)
A simulation with a higher flow of 2800 m 3 /s was carried out for which both Shipyard Lake and the secondary channel became flooded. The results of the two models for the two streams are shown in Figure 12. Again, mixing is evident since the tracer concentrations are increasing in the right stream and decreasing in the left stream, both in the flow direction. Mixing is more rapid for the 2800 m 3 /s scenario since the tracer concentrations

WASP Mixing Validation (Q = 2800 m 3 /s)
A simulation with a higher flow of 2800 m 3 /s was carried out for which both Shipyard Lake and the secondary channel became flooded. The results of the two models for the two streams are shown in Figure 12. Again, mixing is evident since the tracer concentrations are increasing in the right stream and decreasing in the left stream, both in the flow direction. Mixing is more rapid for the 2800 m 3 /s scenario since the tracer concentrations drop to 40% of the original concentration at the end of the 15 km stretch, whereas the concentrations only drop to 55% in the 600 m 3 /s scenario ( Figure 10). Contours of the tracer concentrations simulated in ARM and WASP are provided for the 2800 m 3 /s flow case in Figure 13. The same trend for the mixing was observed for both flow rates studied herein. For the higher flow rate, faster mixing was observed, which can be explained by higher turbulence and more chaotic eddies which follow observations by Dimotakis [33] and Sreenivasan [34].  Figure 13. The same trend for the mixing was observed for both flow rates studied herein. For the higher flow rate, faster mixing was observed, which can be explained by higher turbulence and more chaotic eddies which follow observations by Dimotakis [33] and Sreenivasan [34].

Sediment Transport Modelling
Suspended solids and benthic sediments are key components of the water quality of rivers and lakes. When the flow velocity decreases, as is the case when, for example, a river flows into a lake or reservoir, suspended sediments will settle out. Key constituents in oil sands process water, such as selenium and polycyclic aromatic hydrocarbons, tend to sorb to fine sediments. For the current study, separate size classes were represented using the descriptive solids transport option in WASP. For the descriptive solids transport, constant settling and resuspension velocities were defined for each segment.

Sediment Transport Modelling
Suspended solids and benthic sediments are key components of the water quality of rivers and lakes. When the flow velocity decreases, as is the case when, for example, a river flows into a lake or reservoir, suspended sediments will settle out. Key constituents in oil sands process water, such as selenium and polycyclic aromatic hydrocarbons, tend to sorb to fine sediments. For the current study, separate size classes were represented using the descriptive solids transport option in WASP. For the descriptive solids transport, constant settling and resuspension velocities were defined for each segment. Settling velocity of silt and clay were calculated from Stokes Law, and the settling velocity of sand was estimated based on the equation from Ferguson and Church [35].
Size classes for suspended solids were estimated based on suspended sediment samples collected from the Lower Athabasca River using continuous flow centrifugation in 2012 and 2013 as part of the Joint Oil Sands Monitoring Program [1].
Measured TSS data were fitted versus the flow rate for the sediment sampling stations. Power functions were used to fit the data yielding acceptable r 2 values (r 2 > 0.6) for all correlations. This provided equations to estimate TSS based on the given flow rates. For the sediment transport simulation, we used the flow scenario of 2800 m 3 /s because at that flow, both the secondary channel and Shipyard Lake are flooded by Athabasca River water. A flow of 2800 m 3 /s corresponds to a flow of 250 m 3 /s for the Clearwater River (from Figure 5), and a flow of 2550 m 3 /s along the upper Athabasca River reach, which represents a flow ratio of approximately 10-90%, respectively. Although the model domain is situated downstream of the Athabasca/Clearwater river confluence, for the purpose of showing proof-of-concept in this study, the 10/90 percentage ratio was maintained for the right/left streams.
After obtaining TSS values for the flow scenario for the water quality stations, an interpolation was done to obtain initial sediment concentrations for each segment along the study area (see Figure 14). For segments on the right stream (segments with odd numbers), linear interpolation between Clearwater River at Draper and Athabasca/Firebag confluence was carried out. For left stream segments (segments with even numbers), values were interpolated between the station upstream of the Athabasca/Clearwater confluence and Athabasca/Firebag confluence. After obtaining TSS amounts for each segment, the average ratio of long-term means was used to estimate the concentration of silt, clay and sand in each segment. Initial concentrations for benthic segments were estimated from available bed sediment data for the Lower Athabasca River collected as part of the Regional Aquatics Monitoring Program [36].
confluence was carried out. For left stream segments (segments with even numbers), ues were interpolated between the station upstream of the Athabasca/Clearwater con ence and Athabasca/Firebag confluence. After obtaining TSS amounts for each segm the average ratio of long-term means was used to estimate the concentration of silt, and sand in each segment. Initial concentrations for benthic segments were estim from available bed sediment data for the Lower Athabasca River collected as part o Regional Aquatics Monitoring Program [36]. The model was run until it reached an equilibrium state, after which no signifi changes in the sediment concentrations in the surface water segments were observe was assumed that an equilibrium state has been reached when there was no signifi change in sediment concentration in the water body for more than a simulation time o h. For the current study, clay, silt and sand were modelled, and the results of these both the left and right streams of the Athabasca River are shown in Figure 15.
The mechanics of solid transport in WASP are based on advection, and suspen sediments are only transported transversely when there is a flow routed between ments of the right and left streams. It is evident that the left stream segments have hi The model was run until it reached an equilibrium state, after which no significant changes in the sediment concentrations in the surface water segments were observed. It was assumed that an equilibrium state has been reached when there was no significant change in sediment concentration in the water body for more than a simulation time of 24 h. For the current study, clay, silt and sand were modelled, and the results of these for both the left and right streams of the Athabasca River are shown in Figure 15.
Water 2021, 13, x FOR PEER REVIEW 15 of 20 sediment concentrations than the right stream segments. In addition, sediment concentrations tend to drop in the left stream and increase in the right stream, both with flow direction, which is attributed to the transverse mixing, as is the case for the tracer simulations. For all segments, including left and right streams, it was observed that, for the water columns, the amount of sediments decreases in order from silt, to clay and then sand. Due to the higher settling velocity for sand, it can be expected that there would be more sand deposits in the river compared to silt and clay. Sand is typically transported as a bedload or during higher velocity flows due to its roughness and higher density [37][38][39]. Due to its low settling velocity, clay is the dominant sediment texture in the water column. At Shipyard Lake, a considerable deposition of sediment in the lake was simulated compared to the sediment in the right stream of the Athabasca River (see Figure 16). The flow through the lake is in the direction of: right stream → deep part → shallow part → right stream (further downstream). The higher sedimentation in the lake is due to the drop in flow velocity of the water entering the lake. Much of the sediment of all three textures is deposited in the deeper part of the lake with some additional deposition in the shallow lake area. The amount of solids increases in the surface benthic segments because settling dominates over resuspension in the lake. The mechanics of solid transport in WASP are based on advection, and suspended sediments are only transported transversely when there is a flow routed between segments of the right and left streams. It is evident that the left stream segments have higher sediment concentrations than the right stream segments. In addition, sediment concentrations tend to drop in the left stream and increase in the right stream, both with flow direction, which is attributed to the transverse mixing, as is the case for the tracer simulations.
For all segments, including left and right streams, it was observed that, for the water columns, the amount of sediments decreases in order from silt, to clay and then sand. Due to the higher settling velocity for sand, it can be expected that there would be more sand deposits in the river compared to silt and clay. Sand is typically transported as a bedload or during higher velocity flows due to its roughness and higher density [37][38][39]. Due to its low settling velocity, clay is the dominant sediment texture in the water column.
At Shipyard Lake, a considerable deposition of sediment in the lake was simulated compared to the sediment in the right stream of the Athabasca River (see Figure 16). The flow through the lake is in the direction of: right stream → deep part → shallow part → right stream (further downstream). The higher sedimentation in the lake is due to the drop in flow velocity of the water entering the lake. Much of the sediment of all three textures is deposited in the deeper part of the lake with some additional deposition in the shallow lake area. The amount of solids increases in the surface benthic segments because settling dominates over resuspension in the lake.

Summary and Conclusions
A novel methodology was proposed to simulate transverse mixing along a long stream of a river network with a secondary channel and side lake. A steady-state water quality model based on two flow scenarios (600 and 2800 m 3 /s) was developed to study transverse mixing and sediment deposition in the modelling domain.
The first objective of this study was to model transverse mixing in the river-lake system. Good agreement with the previously developed ARM model was observed. It was established that a quasi-2D approach is a reasonable and accurate means of modelling transverse mixing of the system with a minimal processing time compared to more complex 2D models. The approach is also well suited for scaling up to a larger river network in future work by extending the modelling domain from Fort McMurray to the Athabasca Delta.
The second objective was to quantify deposition of sediments along the system. As expected, much of the sediment load deposits are found within the secondary channel and the side lake, with less deposition in the main channel. It was observed that coarser materials tend to deposit in the deep part of the lake, where the Athabasca River water first enters the lake. Due to the longer settling time for finer sediments, they were the main sediments suspended in the main stem.
The last objective was to emphasize the importance of drawing on the output from additional models to complement the water-quality modelling exercise. To support the water-quality modelling with WASP, the hydrodynamics of the water-quality modelling domain was analyzed using HEC-RAS and ONE-D. The ONE-D model has an extended domain that incorporates recordings from gauge stations. The calibrated and validated ARM model was also useful to calibrate the transverse mixing within the WASP model. Should we scale up to extend the modelling stretch from Fort McMurray to the Athabasca Delta, additional models, in particular a hydrological model of the Athabasca River basin, The simulations with WASP were very fast, taking only about 30 s for the model to reach steady state. This is much faster than a two-dimensional approach that was carried out in which a "relatively coarser grid was used for all of the main chemical constituent runs due to their higher computational requirements" [40]. Our approach presented here would be an alternative to creating a second grid.

Summary and Conclusions
A novel methodology was proposed to simulate transverse mixing along a long stream of a river network with a secondary channel and side lake. A steady-state water quality model based on two flow scenarios (600 and 2800 m 3 /s) was developed to study transverse mixing and sediment deposition in the modelling domain.
The first objective of this study was to model transverse mixing in the river-lake system. Good agreement with the previously developed ARM model was observed. It was established that a quasi-2D approach is a reasonable and accurate means of modelling transverse mixing of the system with a minimal processing time compared to more complex 2D models. The approach is also well suited for scaling up to a larger river network in future work by extending the modelling domain from Fort McMurray to the Athabasca Delta.
The second objective was to quantify deposition of sediments along the system. As expected, much of the sediment load deposits are found within the secondary channel and the side lake, with less deposition in the main channel. It was observed that coarser materials tend to deposit in the deep part of the lake, where the Athabasca River water first enters the lake. Due to the longer settling time for finer sediments, they were the main sediments suspended in the main stem.
The last objective was to emphasize the importance of drawing on the output from additional models to complement the water-quality modelling exercise. To support the water-quality modelling with WASP, the hydrodynamics of the water-quality modelling domain was analyzed using HEC-RAS and ONE-D. The ONE-D model has an extended domain that incorporates recordings from gauge stations. The calibrated and validated ARM model was also useful to calibrate the transverse mixing within the WASP model. Should we scale up to extend the modelling stretch from Fort McMurray to the Athabasca Delta, additional models, in particular a hydrological model of the Athabasca River basin, should also be included. This will be essential to predict the state of the river's water quality under a changing climate scenario in the future. The impacts of an ice cover on the river's water quality should also be considered in future research.