Multiscale Temporal Response of Salt Intrusion to Transient River and Ocean Forcing

Salt intrusion in surface waters endangers freshwater availability, influences water quality


Introduction
Deltas and estuaries occupy a critical interface between riverine and marine environments.They are now exposed to interlinked and sometimes conflicting pressures from climate change and human interventions.One specific example concerns the intrusion of salty ocean waters leading to increased salinity of surface waters.Increased salt intrusion threatens the provision of many high-value estuarine ecosystem services since salinity is an important control on ecosystem functioning.Indeed, saltier waters impact freshwater resourcing from surface waters and aquifers (e.g., White & Kaplan, 2017), crop production (e.g., Smajgl et al., 2015), potentially public health (Shammi et al., 2019), and significantly degrades delicate transitional habitats (Herbert et al., 2015;Tully et al., 2019;White & Kaplan, 2017).Salt dynamics also play an important role on sediment dynamics (e.g., de Nijs & Pietrzak, 2012;Mietta et al., 2009;Postma, 1967), on oxygen depletion (e.g., Bruce et al., 2014;Hong & Shen, 2012) and on the biome (e.g., Crain et al., 2004;Little et al., 2017; M. Zhang et al., 2017).
Abstract Salt intrusion in surface waters endangers freshwater availability, influences water quality, and affects estuarine ecosystem services with high economic and social values.Salt transport and the resulting salinity distributions result from the non-linear interactions between salt and water dynamics.Estuaries are often considered under (quasi)-steady assumption or by focusing on specific timescales.Our understanding of their temporal multiscale response to transient forcing is limited, which hinders the implementation of effective mitigation strategies.We apply wavelet analyses to quantify the variability of salt intrusion from hourly to seasonal timescales and unravel the temporal variability of its response across scales.We focus on an estuary that undergoes significant transient forcing, the Modaomen estuary in the Pearl River Delta, and apply the wavelet analyses to year-long data generated by a coastal ocean numerical model.Our results show that this estuary responds to changes in tidal and riverine forcing throughout the year over interwoven timescales.Our results highlight the temporal variability of the salt intrusion response time both within a given regime and for the transition between regimes.They also suggest that tides control the response time more strongly than river discharge, even though river discharge determines the magnitude of the salt intrusion, and thus modulates the evolution of the salt intrusion response time.We propose a broadly applicable framework to calculate response times with simple data.These results can provide a first-order guidance for design and implementation of estuarine management strategies and mitigation measures that ensure water access and facilitate sustainable development.
Plain Language Summary With more and more people living by coast, more water resources are needed, and the fragile equilibrium of estuaries and coastal ecosystems is at risk.In estuaries and deltas, the salty water from the ocean meets the fresh water from the river and can propagate upstream.We call this salt intrusion.One of the solutions used to control salt intrusion is releasing water to push the salt intrusion back, but for this to be effective it needs to be well timed.We use a computer model to understand how the changes in river and ocean affect how far and when the salty water reaches upstream.We show that the river determines how far it reaches but it is mostly the tide that determines when this will happen.Since the ocean moves at different scales, the salt intrusion also responds at multiple scales with different time delays, and it is difficult to distinguish one from another because they overlap.The time delays are not constant and evolve throughout the year.We need to consider the multiple overlapping scales and the evolving time response to ensure the sustainable management of estuaries.
Salt dynamics and transport in estuaries, and therefore salt intrusion, are determined by a delicate balance between tidal forcing, river discharge (hereafter Q r ), atmospheric effects, and bathymetric control.Even though much of our understanding of estuarine dynamics relies on quasi-steady state assumptions and builds on the early works of Chatwin (1976) and Hansen and Rattray (1965), the importance of time dependence in estuarine dynamics has long been recognised starting with the pioneer work by Kranenburg (1986) in well-mixed estuaries and that of MacCready (1999) in well-mixed and strongly stratified estuaries.Estuarine forcing and salinity distributions are rarely steady and ignoring this fact can have first order effects on salt balances (Banas et al., 2004).Whether the estuary can be assumed to be quasi steady or unsteady fundamentally depends on how the response timescale of the estuary (i.e., the adjustment time) compares with the forcing timescale (Banas et al., 2004;Hetland & Geyer, 2004;Lerczak et al., 2006;MacCready, 2007;Monismith et al., 2002).Nevertheless, the more recent view is that estuaries are in a continual state of adjustment to tidal forcing (Geyer & MacCready, 2014), with unsteadiness an integral part in estuarine dynamics (Banas et al., 2004).
Altogether, tidal forcing, river discharge, atmospheric forcing and bathymetric control operate over a range of intertwined temporal scales.The salt intrusion length (hereafter L s ) is defined as the upstream horizontal distance of an isohaline from the estuary mouth.Commonly used values range from 0.5 psu (e.g., Gong & Shen, 2011) to 2 psu (e.g., Monismith et al., 2002).L s and Q r are typically inversely correlated following a power law relationship   ∝  −  with exponents depending on the estuary and its conditions for example, n = 1/3 (e.g., Hansen & Rattray, 1965;Hetland & Geyer, 2004;MacCready, 1999) or n = 1/7 (Monismith et al., 2002).Estuaries have shown an asymmetry in the response to increasing and decreasing Q r with faster response to a rising Q r than to a falling Q r (Gong & Shen, 2011;Hetland & Geyer, 2004), which has been attributed to non-linear behavior (e.g., S.-N.Chen, 2015;Hetland & Geyer, 2004).The effect of the tide on salt dynamics is complex since it depends on the type of estuary and the relative importance of the different mechanisms at play.Even though L s has been expressed as a function of the tidal current amplitude (e.g., Bowen & Geyer, 2003;Prandle, 2004), the influence of the tide is generally taken into account through horizontal and vertical mixing coefficients.At subtidal temporal scales, L s varies with the spring-neap cycle with a phase that depends on the type of estuary.Larger L s occurs during neap tides for exchange flow dominated estuaries and partially mixed estuaries (e.g., Banas et al., 2004;Lerczak et al., 2006) while for estuaries where tidal dispersion dominates and salt-wedge estuaries larger L s will occur during spring tides (e.g., Ralston et al., 2010;Xue et al., 2009).This tidal influence on L s is modulated by Q r , since L s seems to be particularly sensitive to spring and neap modulation under high Q r conditions (Bowen & Geyer, 2003;Gong & Shen, 2011;Lerczak et al., 2009;Ralston et al., 2010), meaning that temporal scales of varying Q r and varying tides are intertwined.
Valuable information on the temporal response of L s to changes in Q r and tidal mixing (including adjustment time and processes) has been provided using primitive models (e.g., Kranenburg, 1986;MacCready, 1999MacCready, , 2007;;Monismith, 2017;S.-N. Chen, 2015), observational studies (e.g., Banas et al., 2004;Lerczak et al., 2009) and full three-dimensional circulation models (e.g., Hetland & Geyer, 2004;Ralston et al., 2008).These studies have often been limited to a single stratification regime (i.e., well-mixed, partially mixed or salt-wedge) and typically consider an average state or mean state (i.e., subtidal variability/steady state) around which there is some variation (e.g., Kranenburg, 1986;Lerczak et al., 2006;MacCready, 2007).However, estuaries are characterized by a range of forcing conditions, which implies they do not occupy a single point in any estuarine parameter space (Geyer & MacCready, 2014) and may not always remain under a single regime.In the systems that undergo rapid and large changes in external forcing, the estuary is unlikely to remain under the same regime, and the assumption of average state and noise may be problematic (Banas et al., 2004;MacCready, 1999).
In many cases, salt intrusion is intensifying due to for example, sea-level rise (e.g., Hong & Shen, 2012), dredging (e.g., Eslami et al., 2019) or reduction of rainfall and freshflows (e.g., White & Kaplan, 2017), and remains problematic in spite of mitigation efforts (e.g., White & Kaplan, 2017).Since management strategies and mitigation measures will include human interventions ultimately modifying river discharge and tide propagation in the estuary, we need a more complete understanding of the transient response of estuaries to these changing forcings across conditions and regimes.Water management measures for salt intrusion often rely on flow modulation.Questions on the size and duration of the measures needed for effective management and mitigation require a focus on the estuarine response across temporal scales that fully considers the large variability inherent to the river and the ocean forcing.We note that this is important not only for the many estuaries which already have very large river discharge variability, but also more generally because seasonal variability is likely to become 10.1029/2021JC017523 3 of 24 more pronounced due to climate change and associated droughts/floods (e.g., Trenberth et al., 2014;van Vliet et al., 2013).A better understanding of the timescales and transient behavior of estuaries is the starting point toward prognostic modeling tools that would ensure water supply for end users, would support sustainable development and thriving ecosystems, and would allow to cope with future changes.
The aim of this study is to address the gap in characterizing the multiscale transient response of estuaries under large changes in external forcing, to unravel the different scales at play, as well as the evolution of the salt intrusion response time across temporal scales.We will determine the multiscale (i.e., spanning from hours to seasons) response of salt intrusion to river and ocean forcings, when these are undergoing large transient change.
Our approach will apply a widely used state-of-the-art ocean model to a specific case study: the Modaomen Estuary in the Pearl River Delta.This specific estuary is of particular interest in this case as it has been subject to increasingly severe salt intrusion, which is threatening water resourcing for part of the world's largest urban area (World Bank, 2015).Our study covers transient river and ocean forcings by considering a realistic numerical simulation covering a full year and thus including strong seasonal monsoonal variability.Even though wind and atmospheric forcing is known to impact salt transport and intrusion (e.g., Li & Li, 2012;Scully et al., 2005;S.-N. Chen & Sanford, 2009;Valle-Levinson et al., 2019;Xie & Li, 2018), ocean and Q r remain the dominant drivers, particularly so for the Modaomen estuary where tidal range and Q r have the largest effect on the the intraseasonal variability (Lin et al., 2019).We therefore focus on ocean and river forcing and neglect atmospheric forcing.
We discuss the influence of ocean and river discharge on L s variability, we investigate the significant timescales of variability of each one of the time series -tidal forcing, Q r and L s -, how they relate one another and how the relationships evolve overtime.To do this, we use spectral methods: wavelet and wavelet coherence.We use the wavelet to isolate time-scales of variability from hours to seasons, and wavelet coherence to examine the link between L s and the forcings.We first present the case study location, followed by the numerical model and setup used, and the analyses techniques used (Section 2).Results (Section 3) include validation of the model and time-series analyses of the numerical outputs, primarily using wavelet analysis.Section 4 is a discussion of the main findings and we draw the conclusion in Section 5.

Methods
Our approach applies the Finite Volume Community Ocean Model (FVCOM) model to the entire Pearl River Delta.The model is validated against available data for water elevation and salinity at several stations along the Modaomen Estuary to provide confidence in its ability to adequately reproduce the relevant processes.We then use wavelet analysis tools on the numerical results from a year-long simulation (i.e., September 2007to September 2008) to extract information on the multiscale response of the salt intrusion to time-varying forcings.

Case Study: Modaomen Branch of the Pearl River Delta
The region of study is the Pearl River Delta, China, the city of cities and the largest megalopolis in the world both in terms of population and urbanized area since 2010 (World Bank, 2015), which includes two megacities (population >10 M) and several cities with over 1 M habitants.The region has undergone rapid development since the 1980s and the morphology and hydrology of the Pearl River Delta have been heavily modified by human interventions such as dam construction (e.g., Cai et al., 2012;Cai et al., 2019;S. Zhang et al., 2008;W. Zhang et al., 2012), channel dredging (e.g., Cai et al., 2012;Cai et al., 2019;He et al., 2019;Luo et al., 2007;W. Zhang, Cui, et al., 2010;Z.;Wu et al., 2018), and land reclamation of intertidal wetlands (e.g., Cai et al., 2019;He et al., 2019;W. Zhang et al., 2010, 2015, Z. Wu et al., 2018).
Three major tributaries combine to form a complex channel network in the Pearl River Delta with eight outlets discharging into the South China Sea.The West River, the North River, and the East River deliver respectively 77%, 15%, and 8% (C. S. Wu et al., 2016) of the overall freshwater discharge to the delta.The Modaomen estuary is located in the downstream portion of the West River and it is the outlet with the largest discharge (Cai et al., 2012): estimated annual river discharge at Denglongshan Station (Figure 1) of 88.393 billion m 3 (Gong & Shen, 2011).The region presents a monsoonal hydrological regime with one dry season (December-February) and one wet season (April-September) when about 80% of the annual freshwater input is delivered (W.Zhang et al., 2009).Several water treatment plants are located on the Modaomen estuary and provide freshwater to over 10 million people in nearby cities (e.g., Jiangmen, Zhongshan, Macao, Zhuhai).This water supply is now threatened by salt intrusion with a worsening trend in recent years (e.g., Gong & Shen, 2011;Gong et al., 2012;Luo et al., 2007;W. Zhang et al., 2009).
The Modaomen estuary is weakly convergent and it is on average less than 3 km wide.The mean depth of the estuary is 4.3 m and a deep channel meanders through the estuary with depths ranging between 5 and 17 m.Several sand bars result in lateral variability of the bathymetry.
The tide in the Modaomen is mixed slightly semidiurnal (i.e., form factor close to 1, the details are presented later in Section 3.1.Model Validation) with strong ebb dominance both in duration and velocity magnitude, most likely due to the high Q r and the microtidal regime (Gong, Chen, et al., 2018).The tidal range decreases from 1.11 m at the mouth to 0.86 m in Denglongshan and continues to decrease as the tide propagates further upstream.The major constituents are M 2 , S 2 , K 1 , and O 1 (Gong & Shen, 2011).The interaction between constituents generates a marked asymmetry between consecutive spring and neap tides.
The effect of different forcings on salt dynamics and on L s in the Modaomen has been extensively studied with studies acknowledging the role of Q r (e.g., Gong et al., 2012;Xu et al., 2015;W. Zhang et al., 2013;Z. Zhang, Cui, et al., 2010), interacting branches (Gong et al., 2012), tides (e.g., Gong et al., 2012;Gong & Shen, 2011;Yuan et al., 2013), wave current interaction (Gong, Chen, et al., 2018), wind (e.g., Gong et al., 2012;Lin et al., 2019) or sea level rise (Yuan et al., 2013;Yuan & Zhu, 2015).Previous studies on the Modaomen have focused on the dry season when the estuary is partially mixed and the salt intrusion is most acute.These studies identified the contribution of different drivers at subtidal scales (e.g., Gong & Shen, 2011) or intraseasonal and interannual scales (e.g., Lin et al., 2019).However, the multiscale temporal response and its evolution over time have not been systematically studied yet (Liu et al., 2014).Here we take the Modaomen estuary as an example to investigate the variability across scales spanning from hours to seasons of saltwater intrusion and its evolution across a wide range of river discharge conditions.

Numerical Model
We used FVCOM (v3.1.6)(C.Chen et al., 2003) to numerically simulate hydrodynamics, salt dynamics, and salt intrusion in the Pearl River Delta.FVCOM is an unstructured grid, finite volume, three-dimensional hydrodynamical model.Governing equations for momentum, continuity, temperature, salinity, and density are solved under the Boussinnesq approximation and the hydrostatic assumption.Equations are solved in FVCOM using a split-mode method.The external (barotropic) time step is 0.25 s and the ratio of internal (baroclinic) to external time step is five.FVCOM uses a second-order accurate advection scheme, however due to the unstructured grid it can provide similar accuracy to higher-order advection schemes in structured grid models (Huang et al., 2008).We use the General Ocean Turbulent Model (Umlauf & Burchard, 2005) for the vertical turbulence closure model choosing a two equation second-order turbulence model with balance equations for turbulent kinetic energy and dissipation rate (i.e., k-ϵ model).The stability functions are obtained from an explicit algebraic model assuming quasi-equilibrium and using model constants from (Canuto et al., 2001).For the horizontal mixing, we consider a constant Smagorinsky horizontal diffusivity coefficient of 0.2 m/s 2 (Smagorinsky, 1963).The model allows wetting and drying, and we consider a minimum depth of 5 cm with cells becoming inactive if the water level is below that value.Bottom friction is implemented assuming the bottom cell is in the logarithmic layer of the rough bottom boundary layer with a user-defined constant roughness lengthscale z 0 .Due to lack of information on sea bed composition and roughness, we consider a uniform z 0 over the entire domain.z 0 was the only parameter in the calibration.The calibration was based on minimizing errors in tidal amplitude and tidal phase.Sensitivity tests on the value of z 0 have indicated optimal results for a spatially constant roughness length z 0 of 0.001 m.This value is in the range of values used in this type of systems (e.g., Ralston et al., 2010;Ralston et al., 2016).The drag coefficient C d in the quadratic law is set to a minimum value of 0.0025.
The unstructured grid of the model allows for resolutions as high as 20 m in the narrowest channels of the delta (with at least five elements across section) and at the coast but has larger cell sizes (i.e., 20 km) at the offshore boundary (cf. Figure 1a and 1b).In the Modaomen branch resolution varies between 50 and 100-300 m.The overall study area covers 19.5°N-23.5°N and 110°E-118°E and is shown in Figure 1.The domain extends beyond the shelf and includes part of the South China Sea but the focus of the study solely is the Modaomen branch of the Pearl River Delta.The significant extension of the model both upstream and offshore ensures that the model boundaries do not interfere with the area of interest.For the vertical discretization we consider fifteen σ (i.e., terrain following) uniform layers meaning that the thickness of the layers in the region of interest did not exceed 0.7 m.Sensitivity tests on the number of vertical layers (not shown here) confirmed that 15 layers is a good compromise between computational cost and model accuracy.The bathymetry (Figure 1c and 1d) was generated from a set of surveys undertaken in 2008 in order to produce the most accurate representation of the delta at the time of the simulation.The water depth varies from about 700 m deep off the shelf to less than 2 m within the delta.
Upstream boundary conditions (cf. Figure 1 for locations of upstream boundaries) prescribe daily mean river discharge Q r form the three main rivers at the hydrological stations of Gaoyao, Shiiao and Boluo with a salinity of 0 psu and a temperature of 10°C.At the offshore open boundary, the model is forced with tidal elevation using 13 tidal constituents (M 2 , S 2 , N 2 , K 2 , K 1 , P 1 , O 1 , Q 1 , Mf, Mm, M 4 , Ms 4 and Mn 4 ) from the TPXO8 model (Egbert & Erofeeva, 2002).The boundary conditions for the depth-dependent time-varying salinity and temperature rely on the interpolation of HYbrid Coordinate Ocean Model (HYCOM) global reanalysis model (Bleck, 2002;Bleck & Boudra, 1981) along the open boundary nodes.We then use a Blumberg & Khanta implicit open boundary condition for the perturbation (depth-dependent component after removing depth-averaged value) of temperature and salinity with an adimensional nudging timescale coefficient of 0.0014.A sponge layer of 20 km radius (one cell around the boundary) and 0.0001 coefficient is applied at the open boundary to filter high frequency computational noise at the outflow side (C.Chen et al., 2013).

Model Validation
The model was validated against available time series for water level and salinity at several monitoring stations along the Modaomen branch (Figure 1).Validation of the modeled water levels is undertaken on both phase and amplitudes of the tidal constituents extracted from model results and monitoring data.This provides a better assessment of the model veracity than a direct quantitative comparison of the time series, given that atmospheric forcing is not included in the model.We use NOCtide, a MATLAB ® implementation of the long-standing harmonic tidal analyzer used at the National Oceanography Centre, see example Murray (1964), to compute the tidal constituents of both the measurements and the model results.Daily maximum salinity data are available at different stations along the Modaomen branch at a fixed distance above the bed corresponding to 0.8 × h, where h is the mean water depth at the station (Lin et al., 2019).In most cases, these stations match water extraction sites.We interpolated modeled salinity at the elevation above bed of the observations and compared to measured daily maximum salinity.

Wavelet Analysis
Salt intrusion length is extracted from salinity numerical model outputs.It is defined as the distance along the longitudinal section shown in Figure 1 from the river mouth to the bottom 0.5 psu isohaline.The longitudinal section generally follows the deep channel in the Modaomen except when islands are present, in which case the longitudinal section follows the eastern branch where both salinity monitoring and water extraction stations are located.The definition of salt intrusion length is consistent with previous studies in the area (Gong et al., 2012;J.-C. Chen et al., 2004).The 0.5 psu threshold is chosen here as it is the threshold used by water extraction operators in the region: water extraction is stopped when salinity is above this threshold.We also calculate the salt intrusion length based on the near-surface 0.5 psu limit.Differences between the bottom salt intrusion length and the surface salt intrusion length are linked to stratification: in absence of stratification (well-mixed case) the isohalines are almost vertical and salt intrusion length is expected to be the same near-bed and near-surface; when stratification is present isohalines deviate from vertical leading to a difference between salt intrusion length nearbed and near-surface.We thus consider the difference between the bottom salt intrusion length and the surface salt intrusion length as a proxy to the level of stratification in the estuary.In order to remove high-frequency fluctuations due to the primary tidal constituents (semidiurnal and diurnal) and their harmonics (e.g., overtides) and ease interpretation at subtidal scales, both surface and bottom intrusions are low-pass filtered using a Chebyshev type II filter with a stopband for periods up to 26 hr and a passband for periods above 30 hr.Throughout this manuscript L s refers exclusively to the un-filtered signal while the low-pass filtered signal is referred to as LP-L s .
We analyze the multiscale temporal response of salt intrusion to transient forcing by computing the continuous wavelet transforms of the forcing (i.e., Q r , sea surface height and tidal range) and L s , and then calculating the cross-wavelet transforms and wavelet coherence between the various forcings and the salt intrusion length.We refer the reader to Grinsted et al. (2004) and Torrence and Compo (1998) for complete details on wavelet analyses and will only provide a succinct overview below.The wavelet function operates as a band-pass filter to the time series and allows us to obtain the dominant modes of variability and how those modes vary in time.We choose to use the Morlet function, which is non-orthogonal and complex, to calculate the continuous wavelet transform (CWT) because it oscillates with a mean value of zero and is particularly adequate for tidal signals.The wavelet transform at each scale are directly comparable to each other and to the transforms of other time series because the wavelet function at each scale is normalized.To ease the comparison of the different power spectra, the wavelet spectrum is normalized by the variance thus the colorbar in our CWT figures represents the normalized wavelet power spectrum.Statistical significance is estimated against a red noise model using Monte Carlo methods.The thick contour lines in our figures indicate the 95% confidence level thus highlighting the regions in the time-frequency space where the variance is significantly above red noise at the 95% confidence level.The cone of influence (COI) indicates where edge effects due to the finite time series become important and is depicted by the shaded area in all our wavelet figures.
We investigated the relationship between various forcing (either Q r , SSH and tidal range at the mouth) and L s by computing the cross-wavelet coherence (WTC) between them.The wavelet coherence is the square of the cross-spectrum normalized by the individual power spectra and can be thought of as a localized correlation coefficient in the time-frequency space.It varies from 0 to 1, with high values obtained when wavelets are highly coherent.An important advantage of wavelet coherence is that it indicates coherence even when the common power is low.Cross wavelet phase angles are also calculated and provide information on the phase between the two time series analyzed.They are represented by arrows in the figures with arrows pointing right denoting signals in phase, and arrows pointing left denoting anti-phase signals.It is important to note that, when not indicating in-phase or anti-phase signals, the phase angle may not provide a single solution with respect to phase lead or lag, that is, signal 1 leading signal 2 by a phase angle of 90° and signal 1 lagging signal 2 by 270° both correspond to the same phase angle.Other considerations (e.g., cause/consequence) are then required to fully determine phase lead/lag.Statistical significance for the WTC was also determined against a red noise model.

Model Validation
We assess the performance of the model by validating simulated tide and salinity.For all the stations the model shows good agreement with the measurements for tides.Amplitude and phase for diurnal, semidiurnal and shallow water constituents are shown in Figure 2. The best results are obtained for diurnal and semidiurnal constituents, which are well reproduced by the model.The largest discrepancies (underprediction) are observed for Zhuyin, the northernmost station (i.e., 34 km inland from the river mouth), which are likely to be caused by errors in the bathymetry.However the water extraction plants are situated further south, where the agreement is good.Model-data agreement appears less good for the overtides, probably partly because of much lower magnitude of the shallow water constituents.Intriguingly, the best agreement for the overtides appears to be found in Zhuyin whilst M4 and MS4 are underpredicted, and MN4 is overpredicted for the other stations.This different pattern for the overtide may relate to errors canceling out for Zhuyin, for example errors in the bathymetry in the Modaomen could be canceling incorrect incoming overtide that itself could be due to inaccuracies in the bathymetry offshore.For the other stations small differences in diurnal and semidiurnal constituents cascade down to the overtides resulting in larger discrepancies.Aggregated amplitude and phase Root Mean Square Error (RMSE) for the main constituents range between 1 and 2 cm and 5.29°-5.68°for O 1 , K 1 , and S 1 (cf.Table 1).The largest aggregated rmse are obtained for M 2 , the constituent showing the largest amplitude and phase.The form factor F = (K 1 + O 1 )/(M 2 + S 2 ), calculated from the tide gauge data, is 1.04 for Dahengqin, 1.03 for Denglongshan, 1.01 for Sanzao, and 1.05 for Zhuyin characterizing the Modaomen as mixed, with semidiurnal dominance.
The model validation for salinity is presented in Figure 3 by comparing model results with daily salinity maximum at several monitoring stations along the Modaomen.Both observations and model results show marked spring-neap variations, the phasing of which is well captured by the model.Maximum salinity values obtained with the model are approximately within a 18% range of the maximum salinity measured.The observed drop in salinity values linked to a peak in Q r between 8 February and 15 February is also reproduced.RMSE and correlation coefficients (cf.Table 2) vary along the branch with best absolute values obtained for the most onshore station (i.e., Majiao).Worst values for Denglongshan (i.e., RMSE = 3.54 psu and correlation coefficient r = 0.64) may relate to secondary circulation around the island and the bathymetry not being accurately resolved there.
It is important to note that our intention is not to develop close to perfect numerical simulations, but to use a numerical model that does reproduce the key processes we focus on here.To that end, the model does capture adequately the effects of tides and Q r on salinity.Some of the model data discrepancies can be explained by the model not including all processes.For example, the period of poor agreement between 15 February and 7 March corresponds to enhanced SW winds that weaken the landward salt transport (Gong, Lin, et al., 2018) but are not included in the model.In addition we have not considered freshwater extraction from the estuary in our simulations due to a lack of adequate data, which may partly explain the salinity underprediction in the upper estuary.The combination of spurious mixing induced by the numerical scheme and mixing from the turbulence closure 10.1029/2021JC017523 9 of 24 can also result in excessive total mixing of salt in FVCOM (e.g., Ralston et al., 2016).This over-mixing may explain the model's underprediction of upstream salinity (i.e., at Zhuyin, 34 km form the river mouth).

Salt Intrusion Length (L s )
The Q r for all three of the branches of the river is shown in the upper panel in Figure 4.The analysis period includes a relatively isolated high discharge peak in February and is followed by an enhanced Q r regime from April onwards with several discharge pulses each of approximately 2 weeks duration.The largest discharge is reached in June (i.e., total Q r 68,000 m 3 /s) and its value is two orders of magnitude larger than the baseline values of the dry season (i.e., total Q r under 2700 m 3 /s).Two more high discharges are observed before a return toward dry season averages, however the Q r during the tail of the monsoon is still several times (i.e., about 14,000 m 3 /s) that of the dry season.
Both SSH and tidal range are shown in Figure 4b.The tidal range shows a modulation at the fortnightly scales (spring-neap and/or tropical/equatorial cycles).For simplicity, we refer here to the moderate tides in the fortnightly modulation as neap tides, conversely we use spring tides for the large tides in the fortnightly modulation.The tidal range also shows seasonal variability with lower variability around the equinoxes (i.e., mid-March, mid-September).This minima in tidal range variability is concomitant with the moment when Q r starts to significantly increase.
Both top and bottom L s show (Figure 4c) a response (as expected) to changing Q r and tidal forcing.Surface and bottom salt intrusion lengths show a fortnightly modulation (spring-neap and/or tropical/equatorial cycles)  (Figure 4c), and there appears to be a time lag between the fortnightly modulation of the SSH (Figure 4b) and the fortnightly modulation of the L s (Figure 4c).During the first part of the study period before the Q r increases (i.e., until June), neap tides correspond with increased L s , most likely due to reduced tidal mixing (Gong et al., 2014;Gong & Shen, 2011;Ralston et al., 2008).Larger L s at neap tides is characteristic of exchange flow dominated estuaries (Banas et al., 2004) and has been found in other partially mixed estuaries (e.g., Lerczak et al., 2006).Salt intrusion also appears to be increased for the secondary neap tides in the month, the same behavior has been previously shown in the Modaomen for the dry season (Gong & Shen, 2011).This fortnightly behavior due to tides changes from July onwards (concomitantly with a high flow regime).The magnitude of the spring-neap variations in L s is reduced, and L s reaches its maximum during spring tides.This behavior is characteristic of salt wedge estuaries (e.g., Ralston et al., 2010;Xue et al., 2009).Salt intrusion also appears to respond quicker from neap to spring and more gradually from spring to neap prior to June (i.e., during the low Q r season) which is consistent with previous results in the region (Gong & Shen, 2011).This trend seems to reverse from July onwards (i.e., concomitant with a high Q r regime) when the salt intrusion appears to respond quicker from spring to neap than from neap to spring.
Salt intrusion length is also clearly sensitive to Q r as each peak in Q r results in the shortening of the L s salt intrusion length.The highest Q r result in negative L s , corresponding to salt being entirely flushed out of the estuary (0.5 psu isohaline offshore of the estuary mouth).During June, the increased Q r favors the development of a fresh water plume outside the estuary.Salt intrusion responds more quickly to increases than to decreases in Q r (Figure 4c) again matching earlier studies in the Modaomen (Gong & Shen, 2011).Hetland and Geyer (2004) link this asymmetry in the estuary response to nonlinear effects of bottom drag.
Several distinct periods can be distinguished from the time series of L s (Figure 4c).During the first part of the study period before Q r increases (i.e., January to March), bottom and surface salt intrusions are overlapping but   not identical, thus suggesting a partially-mixed estuary.During this period, tides seem to dominate L s variations although probably because of the little change in Q r .Towards the end of this period of low river discharged (April), the surface L s reduces more than the bottom L s with increasing Q r .The resulting periodic divergence between surface and bottom L s suggests increasing (periodic) stratification.However, the increase in Q r only has a limited effect on the bottom L s .This behavior is consistent with increased Q r increasing stratification but not significantly impacting L s due to a negative feedback from increased estuarine circulation (MacCready, 2004;Wei et al., 2017).However tides (i.e., semidiurnal and fortnightly modulation) still appear to be the dominant driver in L s salt intrusion length variability.
Over the main peak in Q r (i.e., June), bottom and surface L s are no longer always overlapping, denoting much increased stratification and the estuary moving from partially mixed to (strongly) stratified.The response of L s is much more sensitive to the spring-neap cycle under high flow conditions with a large fluctuation of L s over the spring-neap cycle.Similar behavior has been shown in the Modaomen (Gong & Shen, 2011) and in other estuaries for example, the Hudson (Bowen & Geyer, 2003;Lerczak et al., 2009) and Merricmack (Ralston et al., 2010).At the same time both L s are reducing corresponding to salt being flushed out of the estuary.
Finally, after the main Q r peak (i.e., from July onwards), L s recovers toward dry season values and bottom and surface L s are only periodically overlapping around spring tides.This last period presents a fortnightly modulation in the stratification with the moment of higher stratification evolving over time from spring tides to neap tides.The response of the estuary during this last period appears to be quicker from spring to neap tides than it is from neap to spring.Similarly to the period of high Q r , the range of the oscillations at the tidal scales is narrower for this period than it was for the period with low Q r .

Multiscale Temporal Analysis
We present in the following section the results of the multiscale temporal analysis.We use CWT and Wavelet Coherence (WTC) to offer a visual representation of the temporal variation of the different forcings and L s , and their relationship in the time-frequency space.The CWT expands time series into time-frequency space and can help find dominant modes of variability and track how those modes vary with time.From the two CWTs of two time series, we construct the WTC that can be thought of as a localized correlation coefficient in time-frequency space between two CWTs.The main strength of the wavelet analysis is that it enables to find localized periodicities (or bands) that can be linked to specific processes and that it offers a graphical representation of how the relationship between the time series (i.e., forcing and L s ) evolves through time.

Wavelet River Discharge (Q r )
Wavelet analysis in hydrology commonly focuses on time scales beyond the capabilities of the present analysis (i.e., seasonal to multi-annual (e.g., Jalón-Rojas et al., 2016;Labat, 2010).The CWT power for the Q r (Figure 5) highlights important variability at short time scales (i.e., days to months here).The high power at seasonal scales (i.e., 3-month band) is representative of a monsoon type hydrological regime.There is also significant power at sub-monthly scales (15-30 days) from January onwards.This high power corresponds to local (in time) variability in the high Q r events -the Q r pulses appear to last about 2 weeks.This variability at sub-monthly scales overlaps with scales associated with fortnightly tidal variability.Between 12 May and 12 July, when high Q r pulses occur in a more persistent manner, small timescales corresponding to the bands of four to 15 days also show high CWT power due to fluctuations at shorter scales in Q r .

Wavelet Sea Surface Height (SSH) and Tidal Range
Even though the wavelet analysis of tidal water elevation may not provide much additional information compared with knowledge of the tidal constituent, we still present these results as they allow a relatively simple description of tidal variability as depicted through the wavelet analysis.One important aspect is that the scales in the wavelet analysis (CWT here but this also applies to the WTC) are not precise enough to exactly identify each tidal constituent.Instead, they are better understood as representing specific tidal and subtidal species or bands: overtide, semidiurnal, diurnal, fortnightly, and monthly bands.We specifically present analyses for both SSH and tidal range as they are expected to provide different and complementary information.The SSH displays clear variability at semi-diurnal and diurnal scales the amplitude of which is modulated at fortnightly and monthly scales (see Figure 4b, for example) and which results in negligible subtidal residual (i.e., for bands above the daily band).As such, the expectation is for the CWT analysis for SSH to show most of the variability at the semi-diurnal and diurnal scales with CWT power modulated according to the tidal amplitude modulation.This is indeed confirmed in Figure 6b.In addition, the CWT analysis shows that the quater-diurnal (i.e., first overtide) variability is not significant.The modulation of the CWT power is altogether explained by lunar and solar changes (see e.g., Pugh (1987), for details).However, very little CWT power directly occurs at subtidal scales, which is linked to a very small subtidal (low-pass) SSH residual.In order to further analyze the variability at subtidal scale, we therefore calculate the CWT power for the tidal range.
By definition, results for the tidal range (Figure 6c) show no variability at the tidal bands (i.e., diurnal, semidurnal, quater-diurnal).The tidal range variability is entirely contained at subtidal bands with CWT power peaking between 8 and 30 days and showing little change in time.The two high CWT power bands correspond to fortnightly and monthly bands.These periodicities are unsurprising and typical of variations in tidal range due to (a) the spring neap cycle occurring twice per synodic month (due to the alignment of Earth Moon Sun and seen as M 2 /S 2 interaction); (b) the cycle between tropic and equatorial tides each occurring twice a tropical month (due to the moon declinational cycle and seen as O 1 /K 1 interaction); and (c) the apogee/perigee cycle occurring once per anomalistic month (due to the varying lunar distance to the Earth and seen as an M 2 /N 2 interaction.)We refer the interested reader to Kvale (2006) and Pugh (1987) for more thorough descriptions of tidal range patterns.However one potential weakness of the wavelet analyses is that they do not enable us to precisely pinpoint the tidal mechanism beyond semidiurnal, diurnal, fortnightly, and monthly variability.In particular, given the known tidal constituents observed and the form factor in the Pearl River Delta, we expect both spring-neap and tropic-equatorial cycles to significantly vary, as confirmed by the CWT power for SSH.The observed seasonal variability for the tidal range (cf. Figure 4b) at seasonal scales cannot be picked by the CWT since the largest scale in the analysis is smaller than the 6 months needed to identify it.

Wavelet Salt Intrusion Length (L s )
The CWT wavelet analysis of L s (Figure 7) reveals at which scales the L s varies and how these scales evolve in time between October 2007 and September 2008.Variability occurs both at tidal and subtidal scales.At tidal scales, variability is mostly only significant at the diurnal scale, with semidurnal variability only being significant during two short events in October and March.Similar to the SSH, CWT power at tidal scales is modulated fortnightly, corresponding to fortnightly modulation in the amplitude of the tidal salinity variations (i.e., seen as gray area in panel b).Such behavior is linked to fortnightly modulation of diurnal tidal amplitude (tropic-equatorial cycle).An evolution in the variability at tidal scales throughout the year is also apparent.During the low flow   Q r season, that is, before early June, the diurnal band is well marked and almost permanent, even though there is a decay from late March on-wards (i.e., onset of the monsoon).In June, concomitant with the main peak in Q r , tidal (both semidurnal and diurnal) variability vanishes and is eclipsed by variability at longer time scales.After June, concomitant with a high Q r season, CWT power at the diurnal time scale is more intermittent than during the low Q r regime, vanishing for each Q r local (in time) peak.Variability at subtidal scales is more persistent than at tidal scales.Most of the energy is localized at the fortnightly scale and the 2-month scale (i.e., 60 days) as indicated by the darkest red, although there is enhanced variability between the beginning of May and end of June at the monthly scale.The 2-month scale variability most likely reflects the influence of Q r since there are no tidal processes at that scale and atmospheric forcing that could play a role at this scale is not included in our model.Even though fortnightly scales are easily explained by tidal processes (e.g., spring neap cycling in a mixed tidal regime), the Q r has variability at similar scales and it is not possible to fully distinguish which (river vs. tide) causes the L s variability based solely on the CWT power.To infer the influence of the different forcings in the temporal patterns at different scales of L s we use the cross-wavelet analysis.

River Discharge (Q r ) Versus Salt Intrusion Length (L s )
The WTC analysis between Q r and L s (Figure 8c) presents a complex pattern.High coherence occurs at different scales for different dates (e.g., different scales for the peak in February vs. the peak in April), which indicates changes in how the response of L s to increased Q r evolves over time.The orientation of the arrows (i.e., the phase angle) shows a gradual variation of about 45° for a given band suggesting locked behavior however it changes greatly between bands.
During the low Q r season until March, high coherence is relatively isolated, occurs across scales and is linked with a single pulse Q r leading the L s response.During the 2 week period that the Q r pulse in February seems to last, the phase angle is fairly constant for a given band however it changes between bands (i.e., 265° for 4 days band, 180° for 8 days band and 260° for the 14 days band).Since the phase angles correspond to temporal lags, we obtain different temporal lags for each one of the bands (i.e., about 3 days for the 4 days band, about 4 days for the 8 days band, and about 10 days for the 14 days band).In March-April, still during the low Q r season, the coherence between Q r and L s is also statistically significant at scales between 3 and 8 days with the phase angle evolving from 225° to 270°, meaning that the response of the L s to the Q r becomes slower as the monsoon onsets.
These phase angles at these periods are equivalent to a temporal lag evolving from [1.8, 5] days in March to [2. 25,6] days in April.This slower response of L s could be linked to the negative feedback of increased stratification delaying the impact on L s .Due to the inverse relationship of between Q r and L s , an instantaneous response of L s to changes in Q r will appear in the WTC as antiphase, that is, 180° and arrows pointing left.The 260° phase angle that we obtain for the fortnightly band is equivalent to a 80° lag relative to the antiphase and would correspond to a time lag of about 3 days relative to the antiphase, that is, between the Q r peak and the decrease in L s .The wavelet analysis also enables us to identify the time lags for all time scales at play, thus providing a disentangled view of the temporal evolution of the salt intrusion response time.
Once the monsoon starts, mid-May to mid-June, L s and Q r show high coherence for the scales of 10-20 days with phase angles decreasing from 225° to 180°, that is, time lags evolving from [6.25, 12.5] days in mid-May to [5,10] days in mid June.For the maximum Q r , Q r and L s show antiphase behavior, and therefore increases in Q r result in almost instantaneous decreases in L s .There is very high coherence around 9 June for the 2-4 days band with a phase angle of 270° meaning that the L s responds to changes in Q r at these periodicities with a delay of between a 1.5 and 3 days, given the inverse relationship of Q r and L s the delay between the increase in Q r and the reduction in L s is actually between 0.5 and 1 day.High coherence is also found for the band of 3-5 days between 15 July and 18 August again with a 270° phase lag (i.e., the equivalent temporal lag at this periodicity is between 2.25 and 3.75 days again, given the inverse relationship between Q r and L s these correspond to temporal lags between 0.75 and 1.25 days relative to the antiphase).Monthly to 2 months scales present high coherence from mid April once the monsoon starts and Q r and L s are mostly out of phase at these scales (i.e., phase lag of 180° and then, between 2 weeks and a month for the response of L s at these scales or instantaneous response with respect to the antiphase).

Sea Surface Height (SSH) Versus Salt Intrusion Length (L s )
The WTC between SSH and L s provides additional information (Figure 9d), in particular for subtidal timescales which do not have sufficiently large variability in comparison with primary tidal scales to pass the significance test for the Cross-wavelet power (not shown here).As expected, the primary tidal scales exhibit high coherence for the entire analysis period (i.e., September to September).Subtidal scales present a less persistent pattern of high coherence indicating that SSH and L s are correlated only at certain times.During the low flow season between October and April, the tidally varying SSH leads the tidal L s variations with a phase angle between 45° and 90° for semidiurnal periodicity and with a phase angle of approximately 45° for the diurnal periodicity.These phase angles would correspond to between 1.5 and 3 hr for semidiurnal variability and to approximately 3 hr (between 1.5 and 5 hr) for diurnal variability.Previous studies for the dry season in the Modaomen have identified a time lag in the response of L s of 5 hr relative to the high slack water (Gong & Shen, 2011) and 4 hr (Liu et al., 2017).From May onwards, SSH and L s are close to be in phase (i.e., phase angle close to zero and arrows pointing right) for both semidiurnal and diurnal signal, which implies that L s responds more quickly and almost instantaneously to changes in SSH at these periodicities.Phase angles increase at the end of the simulation (i.e., 45° for the semidiurnal band and 90° for the diurnal band).These larger phase angles indicate that the L s slows down its response to changes in SSH and it is no longer instantaneous as a lower Q r regime resumes (i.e., response time of 1.5 hr for the semidiurnal band and of 6 h for the diurnal band).For the subtidal bands, the temporal pattern in the coherence between the SSH and L s reflects the transfer of energy between periods.As such between 18 February and 14 April there is some significant coherence at the weekly to fortnightly variability with SSH and L s being out of phase (i.e., phase angle of 180°).This period corresponds to a reduction in the coherence a tidal scales.Between 20 May and 7 July, when the L s decays and moves offshore, the SSH and the L s are out of phase for the fortnightly band (i.e., phase angle of 180° and up to a week for the L s to respond to changes in SSH for the fortnight periodicity).High coherence is also observed at the monthly scales from March onwards, with the SSH and the L s out of phase (i.e., phase angle of 180°indicating a response time of about 2 weeks of L s to changes in SSH for the monthly band).High coherence is also shown between SSH and L s from December to July for the 2-month band with signals being out of phase (i.e., phase angle of 180°) and hence a response time of up to a month of L s to changes in SSH.Since our model does not include atmospheric forcing that could also be responsible for variability at synoptic scales, the subtidal variability is likely reflecting the influence of the river flow on the SSH.The antiphase relationship between SSH and L s at subtidal scales also supports this hypothesis.

Tidal Range Versus Salt Intrusion Length L s
The periods showing highest WTC (cf. Figure 10d) between the tidal range and the L s correspond to those at the fortnightly band (cf. Figure 10d) reflecting the spring-neap and tropical-equatorial modulation of the tidal range and the modulation at the fortnightly scales of the L s .The response of the L s to the tidal range evolves over time.
During the low flow season, from October to March the phase angle between the tidal range and the L s increases from 180° (i.e., tidal range and L s are out of phase) to 225°, meaning that the salt intrusion reacts more slowly to changes in the tidal range in February when the temporal lag is of about 9 days (i.e., phase angle 225°) than it does in October when the temporal lag is of about 7 days (i.e., signals out of phase and phase angle 180°).The results prior to March agree with the results in Figure 4 and with previous results in the Modaomen by Gong and Shen (2011) and provide in addition the temporal evolution of the response the second half of the simulation, from April onwards, the phase lag between the tidal range and the salt intrusion length increases from 135° to 180° (i.e., 5.25 to 7 days for the response of L s to changes in tidal range).The two weeks between the 17th March and the 2nd April have reduced tidal range variability which explains that there is no coherence between tidal range and L s at the fortnightly band.Between 17 March and 14 April, close to the equinox, the only band with high coherence between the tidal range and the L s is that of 4-8 days.The phase angle between tidal range and L s for the 4-8 days band is 270° meaning that it takes 3-6 days for the L s to respond to changes in the tidal range.Additionally high coherence is also found between 24 December and 15 March for the monthly band with a phase angle of 135° (i.e., it takes 11 days for the L s to respond to changes in tidal range at the month scale).
The WTC reflects the quicker response of the L s to variability at fortnightly bandscales under high flow regimes, this same behavior has previously been shown by other studies (e.g., Lerczak et al., 2009;Ralston et al., 2010).
The change in phase angle also points toward a change in the type of estuary and dominant processes.Antiphase means maximum L s for neaps indicative of exchange flow dominated and partially mixed estuaries.On the other hand, in phase would mean maximum L s during springs that is characteristic of salt wedge or tidal dispersion estuaries.Likewise from October to March the phase angle evolves from 180° to 225° reflecting the larger L s during neaps or briefly after neaps and indicating exchange flow dominated or partially mixed estuary.
Conversely from April to September the phase angle increases from 135° to 180° indicating maximum L s after spring tides (i.e., reflecting a salt wedge or tidal dispersion dominated estuary) slowing down until happening at neap tides characteristic of exchange flow dominated estuaries.

Quantitative Summary of Response Times
The

Discussion
We have built and run a FVCOM set up for the Pearl River Delta to analyze the temporal behavior of L s and its response to transient river and ocean forcing.Critically, this study spans both low and high discharge season, whereas earlier studies were limited to the dry season.The model captures the observed estuary behavior and allows us to analyze the complex temporal patterns between Q r , SSH and tidal range, and L s .
The studied estuarine system is under the influence of multiple forcings that operate at multiple intertwined scales.The information that we can gather from the time series analysis (Figure 4) and from the wavelet analysis (i.e., WTC from Figures 8-10) is complementary.Time lags can be extracted from time series (Figure 4) and produce results for precise conditions at a given instant following a perspective that aggregates over forcings and scales.The wavelet analysis provides a different perspective focusing on decomposing into time-frequency domain.This method provides a single framework that can disentangle forcings and scales, and it enables to decompose response times into the time-frequency domain.One advantage is that the WTC automatically identifies the relevant timescales, then determines the response time at those timescales, and how these response times evolve over time (see .The wavelets pick relationships between particular scales.Likewise, if the relevant scales in a forcing evolve, the relevant scales in the response times are also likely to evolve.The limitation on the scales resolved only depends on the frequency and duration of the input data.In our case, this limit is established at daily scales for Q r and hourly scales for tidal forcing.In contrast, most other methods typically require a manual pre-selection of the timescales of interest.For example, calculations of estuarine adjustment time (MacCready, 1999(MacCready, , 2007;;Monismith, 2017;Ralston et al., 2008) rely on averaging over tidal dynamics and therefore implicitly focus on subtidal scales.One limitation of wavelet analyses is the resolution of the frequency domain which do not match exactly tidal constituents and cannot distinguish within harmonics (e.g., M 2 and S 2 ; O 1 and K 1 ).However, the wavelet analysis picks the non-stationary river-tide dynamics not considered by classic harmonic analysis (Guo et al., 2015) and that might cause the fortnightly and monthly variability observed for SSH in Figures 9 and 11b.Another limitation is that decomposition into time-frequency domain is not equivalent to decomposition into specific salt transport processes, especially when processes occur over overlapping scales.Detailed analyses of the salt transport processes that result in the modeled salt intrusion would instead still rely on existing process decomposition (e.g., Lerczak et al., 2006), and their temporal behavior may be analyzed with wavelet analyses again.In the present work, we purposefully limit ourselves to analyses linking salt intrusion to bulk external parameters (SSH, tidal range, Q r ) to provide new and valuable information on the validity of several critical simplifying assumptions of estuarine physics, as discussed in the following paragraphs.Estuaries can transition between regimes depending on external forcing (e.g., tidal forcing, river discharge) as well as on internal estuarine dynamics (e.g., stratification, estuarine mixing) (e.g., Geyer, 2010;Geyer & MacCready, 2014;Hansen & Rattray, 1965).An analysis of the time series (Figure 4) identifies two main regimes linked to river discharge.During the low Q r period before the monsoonal peak, L s values are largest during neaps and the response time from neaps to springs is shorter than that from springs to neaps.Such behavior is characteristic of exchange flow dominated and partially mixed estuaries (Banas et al., 2004;Lerczak et al., 2006).After the monsoonal peak during the high Q r period, the largest L s values occur during spring tides and the response time from neaps to springs is larger than that from springs to neaps (Figures 10 and 11b.3), which is characteristic of tidal dispersion dominated or salt-wedge estuaries (Ralston et al., 2010;Xue et al., 2009).A brief transition between low flow and high flow periods can also be identified and it roughly corresponds to the development of an offshore plume during the monsoonal peak.Our results provide more detailed information on the evolution in time of the temporal response of salt intrusion within and across regimes.In particular, our results reveal a shift from a response to river discharge that slows down in time to a response to river discharge that speeds up in time (see Figure 11b.3 for weekly and fortnightly scales) toward the later stages of the partially mixed regime.A probable explanation is that the effect of increased river discharge on salt intrusion can no longer be mitigated by an increase in estuarine stratification.
The cross-wavelet analyses are particularly useful given the importance of quasi-steady state assumptions in understanding estuarine dynamics.Indeed, the direction of the arrows do provide an easy visual indication of how a response time compares to the forcing timescale.When pointing right for response to tidal forcing, the response is in phase, thus likely instantaneous, and the quasi-steady state assumption is likely valid.When pointing left for response to river forcing, the response is antiphase, thus likely instantaneous, and the quasi-steady assumption is likely valid.During the partially mixed regime there is limited variability associated to the river.Q r and L s show significant coherence at event and seasonal scales, however the quasi-steady assumption is only valid at monthly and longer seasonal scales (i.e., arrows pointing left for long periods but deviating from pointing left for scales shorter than monthly in Figure 8, see also Figure 11b.4).The response across scales to the tidal forcing evolves through time but remains non negligible when compared to the scales of the forcing ).The quasi-steady assumption is hence unlikely to be valid for tidal forcing except for diurnal and semidiurnal scales during the brief transition period when the response time is either instantaneous or minimum.
The leading theory about the dynamical response of the salt intrusion within high/low flow regimes is that of estuarine adjustment theories.Adjustment times in the literature (Lerczak et al., 2009;MacCready, 1999MacCready, , 2007) ) are often given as fraction of the advection timescale T Adv = L s /u o where u o = Q r /a o with a o the average cross-sectional area (15000 m 2 for the Modaomen estuary), L s is the low pass salt intrusion length since the formulations are based on tidally averaged dynamics.For the Modaomen estuary, we take Q r to be the fraction of the discharge entering the Modaomen estuary, that is, 31.85% of the total flow from the West and North rivers (Gong & Shen, 2011).For the dry season, the time lags between the changes in Q r and the response of the L s that we obtain directly from the time series analysis (calculated between the starting point of the increase in Q r and the start in the decrease in L s ) are approximately 3 days.This broadly matches with the information provided by the WTC for fortnightly scales (those comparable to the scales considered by estuarine adjustment theory, Figures 8  and 11b.3 for Q r ).This 3 day lag is comparable to the values found by Liu et al. (2014) in their wavelet analysis (i.e., 3.5 days for 2007/2008 dry season) and by Gong and Shen (2011) in their time-series analysis (i.e., 2.5 days for a similar Q r ).The 3 day lag in the response of L s to changes in the Q r compares best with the analytical theory of MacCready (2007) for estuaries where tidal processes dominate given by the expression T Adj = T Adv /2 = L s /2u o (i.e., 0.7-4.4days).This conflicts with decomposition of salt transport Gong & Shen (2011), which indicated dominance of advective transport and exchange flow driven transport.A similar decomposition applied to our numerical results (not shown) confirms the earlier results of (Gong & Shen, 2011).This therefore seems to indicate that the expression that results in the best fit does so for the wrong reasons.The likely explanation for this discrepancy is that assumptions postulated to derive the time adjustment expression are not satisfied.These simplifying assumptions include simplified geometry and bathymetry, and the "linearization assumption" (i.e., the variations of the salt intrusion length are small compared to the salt intrusion length).We note that the latter is unlikely to hold for the Modaomen estuary (see Figure 4 for example).
Beyond helping to quantify response times (see Figure 11b), the WTC shows lower coherence values for Q r (Figure 8) than it does for SSH (Figure 9) and tidal range (Figure 10).This suggests that the response time is more strongly controlled by the tide, be it the tidal range or the SSH, than by the Q r .Indeed, the complex pattern of the wavelet coherence between the Q r and L s , showing high coherence between Q r and L s for abrupt changes in Q r (early February and also from mid-March onwards), suggests that Q r has a strong influence on the magnitude of L s .This important result aligns with those of Monismith (2017) and their non-linear model in that the timescale of the response does not the depend on Q r but rather on the current value of L s that itself reflects the history of the flow.Since the magnitude of L s itself depends on the flow, the timescale may appear to depend on the flow that is, quicker response for increasing in magnitude Q r and slower response for decreasing in magnitude Q r (Figure 4).Altogether, our wavelet results show that it is not only Q r that will determine the response in L s but also the SSH, the tidal range, and the moment in time, not only in terms of regime but also within the regime due to the temporal evolution of the estuary response (i.e., see Q r pulse in early February during spring tide when the L s recovers toward baseline levels versus Q r pulse in March during neap tide when there is no such recovery, cf.Figures 8 and 11b).
Whilst we use the Pearl River Delta as a case study for a precise period of time, our results can be extrapolated to other years and other areas.The wavelet analysis indicates largest coherence between tide and L s with Q r being significant punctually for weekly and fortnightly scales, and for monthly and larger scales given the bimodal distribution of the monsoon.Furthermore, the year 2007/2008 in the Peal River has been shown to be an average year in terms of Q r (Hong et al., 2020) and so our results are likely to also be relevant for other years.The unique hydrodynamics conditions and bathymetry of different estuaries make that they are likely to show different salt dynamics.For example, the tidal behavior around equinoxes (reduced tidal range here) is not universal.We would therefore expect a different behavior in systems with primarily semi-diurnal tides for which tidal range would be larger at the equinoxes (Pugh, 1987).But other estuaries with mixed tides (e.g., Columbia, San Francisco Bay) may have similar behavior to the one we depict in the Modaomen Estuary.Transient behavior and multiscale temporal response are also likely to be found in other estuaries worldwide with already pronounced seasonal variability (e.g., Gironde, Hudson, Ganges, etc.) and provide further insight into the dynamics of the estuary.
The wavelet approach can also help disentangle overlapping scales in the case of temperate climates with more gradual transitions between seasons and/or discharge events in all seasons and where response time to Q r is likely to not be instantaneous (e.g., British rivers).

Conclusion
We used a 3-dimensional unstructured hydrodynamic model to disentangle the multiscale temporal pattern of salt intrusion in estuaries subject to transient forcing.The model was calibrated and validated against available in situ data and its results analyzed using wavelets analysis to extract time-frequency patterns in the behavior of salt intrusion.
We provide a new study that combines both low and high discharge seasons.We show that the response of salt intrusion to tidal and riverine forcing operates at multiple interwoven scales and evolves through time both intra and inter-regime.The salt intrusion response time is mostly determined by the tide, whilst the influence of the river is noticeable in the modulation of the response times (i.e., the different regimes or monthly scales and longer) and in the magnitude of the L s .During the dry season the estuary is partially mixed and transitions toward a salt wedge state that is attained after the major river discharge in July.During the transition period the estuary responds almost instantaneously to changes in Q r , SSH and tidal range.The time response of L s to changes in SSH and tidal range evolves over time, the response is first slow (until the equinox) and then quick for the partially mixed regime and then slower for the salt wedge regime.The response to spring and neap is also asymmetric for the partially mixed stage, strongest L s are obtained for neap tide and the estuary responds quicker from neap to spring than from spring to neap.Conversely for the salt wedge phase the largest L s are obtained for spring tides and the estuary responds quicker from spring to neap than it does from neap to spring.
With regards to water resources management, mitigation approaches often rely on altering the release of freshwater.In the Modaomen estuary, this occurs via: (a) release of freshwater at the upstream when the salt intrusion is severe; (b) Reservoir construction along the estuary for freshwater storage during wet season and freshwater supply during dry season and (c) pipeline construction for transporting freshwater from the upstream directly to the region where freshwater supply is threatened.Since the L s response evolves over time, understanding the response of the estuary to changes in Q r is vital to adequately time any intervention.Empirical curves such as  ∝    are particularly attractive because of their simplicity.However, the response in L s is not determined solely by river discharge but also by tidal forcing and the moment in time implying control from antecedent conditions.These empirical curves also assume a quasi-steady system.The wavelet method provides a direct assessment of this assumption, which is found to be not valid except for monthly and longer scales for the river discharge.At all other scales for Q r and for all scales for tidal forcing the system is unsteady.When the quasi-steady assumption is not justified, simple relationships to determine the estuarine adjustment time are again particularly attractive.However, our wavelet analyses demonstrate that linear estuarine adjustment theories are not applicable in the case of the Modaomen Estuary.It is finally important to note that we applied WTC to forcing and L s obtained from numerical modeling, but the method can also be applied to other variables that may be easier to measure or model (e.g., salinity).For instance, the response times obtained by Liu et al. (2014) between Q r and salinity at an upstream station are similar to the ones we obtain here but a more precise correspondence would need to be investigated further.

Figure 1 .
Figure 1.Pear River Delta model.(a) Full model domain, colorbar indicates grid resolution (logarithmic scale).(b) Zoom in the Pearl River Delta region, colorbar indicates resolution (logarithmic scale).(c) Full model domain, colorbar indicates bathymetry.(d) Zoom in the Modaomen branch, colorbar indicates bathymetry, orange line indicates the longitudinal profile, yellow lines depict the cross-sections of interest, orange cross is the reference for salt intrusion distance and the point for which sea surface height is analyzed, and squares and diamonds indicate respectively the elevation and salinity stations.

Figure 2 .
Figure 2. Elevation validation.(a) Amplitude of diurnal constituents (in m).(b) Phase of diurnal constituents (in °).(c) Amplitude of semidiurnal constituents (in m).(d) Phase of semidiurnal constituents (in °).(e) Amplitude of overtides (in m).(f) Phase of overtides (in °).The 1:1 line indicate perfect agreement between the measurements and the model.Symbols above the line indicate an overprediction by the model whilst points under the line indicate underprediction.

Figure 3 .
Figure 3. Salinity Validation for the four stations along the Modaomen branch.Hourly model results are presented in blue.Gray lines represent the daily maximum measured data.Station: (a) Majiao, located at about 23.2 km from the river mouth reference for L s ; (b) Lianshiwan, located at about 21.3 km from the river mouth reference for L s ; (c) Denglongshan, located at about 18.8 km from the river mouth reference for L s , and (d) Dachongkou, located at about 15.6 km from the river mouth reference for L s .

Figure 4 .
Figure 4. Time series or river and ocean forcing and salt intrusion length.(a) Daily river discharge for the three rivers, (b) Sea surface height and tidal range both at the mouth of the Modaomen, and (c) Hourly salt intrusion length for the surface layer and the bottom layer (raw values and low-pass filtered ones).

Figure 5 .
Figure 5. Frequency time analysis of the river discharge (a) Daily river discharge for the three rivers and the total river discharge of the delta, (b) Wavelet-power spectrum normalized by the variance of the river discharge (colorbar).The contour designates the 95% confidence level against red noise.The lighter shade designates the COI.

Figure 6 .
Figure 6.Frequency time analysis for the water level (a) tidal water level and tidal range, (b) Continuous wavelet spectrum for the SSH normalized by the variance (colorbar), (c) Continuous wavelet spectrum for the tidal range normalized by the variance.The contour designates the 5% significance level against red noise.The lighter shade designates the COI.

Figure 7 .
Figure 7. Frequency time analysis for the salt intrusion length (a) river discharge, (b) (Bottom layer) salt intrusion length raw model results in gray line and low-pass filtered in black continuous line, and (c) Wavelet spectrum for the salt intrusion length normalized by the variance (colorbar).The contour designates the 5% significance level against red noise.The lighter shade designates the COI.

Figure 8 .
Figure 8. Crosswavelet Analysis River Discharge versus Salt Intrusion Length: (a) River, (b) (Bottom layer) salt intrusion length raw model results in gray line and low-pass filtered in black continuous line, (c) Crosswavelet coherence for sea-level and salt intrusion length (colorbar).The contour designates the 5% significance level against red noise.The lighter shade designates the COI.The arrows provide information on the phase between the two time series analyzed.Arrows pointing right denote signals in phase, and arrows pointing left denote anti-phase signals.

Figure 9 .
Figure 9. Crosswavelet Analysis Sea Surface Height (SSH) versus Salt Intrusion Length (L s ) (a) Daily river discharge.(b) SSH, (c) Salt intrusion length L s raw model results in gray line and low-pass filtered in black continuous line, (d) Crosswavelet coherence for SSH and salt intrusion length L s (colorbar).The contour designates the 5% significance level against red noise.The lighter shade designates the COI.The arrows provide information on the phase between the two time series analyzed.Arrows pointing right denote signals in phase, and arrows pointing left denote anti-phase signals.

Figure 10 .
Figure 10.Crosswavelet Analysis Tidal Range versus Salt Intrusion Length (a) Daily river discharge, (b) Tidal range, (c) (Bottom layer) salt intrusion length raw model results in gray line and low-pass filtered in black continuous line, (d) Crosswavelet coherence for sea-level and salt intrusion length (colorbar).The contour designates the 5% significance level against red noise.The lighter shade designates the COI.The arrows provide information on the phase between the two time series analyzed.Arrows pointing right denote signals in phase, and arrows pointing left denote anti-phase signals.
phase angle arrows in Figures 8-10 provide the response times of salt intrusion length to changes in river discharge and tidal forcing, and we present a quantitative summary of the response times in Figure 11.Values in Figure 11 are split into diurnal and semi-diurnal variability (Figures 11b and 11l); weekly variability corresponding to periods from two to 8 days in Figures 8-10 (Figure 11b.2); fortnightly variability corresponding to periods from 10 to 20 in Figures 8-10 (Figure 11b.3); and monthly variability corresponding to periods from 20 to 30 in Figures 8-10 (Figure11b.4).Values for the response time are only provided for cross-wavelet results above 95% confidence compared to red noise but do not take into account the coherence value (color in Figures8-10).This means that Figure11does not include information about which forcing is more relevant, even though this is included in Figures8-10.We therefore refer the reader to Figures 8-10 to assess the relevance of the forcings relative to each other via the coherence value.The shaded areas correspond to the uncertainty in time response that arises from looking at a range of scales.

Figure 11 .
Figure 11.Summary of response times (a) Regimes identified based on river discharge and estuarine adjustment times.L s (N) and L s (S) are the salt intrusion length for neap and spring tides, respectively.t r (N to S) and t r (S to N) are the response time from neap to spring and from spring to neap, respectively.(b) Quantitative summary of wavelet-derived response times of L s to changes in the sea surface height (SSH), tidal range (TR) and river discharge (Q r ) at the relevant scale.The scales highlighted by the wavelet are from top to bottom (b.1) tidal scales: where D1 is diurnal and D2 semidiurnal.(b.2) weekly scales, (b.3) fortnightly scales and (b.4) monthly scales.For all of the plots the x axis is the time of the simulation.

Table 1
Aggregated Root Mean Square Error for the Main Tidal Constituents

Table 2
Salinity Validation This work is part of the SUPREME (Salt intrusion: Understanding the Pear River Estuary by Modelling and field Experiments) project, supported by a three-way international funding through the Netherlands Organisation for Scientific Research (NWO, Grant ALWSD.2016.015), the National Natural Science Foundation of China (Grant No. 51761135021) and the UK Research Councils (UKRI).The UK funding is through the NEWTON fund under the EPSRC Sustainable Deltas Programme, Grants EP/R024480/1 and EP/R024480/2.Boundary conditions are from the 1/12 deg global HYCOM + NCODA Ocean Reanalysis was funded by the U.S. Navy and the Modeling and Simulation Coordination Office.