ff ects of bottom topography on dynamics of river discharges in tidal regions : case study of twin plumes in Taiwan Strait

Introduction Conclusions References

This paper is focused on the plumes of the Zhuoshui River and the Wu River.These rivers are significant sources of freshwater discharges from the west coast of Taiwan into Taiwan Strait connecting the South China Sea and the East China Sea.While the Zhuoshui River estuary is located at a significant topographic feature called the Chang-Yuen Ridge (hereinafter CYR) and the region is characterized by a very gentle bottom slope and a considerable part of the area is subject to wetting/drying (WAD) during the tidal cycle, the Wu River runs into the ocean at the northern extremity of CYR, where the bottom slope is steeper and no significant WAD takes place.The mouths of the two rivers are situated at a distance of only 30 km from each other, the freshwater discharge rates for the two estuaries are comparable, their plumes are exposed to similar atmospheric and tidal forcing and encounter apparently similar oceanographic conditionsand yet, the two plumes exhibit significant differences in their dynamic behaviour.The Zhuoshui plume is normally smaller and confined to a narrow belt along the shore, while the Wu plume tends to attain a bulge-shaped pattern and is spatially more extensive -even though the Wu River discharge rate is generally smaller than that of the Zhuoshui River.Indeed, analysis of MERIS/ENVISAT visible band satellite imagery (not shown here) since 2002 suggests that the area of the Wu River plume is larger than the area of the Zhuoshui River plume by about 40 % on average.
The structure and evolution of the plume depend on a delicate balance between the input of freshwater and buoyancy from the river and plume dissipation through vertical and horizontal turbulent mixing.Hence, any model of plume evolution is particularly sensitive to how turbulence is parametrized (Luketina and Imberger, 1989;Oey and Mellor, 1993;Ruddick et al., 1995;Garvine, 1995;Kourafalou et al., 1996;Hetland, 2005Hetland, , 2010;;Schiller and Kourafalou, 2010;Wang et al., 2011;Hwang et al., 2011).
One mechanism possibly responsible for the differences between the two plumes could be tidal straining, i.e. the interaction between the cross-shore velocity shear and the crossshore density gradient, and the so-called strain-induced periodic stratification (SIPS) as its consequence.As was shown by Simpson et al. (1990), in regions of freshwater influence, tidal straining generates semidiurnal currents rotating anticyclonically in the surface layer and cyclonically in the bottom layer.These currents interacting with the cross-shore density gradients produce a semidiurnal cycle of stratification.Due to continuity requirements imposed by the coast, the offshore-directed surface currents and onshore-directed bottom currents lead to coastal upwelling tending to displace the plume offshore (de Boer et al., 2009).This effect on river plumes has been investigated in a number of field and model studies (Simpson and Souza, 1995a, b;Souza and James, 1996;de Boer et al., 2006de Boer et al., , 2009;;Burchard and Schuttelaars, 2012).However, energetic tidal currents impinging on the CYR produce strong bottom-generated turbulence penetrating upward which will tend to destroy the classical patterns of tidal straining system existing for a flat or gently sloping bottom.
In this paper, we use the Princeton Ocean Model to investigate the influence of mixing governed by bottom-generated turbulence on the dynamics and evolution of the Zhuoshui and Wu plumes and identify processes primarily responsible for differences in the development of the two plumes.The central hypothesis of this paper is as follows: the dynamic behaviours of the Zhuoshui River and Wu River plumes are different because their evolution occurs under different regimes of bottom-generated turbulent mixing.
In this paper, we attempt to investigate the influence of mixing governed by bottom-generated turbulence on the dynamics and evolution of the Zhuoshui and Wu plumes and show what processes are responsible for differences in development of the two plumes.
Additionally, we apply a Lagrangian model to investigate the role of the tidal WAD cycle in the Zhuoshui mouth area, and show that this effect leads to further decay of the Zhuoshui plume.We use observational data from a field campaign in the study area to tune and validate the model experiments.
The paper is organized as follows.In Sect.2, the study region is described and observational data are presented.Section 3 provides a brief overview of the two models used in this study.In Sect.4, results of numerical simulations of the twin river plumes are given, followed by conclusions in Sect. 5.

Study region and some observational data
Taiwan Strait (TS) connects the East China Sea in the north and the South China Sea in the south, and is sandwiched by the China mainland and Taiwan Island (Zhu et al., 2013).Figure 1 shows the bathymetry in the TS and surrounding areas.The major topographic features are the shallow, sandy Chang-Yuen Ridge north of the Peng-Hu Channel (PHC), the Taiwan Bank (TB) in the southwestern TS, and the deepening topography in the northeastern TS.Most of the strait is shallower than 80 m.The PHC is located between the Peng-Hu Archipelago and Taiwan Island, descending southward sharply from 70 m to deeper than 1000 m.It is generally believed that the PHC is important in conducting the water from the southern opening of TS into the rest of the strait (Chuang, 1985(Chuang, , 1986;;Wang and Chern, 1992;Jan et al., 1994).
In this study, we focus on the central TS with the CYR region.In summer, the general flow here is represented by the Northern South China Current (NSCS) flowing north through TS (Jan et al., 2002(Jan et al., , 2010)).Combined with the tidal current, the NSCS impinging on the CYR is separated so that the surface and bottom waters flow in different directions: the heavier deep waters separate from the Taiwan coast south of the CYR, while the lighter surface waters flow over the CYR and spread along the eastern coast of the TS (Jan et al., 1994).
As was shown by Chuang (1985Chuang ( , 1986) ) and Jan et al. (2002), semidiurnal tidal oscillations play a crucial role in the circulation of TS.Tides penetrate into TS through the northern and the southern openings, then converging in the central part of the strait.Lin et al. (2000) argued that the semidiurnal tide in TS is represented as a standing wave, while Jan et al. (2002) suggested that the propagating and standing semidiurnal tides coexist, with most of the energy of the semidiurnal tide coming from the north.The diurnal tide in TS is less intense than the semidiurnal tide, except in the southernmost sector.Sea level measurements revealed that the diurnal tide propagates from north to south (Hwung et al., 1986).
Another important oceanographic feature in the TS is massive discharge of freshwater from multiple rivers.With its length of 176 km, catchment area of 3160 km 2 and long-term average discharge rate of 210 m 3 s −1 , the Zhuoshui River is the largest of those at the western coast of Taiwan.The Zhuoshui River's ROFI (River Of Freshwater Influence) stretches along the western coast.In the north, it often merges with the ROFI of the Wu River (also known as the Dudu River) whose average discharge rate is about 120 m 3 s −1 .
A map of the field survey conducted in the area on 25-27 June 2013 is shown in Fig. 2. A local fishery boat was used to take measurements at 13 hydrographic stations in the coastal zone.All stations were occupied on 27 June.
Temperature and salinity along the ship's track were recorded continuously using a pump-through CTD system.
In addition, at each of the stations, surface-to-bottom CTD profiling was conducted and water samples were collected from the surface and the near-bottom layer.Standard Sea Bird's SBE19plus instruments were used for the profiling and in the pump-through device.
The water velocity was continuously measured at a single point as 10 min averages by SeaHorse mechanical current metre (Sheremet, 2010) during 25-27 June.The instrument was deployed at a depth of about 5 m at a distance of 4 km to the west from the Zhuoshui River estuary (23.8488 • N, 120.1995 • E).
A portable meteorological station was mounted at a height of about 10 m above the sea level at the top of an old lighthouse at the Wu River mouth, 32 km north of the Zhuoshui River estuary.The station started working at 10:00, 25 June, and the data of wind speed and direction, as well as atmospheric pressure, air temperature and humidity were collected as 20 min averages until 23:30 LT, 28 June.
As can be observed in Fig. 2, during the survey, the Zhuoshui River plume was manifested as a belt of low salinity (< 30 psu) oriented along the coast in a pattern suggesting northward propagation.Yet fresher water (S < 25 psu) is seen in the area adjacent to the river mouth, where the inner part of plume forms a sharp salinity front approximately at the 5 m isobaths.The velocity record for 25-27 June as obtained by the current meter is shown in Fig. 3 (lower panel), indicating dominance of semidiurnal tides.

Models
To investigate the dynamics of the Zhuoshui River and Wu River plumes, hereinafter we use two numerical models: namely, the well-known a finite difference, σ -coordinate Princeton Ocean Model (POM), and the recently developed Lagrangian model STRiPE (Osadchiev and Zavialov, 2013).The two models complement each other for simulating the complex dynamics of the plume.

Implementing POM
The POM is a 3-D primitive equation ocean model that includes full thermodynamics and a level 2.5 Mellor-Yamada turbulence closure (Blumberg and Mellor, 1987;Mellor and Yamada, 1982).The POM model grid area covers the middle part of TS from 116.5 to 121.0 • E and from 23 to 25 • N (Fig. 4).The model domain is divided into 270 × 122 grid cells whose sizes x in longitudinal direction and y in latitudinal direction are both equal to one nautical minute.This corresponds to zonal resolution ranging from 1.723 km at 23 • N to 1.698 km at 25 • N, and meridional resolution of 1.836 km.In the vertical, 21 σ -levels are unevenly distributed, so that higher resolution is set to be near the surface and bottom.The minimum water depth in the model domain was set to 4 m.
The time steps were set to 3 s for the external mode and 120 s for the internal mode.The vertical eddy viscosity and  diffusivity were provided by the Mellor-Yamada turbulence closure scheme with a background value of 10 −5 m 2 s −1 .The horizontal eddy viscosity, K L , was calculated through the embedded Smagorinsky formula with a proportionality parameter, C H , equal to 0.1, and the horizontal eddy diffusivity, K L was obtained based on the inverse Prandtl number of 0.5.Here, U and V are the horizontal components of mean velocity.
The model was forced with tides and river discharge.The river input was imposed at the corresponding grid points by specifying the vertical velocity at the sea surface, w S .Following Oey (1996), we specify it as where Q is river discharge rate (m 3 s −1 ) and N is the number of grid areas around the river mouth where the discharge is to be distributed.We assume that w S injects volume and freshwater (salinity = 0 psu), but no momentum.Note that, when used in combination with the STRiPE, the POM was run without river discharges.One important feature of the region of our interest is the CYR which may significantly affect the structure of the along-shore flow and enhance tidal mixing as the current impinges on the ridge.To describe this effect, it is important to adequately simulate the bottom stress.In the POM, the stress, τ x b and τ y b induced by the bottom friction, is described as where C z is a non-dimensional coefficient, which is a function of the roughness length of the bed, z 0 , and the water depth, H : where κ = 0.4 is the von Karman constant, and σ kb − 1 is the depth of the layer overlying the bottom in sigma coordinates.
In other words, a logarithmic current profile is assumed, and the bottom friction is determined by the current velocity of the layer nearest to the bed.For stability of numerical simulations of flows in shallow regions, such as straits and channels, the Manning formula for bottom friction is often applied (e.g.Liu at al., 2003;Chiou et al., 2010).Instead of Eq. ( 4) we use the following formula for the coefficient C z : where C = n −1 H 1/6 , with n the Manning coefficient and g the acceleration of gravity.The bottom friction is determined with a use of the vertical mean of U and V in Eq. ( 3) at the corresponding location.In this case, therefore, the velocity is much less sensitive to the logarithmic profile assumption used in conventional parametrizations for C z which may lead to large errors and spurious features over steep topography.Following Chiou et al. (2010), for TS, we set the Manning coefficient n to 0.032 s m −1/3 .At the open (northern, N b and southern, S b , see Fig. 4) boundaries of the model domain, two types of boundary conditions for the temperature and salinity were used: namely, the condition for inflow and that for outflow.In the case of the inflow into the domain, the values for T and S at the corresponding open boundary were accepted.In the case of the outflow from the domain, the radiation condition  show little dependence on wind forcing.In order to elucidate this, we also performed a numerical experiment using POM without tide and with constant northeasterly wind at 5 m s −1 acting alone.The modelled wind-driven current at the depth 5 m obtained thereby was below 0.05 cm s −1 , while the tidal current is as fast as up to 0.8 m s −1 .
Hence, in the experiments described in the rest of this paper, we run POM and STRiPE under the tidal forcing only, with no wind forcing.

Implementing the STRiPE model
The STRiPE (Surface-Trapped River Plume Evolution) is a recently developed Lagrangian numerical model designed specifically for simulating river plumes under various forcing conditions.Detailed description of the STRiPE model, as well as its validation and examples of its application to river plumes, are given in Osadchiev and Zavialov (2013).The STRiPE model tracks the motion of river water "particles" discharged into the sea.The particles are the elementary water columns released from the river mouth, assigned with the initial velocity depending on discharge rate.The initial height, density and salinity of the particles are set equal to river mouth depth, density and salinity.While moving, the particles are allowed to mix with the seawater, so their characteristics change in time.Finally they become equal to the characteristics of the ambient ocean -in particular, particle height becomes equal to zero which means that the particle eventually dissipates.In the model, vertical mixing is based on the salinity diffusion equation: where S is particle salinity and D v is the vertical diffusion coefficient parametrized via the Richardson number Ri and scaling coefficient C v as given in Large (1994): Considering the main forces that determine river plume dynamics -namely, the Coriolis force, the force applied from the wind, the friction with the underlayer flow, the lateral friction, and the pressure gradient force -the model then solves the momentum equations applied to the individual particles.Small-scale horizontal turbulent mixing is modelled by the random-walk method (Hunter, 1987) complementing the momentum equations.Horizontal turbulent diffusivity K L used in the random-walk scheme is parametrized by the Smagorinsky diffusion formula (1).At every step of the model integration, the overall set of particles represents the river plume, and hence the temporal evolution of the plume structure is obtained.
The advantage of the STRiPE consists in computational efficiency and relative easiness of simulating a variety of aspects of plume dynamics.Here, we use it to elucidate the role of WAD effects.However, a significant drawback of the model lies in the fact that the circulation of oceanic waters, surrounding the river plume, cannot be simulated by the STRiPE model itself.The background velocity field must be prescribed as input data, obtained either from measurements or from another model.We, therefore, use the STRiPE in combination with POM.
Firstly, we run POM simulations without freshwater inflow to obtain the "background" coastal circulation.Then we use this circulation field as input data for STRiPE and simulate river plumes under WAD and no-WAD conditions using STRiPE.The STRiPE was applied to the study region, including the WAD zone (about 4 km wide) near the Zhuoshui River estuary.The numerical experiments were set up as follows.The model domain was the area between the Zhuoshui and Wu estuaries with the realistic shoreline from 120.0 to 120.65 • E and from 23.7 to 24.5 • N. The freshwater particles were released from the Zhuoshui River and the Wu River estuaries with the initial velocities of 0.1 m s −1 .The initial vertical scale of the released water volumes (i.e. the plume thickness at the mouth) was assumed equal to 5 m.The model integration time step was set to 10 min, while the total simulation period was 2 days, from 01:00 of 26 June 2013, through 01:00 of 28 June 2013.The input data of the ambient velocity fields were imported from POM solved with no freshwater inflow.

Vertical structure of tidal velocity
Figure 5 shows the simulated sea surface height (SSH) and the velocity components at both selected sites near the Wu River (Site 1) and Zhuoshui River (Site 2) mouths.It can be seen from the figure that the along-and cross-shore components of the velocity were tidally forced and exhibited semidiurnal variability.The zonal (cross-shore) velocity at Site 1 exceeded 0.55 m s −1 for the flood flow and flow −0.44 m s −1 for the ebb, while the meridional (along-shore) velocity component was somewhat smaller during both flood and ebb.For the study area, a flood phase of tide occurs when flow turns to the northeast and ebb phase occurs when flow turns to the southwest.A peculiar feature of tidal dynamics in this area is the asymmetry of the SSH curve and current velocity, so that near the Wu River (Site 1), the period of falling tide lasted less (∼ 1 h) than that of the rising tide.Note that here, the flood and ebb peaks perfectly coincide with maximum and minimum of the SSH, respectively.The water slacks occur at the middle of the transition period between the flood and ebb tides.In contrast, a stronger shift between the SSH curve and the tidal current takes place at shallower Site 2 (Zhuoshui River), where the ratio between the falling and the rising tides is the same as above, but the tidal current velocity leads the maximum of SSH by approximately 1.5 h.At the Site 2, the meridional velocity, on the flood flow, exceeded 1.2 and −1.1 m s −1 on the ebb flow, while the zonal velocity component was weaker and its magnitude was less than 0.6 m s −1 .The slight slope of the velocity contours in Fig. 5 suggests that the tidal currents at the surface lag those in the near-bottom layer.This is clearly seen at Site 1 near the Wu River mouth.At the shallower Site 2 tidal mixing partly offsets this effect.It can also be seen from Fig. 5 that the slacks when velocity falls to zero, and the moments of current reversal for the Zhuoshui River area precede those for the Wu River.Such an "asymmetry" may result in periodic mergers and separations of the plumes during the tidal cycle.

Horizontal structure of tidal currents and evolution of the river plumes
In the numerical experiments described in what follows, the discharge rates for both rivers were set to 210 m 3 s −1 , which is the climatic average for the Zhuoshui River.As a matter of fact, the long-term average discharge rate of the Wu River is somewhat lower, but the runoffs are highly variable in time, and for the purposes of this study we set the discharges of both rivers equal to each other in order to emphasize other dynamical factors determining differences in the behaviour of the two plumes.To reach the stabilization of plume sizes, tidal forcing was applied as a seamless perpetual 3-day cycle from 25 to 27 June 2013 -each run started on 25 June from the first water slack event following the ebb, and stopped with the second slack at the end of 27 June.With the chosen discharge rate values, the stabilization of plumes was attained during 3 × 3 computational days.
Figure 6 shows four successive phases of the tidal velocity and the respective Zhuoshui River and Wu River plumes as obtained for 27 June.The velocity fields and the plumes are presented against the position of the 30 m isobath, which approximately delineates the Chang-Yuen Ridge plateau.The tidal cycle was mostly governed by the semidiurnal M2-constituent with a period of approximately 12.4 h. Figure 6b, d clearly illustrate significant enhancement of the flood and ebb flows over the CYR.Here, the velocity peaks exceeded 1.2 and −0.8 m s −1 during the flood and the ebb, respectively.The 30 m isobath helps to determine the mean position of the flow-enhanced area above the CYR during the tidal cycle.However, in the flood phase, this area is notably shifted to the north with respect to the 30 m isobaths (Fig. 6b), while during the ebb the flow-enhanced area remains practically inside the 30 m contour (Fig. 6d).
More details can be seen in the insets at lower panels of Fig. 6b, d presenting the Y -Z structure of the tidal current along 120 • E during the peak flood and ebb.The maximum currents for the flood and ebb are located over the summit of the CYR and decrease towards the rim of the ridge.Asymmetry of the velocity is clearly seen in both cases, so that during the flood a core of meridional velocity > 0.5 m s −1 (yellow shading) propagates with the tide to the north with the frontal flank slightly deepening, while during the ebb, the core of meridional velocity < −0.5 (blue shading) shifts slightly to the south.The asymmetry of the tidal cycle and the shape of the ridge results in asymmetry of the velocity distribution over the ridge.
The numerical experiments highlighted important differences evident in the shape and the behaviour of the two plumes.Although the discharge rates of the Zhuoshui River and the Wu River were set equal, throughout the tidal cycle, the Zhuoshui River plume was somewhat narrower, and less regularly shaped and more eroded than the Wu River plume.The latter attained a "classical" bulge-like shape, which remained well pronounced, albeit its position shifted back and forth with the tidal oscillations.In contrast, the Zhuoshui River plume stretches along the shore but its shape is largely distorted by turbulent mixing.Such differences in the structure and the evolution of the two plumes located so close to each other are likely to be associated with inhomogeneous structure of turbulence caused by the interaction of tidal flow with the bottom topography.In particular, the shallow CYR must play a central role in such interaction.In the next section, we use the so-called Simpson-Hunter parameter (Simpson and Hunter, 1974) to elucidate this effect.This parameter allows us to approximate identification of the limit between tidally well-mixed and stratified zones.Furthermore, we focus on the tidally well-mixed zones near the CYR and the river plumes.

The Simpson-Hunter parameter
Using simple energy budget considerations, Simpson and Hunter (1974) argued that the location of the frontal zone between tidally well-mixed and stratified regions -including those formed by river discharges -can be estimated by a critical value of the tidal mixing parameter K SH : where h is the mean water depth (m), and U is the depthaveraged tidal velocity (m s −1 ).The value of K SH essentially reflects the mean intensity of tidal mixing through the water column at a specific location.As K SH decreases, either because the water depth h decreases or the tidal current u increases, the tidal mixing becomes stronger, and vice versa.This approach has been successfully applied to investigating tidal fronts on the continental shelves in a number of more recent studies (e.g.Zhao, 1985;Garrett and Maas, 1993;Hu et al., 2003;Lü et al., 2007;Holt and Umlauf, 2008; Zhu et al., 2013).For TS, the critical value of the parameter, as was found by Zhu et al. (2013), is equal to K SHc = 1.78.Note that the Simpson-Hunter parameter was originally proposed for thermal stratification resulting from solar radiation producing buoyancy in the surface layer.In the present work, we extend this approach to include lateral contrasts of buoyancy due to river discharges.We computed the distribution of K SH using the velocity obtained from the POM. Figure 7 (left panel) exhibits the contours where the parameter K SH attains this critical value of 1.78, for the study region during the flood flow.As seen, the critical contour of K SH delineates the CYR, indicating that strong tidal mixing occurs over the entire CYR.For comparison, in the right panel of Fig. 7 we present the critical tidal mixing parameter, K SHc , for the entire TS estimated from CTD measurements during the 8 years (2004-2011) by Zhu et al. (2013).These estimates also indicate strong tidal mixing in the area over the CYR in the vicinity of the Zhuoshui and Wu rivers.
As Fig. 7 (left panel) shows, the Zhuoshui River mouth is inside the critical K SH contour, while the Wu River estuary is outside it.This fact may explain the differences in the behaviour of the plumes during a tidal cycle: the Zhuoshui River discharges the river water into the zone of intense tidal mixing; therefore, the plume is exposed to extreme erosion.In contrast, the Wu River delivers river water into a moderate mixing zone, a few miles north from the CYR energetic zone as indicated by the critical K SH contour in Fig. 7, although the ebb phase of the tide advects the river water southwestward towards the enhanced mixing zone.This results in narrowing of the southern extremity of the Wu River plume, evident in Fig. 6d.However, during the entire tidal cycle, most of the Wu River plume lies outside of the intense tidal mixing zone, which may explain the relatively weak destructive effect of mixing on its bulge (Fig. 6a-d).
Thus, the Simpson-Hunter parameter is useful for determining of tidal mixing zones and qualitative elucidation of their potential effect on the fate of the river plumes.However, more in-depth analysis requires more information on processes and turbulent structure in the study region.In the next sections, we analyse some modelled turbulent parameters, giving emphasis to the features that may help to explain the differences in the characteristics of the Zhuoshui River and the Wu River plumes.

Turbulent kinetic energy production rate
In the bottom layer of the ocean, the water flowing above the seabed experiences frictional stress, which extends upwards throughout the water column.The bottom friction is one of the major causes of tidal turbulent kinetic energy (TKE) production in shallow water.The stress due to bottom friction is (see Eq. 3) where C D is the drag coefficient, ρ is the sea water density, and u is the horizontal current velocity vector.It follows then that the work A done by the bottom stress per unit time is expressed as where A is assumed to balance the energy dissipation that occurs in the boundary layer (although this is not the case for non-equilibrium turbulent regimes, as shown by direct measurements of turbulence with ADCP in energetic tidal channels and estuaries; e.g.Rippeth et al., 2002Rippeth et al., , 2003;;Korotenko et al., 2012Korotenko et al., , 2013)).As the tidal current flows across the shoal of small depth, the current velocity increases due to the conservation of mass.TKE production and dissipation rates will, therefore, grow as the third power of the current speed.The  As already illustrated by Simpson-Hunter's parameter above, the P S distributions, shown in Fig. 8a and b indicate that the continental waters discharged from the Zhuoshui River and the Wu River propagate in energetically different ambiances as manifested by the values of TKE near the surface (Fig. 8a, b).While the Zhuoshui River plume spreads through the strong shear-production area over the CYR and is, therefore, subject to faster destruction, the Wu Rivers plume would spread mainly in the weak shear production area, retaining larger area and the classical shape of the plume (Yankovsky and Chapman, 1997).
Dissimilarities of tidal mixing in the areas adjacent to both rivers and their reflection in the salinity gradient in the river plumes can be revealed by the horizontal Richardson number (Monismith et al., 1996;Stacey 1996;Bowen, 2000) defined as Ri x = gβ (∂S/∂x) H 2 u 2 * , where β is the saline contraction coefficient, u * is friction velocity, H is depth.
If Ri x < Ri crit , where the critical value Ri crit is an O(1) constant, mixing is strong and destroys stratification, whereas if Ri x > Ri crit , mixing is not sufficient to prevent the stratification build-up that intensifies with each tidal cycle.Using the POM, Monismith et al. (1996), Stacey (1996) and Bowen (2000) have found that the critical value of Ri crit is close to 0.3.
Figure 8c and d show the normalized ratio Ri x /Ri crit for the phases corresponding to (a) the water slack at 08:00 and (b) the peak flood at 11:00 on 27 June 2013, respectively.As seen, the values of the horizontal Richardson numbers fall below Ri crit (blue shading) in the area of the Zhuoshui River discharge indicating considerable destruction of stratification by tidal mixing, while the Wu River plume is located in area where the effect of the tidal mixing does not prevent the development of stratification.

Time-depth variability of TKE
Generally, in tidal flows, variations of turbulent kinetic energy are tightly connected with the tidal cycle (e.g.Rippeth et al., 2003;Souza and Howarth, 2005;Wiles et al., 2006;Korotenko et al., 2012Korotenko et al., , 2013)).The simulated time series of the vertical structure of TKE for Sites 1 and 2 are depicted in Fig. 9.Note that this particular experiment was performed without taking into account the discharges from both rivers to emphasize the upward propagation of bottom-shear generated turbulence from the seabed due to its generation by bottom stress.As seen, during the tidal cycle, the TKE spans about four orders of magnitude.During the tidal slack, the TKE sharply drops below 10 −4 m 2 s −2 throughout the water column.The TKE attains a maximum in the bottom layer and decreases upwards.The most important finding illustrated by Fig. 9 is the following: while at the Site 2 (Zhuoshui River plume area) TKE values above 10 −3 m 2 s −2 emerge at the surface because of shallower water and higher bottom velocity shear, at the Site 1 (the Wu River plume area), the penetration of TKE generated by the bottom stress is limited to about 25 m above the seabed and, hence, does not reach the surface layer, where it could affect the river plume.To estimate the effect of bottom-generated turbulence on both rivers' plumes, next we examine zonal transects of the term in TKE balance that describes destruction of turbulence caused by buoyancy.

Zonal transects of TKE destruction rate by buoyancy
Following Wang et al. (2011), for convenience, here we define the rate of TKE destruction/production by buoyancy, P b , as −K H N 2 , where K H is vertical diffusivity and N 2 is the buoyancy frequency N 2 = −(g/ρ o )∂ρ/∂z).As can be seen, in the case where N 2 > 0 (P b < 0), i.e. the flow is stable, the stratification suppresses turbulence (the effect of turbulence destruction by buoyancy takes place); whilst in the case where N 2 < 0 (P b > 0), the flow becomes unstable and the instability leads to TKE generation.Figure 10 presents the simulated time series of the vertical structure of TKE destruction rate by buoyancy for zonal transects A and B shown in Fig. 1.To elucidate the intratidal variability of both river plumes, the cross-plume transects A and B were positioned approximately 3 and 5 km north of the Wu and Zhuoshui river's mouth, respectively.The values of −P b varied from 10 −8 to 10 −4 m 2 s −3 .The upper limit corresponds to complete destruction of TKE by buoyancy.
Comparison of the vertical distribution of −P b reveals tremendous differences in balance of TKE for both rivers during a tidal cycle.The plume of Wu River located in deeper water experiences a smaller effect from the bottom-generated turbulence on the plume structure and thickness, although significant mixing does occur during ebb when the stratification in the 20 m bottom layer decays rapidly.During slack periods when the flow reverses to flood (Fig. 10b, f at transect A), the stratification tends to suppress turbulence at the mid-depth.Generally, despite the tidal mixing activity, the thickness of the strong stratified upper part of the Wu River plume remains about 15 m and mixing in the bottom layer that develops during the ebb and flood and propagating upward is too small to compete with the stratification originated by the surface plume.
In contrast to the Wu River plume, the plume of Zhuoshui River (Fig. 10b), even at the sea surface, is unstable and intermittent due to strong tidal mixing above the shallow ridge.As seen, during the entire tidal cycle, the bottom shear-generated turbulence is strong throughout the water column, while the destructive influence of the stratification on TKE is limited to the upper 10 m.It is noticeable that the plume on ebb and slack break apart into several lenses, so that high mixing zones appear at the surface between them (Fig. 10b, d and  e, at the transect B).During flood, tidal current advect's the river water and salinity front northward from the river mouth, creating straining which acts to stratify the water column and competes with the bottom-generated turbulence and mixing.The newly discharged riverine water soon reaches the transect B developing a temporal restoration of Zhuoshui River's plume at the sea surface.

Horizontal diffusivity K L
The discussion above was focused on the vertical turbulent mixing in connection with the bottom topography.However, horizontal mixing can also contribute significantly to the evolution of the plumes.For example, Kourafalou et al. (1996) studied the effect of variability of the horizontal diffusivity on the development and dynamics of a river plume and found that, as the horizontal diffusivity decreased, the along-shore and cross-shore dimensions of the river plume increased while coastal current meandering diminished.Their experiments were performed with constant values of K L , which may be appropriate for experiments with a stable external forcing.However, in reality, the plume dynamics is more complex and the impact of variable horizontal diffusivity on river plume dynamics should be taken into account.Therefore, in this paper, we examine the time-depth variability of the K L (Eq. 1) over the tidal cycle.
Figure 11 exhibits the time-depth variability of the horizontal diffusivity at sites 1 and 2. Comparing the variations of K L at the two sites, we can see considerable differences that are likely to be caused by the flow impinging on the CYR. the ridge creates significant deformation of the velocity field above the CYR.As Fig. 11 indicates, for site 1 near the Wu River during the flood, the maximum of K L occurs at the sea surface, while during the ebb it migrates to the intermediate depths of 5-10 m.In contrast, at site 2 near the Zhuoshui River, the horizontal diffusivity significantly increased during the ebb, i.e. velocity deformation peaked when the SSH over the CYR attained its minimum.Thus, comparing the variability of K L near the Wu River and Zhuoshui River, we again found considerable differences, potentially important for the dynamics of the river plumes.
To estimate the relative role of horizontal diffusion in spreading of both river plumes as compared with that of tidal mixing, we ran a numerical experiment varying the coefficient C H in the formula for K L (Eq. 1).namely, we assumed C H = 0.7, which corresponds to large diffusivity/viscosity.Figure 11c illustrates the result of this experiment for peak flood.The salinity field shown in this figure can be compared with that shown in Fig. 6b.In Fig. 11c, we present the surface salinity distribution without velocity vectors to show the former in detail.As expected, with the elevated value of C H , both river plumes are more disperse and "smooth" than those simulated for the lower C H .It is notable that the bulge-shaped Wu River plume propagates significantly farther northeastward along the shore than the Zhuoshui River plume which is accumulated around the mouth.The general structure of the Zhuoshui plume -including the separation between the "new" branch of the plume north of the discharge point going with the tidal current and the "old" one south of the discharge point remaining from the ebb phase -is conserved throughout the tidal cycle.This indicates that the influence of horizontal diffusion on the formation of Zhuoshui River plume appears to play a secondary role in comparison with vertical mixing associated with bottomgenerated turbulence above the CYR.

The role of WAD
In terms of the tidal wetting and drying (WAD), the areas adjacent to the Zhuoshui River and the Wu River exhibit very different conditions.The shallow and flat bottom near the Zhuoshui River mouth dries during the low tide over an area about 15 km long and 4 km wide (Fig. 12), while the deeper shelf near the Wu River mouth has no significant WAD and constantly remains covered with water.This feature gives us an opportunity to investigate the influence of WAD processes on river plume dynamics.Because the version of POM used in this study was not intended to simulate WAD, for this purpose we used the Lagrangian model STRiPE for this purpose.The model and the experimental setup are described above in Sect.3.2.
We performed two STRiPE simulations, one with WAD and one without it, in order to investigate the influence of WAD on plume dynamics.Both STRiPE simulations were performed without wind forcing, and the ambient sea circulation was adapted from the POM experiments as described  above.Firstly, we simulated the behaviour of the Zhuoshui River and the Wu River plumes with the fixed shoreline corresponding to the high-tide phase of the WAD cycle.Secondly, we ran the model with the WAD cycle and the movement of the boundary taken into account.The movement of the shoreline between its highest and lowest tide positions was set uniform for every tidal cycle.We performed two POM simulations, one with the high-tide shoreline and the other one with the low-tide shoreline configurations (outside the WAD region, the shelf circulations obtained for the two cases exhibited only minor differences).If there was no WAD area in the model domain, the set of intermediate velocity fields between high-and low-tide phases could be computed as transitional averages between the two limiting cases.In the presence of WAD these average velocity fields should be additionally resized in accordance with the contraction of the field domain caused by uniform movement of the shoreline in the WAD region.The resulting set of fields is used as an input for the STRiPE simulation.The boundary condition for the STRiPE in the WAD zone was no-normal flow at the coastline, i.e. once the "particle" hits the shoreline, the cross-shore velocity of the "particle" is set equal to zero at the shoreline, while the along-shore velocity component retains its value without change.
Both STRiPE simulations (with and without WAD) were performed without wind forcing, and the ambient sea circulation was adapted from the POM experiments as described above.
Firstly, we should point out that the results of the non-WAD STRiPE simulation are consistent with the corresponding results of the POM simulation which are described in Sects.4.1-4.2 and illustrated by Fig. 6. Figure 13 shows that the positions, shapes and areas of the STRiPE (green area) and POM (red area) plume contours are in good agreement with each other.Secondly, comparing the results of the two STRiPE simulations, we can obtain some insights into the role of WAD.In Fig. 14, we present the modelled surface salinity distributions outside the WAD region averaged over a 48 h simulation period, namely from 01:00 of 26 June 2013 through 01:00 of 28 June 2013.
It can be seen from Fig. 14 that the presence of WAD significantly changed the Zhuoshui River plume at the WAD zone; plume area decreased considerably compared with the "no WAD" conditions.The salinity anomaly in the plume also decreased by about 3 psu.At the same time, the Wu River plume was affected less significantly.Hence, it is likely that WAD adds significantly to the erosion of the Zhuoshui River plume.The possible physical mechanism is as follows: the emergence of a piece of land during the low-tide phase  alters the along-shore current on the inner shelf and increases shears, which, in turn, enhances mixing down through the water column in the plume area.

Summary and conclusions
A high-resolution Princeton Ocean Model was used to investigate the intratidal variability of currents and turbulent mixing and their impact on the characteristics and evolution of the plumes of two neighbouring rivers, the Zhuoshui River and the Wu River, at the central eastern coast of Taiwan Strait, where the shallow sand Chang-Yuen Ridge is a significant feature of the bottom topography having considerable impacts on the local dynamics.The two estuaries are located close to each other and the conditions they are exposed to are similar in some respects, and yet the two plumes exhibit significantly different behaviour.We attribute this to the differences in the bottom topography.Indeed, as numerical experiments showed, the Wu River plume located in a relatively deep area at the northern extremity of the CYR is characterized by larger spatial extent and tends to maintain its "bulge-like" shape throughout the tidal cycle, whilst the Zhuoshui River plume localized over shallow and flat-bottom terrain is reduced in size, stretched along the shore and is largely destroyed by turbulent mixing.
For preliminary analysis of tidal mixing influence on the stratification development, we applied the Simpson-Hunter criterion (Fig. 7) and also evaluated horizontal Richardson number (Fig. 8c, d) to qualitatively explain the differences between the two plumes, and found out that, according to the criterion, the Zhuoshui River runs into the area that is permanently well mixed, and the Wu River mouth is located in a stratified area outside the region of intense tidal mixing.This is also confirmed by the results of numerical modeling with POM.It was demonstrated that the area of bottom-generated turbulence, which is strong during the ebb and flood flows, coincided with the shallow area above the CYR, thus indicating a dissimilarity of tidal mixing acting on the two plumes.
For more in-depth analysis of tidal mixing on stratification development in the area of both rivers, we studied intratidal variations of TKE shear production rate, P S (Fig. 8a, b), and TKE destruction/generation rate due to buoyancy, P B (Fig. 10).These two terms in the TKE budget both mirrored significant differences in behaviour of the rivers' plumes during a tidal cycle.
We also examined the time-depth variability of TKE (Fig. 9).The predicted bottom-stress-generated TKE was found to vary by 4 orders of magnitude.Near the Zhuoshui River mouth the bottom-generated turbulence penetrated throughout the water column and reached the sea surface during the ebb and flood tide phases, significantly destroying the river plume.On the other hand, in the deeper area near the Wu River plume, the bottom-generated TKE was not strong enough to mix the water column up to the surface and had only marginal influence on the plume.The turbulent kinetic energy production rate was found to vary strongly, attaining maximum values of about 10 −1 W m −3 over the CYR at the seabed during the flood and ebb tide.Away from the bottom layer, the TKE production gradually decreased upwards and was significantly suppressed in the surface layer by strong stratification generated by the river plumes.Inside the plumes, TKE production rate fell below 10 −6 W m −3 .The rates of TKE destruction by buoyancy varied from 10 −8 to 10 −4 m 2 s −3 .The upper limit denotes a complete destruction of TKE by buoyancy.
Horizontal turbulent mixing also contributes to the plume dynamics and demonstrated certain dissimilarities for the two river plumes (Fig. 11).In the Zhuoshui River area during the flood, the maximum of horizontal diffusivity occurred at the surface, while during the ebb it moves to intermediate depths below the plume layer.In contrast, in the Wu River mouth area the horizontal mixing was generally smaller and Experiments with increased horizontal viscosity/diffusivity (C H = 0.7 in Eq. 1) revealed a general smoothing effect on both rivers' plumes, however the shape of the plume remained akin to those predicted for the basic value C H = 0.1.
In addition, we applied a Lagrangian particle-tracking model used in combination with POM to investigate the effect of tidal wetting-and-drying near the Zhuoshui River mouth.It was found that this process adds to the downward mixing of the Zhuoshui River plume.
Altogether, these findings point to the extreme importance of bottom topography features for the dynamics of river discharges in tidal regions.
In conclusion, it should be noted that the environments of both rivers have been seriously threatened in recent years by ongoing construction activities, concrete industry and agriculture.Particularly, the concentration of metals such as copper and iron was elevated over the regulated standard.The elevated level of copper was attributed to the addition of that in the animal feeds for many livestock in the area (Environmental Policy Monthly, 2010).Thus, the present study may be useful for understanding the distribution and fate of various contaminants delivered into the marine environment by the Zhuoshui and Wu rivers.

Figure 1 .
Figure 1.Bathymetry of Taiwan Strait.TB is Taiwan Bank, PH is Peng-Hu Archipelago, CYR is the Chang-Yuen Ridge, PHC is the Peng-Hu Channel (modified fromJan and Chao, 2003).Points 1 and 2 indicate the locations of the sites we refer to in the discussion.The mean water depths at sites 1 and 2 are 33 m and 13 m, respectively.Zonal transects A and B we refer to in the discussion of buoyant destruction of turbulent kinetic energy.

Figure 2 .
Figure 2. Left: bathymetry of the Zhuoshui mouth area (based on ETOPO1).The isobaths are labelled in metres.Right: map of stations occupied during the field survey in June 2013 (bullets).The black box indicates the position of mooring station with current meter near the bottom.Black curves are the estimated surface salinity contours (morning of 27 June).

Figure 3 .
Figure 3. Lower panel: hourly averaged vectors of the current velocity at the depth 5 m obtained at the mooring station.Middle panel: same series as simulated by POM model, see Sect. 3. Upper panel: series of wind velocity as measured for the same period.

Figure 4 .
Figure 4.The POM domain and bottom topography.Locations of Zhuoshui River and Wu River mouths are indicated by white dots.
) was used, where n represents the direction normal to the open boundary.Here, U n is the vertically averaged normal component of the velocity.The barotropic (vertically averaged) velocities on the open boundaries of the POM were estimated using the Flather (1976) formula:

where η 0
is the sea surface elevation at open boundary provided with U 0 n by the barotropic tidal model of Hu et al. (2010); U n is the vertically averaged normal component of the velocity at the open boundary at time t; η is the model sea surface elevation calculated from the continuity equation and located half a grid cell inside of the open boundary in the POM model domain.The tidal model data used corresponded to the period of our field survey mentioned above, i.e. 25-27 June 2013.Figure 3 illustrates current velocities obtained from the POM forced by tide and the northeasterly wind of 5 m s −1 versus the measurements conducted with SeaHorse current metre on 25-27 June 2013.The instrument was deployed at the depth 9 m, about 0.5 m above the bottom.The comparison shows a reasonably good agreement between the modelled and measured velocities, perhaps except some discrepancy from about 00:00 to 06:00 on 27 June.Both the observed and the simulated currents are dominated by tide and appear to www.ocean-sci.net/10/863/2014/

Figure 5 .
Figure 5. Time-depth variability of the sea surface height: zonal (eastward direction is positive) and meridional (northward direction is positive) tidal velocities at Site 1 (upper panels) and Site 2 (lower panels).

Figure 6 .
Figure 6.Successive phases of tidal velocity and salinity at 1 m depth at (a) 08:00, (b) 11:00, (c) 14:00, and (d) 17:00 GMT on 27 June 2013.The phases (b) and (d) correspond to the peak flood (northeastward) and ebb (southwestward) flows, respectively.The 30 m isobath (solid blue line) approximately delineates the Chang-Yuen Ridge.The 32.5 psu salinity contour delineates the river plumes.The insets in (b) and (d) show meridional velocity along 120 • E above the ridge during the peak flood and ebb flows, respectively.

Figure 7 .
Figure 7. Left panel: planar plot of computed SH parameter in the coastal zone over the CYR area including Zhuoshiu and Wu rivers' mouths.The contour of K SHc = 1.78 (red line) and isobath of 30 m delineating the CYR (white line) are shown.The modelled K SH corresponds to flood phase at 11:00 GMT, 27 June 2013.Right panel: MODIS images of Chl a on 9 July 2009.The contours for the critical tidal mixing parameter (K SHc = 1.78) estimated from CTD measurements by Zhu et al. (2013) are shown in red.Locations of the summertime upwelling cores in the TS are shown with blue ellipses.Acronyms are given in Fig. 1.

Figure 8 .
Figure 8. Modelled horizontal distributions of shear production rate of TKE, P S (upper panel) and normalized horizontal Richardson number, Ri X /Ri crit (lower panel) at 1 m depth at (a, c) 08:00 and (b, d) 11:00 GMT on 27 June 2013.

Figure 9 .
Figure 9. Modelled time-depth variability of TKE at Site 1 (Wu River area) and Site 2 (Zhuoshui River area).The ordinate axis is height (in metres) above the bottom (mab).Ocean surface elevation is depicted by black line.For clarity, TKE results are presented in the absence of both plumes.

Figure 10 .
Figure 10.Distribution of TKE destruction by buoyancy, log 10 (−P b ), along transects A and B across the Wu and Zhuoshui rivers, respectively (Fig. 1).Tidal stages and tidal current directions on ebb (a, e) and on flood (c) for both transects are shown in the bottom right corners of the upper panel.

Figure 11 .
Figure 11.Simulated time-depth variability of horizontal diffusivity, K L , at (a) Site 1 (Wu River area) and (b) Site 2 (Zhuoshui River area).Panel (c) presents the distribution of surface salinity under large diffusivity (C H = 0.7 in Eq. 1) for the flood phase corresponding to Fig. 6b.

Figure 12 .
Figure 12.Model domain for numerical simulations with STRiPE.Also shown are the shorelines for low tide (yellow) and high tide (red) conditions.The area between the red and the yellow contours is the WAD region.

Figure 13 .
Figure 13.Comparison of STRiPE and POM simulation results.The green and red colours display the contours of the river plumes obtained by STRiPE and POM respectively at (a) 08:00, (b) 11:00, (c) 14:00, and (d) 17:00 GMT on 27 June 2013.The phases (b) and (d) correspond to the peak flood (northeastward) and ebb (southwestward) flows, respectively.The 32.5 psu salinity contour delineates the river plumes.

Figure 14 .
Figure 14.Modelled surface salinity distributions at the study area without WAD (left) and with WAD (right) under "no wind" conditions averaged for the period from 01:00 of 26 June 2013, through 01:00 of 28 June 2013.