The Response of Lateral Flow to Peak River Discharge in a Macrotidal Estuary

: The Ou River, a medium-sized river in southeastern China, is selected to study the lateral ﬂow response to rapidly varied river discharge, i.e., peak river discharge (PRD). A three-dimensional model based on the Finite-Volume Community Ocean Model is validated by in situ measurements from 15 June to 16 July 2005. PRD, which considers the extra buoyancy and longitudinal momentum in a short time, rebuilds the stratiﬁcation and lateral ﬂow. PRD, compared with low-discharge, generally makes stratiﬁcation stronger and lateral ﬂow weaker. PRD mainly rebuilds lateral ﬂow by changing lateral advection, lateral Coriolis, and lateral-barotropic pressure gradient terms. After PRD, the salinity recovery time is longer than that of the ﬂow because the impact on buoyancy lasts longer than that on longitudinal ﬂow. Longitudinal ﬂow is mostly a ﬀ ected by the momentum transferred during PRD; therefore, the recovery time is close to the ﬂooding duration. However, the lateral ﬂow is a ﬀ ected by the buoyancy, and its recovery time is generally longer than the ﬂooding duration. The lateral ﬂow recovery time depends on transect width, ﬂow velocity and the variation caused by PRD. PRD occurs widely in global small- / medium-sized river estuaries, and the result of this work can be extended to other estuaries.


Introduction
River discharge plays a crucial role in estuarine circulation and material transport [1]. River discharge has been studied in many large river estuaries, such as the Amazon River Estuary [2,3], the Nile River Estuary [4], the Mississippi River Estuary [5,6], and the Yangtze River Estuary [7][8][9]. Due to a lack of connected lakes and tributaries (small drainage basins) acting as sponges, the fluvial flow of small-/medium-sized rivers is likely to form peak river discharge (PRD), i.e., a rapidly varied large-volume flux during a short period under natural conditions (e.g., heavy rainfall or snowmelt) or anthropogenic conditions (e.g., sluice floodwater) [10]. This discharge is different from the steadily increased discharge that occurs during rainy seasons in large rivers (hereafter known as steady flooding) [10]. For example, the maximum discharge of the Amazon River is approximately twice as high as its average river discharge; however, the maximum discharge of small-/medium-sized rivers can be several orders greater than the average discharge with extremely rapid variation [11]. However, according to Milliman and Syvitski [12], the sediment flux from small-/medium-sized rivers has been significantly underestimated, which is also essential to marine geomorphology evolution [13] and biogeochemical processes [14], such as the Trinity River in North America [15], the Tiber River in The ORE is split into the North and South Channels by Ling-Kun Island (LKI). To meet the demand of land use by socioeconomic development, a submerged dike (South Channel Submerged Dike, SCSD) was built in 1979 in the South Channel to silt the South Channel and naturally and gradually form land. The North Channel, as a main channel, is further divided into the northeastward Sha-Tou Passage (STP) and the Middle Passage (MP); it stretches along the North Submerged Dike (NSD), Ling-Ni Dike (LND), and several islands, including Da-Men Island (DMI), Ni-Yu Island (NYI) and Zhuang-Yuan-Ao Island (ZYAI) (see Figure 1b for details).
Fieldwork was carried out in the ORE in the plum-rainy/typhoon season of 2005 ( Figure 1). Tidal elevations (blue triangles in Figure 1) were recorded hourly at three long-term tidal stations (HH, WXT, and ZYA) and several short-term tidal gauges obtained by Temperature Depths (Compact-TDs). Compact-TDs were set to work in a burst mode with a rate of 1 Hz, which records 60 samples within one burst for each hour, and the burst-averaged data was used as the hourly record. Meanwhile, electromagnetic current meters (ECMs) and Conductivity Temperature Depths (CTDs) were used to record the current velocity and salinity (red stars in Figure 1) during the spring tide from 09:00 June 23 to 12:00 June 24 (UTC +0800, day of the year 174. 38-175.50) and the neap tide from 15:00 June 30 to 19:00 July 1, 2005 (UTC +0800, day of the year 181.63-182.79) with one-hour intervals. ECMs and CTDs were also set to work at a sampling rate of 1 Hz, and for each hour they were deployed to measure the profiles of current, temperature, and salinity. The daily average river discharge was obtained from the Xuren hydrologic station in the OR and the Shizhu hydrologic station in the Nanxi River (NR) (red triangles in Figure 1d). A PRD event was captured during the field measurements ( Figure 2b). Fieldwork was carried out in the ORE in the plum-rainy/typhoon season of 2005 ( Figure 1). Tidal elevations (blue triangles in Figure 1) were recorded hourly at three long-term tidal stations (HH, WXT, and ZYA) and several short-term tidal gauges obtained by Temperature Depths (Compact-TDs). Compact-TDs were set to work in a burst mode with a rate of 1 Hz, which records 60 samples within one burst for each hour, and the burst-averaged data was used as the hourly record. Meanwhile, electromagnetic current meters (ECMs) and Conductivity Temperature Depths (CTDs) were used to record the current velocity and salinity (red stars in Figure 1) during the spring tide from 09:00 23 June to 12:00 24 June (UTC +0800, day of the year 174. 38-175.50) and the neap tide from 15:00 30 June to 19:00 1 July 2005 (UTC +0800, day of the year 181.63-182.79) with one-hour intervals. ECMs and CTDs were also set to work at a sampling rate of 1 Hz, and for each hour they were deployed to measure the profiles of current, temperature, and salinity. The daily average river discharge was obtained from the Xuren hydrologic station in the OR and the Shizhu hydrologic station in the Nanxi River (NR) (red triangles in Figure 1d). A PRD event was captured during the field measurements ( Figure 2b).
Two transects along the MP (the red cross-channel lines in Figure 1b,c) are selected for flow and stratification analysis. Transect I is near the observed stations. Meanwhile, Transects I and II represent two typical bathymetric features in the ORE: Transect I is located in the narrower channel with a shallow sill along the channel, and Transect II has a wider width and a submerged dike (NSD), approximately 0.6 m below the mean sea level, on the north. Two transects along the MP (the red cross-channel lines in Figure 1b,c) are selected for flow and stratification analysis. Transect I is near the observed stations. Meanwhile, Transects I and II represent two typical bathymetric features in the ORE: Transect I is located in the narrower channel with a shallow sill along the channel, and Transect II has a wider width and a submerged dike (NSD), approximately 0.6 m below the mean sea level, on the north.
The observed longitudinal flow, lateral flow and salinity near Transect I were considered at four different tidal stages in Figure 3. Note that all the cross-channel plots in this study are landward looking, and the local horizontal coordinates of transects are defined to be positive seaward (perpendicular to the transects downstream, x-axis) and positive to the right (landward looking, northward along the transects, y-axis). Figure 3a,i illustrate that the ebb flow and freshwater began from the south in the North Channel, and longitudinal flow reached the maximum in the Deep Navigation Channel (DNC, labeled in Figure 3j) at maximum ebb tide (Figure 3c). Upstream longitudinal flow also began from the south at low water slack (Figure 3e), and the magnitude was dependent on the bathymetry at maximum flood tide (Figure 3g). The lateral flow seemed discontinuous because the observation stations were not strictly in a straight line (Figure 3b,d,f,h). At high water slack, the southward flow was mostly located at the shoal bottom and the DNC surface (labeled in Figure 3j), and the rest was mostly northward (Figure 3b). Lateral flow was mostly southward at maximum ebb tide, and the water flowed northward near the south head of the transect and the DNC bottom (Figure 3d). At low water slack, the southward flow covered the shoal surface and the DNC, and the other part was northward flow (Figure 3f). Lateral flow was mostly northward at maximum flood tide, and the water near both the north and south heads of the transect flowed southward, although southward flow was weak at the DNC (less than 0.1 m·s −1 ) ( Figure 3h). As the salt wedge was pushed downstream by flooding at some moment during the fieldwork, these  The observed longitudinal flow, lateral flow and salinity near Transect I were considered at four different tidal stages in Figure 3. Note that all the cross-channel plots in this study are landward looking, and the local horizontal coordinates of transects are defined to be positive seaward (perpendicular to the transects downstream, x-axis) and positive to the right (landward looking, northward along the transects, y-axis). Figure 3a,i illustrate that the ebb flow and freshwater began from the south in the North Channel, and longitudinal flow reached the maximum in the Deep Navigation Channel (DNC, labeled in Figure 3j) at maximum ebb tide ( Figure 3c). Upstream longitudinal flow also began from the south at low water slack (Figure 3e), and the magnitude was dependent on the bathymetry at maximum flood tide (Figure 3g). The lateral flow seemed discontinuous because the observation stations were not strictly in a straight line (Figure 3b,d,f,h). At high water slack, the southward flow was mostly located at the shoal bottom and the DNC surface (labeled in Figure 3j), and the rest was mostly northward (Figure 3b). Lateral flow was mostly southward at maximum ebb tide, and the water flowed northward near the south head of the transect and the DNC bottom (Figure 3d). At low water slack, the southward flow covered the shoal surface and the DNC, and the other part was northward flow (Figure 3f). Lateral flow was mostly northward at maximum flood tide, and the water near both the north and south heads of the transect flowed southward, although southward flow was weak at the DNC (less than 0.1 m·s −1 ) ( Figure 3h). As the salt wedge was pushed downstream by flooding at some moment during the fieldwork, these observations cannot reveal the intratidal salinity variation and show only the salinity structure at high water slack (Figure 3i).

Model Configuration
The numerical model used here is the three-dimensional unstructured grid primitive equation Finite-Volume Community Ocean Model (FVCOM) [48,49].  (Figure 1d), and 18 uniformed sigma layers in the vertical direction. The daily average river discharge at the Xuren and Shizhu hydrologic stations was used as the OR and NR discharges (red triangles in Figure 1d). Eight primary tidal constituents (M2, S2, N2, K2, K1, O1, P1, and Q1) were obtained from the NAO.99b tidal prediction system (http://www.miz.nao.ac.jp/staffs/nao99/index_En.html), and three quarter-diurnal tides (M4, MS4, and MN4) were derived from an ECS numerical model [50]. The water level calculated by the above tidal constituents was used as the open boundary forcing (the red curves in Figure 1a,d indicate open boundary). Sea surface wind (hourly and 0.312° × 0.312° grid) was obtained from NCEP Climate Forecast System Reanalysis (CFSR) Selected Hourly Time-Series Products (http://rda.ucar.edu/datasets/ds093.1/). The initial salinity and monthly average salinity at the open boundary were derived from the same ECS model [50]. The summer-observed surface-bottom temperature difference is less than 1 °C in the ORE and has less of an effect on the variation in density than salinity [51]. Thus, the numerical model was simplified with a constant temperature (20 °C). The model was initiated with zero water level and velocity and was first run for 120 days with a fixed river discharge, including both OR and NR (the daily average on 1 June 2005) and the open boundary

Model Configuration
The numerical model used here is the three-dimensional unstructured grid primitive equation Finite-Volume Community Ocean Model (FVCOM) [48,49]. The daily average river discharge at the Xuren and Shizhu hydrologic stations was used as the OR and NR discharges (red triangles in Figure 1d). Eight primary tidal constituents (M 2 , S 2 , N 2 , K 2 , K 1 , O 1 , P 1 , and Q 1 ) were obtained from the NAO.99b tidal prediction system (http://www.miz.nao.ac.jp/staffs/nao99/index_En.html), and three quarter-diurnal tides (M 4 , MS 4 , and MN 4 ) were derived from an ECS numerical model [50]. The water level calculated by the above tidal constituents was used as the open boundary forcing (the red curves in Figure 1a,d indicate open boundary). Sea surface wind (hourly and 0.312 • × 0.312 • grid) was obtained from NCEP Climate Forecast System Reanalysis (CFSR) Selected Hourly Time-Series Products (http://rda.ucar.edu/datasets/ds093.1/). The initial salinity and monthly average salinity at the open boundary were derived from the same ECS model [50]. The summer-observed surface-bottom temperature difference is less than 1 • C in the ORE and has less of an effect on the variation in density than salinity [51]. Thus, the numerical model was simplified with a constant temperature (20 • C). The model was initiated with zero water level and velocity and was first run for 120 days with a fixed river discharge, including both OR and NR (the daily average on 1 June 2005) and the open boundary salinity (the monthly average in June 2005) to reach a quasi-equilibrium state. Then, the model was run with daily average river discharge for another 45 days from 1 June to 15 July 2005, and the results Water 2020, 12, 3571 6 of 20 from the last 35 days were analyzed. The wind wave was nonsignificant in the ORE due to the shelter effect of the island [52] and thus was not considered in this study. This case was marked as the Control Run (Run 0).

Validation
To validate the numerical model, the skill score (SS) is used here, which is the ratio of the root mean square error (RMSE) normalized by the standard deviation of the observation [53]: where X is the variable being evaluated and X is the temporal average. Referring to previous studies [8,54], an SS greater than 0.65 can be acknowledged as an excellent simulation, and an SS less than 0.20 is considered poor. The SS indicates the accuracy of the model result. Furthermore, the correlation coefficient (CC) between the model result and the observation was used to quantify the variation trend: The simulated tidal elevation is compared with the observations at 14 stations in the model domain. The SSs and CCs at each station are given in Table 1. All the values exceed 0.90, indicating that the model can well reproduce the tidal elevation. Comparisons between the observed and simulated water levels at the three long-term tidal stations are shown in Figure 4a-c. During spring tide, most of the SSs for the current velocity are over 0.65, and the CCs are over 0.90 (Table 2). During neap tide, the SSs and CCs are also excellent at most stations, except those in the nonmainstream direction. Taking No. 8 station as an example, visual comparisons between the simulated and observed current velocity/direction during spring tide and neap tide are shown in Figure 4d-g. Table 1. Skill score (SS) and correlation coefficient (CC) between the observed and simulated tidal elevation at different stations and their observation duration in 2005 (UTC +0800).

Station
Observation  A comparison of SSs and CCs between the simulated and observed salinities is summarized in Table 3, which indicates that the model could reproduce the salinity variation. Note that an episode of peak river discharge occurred during the observed spring tide, and the daily average discharge measurement could not accurately reproduce the hourly salt-wedge changes, which has a greater impact on the surface salinity simulation than the bottom (Table 3). Taking No. 8 station as an example, visual comparisons between the simulated and observed salinities during spring tide and  A comparison of SSs and CCs between the simulated and observed salinities is summarized in Table 3, which indicates that the model could reproduce the salinity variation. Note that an episode of peak river discharge occurred during the observed spring tide, and the daily average discharge measurement could not accurately reproduce the hourly salt-wedge changes, which has a greater impact on the surface salinity simulation than the bottom (Table 3) visual comparisons between the simulated and observed salinities during spring tide and neap tide are shown in Figure 5. Overall, this numerical model is able to describe the hydrodynamics in the ORE. neap tide are shown in Figure 5. Overall, this numerical model is able to describe the hydrodynamics in the ORE.

Experimental Design
Based on the ORE model, several experiments are designed to study the response of the lateral flow to a rapid-change discharge during spring tide, which are summarized in Table 4. The daily average river discharge from 1950 to 2008 at the Xuren hydrologic station shows that the duration of a PRD event is usually less than 48 h in the ORE [46]. Thus, the flooding duration is set to 48 h in this study. In contrast to Run 0, wind and NR are ignored in the other cases to focus on the PRD. The OR discharge was given a PRD type with a peak discharge of 5000 m 3 ·s −1 in Run 1 (the solid blue line with squares in Figure 6a), which represents the PRD that occurred during the observation. Run 2 represents a low-discharge condition with 500 m 3 ·s −1 (the solid blue line in Figure 6a), referring to the multiyear average value of 442 m 3 ·s −1 [46]. The maximum PRD of 23,000 m 3 ·s −1 was recorded in July 1952 [46]; however, the maximum discharge from 2005 to 2015 was only 10,900 m 3 ·s −1 , which was recorded in August 2014. Therefore, from Runs 3 to 13, the peak discharge is set to increase from 2500 to 15,000 m 3 ·s −1 with 1250 m 3 ·s −1 intervals to examine the impact of discharge volume on estuarine circulation. Peak river discharge 5000 2 Low-discharge 500 3 Peak river discharge 1250 4 Peak river discharge 2500 5 Peak river discharge 3750 6 Peak river discharge 6250 7 Peak river discharge 7500 8 Peak river discharge 8750 9 Peak river discharge 10,000 10 Peak river discharge 11,250 11 Peak river discharge 12,500 12 Peak river discharge 13,750 13 Peak river discharge 15,000 Water 2020, 12, x FOR PEER REVIEW 10 of 20 The structure of salinity shows distinctive patterns between two cross-sections, and the PRD changes Transect I more than Transect II. Figure 7 shows the salinity structure in Transects I and II in Runs 1 and 2 at four different tidal stages. For Transect I, under the low-discharge condition (Run 2), the freshwater first presents on the shoal water surface at the beginning of the ebb tide ( Figure 7a) and then occupies the entire transect until the maximum ebb tide, while the DNC surface has

Salinity and Stratification
Without a loss of generality, a typical ebb-flood cycle during spring tide (the first ebb-flood cycle in Figure 6a shows as a red line) is described in this section to show how the PRD changes the stratification in the ORE. Under the low-discharge condition (Run 2), the tidally averaged salinity ranges from 16 psu in Transect II to 8 psu in Transect I. The ebb/flood tide average salinity along the section (purple line in Figure 1b along the DNC) is shown in Figure 6b-e (black lines). The isolines are distributed asymmetrically between the ebb (Figure 6b) and flood (Figure 6c) tides at the DNC under the low-discharge condition (Run 2), indicating difference in mixing between ebb and flood tides; the ebb-tidal-averaged isolines are tilted more than the flood-tidal-averaged isolines, which are nearly vertical in the lower layer. In Run 1, the ebb-tidal-averaged isolines become more tilted (Figure 6d), and the flood-tidal-averaged isolines remain vertical in the lower layer (Figure 6e). Regarding the surface layer, the isolines become more tilted after PRD (Figure 6e), indicating that the PRD can enhance only the stratification on ebb tides and the surface layer (above −4 m) on flood tides. The lower layer on flood tides is still well-mixed even under PRD.
The structure of salinity shows distinctive patterns between two cross-sections, and the PRD changes Transect I more than Transect II. Figure 7 shows the salinity structure in Transects I and II in Runs 1 and 2 at four different tidal stages. For Transect I, under the low-discharge condition (Run 2), the freshwater first presents on the shoal water surface at the beginning of the ebb tide ( Figure 7a) and then occupies the entire transect until the maximum ebb tide, while the DNC surface has significant stratification (tilted isolines) (Figure 7b). At low water slack, saltwater intrudes from the shoal bottom (Figure 7c). The intratidal variation in salinity indicates a significant difference in mixing between ebb and flood tides, and the transect is well-mixed at maximum flood tide (Figure 7d). In Run 1, the structure of Transect I salinity has a similar pattern as that in Run 2 at four different tidal stages, but with a smaller magnitude (Figure 7e-h). For Transect II, the saltiest water is mostly located on the part near the north head under the low-discharge condition (Run 2). The freshwater first presents on two heads of the transect surface at high water slack (Figure 7i). At maximum ebb tide, the isolines are mostly vertical, indicating that stratification is not significant, and the saltiest water is mostly located on the north head of the transect (Figure 7j). During flood tide, the saltwater intrudes from two heads of the transect bottom at low water slack (Figure 7k), the water column is well-mixed, and the saltiest water is mostly located on the north head of the transect at maximum flood tide (Figure 7l). In Run 1, the structure of Transect II salinity has a similar pattern as that in Run 2 at four different tidal stages but with a smaller magnitude (Figure 7m-p). For Transect I, the transect-averaged salinity is between 2.2-16.4 psu in Run 2 and between 0.6-15.7 psu in Run 1, and the root mean square (RMS) of the difference between Runs 1 and 2 (Run 1 minus Run 2) is 2.6 psu during the ebb-flood cycle. For Transect II, the transect-averaged salinity ranges between 10.7-23.7 psu in Run 2 and ranges between 8.6-23.5 psu in Run 1, and the RMS of the difference between these two runs is 1.6 psu during the ebb-flood cycle, indicating that PRD has a larger impact on Transect I than on Transect II. transect-averaged salinity is between 2.2-16.4 psu in Run 2 and between 0.6-15.7 psu in Run 1, and the root mean square (RMS) of the difference between Runs 1 and 2 (Run 1 minus Run 2) is 2.6 psu during the ebb-flood cycle. For Transect II, the transect-averaged salinity ranges between 10.7-23.7 psu in Run 2 and ranges between 8.6-23.5 psu in Run 1, and the RMS of the difference between these two runs is 1.6 psu during the ebb-flood cycle, indicating that PRD has a larger impact on Transect I than on Transect II.     (Figure 8b). At low water slack, the lateral flow is southward in the upper layer and northward in the lower layer at the shoal; at the DNC, the lateral flow is mostly northward (Figure 8c). The Transect I shallow sill is a barrier to saltwater cross-channel transport by northward near-bottom lateral flow at low water slack; as a result, the water column is almost well-mixed, and lateral flow is mostly northward at the DNC. At maximum flood tide, lateral flow is mostly northward (Figure 8d). For the difference between Runs 1 and 2 (Run 1 minus Run 2), the pattern of PRD-driven flow is mostly opposite to but smaller than that under the low-discharge condition (Run 2) at high water slack (Figure 8e), and it is mostly southward at maximum ebb tide (Figure 8f). During flood tide, the pattern of PRD-driven flow is also opposite to but smaller than that under the low-discharge condition (Run 2) at low water slack (Figure 8g), and it is mostly southward at maximum flood tide (Figure 8h). DNC, the lateral flow is mostly northward (Figure 8c). The Transect I shallow sill is a barrier to saltwater cross-channel transport by northward near-bottom lateral flow at low water slack; as a result, the water column is almost well-mixed, and lateral flow is mostly northward at the DNC. At maximum flood tide, lateral flow is mostly northward (Figure 8d). For the difference between Runs 1 and 2 (Run 1 minus Run 2), the pattern of PRD-driven flow is mostly opposite to but smaller than that under the low-discharge condition (Run 2) at high water slack (Figure 8e), and it is mostly southward at maximum ebb tide (Figure 8f). During flood tide, the pattern of PRD-driven flow is also opposite to but smaller than that under the low-discharge condition (Run 2) at low water slack ( Figure  8g), and it is mostly southward at maximum flood tide (Figure 8h).     (Figure 8l). This lateral flow structure at high/low water slack has also been reported in studying the influence of submerged dikes in the Yangtze River Estuary by Chen et al. [55]. The Run 1 PRD-driven lateral flow magnitude is mostly less than 0.05 m·s −1 at four different tidal stages in Transect II. For the difference between Runs 1 and 2 (Run 1 minus Run 2), PRD-driven lateral flow is mostly southward in the upper layer (above −3 m) and is mostly northward in the lower layer at high water slack (Figure 8m). At maximum ebb tide, PRD-driven lateral flow is mostly northward, and small magnitude (<0.01 m·s −1 ) southward flows are located on the surface of two transect heads (Figure 8n). PRD-driven lateral flow is southward on the near south transect head surface, and the rest is northward at low water slack (Figure 8o); it is mostly northward, and near the two heads, the bottom of the transect is southward at maximum flood tide (Figure 8p).

Lateral Flow
PRD effect is asymmetric between the discharge rising and falling processes. Figure 9a,b show the difference (color) of the transect-averaged lateral flow between each PRD run and Run 2 (each run minus Run 2) in Transects I and II, which indicate PRD-driven lateral flow. For Transect I, the lateral flow is mostly southward and less than 0.07 m·s −1 in magnitude (Figure 9a). The flooding effect on lateral flow is asymmetric between the discharge rising and falling processes: the maximum speed during the discharge rising process (0.07 m·s −1 ) is slightly larger than that during the discharge falling process (0.06 m·s −1 ) in Transect I (Figure 9b). For Transect II, the lateral flow is mostly northward, and the maximum speed during the discharge falling process (0.08 m·s −1 ) is larger than that during the discharge rising process (0.07 m·s −1 ) (Figure 9b).
Water 2020, 12, x FOR PEER REVIEW 13 of 20 lower at high water slack ( Figure 8i); lateral flow is mostly northward at maximum ebb tide ( Figure  8j). During flood tide, the lateral flow is southward in the upper layer of the southern part and the lower layer of the northern part, and northward flow covers the lower layer of the southern part and the upper layer of the northern part at high water slack in Transect II (Figure 8k), and it is mostly southward at maximum flood tide (Figure 8l). This lateral flow structure at high/low water slack has also been reported in studying the influence of submerged dikes in the Yangtze River Estuary by Chen et al. [55]. The Run 1 PRD-driven lateral flow magnitude is mostly less than 0.05 m·s −1 at four different tidal stages in Transect II. For the difference between Runs 1 and 2 (Run 1 minus Run 2), PRD-driven lateral flow is mostly southward in the upper layer (above −3 m) and is mostly northward in the lower layer at high water slack (Figure 8m). At maximum ebb tide, PRD-driven lateral flow is mostly northward, and small magnitude (<0.01 m·s −1 ) southward flows are located on the surface of two transect heads (Figure 8n). PRD-driven lateral flow is southward on the near south transect head surface, and the rest is northward at low water slack ( Figure 8o); it is mostly northward, and near the two heads, the bottom of the transect is southward at maximum flood tide (Figure 8p). PRD effect is asymmetric between the discharge rising and falling processes. Figure 9a,b show the difference (color) of the transect-averaged lateral flow between each PRD run and Run 2 (each run minus Run 2) in Transects I and II, which indicate PRD-driven lateral flow. For Transect I, the lateral flow is mostly southward and less than 0.07 m·s −1 in magnitude (Figure 9a). The flooding effect on lateral flow is asymmetric between the discharge rising and falling processes: the maximum speed during the discharge rising process (0.07 m·s −1 ) is slightly larger than that during the discharge falling process (0.06 m·s −1 ) in Transect I (Figure 9b). For Transect II, the lateral flow is mostly northward, and the maximum speed during the discharge falling process (0.08 m·s −1 ) is larger than that during the discharge rising process (0.07 m·s −1 ) (Figure 9b). Figure 9. Transect-averaged lateral velocity differences (color, unit: m·s −1 ) between each PRD run and Run 2 (each run minus Run 2) in Transects I (a) and II (b). The black contour line indicates 0.01 m·s −1 . The solid red line is the division of the PRD influence period, and the red dotted line is the peak discharge moment. Figure 9. Transect-averaged lateral velocity differences (color, unit: m·s −1 ) between each PRD run and Run 2 (each run minus Run 2) in Transects I (a) and II (b). The black contour line indicates 0.01 m·s −1 . The solid red line is the division of the PRD influence period, and the red dotted line is the peak discharge moment.

Momentum Balance Analysis
Considering a three-dimensional domain, the lateral Reynolds-averaged momentum difference equation [48] is given as follows: where ∆ indicates the difference between PRD and low-discharge conditions, (u, v) are the along-and cross-channel velocities, respectively, w is the vertical velocity, f is the Coriolis parameter, Pt is the barotropic pressure, Pc is the baroclinic pressure, A m is the horizontal diffusion coefficient, Kv is the vertical eddy viscosity coefficient, x and y are orthogonal in the along-and cross-channel directions, respectively, and z is the vertical axis in the Cartesian coordinate system. In the lateral dynamics equation (Equation (3)), the difference of the lateral advection term (−∆uv x ) and the difference of the lateral Coriolis term (−∆fu) are directly affected by the longitudinal flow, and these two terms together are referred to as the longitudinal flow-related term, where the subscript (x, y, z, t) indicates the partial derivation. PRD drives seaward longitudinal flow and changes the longitudinal flow-related term, and buoyancy taken by PRD also modifies the lateral-baroclinic-pressure-gradient (∆Pc y ); as a result, lateral acceleration changes with them, and the changed lateral flow further alters the remaining terms on the right-hand side of Equation (3). Table 5 shows each transect-averaged lateral momentum term in Runs 1 and 2 after the tidal average, which indicates a PRD effect on lateral momentum; the RMS of the lateral momentum change distribution after the tidal average in the transect is also listed in Table 5 to facilitate the comparison of different term changes in magnitude between Runs 1 and 2. During PRD, compared to Run 2, the lateral acceleration in Transect I decreases from 1.39 × 10 −6 m·s −2 to 0.84 × 10 −6 m·s −2 and that in Transect II increases from −0.69 × 10 −6 m·s −2 to −0.66 × 10 −6 m·s −2 ( Table 5). The RMS of lateral acceleration in Transect I (1.84 × 10 −6 m·s −2 ) is smaller than that in Transect II (7.99 × 10 −6 m·s −2 ), so the impact of PRD on lateral acceleration in Transect II is larger than that in Transect I ( Table 5). The lateral flow in Runs 1 and 2 is controlled by the balance between the longitudinal flow-related term (−uv x −fu), the lateral-barotropic-pressure-gradient term (Pt y ), and the lateral-baroclinic-pressure-gradient (Pc y ) term (Table 5). In Transect I, the difference of the longitudinal flow-related term −∆(uv x +fu), with the RMS of 14.06 × 10 −6 m·s −2 , and the difference in the lateral-barotropic-pressure-gradient (∆Pty), with the RMS of 12.27 × 10 −6 m·s −2 , are two of the most significant terms affected by PRD; the RMS of the remaining terms is less than 10 × 10 −6 m·s −2 (Table 5).

Recovery Time
In this study, taking the low-discharge salinity as a reference, the salinity recovery time is defined in a subtidal manner as the duration from a decrease to 10% to when the salinity increases back to 90% of the reference value. Additionally, taking the subtidal flow under the low-discharge condition (Run 2) as a reference, the recovery time of flow is defined as the duration from when the flow increases to 10% of the largest difference from the reference to when it falls back to 10% of the difference. The transect averaged subtidal salinity and flow are calculated using a 34 h low-pass filter to avoid the tidal signals disturbance [50]. Figure 10a-f show the changes in transect averaged subtidal salinity, subtidal longitudinal flow and subtidal lateral flow after the beginning of PRD; the numbers represent the recovery time and correspond to the line color. The PRD pushes the salt wedge seaward, and the salinity drops; subtidal salinity even becomes 0 during the high discharge period under a large PRD condition (for example, red line in Figure 10a); subsequently, the saltwater intrudes the estuary again, and the salinity rebounds. Different transects have different recovery times under the same PRD. In Run 1, the salinity recovery time shortens downstream from 204.0 h in Transect I to 103.3 h in Transect II (blue lines in Figure 10a,b). The recovery time increases as the peak discharge rises (Table 6), and in Run 13, the time extends to 303.5 h in Transect I and 172.5 h in Transect II (red lines in Figure 10a,b). In particular, the salinity recovery time in Transect I slows down significantly as the peak discharge rises when the minimum salinity is 0 during PRD (Runs 10-13, Table 6). Previous studies [21,22] on the salinity recovery time to steady flooding pointed out that longer salt intrusion length and weaker transect-averaged flow correspond to longer recovery time, and this recovery time is related to steady-state estuarine scales (called the low-discharge condition in this paper). For PRD, a salinity recovery time T salt with an empirical curve (Figure 11a) is given by: where L is the averaged salt intrusion length during recovery time [20], U 0 is the averaged transect-averaged subtidal velocity amplitude during recovery time, and U 0 is averaged amplitude of subtidal velocity during recovery time. U 0 is calculated as [25]: where A is the transect area. R salt , which is the ratio of transect-averaged subtidal salinity variation maximum during PRD to the subtidal salinity under the low-discharge condition (Run 2), is calculated as: where s min is the minimum transect-averaged subtidal salinity during PRD and s low−discharge is the transect-averaged subtidal salinity under the low-discharge condition (Run 2). The recovery time of longitudinal flow is close to the PRD duration, and the recovery time of lateral flow is related to the channel width and the subtidal lateral flow amplitude. For PRD, the effect of the extra momentum by PRD changes with the magnitude of the discharge and is dismissed as soon as the PRD ends; however, the effect of the extra buoyancy continues to increase during the entire PRD and is then slowly reduced. As the subtidal longitudinal flow is more affected by momentum than buoyancy, the recovery time of the subtidal longitudinal flow is close to the duration of PRD in the ORE and slightly changes as the PRD increases ( Table 6). The extra momentum induced by PRD becomes weaker from Transects I to II due to the geometric funneling effect and frictional damping. Thus, in Run 1, the longitudinal flow recovery time decreases slightly from 46.5 h in Transect I to 42.8 h in Transect II (blue lines in Figure 10c,d). As the peak discharge increases, the longitudinal flow recovery time changes slightly (Table 6). For lateral flow, recovery time has some fluctuations with increasing peak discharge; in addition, the Transect I recovery time is longer than that of Transect II in larger peak discharge runs (larger than Run 10) and is slightly shorter (<6.2 h) than that of Transect II in smaller peak discharge runs (Table 6). A lateral flow recovery time empirical curve (Figure 11b) is developed by transect width W and the averaged transect-averaged subtidal lateral flow velocity amplitude during recovery time 0 v , which is similar to Equation (5) but using subtidal lateral flow velocity. Rv, which is the ratio of the transect-averaged subtidal lateral flow velocity amplitude variation maximum during PRD to the subtidal lateral flow velocity amplitude under the low-discharge condition (Run 2), is calculated as: vv v (7) where peak v is the transect-averaged subtidal lateral flow velocity amplitude during PRD and low-discharge v is the transect-averaged subtidal lateral flow velocity amplitude under the lowdischarge condition (Run 2). The lateral flow recovery time Tv is given by:  Table 6. Recovery time (unit: hour) of transect-averaged salinity (S), longitudinal (u) and lateral flow (v) after a PRD with peak discharge (unit: m 3 ·s −1 ) shown in parentheses.

Run 3 (1250)
Run 4 (2500) Run 5 (3750) Run 1 (5000)  The recovery time of longitudinal flow is close to the PRD duration, and the recovery time of lateral flow is related to the channel width and the subtidal lateral flow amplitude. For PRD, the effect of the extra momentum by PRD changes with the magnitude of the discharge and is dismissed as soon as the PRD ends; however, the effect of the extra buoyancy continues to increase during the entire PRD and is then slowly reduced. As the subtidal longitudinal flow is more affected by momentum than buoyancy, the recovery time of the subtidal longitudinal flow is close to the duration of PRD in the ORE and slightly changes as the PRD increases ( Table 6). The extra momentum induced by PRD becomes weaker from Transects I to II due to the geometric funneling effect and frictional damping. Thus, in Run 1, the longitudinal flow recovery time decreases slightly from 46.5 h in Transect I to 42.8 h in Transect II (blue lines in Figure 10c,d). As the peak discharge increases, the longitudinal flow recovery time changes slightly ( Table 6). For lateral flow, recovery time has some fluctuations with increasing peak discharge; in addition, the Transect I recovery time is longer than that of Transect II in larger peak discharge runs (larger than Run 10) and is slightly shorter (<6.2 h) than that of Transect II in smaller peak discharge runs ( Table 6). A lateral flow recovery time empirical curve (Figure 11b) is developed by transect width W and the averaged transect-averaged subtidal lateral flow velocity amplitude during recovery time v 0 , which is similar to Equation (5) but using subtidal lateral flow velocity. R v , which is the ratio of the transect-averaged subtidal lateral flow velocity amplitude variation maximum during PRD to the subtidal lateral flow velocity amplitude under the low-discharge condition (Run 2), is calculated as: where v peak is the transect-averaged subtidal lateral flow velocity amplitude during PRD and v low−discharge is the transect-averaged subtidal lateral flow velocity amplitude under the low-discharge condition (Run 2). The lateral flow recovery time T v is given by:

Conclusions
Small-/medium-sized rivers are more likely to suffer PRD than large rivers, which becomes important as the global warming-induced inundation risks rise and terrestrial sediment sequestration increases. Based on some in situ measurements, numerical experiments were designed to study the lateral flow response to PRD in the ORE. The transient PRD alters the stratification and salinity structure, and changes the Transect I more than Transect II. The lateral flow structure varies significantly over an ebb-flood tidal cycle, and PRD generally makes lateral flow weaker, which has also been reported by Chant et al. [35] in the observation of lateral flow during flooding season and dry season. The extra longitudinal momentum can rebuild the longitudinal flow and impact the lateral flow by lateral advection (−uv x ) and the lateral Coriolis term (−fu) (both are longitudinal flow-related terms). Moreover, the extra buoyancy changes the lateral-baroclinic-pressure-gradient terms from the low-discharge condition. As a result, the longitudinal flow-related term, lateral-barotropic-pressure-gradient and lateral-baroclinic-pressure-gradient are the three most significant terms affected by the PRD effects on lateral flow. In previous studies [25][26][27][28][29][30][31], the lateral Coriolis term and lateral-baroclinic-pressure-gradient term are also found to play essential roles in lateral flow.
We find that the flow rebounds to a quasi-steady state more quickly than the salinity. The recovery time of salinity is positively correlated with the peak discharge and shortens downstream, indicating that the PRD impacts on the saltwater front more. The recovery time of longitudinal flow is close to the PRD duration and essentially shortens from Transect I to Transect II. In general, the longitudinal flow is recovered more quickly than the lateral flow after PRD, as the latter depends more on the salinity recovery. The lateral flow recovery time depends on the transect width, the subtidal lateral flow amplitude, and the variation caused by PRD. This scaling of recovery time may be applicable to other river estuaries.
The role of PRD in estuarine dynamics, especially lateral flow, has been examined in this study. Due to the lack of complete observations of transect lateral flow during neap tide, in this study, we focused only on spring tides. The model results show that the patterns of the lateral flow during neap tides are similar to those during spring tides. In addition, the lateral flow recovery time to PRD is slightly longer than that during spring tide. The morphological evolution of small-/medium-sized river estuaries due to extreme weather conditions is not yet sufficiently understood. Therefore, the impact of the PRD on sediment transport will also be explored in future.