Three-dimensional hydrodynamic model of the Rance estuary (France) influenced by the world’s second largest tidal power plant

ABSTRACT The Rance estuary is a small steep-sided ria, located in the Brittany coast of northern France, with a maximum spring tidal range of 13.5m and an average river discharge of 7m3/s. Taking advantage of this significant tidal range, the Rance tidal power station (RTPS) was built in the 1960s as the world’s first and largest tidal power plant, with peak output capacity of 240 Megawatts. It is currently the second world’s largest tidal power installation after the Sihwa-Lake tidal barrage. The RTPS has two active parts: a barrage of 6 sluice gates and a structure of 24 turbines. Despite a well-known effect of the plant on damping estuarine water levels, little attention has been given to currents vertical distribution and the plant's impact on the dynamics of freshwater-saltwater interface. Therefore, a three-dimensional model of the Rance estuary was developed. Moreover, currents and salinity measurements were carried out to validate the numerical model. Simulated and measured currents showed that (i)the RTPS induces an acceleration of flood currents directly upstream of the sluice gates and (ii)ebb currents are strengthened by the narrowness of the Saint-Hubert-Port. Finally, salinity analyses assessed the dynamics of the freshwater-saltwater interface which is pushed further upstream during summer.


KEYWORDS
Numerical modelling; tidal power station; estuarine hydrodynamics; salinity

Introduction
Estuaries are transitional water bodies found in continuously changing dynamic systems, dominated by cyclically tidal forcing, inputs of freshwater, and other natural constraints. Hydrodynamic and morphological attributes observed in estuaries strongly depend on these driving forces, which in turn might be influenced by the combined effect of multiple human pressures.
The presence of a dam built across the estuary's mouth and/or a lock regulating the natural freshwater discharge can have considerable effects on the natural variability and equilibrium of these systems (Kim et al., 2017). In consequence, static features such as estuarine mean water level as well as dynamic processes such as tidal wave propagation and distortion are expected to be affected as a result of the human pressure on the semi-enclosed body of water (Angeloudis & Falconer, 2017).
Owing to the region's hyper-tidal regime and its narrow and confined characteristics, the Rance estuary, located in Brittany, France, hosts the Rance Tidal Power Station (RTPS). Opened in 1966, this largescale tidal power plant was for 45 years the largest marine energy facility in the world by its installed capacity, surpassed in 2011 by the South Korean Sihwa Lake Tidal Power station. Currently operated by Electricité de France, the dam is 750 m long, with a power plant portion equal to 332.5 m housing 24 bulb turbines (Charlier, 2007). By reaching a peak output of 240 megawatts (MW) and average 57 MW, the RTPS supplies 0.12% of the power demand of France, which is the averaged electric energy consumption equivalent to a medium-size city like Rennes, with a population of approx. 360,000 inhabitants in the urban area.
To better understand the influence of the RTPS on the hydrodynamics and its consequences on sediment dynamics and morphological changes, Rtimi et al. (2021) studied flow patterns and tidal asymmetry for diverse scenarios involving past and present bed elevations, the presence/absence of the dam and the existence of a lock at the estuary's upper limit. The two-dimensional (2D), depth-averaged numerical model implemented by Rtimi et al. (2021) provided detailed analysis on the flood-dominant characteristics of the estuary's dynamics for both spring and neap tides. Nevertheless, the understanding of the flow through the water body, and the effect of the functioning of the sluice gates and turbines of the RTPS, did not account for three-dimensional (3D) processes.
The comprehension of three-dimensional hydrodynamic processes is important for a better prediction of the water circulation, salinity distribution and velocity gradients within the body of water. These hydrodynamic processes influence, for example, the concentration of suspended particulate matter and therefore, the estuarine turbidity maxima (Amoudry et al., 2014;Burchard et al., 2018;Hesse et al., 2019). Due to the natural and anthropogenic-influenced variability of the Rance estuary, these processes are crucial for accurate sediment transport and morphodynamic predictions and may have profound implications on the estuarine ecosystem (Bonnot-Courtois et al., 2002;Kirby & Retière, 2009).
In this work, a 3D Reynolds-averaged Navier-Stokes model of the Rance estuary is implemented. This model accounts for the influence of the RTPS and solves the flow velocity and water salinity fields. It is calibrated and validated with a set of continuous recording probes and high-resolution field survey campaigns over a 10 days period. The 3D flow field and salinity distribution are validated with datasets collected from Acoustic Doppler Current Profiler, STPS and SAMBAT probes, respectively. Field surveys and probes data, in combination with 3D model results, allowed to better understand (a) the complex flow behaviour at specific zones of the Rance estuary; and (b) the mixing processes and salinity distribution within this particular body of water regulated by the presence of the RTPS at the estuary's mouth and a lock at its upper limit.

Site description
The Rance estuary is a small steep-sided ria (Evans & Prego, 2003) located on the Brittany coast of northern France. It is 20 km long and has variable irregular width, with a maximum value of 2 km at 10 km upstream from the mouth. Its maximum spring tidal range reaches 13.5 m at the mouth. Taking advantage of this hyper-tidal regime, the first ever tidal power station in the world was built at 3 km from the estuary mouth ( Figure 1(b)). The plant has been in operation and managed by Electricité de France (EDF) since 1966 and is currently the second largest operational tidal power station in the world (Pelc & Fujita, 2002). The head of the estuary is located at the Chatelier Lock, which delimitates the upstream limit of tide propagation. The Rance River drains a small catchment area, with an average river discharge of 7 m 3 /s, low water flow rate of 0.5 m 3 /s and a decennial flood of 80 m 3 /s.

Currents
Surveys of current velocities were conducted at four locations along the estuary, on 14, 15 and 21 October 2020. At each location, horizontal current velocities were collected through cross-sectional transects with a vessel-towed 1200 kHz RDI Workhorse Acoustic Doppler Current Profiler (ADCP) in 0.25 m bins. At locations of transects 1, 2 and 3 ( Figure 1), currents were collected throughout a semidiurnal tidal cycle (12.42 h), with a total of 30 full transects being completed at each location. At transect 4, close to the TPS, measurements were limited to a period of 3 hours during flood tidal stage and 10 full transects were completed. For each completed transect, measurements were then averaged over 20 ensembles to reduce noise.

Salinity
Salinity data used in this study was obtained through continuous measurements collected at six stations along the estuary from June 2019 to December 2020 ( Figure 1). At stations P1, A2 and A3, a STPS probe by NKE was deployed at fixed height above the bottom. At P1, the probe was fixed at a pier close to the Chatelier Lock. At A2 and A3, probes were fixed at 30 cm above the bed on intertidal mudflats. At these three stations, the probes emerged during part of the tidal cycle. As for S1, S2 and S3, salinity was measured by a SAMBAT probe (by NKE) permanently submerged at 80 cm below the water surface and attached to a floating buoy. At the six stations, salinity was recorded continuously with a time step of 10 min. The probes were deployed continuously for more than 18 months. Raw data was carefully verified in order to select valid shorter time series, discarding malfunctions, battery failure or device drifts.

Mathematical and numerical model
Flow and the salinity fields are computed by solving the conservation equation of fluid mass, the Reynoldsaveraged Navier-Stokes (RANS) momentum equations and a passive scalar equation accounting for the salinity, where the Boussinesq approximation is used. The solution variables are the horizontal velocity u = (u,v), the vertical velocity w, the free-surface elevation H = h + b, with h the water-depth and b the bottom elevation above datum, the pressure p, and the water salinity S. These variables are functions of the Cartesian coordinate space x = (x,y,z) and time t. Initial and boundary conditions as well as the representation of the RTPS are imposed as described by Rtimi et al. (2021). To close the hydrodynamic system, a standard k-eps turbulence model for both horizontal and vertical directions was used. Energy losses due to friction are parameterised with the Strickler roughness closure relationship. The Coriolis force is considered in the model, with the Coriolis parameter equal to 1.0909970 × 10 −4 s −1 . The water density varies with salinity and its average value is equal to 999.972 kg/m 3 . For the active scalar accounting for the salinity, the horizontal and vertical diffusion coefficients are set to 2.70 × 10 −5 m 2 /s. The numerical solution of the governing equations is performed with the module Telemac-3D (T3D), belonging to the TELEMAC-MASCARET modelling system (www.opentelemac.org) (Dutta et al., 2017). Telemac-3D's basic algorithm splits the solution into three steps, namely (i) calculation of the advected velocity components, (ii) update of the velocity field by incorporating the diffusion and source terms, and (iii) the water-depth calculation from the vertical integration of the continuity equation and the momentum equations by only including the pressure-continuity terms. An algorithm that deals with wetting and drying processes at the intertidal zones is incorporated in the numerical solution procedure. The computational domain is discretized with an unstructured triangular mesh in the horizontal directions, extruded along the vertical direction to form triangular prisms, covering the volume delimited by the bottom and the water surface elevation. The horizontal mesh of the Rance estuary and its adjacent coastal sea area consists of 389,766 elements and 199,625 nodes. The mesh size varies from 50 m in the offshore zone, down to 5 m in estuary's basin near the RTPS and between Mordreuc and the Chatelier Lock (Figure 1). The mesh size is 20 m elsewhere. Along the vertical direction, 8 layers were specified at each horizontal mesh, i.e. each prismatic element's thickness depends on the bottom elevation. The nodes of the horizontal mesh contain the bathymetric information, where the digital elevation dataset is identical to the high-resolution data (Lidar and field survey records) used by the model of Rtimi et al. (2021).

Tide propagation
Prior to the analysis of currents distribution and salinity variations, the tide propagation is first validated through field data from three tidal gauges for periods spanning neap and spring tides: July 20-30, 2020 and October 13-23, 2020. The numerical model reproduces well the tidal distortion produced by the RTPS as well as the tidal propagation along the estuary (Figure 2) with Root Mean Square Errors (RMSEs) below 10 cm. Moreover, high frequency oscillations of ~10 cm are observed during high tide, and are amplified at both Chatelier Lock and upstream of the RTPS stations. This seiche-like phenomenon was also noticed by Rtimi et al. (2021) over a fortnight period in 2019.

Currents
This section focuses on the period of October 13-23, 2020 when ADCP measurements were performed. Magnitude of maximum flood currents (respectively ebb currents) were first compared between numerical model and ADCP data for transects 1, 2 and 3 (see, Figure 1 for transects locations) given the availability of measurements within a complete tidal cycle. Afterwards, the magnitude of maximum current was analysed locally at transect 4 (upstream of the RTPS) during the opening period of the sluices.
The numerical model reproduces satisfactorily the vertical and cross-section distributions of maximum flood current along the estuary (Figure 3). Indeed, the outflow at Mordreuc (transect 1) is accelerated in the main channel (200-300 m from left bank) due to the morphology of the transect, then it is decelerated by the centre shoal, which extends between 50 and 150 m from left bank, to finally reach magnitudes below 0.25 m/s in the secondary channel (0-50 m from left bank). However, the maximum current magnitude simulated at transect 1 is slightly greater than the one measured. This could be explained by the lack of reliable measurements at the exact time as the numer- the ADCP maximum flood currents even though it is delayed by ~30 minutes compared to the actual moment of maximum flood current. The alongchannel current is sped up by the narrowness of the cross-section at Saint Hubert Port (transect 2), with a peak current of 1 m/s against 0.8 m/s at transect 1. Interestingly, the complex pattern of flow velocity measured at transect 3 is well reproduced by the numerical model (Figure 3(c, f)). The strongest flood currents of ~1.2 m/s are observed in the eastern side of the main channel (200-650 m from left bank), then the flow is slowed down towards the western channel. It is noteworthy the logarithmic-like vertical velocity distribution in the eastern region (450-550 m) with a maximum (respectively minimum) current magnitude at the surface (respectively near the bed). Finally, maximum flood currents are increased going from the estuary upstream (transect 1) towards its downstream (transect 3).
Similar analysis was performed on maximum ebb currents at the same transects ( Figure 4). The numerical model reproduces well the distribution of ebb velocity magnitude along the estuary. Indeed, an interesting pattern of vertical velocity distribution is observed at Saint Hubert Port (transect 2), with flow increase from the bottom to the surface. Contrastingly to maximum flood currents noticed at the most downstream transect (Saint Suliac transect 3), maximum ebb currents are observed at Saint Hubert Port transect Figure 4(b, e) where the morphology of the estuary impacts considerably currents variation. Furthermore, for all transects, flood currents maxima (at the most 1.2 m/s) are stronger than ebb currents maxima (at the most 0.7 m/s). Indeed, this could be explained by the tidal asymmetry in the Rance estuary assessed by Rtimi et al. (2021), in addition to the low river discharge at the uppermost limit (Chatelier Lock).
The depth-averaged numerical model implemented in Rtimi et al. (2021) assessed that the location of maximum flood current across the section upstream of the RTPS is highly sensitive to the operating modes of the dam. Indeed, ADCP measurements support this modification induced by the plant and is also well reproduced by the three-dimensional numerical model ( Figure 5(a,   Figure 3. Comparison of maximum flood current magnitude between (a; b; c) ADCP measurements and (d; e; f) numerical model T3D at transects 1, 2 and 3 respectively (see Figure 1 for transects locations).

Figure 4.
Comparison of maximum ebb current magnitude between (a; b; c) ADCP measurements and (d; e; f) numerical model T3D at transects 1, 2 and 3 respectively (see Figure 1 for transects locations). b)). The strongest along-channel velocities (at the most 2.5 m/s) occurred in the eastern region, between 350 m and 550 m from the left bank, corresponding to the zone directly upstream of the sluices. In the western region (0-70 m) corresponding to a part of the turbines upstream, weaker currents (at the most 0.75 m/s) are noticed due to the operating mode of the turbines. Finally, the centre shoal (200-300 m) and the presence of a dike (Figure 1) significantly reduce the currents to magnitudes below 0.25 m/s. Analysing simultaneously maximum flood currents distributions along the estuary (Figures 3 and 5), the RTPS generates a jet-like flow during the opening of the sluices with peak outflow of 2.5 m/s, which will be reduced by ~40% 5 km upstream.

Salinity
The validation of simulated salinity was carried out using ten days period's dataset (July 20-30, 2020) provided by six salinity gauges located along the   estuary (black squares in Figure 1). The examination of numerical and measured salinities revealed a good agreement between the model and the observations (Figure 6). At stations S1, S2 and S3, salinity is oscillating between 33 and 35 g/l as the probes are constantly submerged. Therefore, the lower estuary is mainly dominated by salt water (35 g/l). Contrastingly, stations P1, A2 and A3 revealed a temporal variation between 0 during low tide and 34 g/l otherwise. This pattern is explained by the emergence of the probes at low water, noticed from water levels at each station (grey lines in Figure 6(d, e and f)). However, the model tends to overestimate the salinity at the uppermost limit (P1) between the 7th and the 10th day (Figure 6(f)). Actually, the salinity in this region is highly sensitive to the river discharge at the Chatelier Lock. Since this latter is not measured, it made it difficult to predict the salinity decrease at P1. Moreover, this variation could also be explained by the presence of the rain which is not taken into account in the simulations. Nevertheless, the model reproduces satisfactorily the overall salinity dynamics in the Rance estuary.
To analyse the freshwater-saltwater interface in the Rance estuary, simulated salinity during July 2020 is compared with historical measurements from Bonnot-Courtois et al. (2002) before and after the construction of the plant (Figure 7). Before the RTPS, the freshwatersaltwater interface could be flushed up to Saint Suliac, but after the RTPS operation, saltwater penetrates further upstream the estuary, up to Saint Hubert Port. In July 2020, horizontal salinity distribution (Figure 7 (c)) revealed that the freshwater-saltwater interface is indeed pushed by the barrage towards the lower estuary. However, the mixing area between fresh and salt waters does not exceed 2 km downstream of the Chatelier Lock, while it could cover, during winter, a larger portion, up to 5 km from the Chatelier Lock (Bonnot-Courtois et al., 2002). Examination of the vertical salinity distribution along the channel (see Figure 1 for exact location) assessed that the Rance estuary could be divided into three regions (Figure 7(d)): (i) 0-200 m from the Chatelier Lock where salinity is below 20 g/l; (ii) 200-1800 m from the Chatelier Lock where the channel becomes deeper, salinity is increased and varies between 25 and 30 g/l, and (iii) above 1800 m from the Chatelier Lock where salinity is around 35 g/l.

Discussion
The present work aims to validate the 3D structure of currents and salinity simulated along the Rance estuary influenced by the world's second largest tidal power plant.
Comparisons between simulated and measured velocities revealed the ability of the numerical model to reproduce the 3D flow structure from the lower estuary (transect 1) to the upstream of the RTPS (transect 4). Flood currents are mainly modulated by the operation modes of the RTPS, specifically the opening of sluice gates. Locally, the flow is substantially accelerated near the eastern channel upstream of the plant, with a homogeneous distribution along the water column. These findings confirm the predictions of the 2D depth averaged model of Rtimi et al. (2021). Furthermore, the morphology of the estuary significantly influences the propagation of ebb currents, with strongest inflows observed in the narrowing at Saint Hubert Port. Moreover, currents do not seem to considerably vary along the water column except at Saint Hubert Port (transect 2) and Saint Suliac (transect 3) during floods.
After over 50 years of service, the Rance hydroelectric power station induces modification in the natural hydraulic regime which led to a different distribution of salinity compared to that before the dam's construction. Dynamics of freshwater-saltwater interface is sensitive to seasonal river discharge variation, hydraulic flushes occurred at the Chatelier Lock and plants operating cycles. Seasonal variability of this interface ranges between 2 km in summer (5 km in winter (Bonnot-Courtois et al., 2002)) from the Chatelier Lock.
In addition to tidal asymmetry depicted from flood and ebb currents (Friedrichs & Aubrey, 1988), density variation induced by salinity gradient and cross-channel flow could influence the estuarine turbidity maximum as well as sediments transport and deposition in the basin (Amoudry et al., 2014;Burchard et al., 2018;Hesse et al., 2019;Sottolichio et al., 2000). Further analysis of the hydrodynamics coupled with sediment transport processes is therefore necessary to investigate the impact of the RTPS on the sedimentation in the estuary.

Conclusions
Analyses of the three-dimensional flow field and salinity distribution, along the estuary over a 10 days period, collected by acoustic Doppler current profiler and continuous measurements probes respectively, have enabled the validation of a 3D numerical model accounting for the presence of the RTPS. The model correctly reproduces ebb and flood currents distribution along the estuary. Moreover, salinity patterns are also satisfactorily captured, with a freshwatersaltwater interface expected 2 km downstream of the Chatelier Lock during summer. Furthermore, these results highlight the added-value of 3D numerical modelling of the Rance estuarine system with respect to the 2D depth-averaged model (Rtimi et al., 2021). The RTPS modifies tidal patterns and the flow structure within the water body, which in turn might influence the sediment transport and estuarine morphological evolution. This subject has not been studied extensively and warrants further investigation.