Estuarine salinity response to freshwater pulses

,


Introduction
Freshwater pulses, during which the freshwater discharge by rivers exceeds three times its long-yearly average value and which last no longer than one month, are common features in many estuaries around the world.They are mostly the result of strong precipitation in the upstream river catchment area (Tee & Lim, 1987;Valle-Levinson et al., 2002;Gong et al., 2007;Liu et al., 2008;Du & Park, 2019;Du et al., 2019;Guerra-Chanis et al., 2021), opening of a freshwater reservoir (Ingram et al., 1986;Lepage & Ingram, 1988), or a combination of those two (Díez-Minguito et al., 2013).The increased freshwater discharge causes a strong downstream transport of salt, which has a large impact on the ecology in the estuary and on the agriculture of the lands around the estuary (Paerl et al., 2006;McFarland et al., 2022).All the above-cited studies indicate that the adjustment time, here defined as the time during which the salinity in an estuary adjusts to high river discharge, is in the order of 1-2 days.Observational studies report that freshwater pulses can cause the salt intrusion length, which is defined as the distance of the 2 psu isohaline to the estuary mouth (Monismith et al., 2002), to shift by tens of kilometers (Díez-Minguito et al., 2013).An estuary can even become entirely fresh (Du & Park, 2019).After such pulses, the estuary returns to its non-disturbed state.Values of the recovery time, defined as the time it takes for the salt intrusion length to reach its background value again, widely vary, but typically they are considerably larger than values of the preceding adjustment times.For example, Valle-Levinson et al. (2002) found 10 days for the Chesapeake Bay (USA), whilst Gong et al. (2007) reported four months for York River estuary, which is located in the same area.
The overall aim of this study is to gain a more detailed understanding of how an estuary will respond to freshwater pulses with different intensity and duration.For such purposes, it is helpful to employ idealized models, which only represent the most dominant physical processes and assume a simplified geometry.Besides yielding insight into the dynamics, these models are fast, flexible and are thus suitable for extensive sensitivity analysis.Earlier studies on estuarine physics (Hansen and Rattray (1965); Chatwin (1976); MacCready (2004); Geyer and MacCready (2014)) have demonstrated the added value of idealized models with respect to detailed numerical models.
The current knowledge of estuarine adjustment to changes in river discharge originates from both simplified and more sophisticated numerical models.Kranenburg (1986) demonstrated, by using analytical arguments applied to a one-dimensional model, that the response timescale, i.e. the time during which an estuary responds to a decrease or increase in river discharge, is inversely proportional to the river discharge after the change.
This finding explains the difference between adjustment time and recovery time.Hetland and Geyer (2004) used a three-dimensional primitive equation model with idealized geometry and simple turbulence formulations to study response timescales.They found a clear difference between adjustment and recovery time, which is in line with the finding of Kranenburg (1986).They argued that during net upstream transport of salt, the motion of the salt intrusion adds constructively to the (subtidal) bottom layer flow.This means that velocities in the bottom layer are stronger than during net downstream transport, so import of salt will experience stronger resistance from the bottom drag and will thus be slower than net export of salt.Chen (2015) extended the analysis of Kranenburg (1986) by allowing the density-driven flow in his model to be time-dependent.He argued that the difference between adjustment and recovery time is the result of the non-linear response of salt intrusion length to changes in river discharge.Monismith (2017) employed a modified version of the model of Chen (2015) to study the unsteadiness of the salt intrusion length under different time-dependent forcings.His model showed good skill in hindcasting salt-intrusion lengths in the northern part of San Francisco Bay.
These studies yielded important insights into the timescales associated with the response of salt intrusion to changes in river discharge.Important to mention here is that the idealized models for estuarine adjustment assume that creation of salinity stratification by vertically-sheared velocity is balanced by destruction of stratification by vertical mixing.This assumption is based on Pritchard (1954), who analyzed observations in the James River estuary under relatively low river discharge.Hereafter, we will refer to this balance of processes determining the stratification as the Pritchard balance.
Studies by MacCready (2007) and Ralston et al. (2008) demonstrated that this assumption works quite well in cases that they consider, but these cases do not include strong freshwater pulses.Dijkstra and Schuttelaars (2021) showed that in steady state the Pritchard balance does not hold in the high-discharge regime.It may be expected that this is also true for time-dependent cases.Knowledge gaps also exist with regard to the sensitivity of the response of the estuary to freshwater pulses for different environmental settings, e.g.different strengths of tides.
The specific aims of this study are twofold.The first is to show the limitations of the Pritchard balance when investigating strong freshwater pulses.The second is to investigate the sensitivity of the estuarine salinity response to a freshwater pulse to different parameters.We quantify the estuarine salinity response by calculating adjustment timescales, recovery timescales and changes in salt intrusion lengths.There are three research questions associated with this second aim: 1) What is the effect of the background conditions of the estuary, i.e. the background river discharge and the strength of the tides, on the salinity response?2) What is the effect of the strength of the peak river discharge on the salinity response?3) What is the effect of the duration of the pulse on the salinity response?
The remaining of this paper is organized as follows: In Section 2.1, deficiencies, including negative salinity values, are identified when the model of MacCready (2007) (MC07 hereafer), which uses the Pritchard balance, is forced with observed river discharge during a strong freshwater pulse in the Guadalquivir Estuary (Spain).A new model, which does not rely on the Pritchard balance, is presented in Section 2.2.This model does not have the deficiencies of MC07 when used to simulate freshwater pulses in the Guadalquivir Estuary (Section 2.3).Afterwards, a sensitivity analysis is done in a more idealized model setup.The experimental setup is given in Section 2.4, followed by the results and discussion (Section 3) and the conclusions (Section 4).

Limitations of the Pritchard balance
In order to show the limitations of available idealized models for estuarine adjustment, the MC07 model is used to simulate the estuarine response to an observed strong freshwater pulse.This model simulates time-dependent, tidally averaged, width-averaged estuarine flow and salinity, building on Hansen and Rattray (1965).The vertical momentum balance is hydrostatic, while in the horizontal a balance is assumed between the pressure gradient force and internal friction.Furthermore, the Pritchard balance is used to describe the vertical structure of salinity.The MC07 model is here applied to the freshwater pulse in February 2009 in the Guadalquivir estuary (Díez-Minguito et al., 2013;Wang et al., 2014;Losada et al., 2017).This pulse has a maximum discharge (main river + tributaries) of 889 m 3 s −1 , while the river discharge in the month before the pulse has an average value of about 32.3 m 3 s −1 .The model settings are as follows: the estuary is 110 km long and its width increases exponentially from 150 m at the upstream limit to 650 m at the mouth.The thalweg has an average depth of 7.1 m (Díez-Minguito et al., 2013), so this is used as the depth of the estuary.The vertical eddy-viscosity coefficient, vertical eddy-diffusion coefficient and horizontal eddy-diffusion coefficient are chosen as in Guha and Lawrence (2013).This means that they depend on the strength of the tidal current and a turbulent length scale, which is the estuary depth for the vertical coefficients and the estuary width for the horizontal coefficient.The model is forced with the observed river discharge from the Alcala del Rio dam and from the four main tributaries after this dam: Aznalcázar, El Gergal, Guadaíra and the Torre del Águila (Agencia de Medio Ambiente y Agua de Andalucía, see chguadalquivir.es/saih/Inicio.aspx).
The representative tidal current amplitude is based on measurements (Navarro et al., 2011) and set to be 1.15 m s −1 and the ocean salinity is 35 psu.The horizontal grid size is 250 m and a timestep of 15 seconds is used to ensure numerical stability.
Results from this simulation are displayed in Fig. 1a and b.Before the pulse, the salinity field is only slightly disturbed by the variations in river discharge.During the freshwater pulse, surface salinity values drop and within a few days after the start of the pulse they reach values of −4.7 psu close to the mouth.The minimum value for surface salinity is thus below zero.Note that at the same time, bottom salinity values at the estuary mouth are prescribed to be 35 psu, which means that the estuary is strongly stratified during the pulse in this simulation.After the pulse, negative salinity values disappear.Salt intrusion recovers in about three weeks to values comparable to the ones before the pulse.The negative values of surface salinity during the pulse are unphysical and motivated the development of a new model that is presented in the next section.

Domain
The study area consists of two parts: an estuary and the adjacent sea.We use x as the along-channel coordinate, where x = −L e is the upstream limit, x = 0 is the estuary mouth, and x = L s is the boundary of the adjacent sea.The width of the es- where b 0 is the width at the upstream limit and L b is the e-folding length scale which controls the width convergence.The estuary and the open sea have different values of L b ; the depth H is constant throughout the domain.

Hydrodynamic module
The hydrodynamic equations are identical to those in MC07.The equations are averaged over the cross-channel y-coordinate and over the tides.Wind is ignored and the Boussinesq approximation is used, with an equation of state that expresses a linear relation between variations in salinity and variations in density.The only difference with respect to MC07 concerns the boundary condition at the bottom z = −H, which is taken to be a partial slip condition (instead of no-slip), as in Dijkstra and Schuttelaars (2021), and reads Here, S f = 2Av H is the friction coefficient, A v is the vertical eddy-viscosity coefficient (assumed constant, see later), u is the along-channel velocity and z is the vertical coordinate.At the upstream limit a river discharge Q is imposed: The along-channel velocity u and salinity s are split in their respective depth-averaged parts (denoted by a bar) and depth-dependent parts (denoted by primes): The solutions of the equations for along-channel velocity read where α = gβH 3 48Av .Here, g = 9.81 m s −2 is gravitational acceleration and β the isohaline contraction coefficient of water (= 7.6×10 −4 psu −1 ).The vertical eddy-viscosity coefficient is parametrized as , where U T is the amplitude of the tidal current and c v = 7.28×10 −5 is an empirical constant (Ralston et al., 2008).This formulation is based on the assumption that the relevant velocity scale for turbulent mixing in an estuary is the amplitude of the tidal current and the limiting vertical length scale of the turbulent eddies is the depth of the estuary.The physical interpretation of Eq. 5 is that the depth-averaged current is solely due to the river discharge, and that the vertical velocity shear is caused by the river current and the density-driven flow.Hereafter, we will refer to u as the exchange flow (Geyer & MacCready, 2014).The vertical velocity w follows from continuity, which results in

Salt module
The salt conservation equation is where t is time.The horizontal eddy-diffusion coefficient is parametrized as where c h = 0.035 is an empirically determined constant (Banas et al., 2004).A closure relation for the vertical eddy-diffusion coefficient Schmidt number (Ralston et al., 2008).At the upstream limit (x = −L e ) a river salinity s ri imposed.The simulated domain stretches well beyond the limit of salt intrusion, to avoid that this condition affects the salinity dynamics in the estuary.In the part that represents the adjacent sea, width increases strongly with distance to the mouth, so that the river flow will become very weak at the sea boundary x = L s .This allows us to assume that at this downstream boundary of the domain (located seaward of the estuary mouth at x = 0) salinity will be well-mixed over the vertical and we can set salinity to be equal to the ocean salinity s oc over the entire depth.Hence, At the bottom and the free surface the vertical salt flux vanishes: At the transition between the parts at x = 0, both s and the salt transport b(us−K h ∂s ∂x ) have to be continuous.Since u and b are continuous, this last condition implies that ∂s ∂x has to be continuous as well.This can be written as

Solution method
To solve for salinity, Eq. 4 is inserted into Eq.8 and this equation is averaged over the depth, resulting in the depth-averaged salt balance: Here T 1 is the tendency term.Terms T 2 -T 4 contain along-channel variations of three width-integrated and depth-mean salt fluxes: that due to river flow (T 2 ), due to exchange flow (T 3 , which can be split into a contribution by the density-driven current and a contribution induced by the river current) and a diffusive flux (T 4 ).The equation for the evolution of s is found by subtracting Eq. 12 from Eq. 8, yielding ∂s ∂t (13) Term T 5 is the tendency term.Terms T 6 and T 7 represent the horizontal advection of s and T 8 the creation of stratification by vertical velocity shear.Term T 9 is equal to minus T 3 , T 10 represents vertical advection, T 11 vertical diffusion and finally T 12 represents horizontal diffusion.Note that when the Pritchard balance is applied, only terms T 8 and T 11 are taken into account in Eq. 13.
This set of equations is solved for s and s .To deal with the vertical variations, a Galerkin method (see e.g.Canuto et al. (2012)) is used.For this, the depth-dependent salinity is written as a Fourier series where N is the number of vertical modes and s n are the Fourier components, which depend on the horizontal coordinate x and on time t.This expression is substituted in to Eq. 13, and afterwards this equation is projected onto the Fourier modes.Combined with Eq. 12, this yields N +1 equations for s(x, t) and the s n (x, t), n = 1, 2, ...N , which are numerically solved by using central differences on a spatially uniform grid in x, while time integration is performed with the Crank-Nicolson method (Crank & Nicolson, 1947).This results in a system of N +1 algebraic equations at every grid point, containing values of s and s n at the previous and current timestep.This system of equations is solved with the Newton-Raphson method (see e.g.Galántai (2000)).

Performance of the model
The new model is next used to simulate the same freshwater pulse in the Guadalquivir Estuary, as was done with the MC07 model.Model settings are the same as in Section 2.1.Additionally, the number of vertical modes is chosen as N = 10 and the sea domain is modeled as a 25 km long channel with an e-folding length scale of 2.5 km.As the numerical scheme is now implicit, it allows for larger timesteps.A standard value of ∆t = 12 hours is chosen, but to guarantee accuracy a smaller ∆t is chosen when the salt field changes fast.When the change in salinity is large, the Newton-Raphson algorithm may not converge for too large timesteps.When this happens, also a smaller timestep is chosen; a minimum timestep of 15 minutes is used.
Results from this simulation are displayed in Fig. 1c.Before and after the freshwater pulse, salt intrusion is stronger than in the simulation with the MC07 model (Fig. 1b).
However, during the pulse, no negative values of salinity are simulated, which indicates that our model has overcome the problematic behavior of the MC07 model when simulating strong freshwater pulses.To check the numerical accuracy of the solutions, additional simulations were done where the time and spatial resolution were taken twice as small and with the number of vertical modes increased to N = 15.The maximum difference in salinity between the simulations was smaller than 0.01 psu, assuring that the results are sufficiently accurate.
There are three possible explanations for the difference between the MC07 model and the new model: the different boundary condition for momentum at the bottom, the different boundary condition for salinity at the sea boundary (and the inclusion of the sea domain) or the generalized equation for the evolution of s .First, to determine whether the boundary condition for momentum is the cause, an additional simulation was done where the no-slip boundary condition from MC07 was used as a boundary condition for momentum.In this simulation, salt intrusion before and after the pulse is smaller than in the simulation with the partial slip boundary condition, but no negative salinity values were simulated.Second, the effect of the boundary condition for salinity at the sea boundary was investigated by calculating the value of salinity at the bottom at the estuary mouth.This has a minimum value of 32.8 psu, which is a relatively small (≈ 10%) deviation from the prescribed value in the MC07 model.Third, the effect of the additional terms in Eq. 13 was studied.Fig. 2 displays the terms in this equation during the simulation.It is clear that the largest terms during the entire simulation are T 8 and T 11 , which are the only terms taken into account in MC07.However, during the pulse, also other terms become important, in particular T 10 , the vertical (upward) advection of salt.
This explains why negative salinity can occur in MC07: the amount of destruction of stratification by vertical mixing during freshwater pulses is too small compared to the creation by vertical shear.This leads to an overestimation of stratification in MC07, leading to negative salinity in the surface layer.

Set-up of the numerical experiments
Next, the model is used to study the sensitivity of the estuarine salinity response to freshwater pulses.For this a model domain is chosen that is a straight channel with a width of 1000 m and a depth of 10 m.The adjacent sea part is 25 km long and it has a convergence length of 2.5 km.This setting is chosen to mimic an 'average' coastal plain estuary.Sea salinity is 35 psu and river salinity is 0 psu.The horizontal grid size varies between a minimum of 125 m and maximum of 500 m for different simulations and the number of vertical modes ranges from 5 to 15, depending on the strength of the stratification.The timestep has values between 15 minutes and 24 hours, giving sufficiently accurate solutions.To model a freshwater pulse, an initial state is chosen in which the subtidal estuarine salinity is steady and in equilibrium with the background river discharge Q bg .The pulse starts when the river discharge increases suddenly to its peak value Q p .The river discharge remains then for a time T pulse at this peak discharge.After the pulse, the river discharge instantly returns to Q bg .Each simulation is continued until the salt intrusion length has recovered to its initial value.The salt intrusion length X 2 is defined as the distance between the estuary mouth and the most upstream position where the salinity exceeds 2 psu.
To quantify the salinity response to a freshwater pulse, several output quantities are defined: change in salt intrusion length ∆X 2 , adjustment time T adj and recovery time T rec .Change in salt intrusion length ∆X 2 is the difference between the value of the salt intrusion length before the pulse and its minimum value during the pulse.Adjustment time T adj is defined as X 2 (t = T adj ) = X 2 (t = 0) − 0.9∆X 2 , which is the time it takes for the salinity in the estuary to adjust to the peak river discharge.Recovery time T rec is the time after the pulse when X 2 (t = T rec ) = X 2 (t = 0) − 0.1∆X 2 , so it is the time when 90% of the recovery of X 2 after the pulse has taken place.These quantities are scaled to make the resulting dependencies more general.As a scale for ∆X 2 , the background salt intrusion length X 2 (t = 0) is used.For T adj , the adjustment time found by Kranenburg (1986) is used as a scale, which reads T adj,sc = 0.9bH∆X2

Qp
. The factor 0.9 accounts for the fact that here T adj is defined when 90% of the change in X 2 has occurred.Finally, the scale for T rec is T rec,sc = 0.9bHX2(t=0) . This is the timescale that results from the assumption that recovery is primarily due to salt transport by the exchange flow.Classical theory (Chatwin, 1976) is applied, i.e., u s during the recovery is approximately balanced by salt transport due to river flow.The factor 0.9 is added for the same reason as that in the formulation of T adj,sc .
The T adj is the scaled duration of the pulse.
Specifically, for addressing research question 1, a number of simulations is performed where F r R,p is fixed and F r T and F r R,bg are varied, since these two quantities are shown to determine the equilibrium state of an estuary (Geyer & MacCready, 2014).The duration of the pulse is chosen to exceed the adjustment time, so equilibrium with the peak river discharge is reached during the pulse.This set of simulations will be referred to as experiment set Background.For answering research question 2, F r R,bg and F r R,p are varied.The tidal Froude number F r T is fixed at a value of 0.62 (U T = 1 m s −1 ).The duration of the pulse T pulse again exceeds the adjustment time.This set of simulations will be referred to as experiment set Peak.For adressing research question 3, the duration of the pulse is varied, as well as F r R,bg and F r R,p .The tidal Froude number F r T is again fixed at 0.62.The values of F r R,bg and F r R,p are equal to those in set Peak.Two series of simulations are done where Tpulse = 1 2 and 1 4 .These simulations will be referred to as experiment set Short.Table 1 contains the range of values of the dimensional parameters for all the experiments that were conducted.
The range of the parameters is based on the following.The amplitude of the tidal current U T is chosen between 0.75 and 1.5 m s −1 , which results in F r T ranging from 0.46 to 0.93.Weaker tides are not investigated, because the momentum balance relies on the assumption of moderate to strong tidal currents with respect to the river current.Larger tidal currents are not investigated because they are considered to be non-realistic.The range of values of the freshwater Froude numbers is based on daily discharge values from five estuaries where freshwater pulses are identified.The considered estuaries are the Gironde (France), the Guadiana (Spain/Portugal), the Guadalquivir (Spain), San Francisco Bay (USA) and the Tagus (Portugal).Specifics of the river discharge datasets are given in Table 2. Freshwater pulses are identified in the river discharge datasets and displayed in (F r R,p , F r R,bg ) parameter space in Fig. 3. Based on these observations, a value of F r R,p = 0.15 is chosen as a standard value for the simulations of experiment set Background.Since a freshwater pulse is defined here as an event when the river discharge exceeds three times its long-yearly average value, an obvious upper bound for F r R,bg is 0.05 for this set of experiments, one-third of the value of F r R,p .The lower bound is F r R,bg = 0.001.For experiment set Peak, values of F r R,bg range from 0.001 to 0.075 and those of F r R,p range between 0.02 and 0.3.These boundaries are indicated by the black lines in Fig. 3.The majority of the observed freshwater pulses fit within these boundaries, but not all of them.
Observed pulses for which 3 F r R,bg > F r R,p probably started far from a steady state (shortly before the pulse, another pulse occurred) and are thus not considered.The strongest freshwater pulses in the Guadalquivir and Guadiana have F r R,p > 0.3 and are also outside the investigated parameter space.This is done because multiple model assumptions are not valid anymore under such extreme circumstances, in particular the width and depth being constant.Such strong freshwater pulses will increase the water level significantly and flood lands next to the estuary.Moreover, simplifying assumptions regarding the momentum balance, which rely on the estuary being partially to well mixed, do not hold during such extreme events.(Table 1).Clearly, all dimensional response characteristics (contours in Fig. 4) become lower for higher F r T and higher F r R,bg .The scaled quantities in panels 4b-d show a different behavior: they only weakly depend on F r T and for increasing F r R,bg the relative The dependence of background salt intrusion length X 2 (t = 0) on F r T and F r R,bg (Fig. 4a) follows the power-law relationship X 2 (t = 0) ∼ F r R,bg −1/3 F r −1 T , according to classical theory on estuarine salt dynamics, in which a dominant balance is assumed between salt export by river flow and salt import by exchange flow (Hansen & Rattray, 1965;Chatwin, 1976;Geyer & MacCready, 2014).However, for low values of the river flow (F r R,bg < 0.005), horizontal diffusion of salt is important, next to the salt import by exchange flow, and this power-law is not valid.Excluding this regime, a least squares fit to the numerical results yields X 2 (t = 0) ∼ F r R,bg −0.40±0.03F r −1.00±0.16 T , in good agreement with classical theory.This theory also explains the patterns found in Fig. 4b, as it predicts that ∆X 2 ∼ (F r T , while a least squares fit to the data yields ∆X 2 ∼ (F r −0.28±0.00 The patterns shown in Fig. 4c,d can be understood by identifying and analyzing the processes that act during adjustment and recovery.Fig. 4c shows that for F r R,bg < 0.015 the adjustment time T adj T adj,sc .In this 'high-pulse regime', where the peak river discharge is relatively large compared to the background river discharge, the dominant process for adjustment is the export of salt by river flow during the pulse.In the 'moderate-to-low pulse regime' (the upper part of the panel) T adj is considerably larger than T adj,sc .During the adjustment, other salt transport mechanisms are then effective as well, viz.import of salt by both the exchange flow and by horizontal diffusion.As they oppose the salt export by river flow, the adjustment time is larger than that would result from river flow alone.The fact that the value of F r T does affect the dimensional adjustment time but not the scaled adjustment time indicates that its effect is mostly through a larger change in salt intrusion length (see panel b), but that the celerity of the adjustment is not sensitive to F r T .
A similar reasoning applies to the recovery time: it will be close to the scaled value T rec,sc if the recovery process is controlled by salt transport due to the exchange flow, as described by the classical theory.Fig. 4d shows that this only approximately holds in the 'weak pulse regime', i.e., in the upper part of the diagram.For moderate to stronger pulses, values of the recovery time are approximately half of T rec,sc .This deviation from quasi-steady classical theory exists because immediately after the pulse, the landward salt transport due to exchange flow is substantially larger than the seaward transport by river flow.A larger value of F r T , i.e. stronger tidal mixing, will result in slower recovery, because the magnitude of the exchange flow is inversely proportional to the value of U T .Yet the recovery time is not very sensitive to the value of U T because this effect is compensated by the fact that the change in salt intrusion length also decreases approximately linearly for higher U T .
To look at this in more detail, we present results of the change in salt content of an estuary for different values of the background river discharge and of the tidal current amplitude.The integrated salt balance is obtained by integrating Eq. 12 over the volume of the estuary: Here, it is assumed that salt transport vanishes at the upstream limit.Term S 1 represents time rate of change of salt content in the estuary, and S 2 -S 4 are depth-averaged salt fluxes at the estuary mouth due to river flow, exchange flow and horizontal diffusion, respectively.Fig. 5 shows time series of the river discharge, salt intrusion lengths and terms S 2 , S 3 and S 4 .Regarding the recovery time, we see that a substantial part of the recovery takes place in the first few days after the pulse (Fig. 5a).This means that just after the pulse the salt transport due to exchange flow is very important for the total recovery time.The value of T rec,sc is calculated by assuming this transport scales with the transport of salt by the background river flow, which is not a good estimate during this period, especially for strong pulses.Thus T rec will be shorter than T rec,sc for strong pulses, which is indeed found.The effect of U T is clearly illustrated in Fig. 5c and d  Results of experiment set Peak are shown in Fig. 6.The same quantities as in Fig. 4 are displayed, except now for different values of F r R,bg and F r R,p , while the amplitude of the tidal current and duration of the freshwater pulse are kept fixed (Table 1).Panel The patterns in parameter space in experiment set Background are mostly explained by whether a pulse is 'weak' or 'strong', i.e. from the ratio between F r R,bg and F r R,p .
Regarding the change in salt intrusion length, its dependence on F r R,p follows again from . This behaviour is visible in Fig. 6b.
A least-squares fit to this data yields an exponent of 0.43 ± 0.01 in this relation, displaying the validity of classical theory.The linear dependence of ∆X 2 on T pulse is a consequence of the fact that the time rate of change of X 2 is linear in time during most of the adjustment (Figs. 5 and 7).Thus the change in salt intrusion length can be estimated from multiplying the downstream velocity of the salt intrusion with the duration of the pulse.Because of this, ∆X 2 will indeed depend linearly on T pulse when no equilibrium is reached.The finding that T rec does not depend on T pulse is due to the fact that the exchange flow is not fully developed before equilibrium with the peak discharge is reached (see Figs. 5     and 7).When the pulse ends at this stage, salt import will thus be weaker compared to when the pulse would have reached equilibrium.So during the recovery after a shorter pulse, the upstream velocity of the salt intrusion will be smaller.At the same time, ∆X 2 is also smaller for shorter pulses.These two effects have similar magnitude and will compensate each other.

Specific application
In this section the model performance is assessed by applying it to observed freshwater pulses in the Guadalquivir Estuary.For this purpose, the model was slightly extended to a new geometry that consists of multiple (instead of one) estuarine parts, as is shown in Fig. 9.In each of these parts the equations as presented in Section 2.2 are solved.For salinity the matching conditions shown in Eq. 11 are used at the boundaries of the parts.Furthermore, additional river discharge of four tributaries are added at the beginning of each part.The other model settings are equal to those used in Section 2.3, with one exception: for salinity at the river boundary a value of 0.5 psu was used, based on observations.Details about the observations are given in Navarro et al. (2011Navarro et al. ( , 2012)).To determine the subtidal salt intrusion length, first a Gaussian filter with a half-amplitude of 12 hours is applied to the raw salinity measurements to average over the tides.Afterwards, the observed salinity (observations are done at the surface) is linearly interpolated between the measurement points and the most upstream point where the salinity exceeds 2 psu is identified.During the observational period, several freshwater pulses occurred: one in February 2009 and a series of three pulses in 2010.
Simulations are done in order to capture the effects of these pulses.Fig. 10 displays the results of these simulations and the observations in the Guadalquivir.To quantify the differences, the root-mean-square error of the observed and simulated salt intrusion length is calculated, which will be noted as RMSE(X 2 ).For the 2009 case, RMSE(X 2 )= 9.6 km.However, this number does not reflect the temporal differences: before day 50 of the year 2009, RMSE(X 2 ) = 3.7 km and after this date it is 15.7 km.For the simulations of the pulses in 2010, we have RMSE(X 2 )= 5.4 km.These values indicate that the model is capable of simulating the temporal behavior of the salt intrusion length in the Guadalquivir Estuary during freshwater pulses.Clearly, there are differences between simulated and observed salt intrusion length, which could be reduced by applying detailed model tuning.For example, Wang et al. (2014) and Losada et al. (2017) argued that the 2009 freshwater pulse in the Guadalquivir created a mud layer on the bottom of the estuary, which decreased the hydraulic drag.This effect could be taken into account by adjusting the value of the partial slip parameter after the pulse.

Other remarks
An interesting difference between the results presented here and existing literature concerns the recovery time.Here, we find that this quantity depends only on the river discharge during the recovery, whereas in previous studies (e.g.Kranenburg (1986); Hetland and Geyer (2004); Chen (2015); Monismith (2017)) it is stated that it depends on the change in salt intrusion length.The reason for this difference is that Kranenburg (1986) assumes that during the recovery the exchange flow does not vary in time.However, here we show that the evolution in time of the exchange flow during the recovery is important for the recovery time (Fig. 5 and Fig. 7).Chen (2015) accounts of time-varying exchange flow, but he estimates recovery time from linearized equations, thereby assuming small changes in exchange flow.Our study, on the other hand, clearly shows that these changes are large.Finally, Monismith (2017) accounts for large changes in exchange flow during recovery, but he assumes that depth-averaged salinity at the estuary mouth equals ocean salinity.Certainly, during strong freshwater pulses that condition is too restrictive.
Finally, we remark that the estuarine salt response to pulses with a duration that is shorter than the adjustment timescale of the system is not considered in the existing literature.We find that the change in salt intrusion length is linearly related to the duration of the pulse (Fig. 8).This is relevant in real estuaries.The duration of freshwater pulses in the observational datasets can be compared to the theoretical adjustment time given by our model.In the smaller estuaries analysed, less than half of the pulses do not reach the equilibrium state (Guadiana: 0.49; Guadalquivir: 0.39; Tagus: 0.38) but in the larger estuaries this portion is even larger (Gironde: 0.76 ; San Francisco Bay: 0.60).
So the duration of the pulse is often the limiting factor for the change in salt intrusion length.

Conclusions
The aim of this study was to quantify the dependence of the estuarine salinity response to freshwater pulses to the background conditions, the intensity and the duration of the pulse.Application of the MacCready (2007) model, which relies on the Pritchard balance, to observed freshwater pulses in the Guadalquivir Estuary showed that use of this balance results in negative salinity values.We therefore developed a new model, which uses a more detailed description of the vertical salinity structure.Simulations with this model did not show negative salinity and moreover, the model performs well when applied to observed freshwater pulses.
Model simulations revealed that the influence of the background conditions on the salinity response to a given freshwater pulse is mainly through the background river discharge; the strength of the tides is of minor importance.Changes in salt intrusion length ∆X 2 can be estimated successfully from classical theory, but this theory is incorrect regarding adjustment time T adj for weak pulses and recovery time T rec for strong pulses.
Simulations with different strengths of the peak river discharge revealed that for ∆X 2 the ratio of peak to background river discharge determines the response.Interestingly, the peak river discharge is the most important control for T adj while for T rec its value is not important.When the duration of the freshwater pulse is too small to reach equilibrium, ∆X 2 will be linearly related to the duration of the pulse, but T rec is not affected.
Observed freshwater pulse characteristics indicate that this control on ∆X 2 is important in real estuaries.

Open Research
Software used to generate the data and create the figures used in this study can be found at git.science.uu.nl/w.y.biemond/code-and-data-freshwater-pulses .git,as well as the river discharge datasets used.Observational data of the Guadalquivir Estuary can be found at zenodo.org/record/3459610.

Figure 1 .
Figure 1.Application of two different models for the case of an observed freshwater pulse in the Guadalquivir Estuary in February 2009.(a) Time series of observed river discharge.(b) Simulated surface salinity ssur with the MacCready (2007) model, which uses the Pritchard balance, versus time and along-channel coordinate x, where x = 0 is the estuary mouth.The white area indicates where ssur is negative.(c) As b, except for the simulation with the model presented in Section 2.2.

Figure 2 .
Figure 2. The magnitudes of the different terms in Eq. 13 versus time and depth in the first grid cell upstream of x = 0 during the simulation of the February 2009 freshwater pulse in the Guadalquivir Estuary with our model.a) T5, b) T6, c) T7, d) T8, e) T9, f) T10, g) T11 and h) T12.
research questions as formulated in the introduction separated the quantification of the estuarine salinity response to freshwater pulses into three parts: the sensitivity to the background state (research question 1), the sensitivity to the peak river discharge (research question 2) and the sensitivity to the duration of the pulse (research question 3).These different research questions motivate the variation of four dimensional parameters: U T , Q bg , Q p and T pulse .These are converted into four dimensionless parameters, which are F r T , F r R,bg , F r R,p and Tpulse .Here, F r T = U T c is the tidal Froude number, with c = √ gβHs oc an internal velocity scale that equals twice the maximum internal wave speed.Furthermore, F r R,bg = Q bg bHc is defined as the background freshwater Froude number and F r R,p = Qp bHc the peak freshwater Froude number.Finally, Tpulse = T pulse

3. 1
Sensitivity analysis Results of experiment set Background are displayed in Fig. 4. Panels show the dependence of background salt intrusion length, change in salt intrusion length, adjustment time and recovery time on the tidal Froude number F r T and background freshwater Froude number F r R,bg .Note that intensity and duration of the freshwater pulse are kept fixed

Figure 3 . 378 Figure 4 .
Figure 3. Observed freshwater pulses in the F rR,p, F r R,bg parameter space.The cross-shaped markers are freshwater pulses where the peak river discharge is less than three times the background river discharge, the diamond-shaped markers indicate events where the peak river discharge exceeds this value and the circular markers indicate freshwater pulses where F rR,p > 0.3.The grey line is where F rR,p = 3 F r R,bg .The black box indicates the part of the parameter space that was investigated by experiments Peak and Short.

Fig. 5a reveals
Fig.5areveals that the adjustment of the salt intrusion length to a freshwater pulse is to a good approximation linear in time.Panel b shows that the magnitude of the salt : a lower value of U T means that the change in salt intrusion length is larger (panel c), but also the salt flux due to exchange flow S 3 is stronger (panel d) and these effects compensate each other.

Figure 5 .
Figure 5. (a) Time series of river discharge (black line) and salt intrusion length (coloured lines) for different background river discharges Q bg .Only the discharge for the case Q bg = 50 m 3 s −1 is plotted.These simulations are from experiment set Background (i.e.Qp = 2423 m 3 s −1 ).(b) Time series of the different terms at the right-hand side of Eq. 15.The colours refer to the same simulations as in panel (a).(c)-(d) As (a)-(b), except for different values of UT .

a
shows background salt intrusion lengths for reference purposes.Panels b and c show the same patterns as in experiment set Background : with increasing strength of the pulses the change in salt intrusion length becomes larger and the adjustment time becomes smaller.The recovery time barely depends on the value of F r R,p (panel d).

Figure 6 .
Figure 6.Results of experiment set Peak.(a) Contour plot of background salt intrusion length X2(t = 0) (values in km) as a function of peak freshwater Froude number F rR,p and background freshwater Froude number F r R,bg .(b) As panel a, except for the change in salt intrusion length ∆X2 (contours, values in km) and the scaled change in salt intrusion length ∆X2/X2(t = 0) (colors).(c).As panel b, except for the adjustment time T adj (contours, values in days) and the scaled adjustment time T adj /T adj,sc (colors).(d) As panel c, except for the recovery time Trec and the scaled recovery time Trec/Trec,sc.The black area indicates where F rR,p < 3 F r R,bg .

Figure 7 .
Figure 7.As Fig. 5a-b, except for different values of peak river discharge Qp from experiment set Peak.In panel a only the discharge for the case Qp = 1000 m 3 s −1 is plotted.

Figure 8 .Figure 9 .
Figure 8. Results of experiment Short.(a) Contour plot of change in salt intrusion length ∆X2 (contours, values in km) and scaled change in salt intrusion length ∆X2/X2(t = 0) (colors) a function of peak freshwater Froude number F rR,p and background freshwater Froude number F r R,bg for T pulse = 1 2 T adj .(b) As panel a, except for T pulse = 1 4 T adj .(c) As panel a, except for the recovery time Trec (contours, values in days) and the scaled recovery time T adj /T adj,sc (col-

Figure 10 .
Figure 10.Time series of observed river discharge Q and observed and simulated salt intrusion length X2 for the Guadalquivir Estuary.The discharge (the red line) is the sum of the main river plus the four tributaries.The dark blue line is the observed X2 and the light blue line is the simulated X2.(a) For the freshwater pulse in 2009.(b) For the series of freshwater pulses in 2010.

Table 1 .
Amplitude of tidal current UT , background river discharge Q bg , peak river discharge Qp and duration of the pulse T pulse for the different sets of experiments.

Table 2 .
Specifications of river discharge datasets for five estuaries where freshwater pulses occur.For the Gironde, river discharge from the Garonne and Dordogne are added.For the San Francisco Bay, the dataset combines multiple sources.