Vertical water renewal in a large estuary and implications for water quality

• A 3Dagemodel is built to quantify vertical water renewal in a micro-tidal estu-


H I G H L I G H T S
• A 3D age model is built to quantify vertical water renewal in a micro-tidal estuary. • The bottom age can reach up to 8 days at the lower-PRE during the wet season. • During the dry season, the maximum bottom age is~1 day and occurs at the upper-PRE. • Saline wedge has impeding effect on vertical water renewal, resulting in high age. • Oxygen depletion in the estuarine is closely related to the vertical water age.

G R A P H I C A L A B S T R A C T
a b s t r a c t a r t i c l e i n f o

Introduction
Estuaries are commonly characterized by the interaction between rivers and the sea, where heterogeneities in environmental parameters are key for producing abundant biological resources and maintaining species diversity (Lohrenz et al., 1990;Sun et al., 2014a). In recent decades, environmental issues in estuaries have become an increasing cause of concerns, especially regarding water pollution and ecological degradation, which are closely related to industrial development, population growth and agricultural intensification (Boesch et al., 2001;Brye et al., 2012;Kemp et al., 2005;Liu et al., 2011;Sheldon and Alber, 2002;Wang et al., 2013). Pollution, eutrophication, and hypoxia occur frequently in water bodies with a limited water exchange capacity, when there is excessive amount of nutrients and organic matters (Curtis et al., 2011;He et al., 2014;Sun et al., 2014a). Water exchanges in both vertical and horizontal directions in estuaries play an important role in estuarine self-purification and self-restoration, and have therefore gained increasing attention over the past decades (Pfeiffer-Herbert et al., 2015;Valle-Levinson et al., 1996;Wei et al., 2007;Zhu et al., 2015).
Previous studies have focused more on the horizontal renewal of water, especially with respect to the water exchange between freshwater and seawater (Kärnä and Baptista, 2016;Ren et al., 2014;Wang et al., 2013). However, vertical water renewal in estuaries is still not fully understood despite its significance to estuarine environments being widely known. Vertical mixing and circulation patterns are attributable to the multiple dynamic conditions within estuaries, including tidal currents, freshwater discharge and density gradient (Codiga, 2012;Geyer and Maccready, 2014). An insufficient vertical renewal can lead to inadequate DO reaeration in the water column, and when excessive oxygen consumption surpasses the oxygen reaeration, low DO near the bottom of an estuary can occur.
In estuaries, stratification is an important factor impeding the vertical exchange of water and dissolved materials, whereas tides promote turbulence and facilitate reaeration. Hypoxia is more likely to occur in micro-tidal estuaries, such as the PRE (Li et al., 2018;Qian et al., 2018;Zhang and Li, 2010) and the Chesapeake Bay (Li and Zhong, 2009;Muller et al., 2016;Stow and Scavia, 2009), and in meso-tidal estuaries, such as the Yangtze Estuary Zhang et al., 2017). In these micro-and meso-tidal estuaries, the tide-induced turbulence is not strong enough to eliminate stratification. In contrast, macro-tidal estuaries such as the Scheldt Estuary (de Brye et al., 2013;Maes et al., 2007;Mialet et al., 2011) and the Chikugo Estuary (Azhikodan and Yokoyama, 2016;Shahidul Islam and Tanaka, 2007) are characterized by strong tidal mixing, which leads to active oxygen replenishment from the atmosphere. Nevertheless, in some macro-tidal estuaries, for instance the Humber Estuary (Uncles et al., 2000) and the Gironde Estuary (Lajaunie-Salla et al., 2017), hypoxia still occurs especially in the turbidity maximum zone due to suspended particulate organic. The present study focuses on micro-and meso-tidal estuaries, where the threat of lower oxygen is more severe, and vertical water renewal is an important factor controlling, or influencing, the occurrence of hypoxia.
The PRE, south China, is a typical micro-tidal estuary, which, along with the wide catchment and coastal region, faces intensive anthropogenic pressure, especially in relation to rapid urban development and a growing population. Incomplete treatment of domestic sewage and industrial wastewater in addition to increased fertilizer usage have seriously threatened the regional environment (Su et al., 2017;Sun et al., 2014a), with hypoxia being one of the main environmental issues in the PRE (Harrison et al., 2008;He et al., 2014;Qian et al., 2018;Yin et al., 2004).
The concept of age can provide useful measurement for the investigation of flow and mass transport processes in relation to water quality (Brye et al., 2012;Gourgue et al., 2007;Kenov et al., 2012;Viero and Defina, 2016). The age of a water parcel is defined as the time interval that has elapsed since the parcel left a specified region, where the age is set to zero. The water age distribution can be calculated by solving a set of partial differential equations based on the Eulerian theory (Deleersnijder et al., 2001;Delhez and Deleersnijder, 2002;Gourgue et al., 2007). One of the distinct advantages of this theory is the simple form of the governing equations, which are readily implemented into the widely used Eulerian framework of numerical simulation (Brye et al., 2012;Kärnä and Baptista, 2016).
The measure of water age has been applied to quantify the horizontal hydrodynamic processes in estuaries (Kärnä and Baptista, 2016;Ren et al., 2014;Wang et al., 2013). Wang et al. (2013) calculated water age via a 3D advection-diffusion hydrodynamic model to investigate water transport processes in the Daliaohe Estuary, north China. Ren et al. (2014) obtained water age distributions in the PRE for various scenarios and compared the effects of barotropic and baroclinic conditions on the age distributions. Kärnä and Baptista (2016) compared different timescales for the Columbia Estuary and developed a regression model for daily maximum and minimum water age predictions. The aforementioned studies were mainly conducted to investigate the water exchange between the freshwater and seawater (Kärnä and Baptista, 2016), or water age as a measure of the horizontal transport process of water, nutrients and pollutants discharged from upstream rivers (Ren et al., 2014;Wang et al., 2013). However, the timescale of vertical renewal is poorly investigated, which means that values and distributions of water age in this sense are generally unknown.
In the present study, the vertical water renewal process in the PRE is investigated using a 3D age model. We focus on the characteristics of water age, factors influencing the timescale, and the impacts on the DO. The age model was built by incorporating a 3D age sub-model into a hydrodynamic model. The model was validated against measured water levels, velocities and salinity data. The vertical and horizontal distributions of water age were obtained, and their seasonal variations were investigated. An attribution analysis of the factors affecting the age was conducted, including the salinity distributions and salinewedge-induced water circulations. Based on the water age, we propose an equation for estimating the DO level in estuarine waters, and apply it to the PRE. In Section 2, a brief introduction of the study area is given before describing the model structure. Section 3 presents the results of the model-predicted salinity and water age distributions. In Section 4, the effects of circulations and the estimation method for DO based on water age are discussed. Finally, the main conclusions are given in Section 5.

Study area
The study area includes the PRE and its adjacent water bodies (Fig. 1,21°30′ N,112°30′ E to 23°00′ N,115°30′ E) and covers 30,000 km 2 . The water depth varies from b20 m in most areas of the PRE and the coastal bays to 90 m at the estuary-sea boundary.
Freshwater from the Pearl River catchment is discharged into the estuary through eight main inlets (Hu-men, Jiao-men, Hongqi-men, Heng-men, Modao-men, Jiti-men, Hutiao-men, and Ya-men, from northeast to southwest, as shown in Fig. 1). The total annual discharge is~336 × 10 9 m 3 /y (Cai and Li, 2011;Zhang et al., 2008). In monsoon conditions, the precipitation varies seasonally over the catchment. The period from April to September is considered to be the wet season, and the dry season is from October to March. Approximately 80% of the annual total discharge into the estuary occurs during the wet season (Cai and Li, 2011;Dong et al., 2004;Zhang et al., 2008).
In the PRE, the tidal system is characterized by a semi-diurnal mixed tidal regime with daily inequalities. The tidal range is relatively small and varies between 0.8 m and 1.7 m (Mao et al., 2004;Wong et al., 2003); thus, the PRE is deemed a micro-tidal estuary. In general, the four major tidal components (M2, S2, K1, and O1) propagate from east to west of the domain, with increasing amplitudes. Among them, the semi-diurnal lunar tide (M2) and the diurnal solar tide (K1) are the first and second largest components, respectively (Sun et al., 2014b).

Data collection
Hydrological data, including water level, velocity, salinity and water quality indicators, including DO and biochemical oxygen demands, are collected. Field surveys are undertaken in 2006 (July and October) and 2007 (March and September). March and July are in the dry and wet season, respectively, while September and October stand for the midflow condition. The velocity was recorded using an Aanderaa current meter at E41, and the salinity was measured at E38 and E42 using a multi-parameter water-quality monitor sonde (YSI). The water level was collected from an existing hydrological station H2. Temperature, carbonaceous biochemical oxygen demand (CBOD) and nitrogenous biochemical oxygen demand (NBOD) were measured from the cruise surveys. The data of DO was obtained from existing literature (Ambrose et al., 1993;He et al., 2014;Hu and Li, 2009;Li et al., 2018;Zhang and Li, 2010). These data are used for model verification, and for estimation of the bottom DO in an age-based approach.
Freshwater discharge at the eight river inlets was recorded by the Guangdong Research Institute of Water Research and Hydropower (Sun et al., 2014a); the amplitudes and phase-lags of tidal components were obtained from the global tidal model TPXO7.2 (Egbert and Erofeeva, 2002). Meteorological parameters were extracted from QSCAT/NCEP blended data, in which the 6-hour time series were used. These data are used to set up the hydrodynamic model.

3D hydrodynamic model
The Delft3D-FLOW numerical model was set up to simulate the hydrodynamic processes in the PRE. The model was developed by Deltares (2014), and has been widely applied to estuaries and coasts for investigating hydrodynamical problems (Lin et al., 2017;Lin et al., 2015;Luijendijk et al., 2019;Zhang et al., 2018). In the model, the Reynolds-averaged Navier-Stokes (RANS) equations are solved for incompressible fluids under the hydrostatic and Boussinesq assumptions, and the convection-diffusion equation is solved for water and material transport. The baroclinic effect on the coastal flow, which is a key hydrodynamic factor in estuarine studies, has been included in this model. The k-ε turbulence closure model is used to predict the vertical eddy viscosity and diffusivity. Horizontally, an orthogonal curvilinear mesh is applied to represent the estuarine and coastal bathymetric features. In the vertical direction, the σ-coordinate system is used to fit the varying bathymetry. The ADI (alternative direction implicit) scheme is adopted to solve the equations with a second-order accuracy in time and space. The water level, velocity, salinity, diffusivity, and Richardson number distributions can be obtained by running the model. The model-suite is open-source, and thus it is straightforward to integrate the age model with the basic hydrodynamic model.

Water age model
The focus of the present study was the process of water masses leaving water surface of an estuary. The water age is defined here as the time interval that has elapsed since a water parcel left the free surface, and represents the timescale of vertical transport of the water parcel in the estuary. Fig. 2 shows a schematic diagram of a typical dynamic process of water age evaluation. For clearness and simplification, the trajectory of a water parcel is shown in the Lagrangian description. The age of the water parcel is regarded as zero when located at the free surface, i.e. A or E in Fig. 2. If the parcel leaves the surface, its age starts to  Schematic diagram of water age for vertical water renewal, whereby the age of a water parcel accumulates until re-touching the surface. At surface points A and E the age is zero (a = 0), while at points B, C, and D, the age values are δ 1 , δ 1 + δ 2 , and δ 1 + δ 2 + δ 3 , respectively. The movement of the water is driven by the river flow and tides, and influenced by density-induced currents.
increase and the age value accumulates over time (e.g. from B to D). Even in the case that the parcel touches seabed, the parcel age still increases (from C to D), because the seabed is assumed to be impermeable. If the parcel moves back to the surface, its age is reset to zero. When it re-enters the water body, the timing of age re-starts from zero again rather than from its past age value. It should be noticed that a water parcel do not necessarily touch the seabed, i.e. point C is not always at the seabed, and in such conditions, the age value can still be calculated (Fig. 2).
In the Eulerian approach, the general kinetic process of age can be described by two partial differential and one algebraic equations (Deleersnijder et al., 2001), as expressed by Eqs. (1a), (1b) and (1c): where a is the water age; α is the age concentration; C is the tracer concentration referring to the proportion of traced water that has touched the surface; u * ¼ ðu; vÞ and w are the horizontal and vertical velocity components, respectively; D h and D v are the horizontal and vertical diffusivities, respectively; ∇ h = (∂/∂x,∂/∂y) is the horizontal vector differential operator (h-Del operator).
To calculate the water age based on Eqs. (1a), (1b) and (1c), four kinds of boundary conditions are specified, i.e. the water surface, seabed, landward and seaward boundaries. At the water surface, where water parcels have touched the atmosphere, the concentration is set to unity, C = 1, while the age concentration is set to zero, i.e. α = 0. On the seabed, a gradient-free Neumann boundary condition is adopted for both C and α with the assumption that the bed is impermeable. The vertical water age in the freshwater-seawater mixing region is not sensitive to the landward and seaward boundaries.
The age model was built by modifying the standard advectiondiffusion component in the Delft3D-FLOW model according to Eqs. (1a) and (1b) and imposing the corresponding boundary conditions, as mentioned above. The basic hydrodynamic parameters obtained from the hydrodynamic model, such as water level, velocity and diffusivity, were used to drive the water age evolution. The age model was then integrated with the hydrodynamic model to simulate the kinetic process of water age.

Model set-up
To better reflect the features of bathymetry in the PRE, the orthogonal curvilinear mesh of the model was refined locally in the estuary: the size of the mesh ranged from 200 m inside the PRE to 1.5 km at the seaward boundary. Main islands could be identified within this mesh system, where an impermeable wall condition was applied. With regard to the σ coordinate system, 10 layers of equal thickness were used. The time step was set to 5 min.
In estuaries, the water age and environmental indicators are influenced by both regular driving forces, including astronomical tides, riverine discharge and wind, and extreme conditions, such as flooding and storm. This study focuses on the 'normal' conditions in relation to the long-term aquatic and hydrodynamic characteristics of estuaries. The effect of tides and riverine discharge are mainly investigated, while the wind is taken into account as a background force in the present study. Issues on flood and storm are not the focus of this study. Because these extreme events would significantly increase the vertical mixing, leading to a very small water age. In such conditions, the vertical renewal is not the control factor for hypoxia occurrence.
At the eight river inlets, the monthly-averaged flow rates were specified. Along the seaward boundary, the time-varying water level distributions were specified according to eight tidal components, including four major components (M2, S2, K1, and O1) and four minor components (N2, P2, K2, and Q2). The amplitudes and phase-lag degrees of these tidal components were interpolated to each grid point along the seaward boundary. A time series of meteorological parameters with an interval of 6 h was applied at the surface of the study domain. The salinity was set to zero at the river inlet, whereas at the seaward boundary the salinity was set to 35 according to the measured data (Qiu et al., 2019;Wei et al., 2016).
The numerical simulations were started from the rest, i.e. zero velocity. Through trial and error, it was determined that the tidal range reached a stationary status after 10 tidal cycles, while the salinity distribution required around 3 months of simulation to reach equilibrium from an initial condition of uniform salinity. Therefore, a 3-month spin-up time was set for every hydrodynamic model simulation. The initial water age was set to zero over the entire computational domain.

Bottom DO estimation
We propose an estimation method based on the water age to investigate DO levels in estuarine waters. Eq. (2) (described fully in Appendix as A3) is used to link DO with water age: where a b is the bottom water age (day), DO s and DO b are surface and bottom DO concentrations (mg O 2 L −1 ), respectively, CBOD and NBOD are carbonaceous and nitrogenous biochemical oxygen demands (mg O 2 L −1 ), respectively, and K C and K N are synthesis oxidation and nitrification rates (day −1 ), respectively. It can be seen from Eq. (2) that the surface-bottom DO depletion is highly related to both the water age and oxygen demand. The novelty of Eq.
(2) is that it explicitly includes the hydrodynamic characteristics of vertical renewal, in the form of water age, and that it is linked to the evaluation of DO. The surface DO was calculated using the saturated dissolved oxygen (DO sa ) and the surface DO saturation rate (r s ) as expressed by Eq. (3): where S is the salinity and T is the water temperature at the surface layer. The equation used to calculate DO sa followed the latest technical guidelines for surface water environmental impact assessment (HJ2.3, 2018). The range of r s was calculated based on the measurements of He et al. (2014). More details regarding parameter definitions and calculation approaches are provided in the Appendix. Parameters and variables for the kinetic processes of DO were collected from existing literature (Ambrose et al., 1993;He et al., 2014;Hu and Li, 2009;Li et al., 2018;Sun et al., 2014a;Zhang and Li, 2010) and survey data (Sun et al., 2014a), which are listed in Table 1. In using the simplified model, based on Eq.
(2), the DO value varies with depth and time, while the levels of CBOD and NBOD are regarded as background environmental conditions that do not change in this case. The measured data indicate that the surface-to-bottom differences in the CBOD and NBOD values are on average 10% and 20% of their maximum ranges, respectively. Thus, the error introduced caused by this assumption can be regarded as small and acceptable. The salinity and bottom age values were taken from the results of the present study (Section 3).

Model verification
The hydrodynamic model predictions were verified against field measured water levels, velocities and salinities, as shown in Fig. 3 and Table 2. The comparison between the model predicted and field measured water levels are shown in Fig. 3a and b for spring and neap tides, respectively. It can be seen that the tidal amplitude was~1.5 m during spring tides, whereas it was b1 m during neap tides, which agrees with the previous study (Mao et al., 2004). The RMSE (root mean square error) was 0.16 m and the NSE (Nash-Sutcliff efficiency) was 0.93. Hence, the predicted water levels generally agreed with measured water levels in both spring and neap tides. Therefore, the present hydrodynamic model can be used to simulate the tidal propagation and fluctuation processes in the PRE. Fig. 3c and d show the tidal current speed and direction, respectively. Comparisons are made for both depth-averaged velocity and threelayer velocities (surface-, middle-, and bottom-layers). Fig. 3c shows a positive current speed for flood flow and a negative value for ebb flow. For both model-predicted and filed-measured results, the ebb current speed (−1 m/s) was larger than the flood current speed (0.6 m/s) which was caused by the combination of tides and riverine discharge. The current direction was a north-based geographical azimuth (Fig. 3d). The flood flow was north-directed whereas the ebb flow is south-directed, and the transition between flood and ebb flows was rapid. The bidirectional property in this reciprocating flow was mainly caused by the shape of the coastline and the topography of the PRE. We found a good agreement between the model-predicted and fieldmeasured velocity values. The RMSE values of flow speed and direction were 0.18 m/s and 30.9°, respectively, and the corresponding NSE values were 0.84 and 0.89.
The time-varying salinity values at the upper-and lower-PRE during the spring and neap tides are shown in Fig. 3e and f, respectively. Both surface and bottom salinities are plotted to show the vertical salinity differences. The salinity value varied from 10 to 30 at the lower-PRE station (E38, near the outer sea), whereas at the upper-PRE station (E42) the salinity values were generally b10. The vertical salinity difference at station E38 was larger than of that at station E42. The salinity varied with the tides at both stations. During the spring tides, the vertical salinity difference disappeared intermittently over the tidal cycles (Fig. 3e), whereas in the neap tides, the difference increased and lasted over a longer period, especially at station E38. The RMSE and NSE of salinity were 2.71 and 0.95, respectively. It can be seen that the model predictions positively captured the measured salinity variation patterns, which indicates that it is capable of reproducing the major features of hydrodynamic and mass transport processes in both the vertical and horizontal directions. Fig. 4 shows the salinity distributions at the surface and bottom layers of the PRE. A large discharge of freshwater significantly affected salinity gradients in both wet and dry seasons. Since the major riverine inlets are located along the west coast, the salinity level increased from the northwest to the southeast, ranging from almost zero to~35. The salinity spatial distributions varied considerably between the wet and dry seasons. During the dry season, a sharp salinity gradient occurred to the northwest of the Neilingding (NLD) Island ( Fig. 4a and c), while during the wet season the gradient front moved seaward, generally southeast towards the island ( Fig. 4b and d). As a result, in the northern area of NLD Island the wet-season salinity value was almost below 5, especially at the surface layer. These salinity distributions are consistent with the field observations reported by Ji et al. (2011). Between the surface and bottom layers, the salinity changed only slightly during the dry season ( Fig. 4a and c). However, during the wet season, there was a significant variation in salinity between the two layers. As shown in Fig. 4b and d, the surface salinity was generally lower than the bottom salinity, especially in the region between NLD Island and the outer sea, and there was a landward salt-wedge along the deep trench. The maximum bottomto-surface salinity difference was~20 during the wet season. The spring-neap tide variation can affect the vertical salinity distribution, while the effect is not very significant. The salinity spatial distributions and seasonal variations were the dominant phenomena in the PRE in the present study, and this can be attributed to the horizontal and vertical mixing between freshwater and seawater.

Water age distribution
The bottom-layer water ages during both dry and wet seasons are shown in Fig. 5a and b, respectively. It can be seen that significant variations in both magnitude and distribution patterns occurred during the two seasons. Firstly, the maximum bottom age varied considerably between the dry and wet seasons. The bottom age was generally b1 day during the dry season, whereas it reached up to 8 days during the wet season. Furthermore, the bottom age distributions in the two seasons were almost opposite to each other. During the dry season, the maximum age appeared at the upstream end of the PRE (~1 day near Humen inlet), and decreased to~0.25 days near NLD Island, before decreasing to~0.1 days in the area adjacent to the open sea. In comparison, during the wet season, the bottom age increased from the upstream end of the PRE to the open sea, with the maximum age being up to 8 days south of the Lantau Island. This indicates that the wet-season vertical renewal at the lower reach of the PRE was weaker than of that at the upper reach. In addition, a common feature in both seasons was that there were significant bulges in the age-contours despite the opposite directions during the dry season (seawards) and wet season (landwards). Fig. 5c and d show the longitudinal and vertical age distributions along the middle deep trench of the PRE. The longitudinal profile was extracted from the estuarine upstream end to the outer sea, as shown in Fig. 5a and c. The value of water age generally increased from the free surface to the seabed. As a result, the maximum age occurred at the bottom of water column. During the wet season, the vertical variation of water age was significant in the downstream region of the PRE and outer sea. At the upper layers, the age contours were generally parallel to the free surface, which may have been due to vertical mixing. However, at the bottom layer, the age contours followed the shape of the undulating bed. The age value was up to 8 days near the bottom regardless of whether the depth was 20 m or 30 m. This indicates that vertical mixing was not the dominant factor controlling the vertical age distribution. Fig. 5b and d show that the larger values of bottom age in the PRE were distributed along the main trench from the outer sea towards the inner estuary. This reveals that the bottom age was influenced heavily by the horizontal movement at the bottom layer. During the  dry season (Fig. 5c), the pattern was similar with the exception that the maximum bottom age was at the upper reach of the PRE and was onlỹ 1 day. Therefore, we interpret that both the vertical and horizontal hydrodynamic processes contributed to the distributions of water age.

Vertical mixing and effect of density stratification
Since we define the age of a water parcel as the time interval that has elapsed since the water parcel left the surface, which represents its vertical transport in the estuary. The vertical mixing should have a direct impact on water age. Fig. 6a shows the vertical eddy diffusivity and bottom along the longitudinal profile of the PRE. This shows that a large diffusivity value coexisted with a small age along the estuary, and vice versa (i.e. there was an antiphase correlation between these two variables). During the wet season (red lines in Fig. 6a), the vertical eddy diffusivity, D v , was in the order of magnitude of~10 −2 m 2 /s near the riverine inlet, whereas the bottom age was~10 −1 day. In contrast, at the estuary mouth, D v decreased to~10 −3 m 2 /s, and the age increased to~10 days. Similarly, the antiphase pattern can be seen during the dry season as well, although the D v value was smaller and the age was larger in the upper reach than of those in the lower reach. This phenomenon demonstrates that the vertical eddy diffusivity, as an indicator of vertical mixing, had an inverse impact on age value. Fig. 6b and c show the relationship between the vertical eddy diffusivity, D v , and the Richardson Number, R i . The Richardson number (dimensionless) is defined as the ratio of the gravitational buoyancy force to the turbulent shear production (Deltares, 2014;Miles, 2010;Richardson, 1920;Taylor, 1931), and can be expressed by Eq. (4): A positive R i indicates a steady density gradient, which is the case in the present study. When the Richardson number is approximately equal to or larger than unity, the buoyancy force is significant and density stratification usually takes place. If R i is small and positive, the shearinduced turbulence can overcome the buoyancy force and the mixing process would dominate. It can be seen from Fig. 6b and c that the relationships between the D v and the R i , shared a similar pattern in both dry and wet seasons. Within the turbulence-dominated range (i.e. R i ≪ 1) the D v did not change noticeably with the R i , and the value of diffusivity was in the order of magnitude of~10 −2 m 2 /s. However, within the buoyancy-dominated range (i.e. R i N 1) the vertical eddy diffusivity varied inversely with R i . The D v dropped to~10 −3 m 2 /s and~10 −4 m 2 /s during the dry and wet seasons, respectively, and the R i was~1 and 10, respectively. This suggests that vertical mixing was significantly hindered where the water was strongly stratified.
When comparing Fig. 6b with Fig. 6c, it is evident that the turbulence-dominated flow and buoyancy-dominated flow took place in different regions of the PRE during the dry and wet seasons. During the dry season, the buoyancy-dominated flow took place in the upper reach (0-60 km and the black, green, and blue datapoints in Fig. 6b), whereas the turbulence-dominated flow was in the lower reach (60-100 km, the yellow and red data-points in Fig. 6b). The density gradient was derived from the saline wedge, which was located in the upper reach during the dry season due to a small freshwater discharge (Fig. 6d). The salinity difference was about~5 between the surface and bottom layers at the upper reach and R i was~1. At the lower reach, however, there was almost no vertical salinity difference (R i ≪ 1). As a result, vertical mixing was weak in the buoyancy-dominated upper reach, with the D v being 10 −3 m 2 /s, whereas vertical mixing was strong in the turbulencedominated lower reach (D v being~10 −2 m 2 /s) according to the relationship in Fig. 6b.
During the wet season, the riverine freshwater discharge was large, which pushed the saltwater towards the downstream area, and the saline wedge occurred in the middle and lower reaches of the PRE (Fig. 6e). Especially in the middle reach (40-80 km), the surface-to-bottom salinity difference was as large as 20, which resulted in a high Ri of~10 (indicated by the blue and yellow points in Fig. 6b). Consequently, turbulence mixing was significantly suppressed by the salinity/density stratification, and thus the D v reduced to~10 −4 m 2 /s, which led to the bottom age being as high as 8 days (Fig. 6b). In contrast, Fig. 6e shows that the upper reach of the PRE was mainly occupied by freshwater during the wet season due to a large riverine discharge. Therefore, the vertical density gradient was rather low, and the D v was as high as~10 −2 m 2 /s (the black and green points in Fig. 6c), with the bottom age being b0.1 days (Fig. 6a).
To summarize, the saline-wedge-induced density stratification had a direct impeding impact on turbulent mixing, and consequently hindered vertical diffusion, which led to a high water age, especially in the bottom layer. In comparison to the dry season, a large freshwater discharge during the wet season resulted in a significant saline wedge and vertical density stratification, which further led to a higher R i , a lower D v , and hence a large water age. In addition, the buoyancydominated position shifted towards the downstream area due to the effect of relatively large riverine flow during the wet season.   Fig. 7 shows the water circulations and age distributions in the vertical-longitudinal section of the PRE. The results showed that the residual currents were stratified into two layers and that the circulations were highly related to the water age distributions. This is because both were caused by freshwater-seawater interaction, especially the saline wedge. The saline-wedge-induced circulation-loops were located at the upper and lower reaches during the dry and wet seasons, respectively, (Fig. 7a and b), which is consistent with previous investigations (Ji et al., 2011;Wong et al., 2003;Zhang and Li, 2010). In the present study, further analysis of the relationship between the circulation and water age was conducted to investigate the vertical renewal hydrodynamic process. In each circulation-loop, the water flow was seaward at the surface, whereas it was landward at the bottom. The circulation was stronger during the wet season than of that during the dry season, which was related to a larger freshwater discharge during the wet season.

Effect of circulation on age distributions
During the wet season, the bottom landward flow carried the water with a high age of 8 days upstream along the seabed. This indicates how the circulation transported high-age water along the PRE and subsequently increased the vertical renewal time in shallow zones. With regard to the dry season, the situation was similar except that the highest-age water was located in the upper reach and the bottom water age was usually b1 day. In comparison to the seasonally variable discharge, the spring-neap tide variation played a minor role in the circulations and age distributions during both wet and dry seasons in this study.

Age-based estimation of DO
DO is one of the major environmental indicators, and provides a representation of ecosystem health in a water body. Sufficient DO is essential for sustaining the functioning of estuarine and coastal ecosystem. However, low DO or hypoxia has been widely reported in estuaries and coasts worldwide (He et al., 2014;Qian et al., 2018;Wei et al., 2016;Zhang and Li, 2010). Excessive oxygen consumption and poor water exchange for re-oxygenation are the two key processes in forming hypoxia. In existing literature, research attention mainly focused on the biochemical processes (He et al., 2014;Li et al., 2018), while few quantitative investigations have been carried out for the fundamental hydrodynamic processes associated with DO. In the present study, the impact of the vertical renewal process on DO in the PRE was investigated using the concept of age.
The bottom dissolved oxygen, DO b , estimated based on Eqs.
(2) and (3) are listed in Table 3 along with the observed DO b values. During the dry season, the modelled DO b was~7 mg O 2 L −1 in both the upper-and lower-reaches of the PRE and 6.6 mg O 2 L −1 in the lower reach, which agrees well with observed DO b . The oxygen depletion from the water surface to the bottom, ΔDO, was not very large (b1 mg O 2 L −1 ) due to the small bottom age (b1 day). During the wet season, the DO b value at the upper reach of the PRE was as high as 6.8 mg O 2 L −1 . However, at the lower reach of the PRE, DO b decreased to 3.5 mg O 2 L −1 , which we interpret as having being related to the large bottom age of 8 days. This suggests that hypoxia is more likely to form in the lower reach of the PRE during the wet season. Generally, our results agree well with the measured data ( Table 3). The estimated surface-to-bottom oxygen depletion, ΔDO, was 2.7 mg O 2 L −1 in the lower reach during the wet season, and accounted for approximately 3/4 of the observed depletion (3.5 mg O 2 L −1 ). The minor differences between the estimated and observed DO b can be mainly attributed to other oxygen-consumption factors, such as sediment oxygen demand. Therefore, we consider that our estimation method can represent a large proportion of the estuarine oxygen consumption, and can be used to detect possible occurrences of hypoxia, in which the water age is a key parameter for quantifying the vertical water renewal.
(2) shows that the surface-to-bottom DO change is related to the product of biochemical oxygen demand (BOD) and water age (a b ). With regards to micro-and meso-tidal estuaries, such as the PRE, the Chesapeake Bay and the Yangtze Estuary, bottom hypoxia events occur frequently Muller et al., 2016;Qian et al., 2018;Stow and Scavia, 2009;Zhang and Li, 2010;Zhang et al., 2017). In this type of estuaries, salinity stratification usually impedes the vertical mixing and slowdown the vertical exchange process, thus leading to a large water age. In this situation, hypoxia is more likely to occur, even though the values of CBOD and NBOD are not necessarily very high. Therefore, it is important to investigate and quantify the vertical water renewal in these estuaries such that its information can be used for environmental management and eco-biological restoration purposes and hence in related decision-making and policies.
Our equations are still valid for macro-tidal estuaries, however, the vertical age in macro-tidal situations is not key to the environmental issue of hypoxia. For instance, in the Scheldt and Chikugo estuaries, low oxygen events rarely occur (Azhikodan and Yokoyama, 2016;Maes et al., 2007;Mialet et al., 2011), because strong vertical mixing produces a very small water age value. Even in the scenario when BOD is not low, there would be no hypoxia as long as the product of the BOD and the water age is still limited. However, if oxygen consumption is very high (e.g. if there is a large amount suspended particulate organic matter in the turbidity maximum zone), hypoxia could occur, as has been observed in the Humber and the Gironde estuaries (Lajaunie-Salla et al., 2017;Uncles et al., 2000). In such cases, the vertical age was not the limiting factor, and more attention should be given to bed organic materials and their suspension.

Conclusions
Water renewal duration can provide insight to environmental and water quality issues, that can have a significant impact on ecosystem functioning. In this study, the vertical water renewal process, which is significant to DO and thus estuarine health, was investigated in the PRE by building a 3D water age model. The age of a water parcel was defined as the time interval that had elapsed since it left the water surface. The water age distributions were obtained for wet and dry seasons, and the related mechanisms and influencing factors were analysed and discussed. Based on the water age, we developed an equation for estimating the DO level in estuaries. The key findings from this study are summarized as follows: 1) The water age distributions in the PRE varied both spatially and temporally between the dry and wet seasons. During the wet season, the maximum bottom age value reached up to 8 days at the lower reach of the PRE. However, during the dry season the bottom age was rather small, with the maximum age being~1 day at the upper reach of the PRE. 2) Density stratification had an impeding effect on vertical water renewal within the PRE, thus leading to an increase in water age. The age was correlated in antiphase to the vertical diffusivity, while the vertical diffusivity was related negatively to the Richardson number. 3) Within the turbulence-dominated range, the vertical diffusivity remained at a high level, and the water age was very small. However, within the buoyancy-dominated range, the vertical diffusivity varied inversely with the Richardson number. During the wet season, the vertical diffusivity reduced to~10 −4 m 2 /s when the Richardson number was~10. 4) The water circulation pattern played an important role in vertical water mixing. In the saline-wedge-induced circulation zone, water flowed seaward at the surface layer, whereas it flowed landward at the bottom layer. Therefore, the high-age bottom water moved landward along the bed of the PRE. 5) Our findings suggest that the estuary is threatened by a low level of DO at the lower reach during the wet season, with the main cause being the large water age of up to 8 days. During the dry season, the small water age (b1 day) leads to a high DO level. 6) We propose an estimation method based on the water age to investigate DO levels in estuarine waters. Our findings suggest that oxygen depletion from the surface to bottom layers can generally be predicted correctly.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
where CBOD and NBOD are the carbonaceous and nitrogenous biochemical oxygen demands (mg O 2 L −1 ), respectively; k C and k N are the oxidation and nitrification rates at 20°C (day −1 ), respectively, with θ C and θ N being their temperature coefficients; k CBOD and k NBOD are half saturation constants for the oxygen limit of oxidation and nitrification (mg O 2 L −1 ), respectively. Integrating Eq. (A1) from the surface (t = 0) to the bottom (t = a b ) yields Eq. (A2): Eq. (A2) can be rearranged as Eq. (A3):  Zhang and Li (2010) where DO s and DO b are the surface and bottom DO concentrations, respectively, and K C and K N are synthesis oxidation and nitrification rates, respectively, expressed as Eqs. (A4) and (A5): In Eqs. (A4) and (A5), a representative dissolved oxygen, DO, is used to calculate the synthesis reaction rates, and a linear approximation is applied, as shown by Eq. (A6): A predictor-corrector estimation approach was used to calculate the DO b , based on Eqs. (A3)-(A6). The results were compared with a numerical solution for Eq. (A1), with the error of the estimation being b1%. Therefore, the present age-based method was used to estimate the current DO b in this study.