Analytical and numerical study of the salinity intrusion in the Sebou river estuary (Morocco). Effect of the ‘Super Blood Moon’ (total lunar eclipse) of 2015

The longitudinal variation of salinity and the maximum salinity intrusion length in an alluvial estuary are important environmental concerns for policy makers and managers since they influence water quality, water utilization and 10 agricultural development in estuarine environments and the potential use of water resources in general. Super-Moon total lunar eclipse is rare event. According to NASA, they have only occurred five times in the 1900s-in 1910, 1928, 1946, 1964 and 1982. After the September 28 th , 2015 Total lunar eclipse, a Super-Blood-moon eclipse will not recur before October 8 th , 2033. In this paper, for the first time, the impact of the combination of a Super-Moon and a total lunar eclipse on the salinity intrusion along an estuary is studied. The 28 th September 2015 Super-Moon total lunar eclipse is focused by the study and 15 the Sebou river estuary (Morocco) is taking as an application area. The Sebou estuary is an area with high agricultural potential, is becoming one of the most important industrial zones in Morocco and it is experiencing a salt intrusion problem. Hydrodynamic equations for tidal wave propagation coupled with Savenije theory, and a numerical salinity transport model (HEC-RAS) are applied to study the impact of the Super-Moon total lunar eclipse on the salinity intrusion. Intensive salinity measurements during this extreme event were recorded along the Sebou estuary. Measurements showed a modification of 20 the shape of axial salinity profiles and a notable water elevation rise, compared with normal situations. The two optimization parameters (Van Der Burgh‟s and dispersion coefficients) of the analytical model are estimated based on the LevenbergMarquardt ́s algorithm (i.e. solving non-linear least-squares problems). The salinity transport model was calibrated and validated using field data. The results show that the two models described very well salt intrusion during the Super-Moon total lunar eclipse day. A good-fit between computed salinity and measurements is obtained, as verified by statistical 25 performance tests. These two models can give a rapid assessment of salinity distribution and consequently help to ensure the safety of water supply, even during such infrequent astronomical phenomenon. Key-words: 2015‟s Super-Moon total lunar eclipse; salinity intrusion; hybrid model; Savenije theory; HEC-RAS; Sebou estuary. 30


Introduction
A total lunar eclipse is one of Nature's loveliest celestial events (Espenak, 2000). Any total lunar eclipse offers unique and breathtaking wonders; the full moon becomes not just a red "blood moon" but also is a Harvest Moon and a "super moon".
Total lunar eclipses appear red and for this reason are often called blood moons (Fig.1). 5 Tidal motions arise as a response to forces associated with the interaction of the earth-moon-sun system, the effect of the moon being about twice that of the sun (Stronach, 1989). A total lunar eclipse occurs when the moon enters the cone of shadow cast by the sunlit earth (i.e. the sun and the moon are in line with the earth) (Fig.2). This shadow covers the entire moon and causes a total lunar eclipse. The high tide is at its highest point and the low tide at its lowest point, the gravitational pull on the ocean is then strong. Lunar eclipses are relatively infrequent events and on average only two of 10 them may be contemplated each year (Mallama, 1996 , Muñoz andPallé, 2011).
Estuaries form essential parts of the human-earth system (Savenije, 2015). As the connecting element between marine water and river, estuaries have properties of both: they contain both fresh and saline water; they experience tides, but also river floods; and they host both saline and fresh ecosystems (Savenije, 2015). The diverse estuarine environment plays an important role in the life cycle of many species but also serves as a site for many human activities. In short, estuaries are 15 important water bodies where many dynamic factors interact and unfold (Xu et al., 2015). For decades, explosive increases in industrial and agricultural productivity, as well as the growing population in estuary regions, have led to numerous environmental concerns (Mai et al., 2002). Salinity intrusion is an important phenomenon in an estuary, and can constitute a serious problem. It influences the water quality and threatens potential water resource use. Intake of fresh water for consumption, agricultural purposes or use by industries may take place within a region not far landward of the limit of salt 20 intrusion. To support policy and managerial decisions, a profound knowledge of processes associated with the salinity structure in estuaries is required (Kuijper and Van Rijn, 2011). Models have been widely used to research salinity intrusion.
One-dimensional mathematical models (analytical or numerical) can constitute the appropriate tools for quick-scan actions in a pre-phase of a project or for instructive purposes. Also, it is methodologically correct to start with the simplest description of the phenomena under study and to evaluate the limits of this approximation before investigating more complications . 5 Our previous studies on the Sebou estuary have shown that the one-dimensional (analytical or numerical) methods compute properly salt intrusion (Haddout et al., 2015a, Haddout et al., 2015b. In addition, Haddout et al., 2015b showed that salinity profiles of Sebou estuary show steep decrease. This characteristic is specific to narrow (recession shape) estuaries (i.e. having a near-prismatic shape and significant freshwater discharge). Such estuaries are called positive estuaries.
The aims of this paper it to investigate the applicability of these two methods (analytical or numerical) during the Super-10 moon day of 28 th September 2015 (total lunar eclipse). Measurements showed a modification of the form of the salinity profiles along the estuary and a notable water level increase, compared with normal situations studied in our earlier works. In addition, calculations during lunar eclipse using the coupled analytical hydrodynamic-salt intrusion model required the recalculation of the geometric parameters of the estuary i.e., cross-sectional areas A 0 , convergence lengths a.
Even in these extreme conditions, a good fit was obtained between the computed and observed salinity distribution for the 15 two models. These models constitutes powerful tools for evaluating salinity intrusion pattern in the Sebou river estuary, even during extreme events inducing sea level rise like lunar eclipse or climate change.

One-dimensional salt intrusion model 20
The analytical salinity intrusion model of Savenije (2005) has been adopted to predict the salinity distribution and salinity intrusion length in alluvial estuaries. This method is fully analytical, although it makes use of certain assumptions, the most important being: the exponential shape of the estuary, the longitudinal variation of the dispersion according to Van der Burgh (1972), and the predictive equations for the boundary condition and the Van der Burgh coefficient. The equations are based on the 1-D cross-sectionally averaged, and tidally averaged, steady-state salt balance equation, in which the advective 25 salt transport is caused by the downstream fresh water discharge, counteracted by the landward dispersive salt transport induced by the different mixing processes (Cai et al., 2015b). In a convergent estuary, where the main geometric parameters (cross-sectional area A , width B and depth h ) can be described by exponential functions (See equations (2.38)-(2.40) in Savenije, 2012).
In a steady-state situation, the partial temporal derivative in the salt balance equation is zero (Gisen et al., 2015b). 30 Considering constant fresh water discharge Q f and tidally averaged cross-sectional area A, the salt balance equations for High Water Slack (HWS), Low Water Slack (LWS) and Tidal Average (TA) situation can be rearranged as: Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2016-276, 2016 Manuscript under review for journal Hydrol. Earth Syst. Sci. Published: 13 June 2016 c Author(s) 2016. CC-BY 3.0 License.
If one assumes that D is constant along the length of an estuary (i.e., D=D 0 , where D 0 is the dispersion coefficient at the estuary mouth), simple analytical solutions for axial distribution of salinity can be derived. 5 Since dispersion depends both on the river discharge and the salinity distribution itself, the constant dispersion is not a correct solution. An efficient and accurate approach to simulate the longitudinal variation of dispersion is presented by Van der Burgh (1972), and also adopted by Savenije (1986Savenije ( , 1989Savenije ( , 1993aSavenije ( , 2005Savenije ( , 2012: Which, using equation (1), can be demonstrated to be the same as the following equation (Savenije, 2005(Savenije, , 2012: where S 0 (g/l) is the boundary salinity at the estuary mouth, D 0 (m 2 /s) is the longitudinal dispersion at the estuary mouth for HWS, LWS or TA condition, K is the Van der Burgh"s dimensionless coefficient, which has a value between 0 and 1. If K 15 = 0, equation (3) reduces to the case with constant dispersion D=D 0 . For the case where K =1, we see that the curves D/D 0 and S/S 0 coincide. If K is small, then tide-driven mixing is dominant near the toe of the intrusion curve; if K approaches unity, then gravitational mixing is dominant (Savenije, 2006;Shaha and Cho, 2011).

5
Making use of the dimensionless parameters (Cai et al., 2015a), equation (5) can be scaled as: Where: S * is dimensionless salinity that is normalized by the salinity at the estuary mouth,  is the estuary shape number representing the convergence of an estuary, D * is the dimensionless dispersion at downstream boundary condition and x * is the dimensionless longitudinal coordinate that is scaled by the frictionless wave length in prismatic channels (Cai et al., 2015a).
This one-dimensional steady advection-diffusion model has been applied to describe the salinity distribution along numerous 10 well-mixed and partially-mixed estuaries (i.e. the adopted salt intrusion model assumes a partially to well mixed situation, which under low flow is the dominant process) (Nguyen and Savenije, 2006) for HWS, LWS or TA condition.
The salt intrusion length * L , defined as the distance from the estuary mouth to the location with fresh water salinity (assumed to be <1g/l isohaline) for HWS, LWS or TA condition, can be determined by setting * S =0 in equation (6): 15 On other hand, in most estuaries, there usually exists an inflection point near the mouth, where the geometry changes (e.g., Gisen et al., 2015b). This inflection point is associated with the transition of a wave-dominated regime to a tide-dominated regime (Gisen et al., 2015b). Making use of this phenomenon, Gisen et al. (2015b) recently expanded the underlying data base and reanalyzed these equations, resulting in: 20 2 0.21 0.57 1 (Revised) 6 Where: f B is the river regime width (typical width in the river dominated region), and C is the Chezy roughness coefficient.
The symbols 1 H (m); 1 B (m); 1 v (m/s); 1 b (m); 1 h (m); 1 D (m 2 /s); 1 E (m) represent the tidal range, stream width, velocity amplitude, width convergence length, depth, dispersion coefficient, and tidal excursion at the inflection point, respectively. If 5 there is no inflection point near the estuary mouth, then these parameters refer to the situation at the mouth itself.
The roughness coefficient can be estimated by  (-) which is defined as the ratio of potential energy of the buoyant fresh water to the kinetic energy of the tide (Fischer et al., 1979) is given by: Where  (kg/m 3 ) is the water density,   is the density difference of ocean and river water over the salt intrusion length (in estuaries, the ratio /   is about 0.025).

Analytical hybrid (hydrodynamic) model 15
Since 1960-s there exists a long tradition of one-dimensional analytical solutions for tidal dynamics in estuaries (e.g., Dronkers, 1964;Ippen, 1966;Prandle and Rahman, 1980;Leblond, 1978;Godin, 1985Godin, , 1999Jay, 1991;Friedrichs and Aubrey, 1994;Lanzoni and Seminara, 1998;Kukulka and Jay, 2003;Horrevoets et al., 2004;Jay et al., 2011;Cai et al., 2012a). These analytical solutions usually made assumptions to simplify or linearize the non-linear set of equations . Of these, most authors used perturbation analysis, where scaled equations are simplified by neglecting higher 20 order terms, generally neglecting the advective acceleration term and linearizing the friction term, higher-order terms, whereas Savenije (2005) uses a simple harmonic solution without simplifying the equations (Cai et al., 2013). Others used a regression model to determine the relationship between river discharge and tide. Exceptions are the approaches by Horrevoets et al. (2004) and Cai et al. (2012b), who provided analytical solutions accounting for river discharge, based on the envelope method originally developed by Savenije (1998) (Cai et., 2013). Recently, Cai et al. (2012a) proposed a new 25 analytical framework for understanding the tidal damping in estuaries. They concluded that the main differences between the examined models (e.g., Savenije et al., 2008;Toffolon and Savenije, 2011;Van Rijn, 2011) lies in the treatment of the friction term in the momentum equation. Furthermore, Cai et al. (2012a) presented a new "hybrid" expression for tidal damping as a weighted average of the linearized and fully non-linear friction term (Cai et al., 2013). Additionally, Cai et al. (2014b) included for the first time the effect of river discharge in a hybrid model that performs better . 30 Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2016-276, 2016 Manuscript under review for journal Hydrol. Earth Syst. Sci. It can be demonstrated that tidal hydrodynamics is controlled by three dimensionless parameters that depend on localized geometry and external forcing (e.g., Toffolon et al., 2006;Savenije et al., 2008), i.e.: ζ the dimensionless tidal amplitude (indicating the seaward boundary condition), γ the estuary shape number (representing the effect of cross-sectional area convergence) and χ the friction number (describing the role of the frictional dissipation). These parameters are defined in 5 Table 1, where η is the tidal amplitude and s K is the Manning-Strickler friction coefficient.
Note that the friction number reflects the non-linear effect of the varying depth (Savenije, 2012). The tidal hydrodynamics analytical solution can be obtained by solving a set of four analytical equations, i.e., the phase lag equation, the scaling equation, the damping equation and the celerity equation (Cai et al., 2013). In Table 2, we present these equations for the 10 general case as well as the special case of the ideal estuary (δ = 0).

Coupled model for salt intrusion
Since tidal dynamics in convergent alluvial estuaries can be reproduced reasonably well by one-dimensional analytical solutions, in principle the output of such model can be used to predict the longitudinal tidal excursion * E (Cai et al., 2015a): 15 (i.e. * E is the dimensionless tidal excursion scaled by the frictionless tidal wave length), defined as: Where: v is the velocity amplitude and ω is the tidal frequency, 0 c is the classical wave celerity of a frictionless progressive wave defined as: In which g is the acceleration due to gravity and s r the storage width ratio (e.g., Savenije et al., 2008).
On other hand, Van der Burgh (1972) assumed that the salinity curves for the HWS and LWS situations can be obtained by applying a horizontal translation over half the tidal excursion in the landward and seaward direction from TA situation and subsequently was demonstrated by Savenije (1986Savenije ( , 1989Savenije ( , 2005Savenije ( , 2012. Thus Eq. (6) can be used to describe the two 25 situations of: Here the asterisk denotes a dimensionless variable. 8 The proposed analytical model by Cai et al., (2012a) for tidal hydrodynamics can be used to predict a variable velocity amplitude v (and hence tidal excursion E * ) for given tidal amplitude at the seaward boundary, estuary shape and friction. 5

Hydrodynamic model
Salinity distribution is influenced by the hydrodynamic regime, which in turn depends highly on the river estuary morphology. In the hydrodynamics module, HEC-RAS solves the following one-dimensional equations of continuity and 10 momentum, known as the Saint-Venant equations (Brunner, 2010): where Q is the discharge (m 3 /s), A is the cross sectional area (m 2 ), x is the distance along the channel (m), t is the time (s), ql 15 is the lateral inflow per unit length (m 2 /s), g is the acceleration due to gravity (m/s 2 ), h is the flow depth (m), 0  (-) is the bottom slope, n is the Manning"s roughness coefficient (m -1/3 /s) and R is the hydrodynamic radius (m).
Manning coefficient used in the momentum equation is evaluated initially by the empirical formula equation (16) proposed by Cowan (1956) and Chow (1973): where n 0 is a basic, n is the value for a straight, uniform, smooth channel, n 1 is the adjustment for the effect of surface irregularity; n 2 is the adjustment for the effect of variation in shape and size of the channel cross section; n 3 is the adjustment for obstruction; n4 is the adjustment for vegetation; and m 5 is a correction factor for meandering channels.
The equations (14) and (15) are solved using the well known four-point implicit box finite difference scheme (Brunner, 25 2010).

9
This numerical scheme has been shown to be completely non-dissipative but marginally stable when run in a semi-implicit form, which corresponds to weighting factor (θ) of 0.5 for the unsteady solution. This value represents a half weighting explicit to the previous time step"s known solution, and a half weighting implicit to the current time step"s unknown 5 solution. However, practically speaking, due to its marginal stability for the semi-implicit formulation, a θ weighting factor of 0.6 or more is necessary, since the scheme is diffusive only at values of θ greater than 0.5. In HEC-RAS, the default value of θ is 1. However, the user can specify any value between 0.6 to1 (Billah et al., 2015).

Mass transport model 10
In the advection-dispersion module, the basic equation is the one-dimensional advection-dispersion one of a conservative constituent (Brunner, 2010): Where C is the salinity concentration (g/l), A is the cross-sectional area of the river (m 2 ), Q is the discharge (m 3 /s) and D x is 15 the longitudinal dispersion coefficient (m 2 /s). This module requires output from the hydrodynamics module in terms of discharge, water level, cross-sectional area and hydraulic radius. The advection-dispersion equation is solved using the ULTIMATE QUICKEST explicit upwind scheme (Brunner, 2010). The resultant finite difference solution for equation (17) is as follows: 11 ..
where C n+1 is the concentration of a constituent at present time step (g/m 3 ), C n is the concentration of a constituent at previous Inputs of the transport model are initial and boundary salinity concentrations and the dispersion coefficient (parameter D in equation (17)). The dispersion coefficient was estimated using Kashefipour and Falconer (2002) formula: Where U * is the shear velocity,U is the cross-sectional average velocity, h is the depth of flow.

Overview of the Sebou estuary
The Sebou is the largest Moroccan river, draining approximately 40.000 km 2 , stretching about 614 km from its source in the Middle Atlas Mountains to the Atlantic Ocean, which represents 6% of Morocco's total land area (Fig. 3). Kenitra harbor, about 17 km from the ocean, has commercial traffic, while Mehdia harbor at only 2 km from the mouth is busy with fishing 10 activities. The flow regime at the level of the Sebou estuary is marked by considerable seasonal and inter-annual variations.
It is under the influence of the tide regime and under the control of many dams . During low flow periods, hydrodynamic regime is controlled by the Lalla Aïcha dam situated 62 km upstream. This dam has been constructed to preserve water for agricultural pumping stations and to avoid that salty waters rise towards these stations. Before the dam construction, excessive salinity reached up to 85 km upstream (Combe, 1969). The tidal height varies from 0.9 to 3.10 m depending on the condition of the tide and the average flow is about 200 m 3 /s at the river mouth Combe (1966). In addition, the tide near the estuary mouth is mainly semi-diurnal with a 44820-s tidal cycle (Haddout et al., 2014) and is a meso/micro-tidal estuary. The topography of the Sebou estuary is presented in Fig. 5. Fig. 5 shows the shapes of cross-sectional area, channel width, and depth during neap-spring and total lunar eclipse tides. The 20 cross-sectional area and width are plotted directly from bathymetric data, and the water depth represents the ratio of these two geometric values.
The entire estuary can be divided into two reaches with different convergence length for cross-sectional area, width, and depth. The first inflection point is located at x =5km (between the river mouth (Mehdia) and Kenitra). The second inflection 25 point is located at x=35km (between Oulad Salma and M"Rabeh), where the convergence length switches and the estuary becomes more riverine. Geometrical characteristics along the Sebou estuary are summarized in Table 3, with the crosssectional area; width and depth are well described by exponential functions (Fig.5).

11
The convergence length is shorter in the seaward reaches (x = 0-5km and x=5-35km) where the tidal influence is dominant, regardless of the river flow upstream Lalla Aïcha dam. In the landward reach (x=35-62km), the influence of river flow becomes important and the convergence length much greater. In addition, we see a gradual depth increase along the estuary 5 in the seaward reaches. Upstream the second inflection point, the depth gradually increases, while the cross-sectional area and width remains roughly constant.

Field data interpretation
In this study, five locations (Mehdia, Oulad Berjel, Kenitra, Oulad Salma, and M'Rabeh) along the Sebou estuary are chosen for salt measurements (Fig. 6) during 12 hours of 28 th September 2015 (Super-moon day). A W-P600 conductivity meter is 10 used in each location. The salinity can be expressed in parts per thousand (ppt or g/l) and the average value in the ocean is 35g/l. The period of measurements included LWS, TA and HWS situations (HWS and LWS correspond to periods when water velocity changes its directions and becomes nearly null).
The maximum and minimum salinity curves at HWS and LWS were thus observed, representing the envelopes of the salinity 15 variation during tidal cycle. Fig. 6 shows vertically averaged salinity concentration measurements conducted at the five locations, during 12-h with 12-min interval. A maximum salinity concentration is recorded during high tides with 35.5 g/l at the river mouth, 32.7g/l at Oulad Berjel, 30g/l at Kenitra, 1.2g/l at Oulad Salma and 0.9g/l at M"Rabeh. A minimum salinity concentration is recorded during low tides with 17.5g/l at the river mouth and less than 1g/l in the other four locations.

20
On other hand, figure 7 shows vertically averaged salinity concentration measurements at normal situation during (spring tide) in three locations (Oulad Berjel, Kenitra, Oulad Salma). A maximum salinity concentration is recorded during high tides with 28.5g/l at Oulad Berjel, 19.8g/l at Kenitra and 0.9g/l at Oulad Salma. A minimum salinity concentration is recorded during low tides with 5g/l at Oulad Berjel and less than 3g/l in the other two locations.

25
The mixing mechanisms in estuaries are guided by tidal dynamics, the dispersion mechanisms and the amount of fresh water discharge from the river estuary. Dispersion includes longitudinal mixing, that takes place by mass travelling in streamlines at different velocities (Nylén and Ramel, 2012)

12
The water level drops in the low tide and then rises and peaks with the high tide to 4 m at total eclipse day (less than 3.2 m in normal situation). This indicates the total lunar eclipse influence on the estuary hydrodynamic regime.

Results and analysis 5
In this paper, the analytical (tidal propagation and salt intrusion) and numerical models introduced in the previous sections are applied to the Sebou estuary to evaluate the total lunar eclipse effect on salinity distribution.

Analytical hydrodynamic hybrid and salt intrusion method application
To turn the steady state model into a predictive model, semi-empirical relations (equation 6) are required that relate the two 10 optimization parameters K and D 1 to hydrodynamic and geometrical bulk parameters. These bulk parameters are dimensionless numbers composed of geometrical (a, b, A 0 , B 0 , h 0 ), hydrological (Q f ) and hydraulic (H, E, C, υ, η) parameters that influence the process of mixing and advection. An analytical hydrodynamic model for tidal wave propagation is used to reproduce the main tidal dynamics along the estuary axis and subsequently for predicting the main parameters (Van der Burgh"s coefficient K, dispersion coefficient D 1 and tidal excursion E) that influence the salt intrusion process during total 15 lunar eclipse day of the Sebou estuary.
Field measurements of salt intrusion along the Sebou estuary axis, which were conducted at the 07 th -18 th -27 th May 2015 (covering a spring-neap cycle) and at the 28 th total lunar eclipse day, are used to test the predictive method. Each tidal cycle consists of two HWS and two LWS salinity observations, which corresponds to the tidal wave periods. The hybrid 20 (hydrodynamic) model presented in Sect.1 was calibrated during (07 th -18 th -27 th May 2015) and at total lunar eclipse day. The calibrated parameters, the Manning-Strickler friction coefficient s K and the storage width ratio r s , are presented in Table 4. It is worth noting that the storage width ratio r s is different in the seaward reaches for various tidal situations. It is important to point out that the model uses a variable depth in order to account for along-channel variations of the estuarine sections. Fig.9 shows the longitudinal computations of tidal amplitude and travelling time along the Sebou estuary for the selected periods 25 (representing the spring tide, the moderate tide, the neap tide, and total lunar eclipse day). The agreement between analytically computed and observed tidal amplitude and travelling time for HW and LW is good. On other hand, these results shows, the difference between eclipse day and spring-neap situation is very remarkable. Additionally, a velocity amplitude and dumping number at eclipse day are show in the same figure, which suggests that the hybrid model proposed by Cai et al., (2012aCai et al., ( , 2014a can well reproduce the tidal dynamics and velocity amplitude with a significant range of dam discharges. Based on the hydrodynamic parameters (e.g., tidal amplitude and velocity amplitude) from hybrid (hydrodynamic) model, it is possible to estimate the main parameters that determine the salinity intrusion from predictive Eqs. (8)-(9) at the inflection point. In most estuaries, there usually exists an inflection point near the mouth, where the geometry changes (e.g., Gisen et al., 2015b). On other hand, the Van Der Burgh"s (K) and dispersion (D 1 ) coefficients are initially estimated by equations 8 5 and 9. However, due to the large uncertainty of these predictive equations, the K and D 1 estimations should be refined on the basis of salinity measurements. The optimization process of K and D 1 has been carried out using the Levenberg-Marquardt non-linear minimization method (Marquardt, 1963). In this method, the following objective function  is minimized during the parameters optimization process. Because the predictive uncertainty of these equations, an optimization is applied to K and D 1 so as to obtain the best fit 15 between computed salt-intrusion curve and in-situ data (salinity). Values of the optimized parameters are summarized in Table 5 (at HWS), where the dispersion coefficient at the estuary mouth D 0 can be obtained by substituting D 1 ; x 1 and K into Eq. (4). It can be shown that the estimated values of the parameters D 1 and K are very close to their reference value  Additionally, salt-intrusion exhibit three distinct tendencies: a dome shape at HWS and TA, a recession and bell shapes at LWS (see Appendix A). These three salinity shapes were observed during the eclipse day.
On other hand, in a positive estuary, the salinity gradient is always negative due to the decreasing salinity in landward 5 direction as it is the case for the Sebou river estuary (Fig. 11a). The dome shaped intrusion curves have a negative curvature in the seaward part of the estuary. As water level increases, the position of the maximum salinity gradient moves towards to the estuary downstream, after which it has a dome shape with monotonous increasing of the salinity gradient . In addition,

The predictive model compared to other methods
This salinity intrusion model has been applied in different estuaries all over the world and by several methods (e.g. Van den Burgh, 1972;Rigter, 1973;Fischer, 1974;Van Os and Abraham, 1990and Savenije, (1993a, 2005

see Appendix B). 15
In Fig. 13, the observed vs computed intrusion lengths are plotted together with data specific to Sebou estuary. As one can see in the same figure, the result of Savenije model is considered most accurate.

Hydrodynamic model
The hydrodynamic regime was first studied and modeled in HEC-RAS. Outputs from the hydrodynamic model (velocity and water level evolution) were used in the salt transport study. The final resolution of hydrodynamic model equations (14 and 15) requires spatial discretization of the study area. The river reach (62 km) was discretized into 203 grids with a length varying between 58 and 996 m (Haddout et al., 2015a). Data on cross sectional (bathymetry) areas from the ANP (National 25 Agency of Ports) and other sources were used. The upstream boundary (at the Lalla Aïcha dam) was given values of discharge as a function of time from 27-09-2015 to 29-09-2015. Also, the downstream boundary (at the mouth) was given values of the water level as a function of time.
The factor n 0 equation (16) is evaluated from granulometric measurements that were carried out from upstream to downstream in the studied reach. The others coefficients were evaluated from observations of the river in aerial photos, from 30 the cross-sectional areas and available photos, and from field visits.

Salt transport model
The salinity model has been calibrated by systematically adjusting the values of a selected system parameter to achieve an acceptable match between the measured salinity and the corresponding values computed by the one dimensional advection dispersion model. The calibration parameter is the dispersion coefficient in the river which is adjusted by trial and error. Fig.  15 15 shows the comparison of the observed and simulated salinity at three locations (Oulad Berejel, Kenitra, Oulad salma) during the total eclipse day.
In the calibration procedure, the dispersion coefficient was modified to the same degree along the studied reach because we assumed that the sources of errors involved in its evaluation are identical for all the grids. The calibrated value of coefficient "D" ranges from 500 to 900 m 2 /s along the river. The results show that the simulated salinity concentration follows a similar 20 trend to the observed data.

Models Performance Verification
The statistical indicators used for evaluating the performance of the numerical and analytical models are: root mean squared error (RMSE); mean absolute error (ABSERR); the Nash Sutcliffe modelling efficiency index (EF); the goodness-of-fit (R 2 ) and the % of deviation from observed streamflow (PBIAS). The statistical parameters were defined as follows (Moriasi et 25 al., 2007;Stehr et al., 2008;Conversa et al., 2015): Where O meas,i is the observed value and S pred,i the computed value of salinity or water level. meas O is the mean observed salinity or water level data and pred S is the mean computed salinity or water level.
The closer the values of RMSE and ABSERR to zero, and the R 2 to unity, the better the model performance is evaluated 10 (Abu El-Nasr et al., 2005). For Percent bias (PBIAS) measures the average tendency PBIAS, expressed as a percentage, of the simulated data to be larger or smaller than their observed counterparts (Gupta et al., 1999). The optimal value of PBIAS is 0, with low-magnitude values indicating accurate model simulation (Moriasi et al., 2007). Positive values indicate model underestimation bias and negative values indicate model overestimation bias (Gupta et al., 1999). The Nash-Sutcliffe efficiency (EF) (Nash and Sutcliffe, 1970) is a normalized statistic that determines the relative magnitude of the residual 15 variance (noise) compared to the measured data variance (information). EF ranges between −∞ and 1 (1 inclusive), with EF = 1, the closer the model EF efficiency is to 1, the more accurate is the model. Values between 0 and 1 are generally viewed as acceptable levels of performance, whereas values ≤ 0 indicate unacceptable performance (Moriasi et al., 2007). 20 Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2016-276, 2016 Manuscript under review for journal Hydrol. Earth Syst. Sci. Published: 13 June 2016 c Author(s) 2016. CC-BY 3.0 License.
The indicators of hydrodynamic-salinity intrusion model are summarized in Tables 6 and 7. The two models are the EF and R 2 coefficients are very near to unity. This result demonstrates the good performance of the analytical model. Also, this shows that the proposed coupled analytical model by Cai et al., (2015a) is applicable and useful. 5 The statistical performances of the numerical model uses water level for comparison. Values of statistical parameters indicated in Table 8 show good correlation model calculations and measurements during calibration and validation. These indicate that the model can estimate the water level at Kenitra fairly well.

10
The statistical indicators for transport model are summarized in Table 9. The results show that the computed salinity concentration follows observed data, which suggest that the presented mass transport model is a reasonably efficient tool for predicting the impact of total lunar eclipse on salt intrusion in alluvial estuaries.

Total lunar eclipse impact
The impact of eclipse on river hydrodynamics is mainly caused by a moon closest to the earth, causing strongest 15 gravitational pull. The maximum salinity at high water along the Sebou estuary has been described in section 4. Total lunar eclipse impact on the maximum salinity at different locations compared to the normal situation is given in Table 10. The results clearly show that the lunar eclipse impact on salinity intrusion is highly significant.
On other hand, Fig. 16 shows the profiles of salinity during total lunar eclipse compared to the spring-neap tides. It appears 20 that the salt intrusion curve computed in the neap-spring tides are recession-type, while it becomes a dome-type shape at eclipse day. According to Nguyen et al. (2012), this is subjected to changes in the degree of convergence of the crosssectional profile, and the effect of the mixing due to freshwater discharge (i.e. that increasing the tidal amplitude at the mouth tends to produce shorter convergence lengths of the cross-sectional area and width). The convergence or divergence of the channel can dramatically change the shape of the salt intrusion curve (Gay and O"Donnell, 2007;Cai et al., 2015b). In 25 addition, Savenije (2005), shows that the recession-type curve occurs in narrow estuaries having a near-prismatic shape and high river discharge and dome-type shape which occurs in strong funnel-shaped estuaries (with a short convergence length) (see Appendix A). At eclipse day, when the channel converges strongly, the mixed water retains relatively higher salinity from the estuary mouth. However, salinity profiles under all spring-neap tides show a gradual decrease from the mouth to the upstream reach . 30 Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2016-276, 2016 Manuscript under review for journal Hydrol. Earth Syst. Sci. Published: 13 June 2016 c Author(s) 2016. CC-BY 3.0 License.
Additionally, it can be shown that the part of the Sebou estuary that is affected by the total lunar eclipse is from 20 to 25 km upstream of the river mouth. A water level rise as showed above during this exceptional event moves the excessive salinity (25g/l) until 20 km upstream. On other hand, water level rise causes large augmentation of salinity in the mesohaline zone of 5 the Sebou estuary. Also, if we considerer Kenitra location, an increase of water level for 0.8m causes an increase of 9.4 g/l in salinity witch correspond to a salinity augmentation of 11.75 g/l per metre of increased water level. At the Oulad salma drinking water station salinity increased to a value of 1.2 g/l that exceeds the limit value of 500 mg/l recommended by the World Health Organization (WHO) for drinking water.

10
Computations using hybrid (hydrodynamic) and salt intrusion models during total lunar eclipse required the recalculation of the geometric parameters of the estuary i.e., cross-sectional areas A0, convergence lengths a and the optimization of the dispersion coefficients D. Geometry is one of the most important parameters in the hydrodynamic and salt intrusion models.
It affects the character of the salinity distribution and appears prominently in the shape of salt intrusion curves during extreme events. Computed results reveal that variations in the sensitivity of these parameters are likely to depend on changes 15 in geometric characteristics.

Conclusions
The purpose of this paper was to study the impact of the total lunar eclipse of 28th September 2015 on salt intrusion in Sebou river estuary. It is, to our knowledge, the first time that this infrequent phenomenon is studied in terms of its influence on water quality. Field measurements showed a change of the salinity profiles form along the estuary axis and a notable 20 water level rise, compared with normal situations studied in our earlier works. In addition, results show that the average salt content increased in the reach between 0-25 km, as a result of water volume rise at mouth. An hybrid proposed bu Cai et al. (2014) coupled to analytical salt intrusion model in alluvial estuaries (Savenije, 2005) and a numerical model (HEC-RAS) have been applied in the Sebou river estuary. Calculations during lunar eclipse using the coupled hybrid-salt intrusion model required the recalculation of the geometric parameters of the estuary i.e., cross-sectional areas A 0 , convergence lengths a. A 25 good fit was obtained between computed and measured salinity during this extreme event. These models reproduce very well the salinity rise. Statistical indicators show that these models fit adequately salinity observations in the Sebou estuary.
A comparison between the two applied models is not the objective of this study since each one can be applied for specific management purposes. The analytical models are helpful for situations where a quick longitudinal salinity profiles is needed.
On other hand, the numerical 1-D model is powerful where a temporal salinity variation is carried out in a specific location, 5 but it needs more data and time for its implementation. Hence, these tools can be very helpful for water managers and engineering to make preliminary estimates on the salt intrusion along the estuary axis even during extreme events. These extreme events can concern similar total lunar eclipse, see level rise due to climate change, a sea tsunami.
Finally, the impact of extreme events on the water quality of Sebou estuary should be considered by managers. Rapid interventions, based on the predictions of our mathematical models can thus be taken. These interventions may involve the 10 pumping station closure for example.

Acknowledgement
The authors would like to express their gratitude to: H. Qanza, M. Hachimi, O. khabali and O. El Mountassir for the efforts in the field measurements during the total lunar eclipse day. The authors would also like to acknowledge the technicians the 15 water services of Kenitra; and the engineers of the National Agency of Ports for their availability and collaboration.

Appendix A. Types of salt intrusion and shape of salt intrusion curves
Salinity distribution is a veritable fingerprint of each estuary and in direct relation to both its geometric form and hydrology. 20 For partially mixed and well-mixed estuaries, a number of designations are used to classify salinity profiles into three types depending on their shape. The following types are distinguished (Savenije, 2005(Savenije, , 2012) (see Fig. 17): -Recession shape, which occurs in narrow estuaries with a near-prismatic shape and a high river discharge (Savenije, 2005(Savenije, , 2012. 25 -Bell shape, which occurs in estuaries that have a trumpet shape, i.e. a long convergence length in the upstream part, but a short convergence length near the mouth (Savenije, 2005(Savenije, , 2012).
-Dome shape, which occurs in strong funnel-shaped estuaries (with a short convergence length) (Savenije, 2005(Savenije, , 2012 Where h 0 is the tidal average depth, v 0 is the maximum tidal velocity, u 0 is the velocity of freshwater, N is the Canter Cremers" estuary number, K is the Van der Burgh"s coefficient, f is the Darcy-Weisbach"s coefficient, F is the Froude 15 number and F D is the densimetric Froude number.  Table 2. Hybrid solution of tidal wave propagation in convergent estuaries (Cai et al., 2015b).  Table 5. Salinity distribution data showing the salinity at the mouth (S 0 ), tidal excursion (E), Richardson number (N R ), dispersion coefficient at High Water Slack, Van Der Burgh"s coefficient (K) and salt intrusion length at High Water Slack 5 (L).         Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2016-276, 2016 Manuscript under review for journal Hydrol. Earth Syst. Sci.