PUBLICATIONSJournal of Geophysical Research: Oceans
RESEARCH ARTICLE
10.1002/2014JC009855

Special Section:
Pacific-Asian Marginal Seas
Key Points:
 Thirty three repeat shipboard ADCP
surveys were conducted in the
Qiongzhou Strait
 Tidal constituents were estimated bya modiﬁed tidal harmonic analysis
method Tidal rectiﬁcation and sea level are
important for residual current
formation

Correspondence to:
X.-H. Zhu,
xhzhu@sio.org.cn

Citation:
Zhu, X.-H., Y.-L. Ma, X. Guo, X. Fan,
Y. Long, Y. Yuan, J.-L. Xuan, and D.
Huang (2014), Tidal and residual
currents in the Qiongzhou Strait
estimated from shipboard ADCP
data using a modiﬁed tidal harmonic
analysis method, J. Geophys. Res.
Oceans, 119, 8039–8060, doi:10.1002/
2014JC009855.
Received 26 JAN 2014
Accepted 30 OCT 2014
Accepted article online 4 NOV 2014
Published online 25 NOV 2014

Tidal and residual currents in the Qiongzhou Strait estimatedfrom shipboard ADCP data using a modified tidal harmonic
analysis method
Xiao-Hua Zhu1,2, Yun-Long Ma1, Xinyu Guo1,3, Xiaopeng Fan1, Yu Long1, Yaochu Yuan1,
Ji-Liang Xuan1, and Daji Huang1,21State Key Laboratory of Satellite Ocean Environment Dynamics, Second Institute of Oceanography, State Oceanic
Administration, Hangzhou, China, 2Ocean College, Zhejiang University, Hangzhou, China, 3Center for Marine
Environmental Study, Ehime University, Matsuyama, Japan

Abstract In spring 2013, 33 repeat shipboard Acoustic Doppler Current Proﬁle (ADCP) surveys were conducted to measure the tidal current in the Qiongzhou Strait (QS). The major tidal currents and the residual
current along a section across the QS were estimated using a modiﬁed tidal harmonic analysis method
based on the inverse technique. A simple simulation and comparisons with previous observations demonstrated that the tidal currents estimated using the modiﬁed tidal harmonic analysis method are reasonable,
and this method was able to control the magnitude and deviation of the estimation error. The direction of
the major axis of tidal current ellipses is generally along the strait. Diurnal tidal constituents are dominantamong the ﬁve tidal current constituents (K1, O1, M2, S2, and MSf). The ratio of the amplitudes of O1, K1, M2,
S2, and MSf, averaged along the section across the QS is 1:0.79:0.42:0.27:0.29. The residual current along theentire section is all westward; the averaged velocity over the section is 6.062.1 cm s21; the associated volume transport through the section is 20.065 6 0.046 Sv (Sv 5 106 3 m3 s21), in which the second value
denotes the uncertainty of ﬁrst value. Dynamic analysis indicates that tidal current activity is more dominant
than mean current and eddy activity, and tidal rectiﬁcation and sea level difference between two entrancesof the QS are important in maintaining the residual current through the strait.

1. Introduction
The Qiongzhou Strait (QS) is a channel connecting the northern part of the South China Sea (SCS) andthe Beibu Gulf (Gulf of Tonkin) (Figure 1). The average width, length, and mean water depth are approximately 25 km, 70 km, and 70 m, respectively. The deepest location in the strait lies below 100 m. The
QS is a key passage for water exchange between the northern SCS and the Beibu Gulf. Quantifying thetidal current and volume transport through the QS is vital to understanding the continental shelf circulation in the northern SCS [Wang et al., 2010] and the cyclonic circulation in the Beibu Gulf [Wu et al.,2008].
The tidal motion in the SCS is maintained mainly by the energy ﬂuxes from the western Paciﬁc into the SCS
through the Luzon Strait [Fang et al., 1999; Zu et al., 2008]. The phase propagations of diurnal and semidiurnal tidal waves are southwestward across the whole SCS. The northern part of the tidal wave propagatesalong the southern coast of China and directly reaches the QS from its eastern side. A component of thesouthwestward tidal wave bifurcates near Hainan Island (Figure 1), propagates into the Beibu Gulf from its
southern boundary, and then approaches the QS from its western side. When tidal waves enter the QS, thediurnal and semidiurnal waves propagate in opposite directions [Cao and Fang, 1990], and nonlinear interaction is expected after they meet each other in the QS [Shi et al., 2002].
Since the early 1960s, some investigators have focused on the circulation around the QS and its effects on
the northern SCS and the Beibu Gulf. Wyrtki [1961] ﬁrst showed the possibility of a westward surface ﬂow
through the QS in October. However, this westward current disappeared in June. A similar current patternin the QS (with a westward surface ﬂow in fall that disappears in summer) was also suggested based on thewind-driven circulation theory [Ke, 1983; Yu and Liu, 1993].

ZHU ET AL.

C 2014. American Geophysical Union. All Rights Reserved.
V

8039

Journal of Geophysical Research: Oceans

10.1002/2014JC009855

o

(a)

CHINA
o

o

(b)

Guangdong

o

21 N

th

VIETNAM

u
So

15oN

22 N

0

Luzon

−20

20 N

10 N

Qiongzh

o

20 N

a
Se

o

ina

−4

00

Ch

0

Gulf of Beibu

o

19 N

5oN

18oN o
108 E

o

0 o
100 E

o

o

105 E

110 E

o

115 E

o

120 E

o

ou Strai

t

H
ai
na
n

25 N

o

o

109 E

110 E

o

o

111 E

112 E

125 E

(c)
24’

North

South

0

Leizhou Peninsula
18’

10

12’

40

30

50

90

70

Dep(m)

20

6’

60

50

30

80

o

(d)
0

20 N
2

4

6

8
10
12
Distance(km)1416Hainan Island

18
54’

48’

o

110 E

12’

24’

36’

Figure 1. (a) Bathymetry of South China Sea, (b) location of Qiongzhou Strait and Hainan Island, (c) Bathymetry of Qiongzhou Strait, and(d) vertical view (looking toward the east) of the bottom topography along a section (the black heavy line in Figure 1c), to which the velocity data collected along the different sections (blue lines in Figure 1c) are converted along the gray lines in Figure 1 c. The bottom topography obtained by the shipboard ADCP measurements using bottom track mode and by the ETOP1 data set is indicated by blue and redlines, respectively. The ‘‘star’’ indicates the position of the Haikou tidal gauge station.

Shi et al. [1998] reported the results of the ﬁrst synchronous current measurement using 15 ships anchoredalong a section across the QS and found a maximum westward current speed of 172 cm s21 during ﬂood
tides and a maximum eastward current speed of 142 cm s21 during ebb tides. They also showed a westward residual current with a speed of 10–18 cm s21 in the central and northern parts of the QS. However,
due to observational limitations (as only three stations were located at the center of the QS and the measurement only extended from the surface to a depth of 20 m), they could not determine the spatial structure of the tidal current across the QS.
By analyzing 37 years of current meter and tide gauge data from 1963 to 1999, Shi et al. [2002] showed thatthe residual current through the QS is westward throughout the year. They also provided a very rough estimation of the associated volume transport: approximately 0.2–0.4 Sv in winter and spring (using a meancurrent speed of 20–40 cm s21 at the surface and 10–15 cm s21 at the bottom); and approximately 0.1–0.2
Sv in summer and autumn (using a mean current speed of 5–10 cm s21 at the surface and 10–30 cm s21 at
the bottom). Similar results were obtained later using nonsynchronous current mooring data collected atseveral stations [Chen et al., 2006; Chen et al., 2007; Yan et al., 2008], drifting bottle data [Yang et al., 2003;
Bao et al., 2005], hydrographic calculation [Yang et al., 2006], and numerical studies [Chen et al., 2009].
Measurements to show the spatial structure of the tidal current (including its major tidal constituents andresidual current covering the whole strait), are, however, lacking. There is, therefore, no estimation of thevolume transport through the QS using synchronous data.
To provide these important data, we measured the tidal current by carrying out 33 repeat shipboard ADCP
surveys across the QS. We estimated the harmonic constants of ﬁve major tidal constituents (K1, O1, M2, S2,
and MSf) and the residual current using a modiﬁed tidal harmonic analysis method, evaluated our results bycomparing them with previous observations, and calculated the volume transport through the QS using theZHU ET AL.

C 2014. American Geophysical Union. All Rights Reserved.
V

8040

Journal of Geophysical Research: Oceans

10.1002/2014JC009855

m

1.5
0
−1.5

(a)
4

(b)
03/07

1

2

3

03/09

03/11

03/13

03/15

03/17

6

7

5

8

03/19

9
10

03/21

12

14

16

11

1315

03/23

18
20

22
24

26 28
30

17 19

21
23

25
27 29

03/25

03/27

32
31

03/29

33

03/31

04/02

Date
32

33
31

7

29

30

28

1

26
16

25

5

22

8

21
6

(c)
0

2
27
17

11
8
1

18

4

6

8

27

3

10

17
24

16

18

6

8

10

12

14
18
16
4

(e)
2

4

33
29
21

25

18

10

8

1

12

32

14
9
4

2

3

4

28
21
16
13

4

3

19
25

22
32

4

5

24
12

7

20

27

6
24
12

7

18

6
26

20

22

26
2

20
32
17

10

11
8
17

15

22

24

12
28

31
15
5

7

8
10

9

8

11

12
20

7

9

11

12

27

24

26
5
2

7

10

30

23

33

6
Hour

30
5
23
10

20

14

9

14
2

10

9

16

30

23
33

2

6
Hour

31

29

18
3

(g)
2

5

14
33

24
5

25
19

20

25
19

22
3

22
31

18

27
24

13

16

14

26

3

1

1

23

19

22

6

(f)

0

8
30

7
1

0

6

32
22
9

3

31
29

13

0

29

16

21
28
1

20

28
21
18
16
13

(d)
4

4
15

14

14

4

2

10

20
13

12

11

2

15

0

12

9
23
19

8
17

15

11

6

10

11

12

Figure 2. (a) Tidal level anomaly (black lines) observed at the north side of the Qiongzhou Strait and (b) the time schedule of the shipboard ADCP survey, and the phase distribution of ADCP survey time for a (c) O1, (d) K1, (e) P1, (f) M2, and (g) S2 tidal period. The cruise number is indicated above or under the bars.

residual current obtained from this study. We validated the performance of the modiﬁed tidal harmonic
analysis method through a simple simulation using tidal level data measured in the QS during the shipboard ADCP surveys. We carried out dynamic analysis to conﬁrm the previously suggested mechanisms forthe formation of residual currents in the QS.

2. Data Collection and Processing
A shipboard ADCP was used to measure the currents in the QS (Figure 1). The ADCP survey was performed
when maintenance was done on the acoustic instruments anchored on both side of the QS. Because there
were many ﬁshing nets in the QS, the shipboard ADCP survey was conducted only during the daytime. During the observation period (from 8 March to 2 April 2013), which included a period from a spring tide (17March) to next spring tide (31 March), 33 shipboard ADCP tracks were performed over the QS (Figure 1c).
Due to bad weather conditions, the ADCP survey was stopped for 6 days from 11 to 16 March (Figure 2).
The 33 shipboard ADCP sampling times (Figure 2) were almost distributed uniformly over all the phases ofO1 tide (Figure 2c) and two major semidiurnal tides (Figures 2f and 2g). The shipboard ADCP sampling could
not cover the entire phase range of K1 and P1 tidal constituent (Figures 2d and 2e).
A 300 kHz ADCP (Workhorse Self-Contained ADCP) was used to measure the velocity with bottom track
mode. The ADCP was set up at the side of a wooden ﬁshing ship using a stainless steel frame and its transducers were maintained at 1 m depth below the sea surface. The ship speed during the survey was kept at
3–5 m s21. It took approximately 1–2 h for the ADCP survey to be completed from one side of the strait tothe other side along a section in the QS.
The ADCP data-sampling rate was set to one ping per 2 s, which is close to the maximum data sampling
rate for this type of ADCP. The bin length and bin number were set to 2 m and 60, respectively. For theseobservation conditions, the error in raw velocity (standard deviation) with a single ping was less than
7.0 cm s21 [Teledyne, 2013]. The depth of the ﬁrst ADCP bin is approximately 3m (5 sensor depth

ZHU ET AL.

C 2014. American Geophysical Union. All Rights Reserved.
V

8041

Journal of Geophysical Research: Oceans

10.1002/2014JC009855(1 m) 1 blank distance (0.14 m) 1 bin length (2 m)), and the blank distance near the bottom is approximately 4 m (5H3ð12cos20 Þ, 20 is the ADCP beam angle, H is water depth with a typical value as 70 m).
The velocity in the top 3 m layer and the velocity 4 m from the bottom could not be measured and were
given by extrapolation following the method proposed by Gordon [1989]. A Global Positioning System(GPS) receiver was used to measure the geographical position of the ship. The GPS data-sampling interval
was 1 s. The ADCP velocity and GPS position data were recorded by the onboard computer using ‘‘VmDas’’software provided by Teledyne RD Instruments, Inc.
The ADCP data-processing procedures were as follows. First, the absolute velocity proﬁles along a section
were obtained using the bottom-tracking mode. The compass biases, which were mainly due to the magnetic declination, were carefully corrected by comparing the ship track determined by the GPS position and
the ADCP bottom track mode; and the corrected absolute velocity at each depth bin and each time point
are obtained after rotating the direction of the bottom track velocity vectors [Kaneko et al., 1992]. After compass correction, the absolute velocity proﬁles were averaged within 100 m horizontal range (approximately
15 proﬁles mean) to produce a standard data set with a spatial resolution of 100 m.
The ADCP survey was not along a ﬁxed section, because it was necessary to make other observations concurrently. However, the ship tracks were close enough to each other (Figure 1c) that we were able to project
the data to a common section normal to the coastal line of QS (the heavy line in Figure 1c, with a length of
19.1 km and an averaged depth of 51 m). The projection of ADCP data from different sections to the common section is based on the assumption that there is no change of velocity along the direction of the QS(the thin lines in Figure 1c). The northward and eastward components of velocity projected to the common
section were then converted into the along-strait (u) and cross-strait (v) velocity components (along-strait
positive direction is set to 8 anticlockwise rotated from the east and across-strait positive direction is set to98 anticlockwise rotated from the east). Finally, u and v were averaged over a horizontal distance of 200 m
to make a 200 m 3 2 m grid data set for further analysis.
Two water-level meters (HOBO U20-001-02) were deployed at the bottom on both sides of the strait. However, only the water-level meter on the northern side of the QS was recovered. During the bad weather
period, we recovered the water-level meter for safety, thus, the water level data have a gap from 10 to 17
March. The data-sampling rate was set to 6 per hour for the water level meter. The raw data were smoothedthrough a running mean of 1 h (Figure 2a). The hourly water level anomaly data are ﬁnally used in the harmonic analysis. The accuracy of the hourly water level data are expected to be higher than 0.57 cm.

3. Modified Tidal Harmonic Analysis
In a coastal region, an observed current velocity v includes both tidal and subtidal components and can be
written in the following form:
vðtÞ5v0 ðtÞ1

MX
ðxj cosrj t1yj sinrj tÞ

(1)

j51

where j denotes tidal constituents, xj ; yj are harmonic constants, and v0 is the residual current. The time t is
the Greenwich Mean Time (8 h behind the local time), and rj is the angular frequency. M is total number oftidal constituents. Assuming v0 is constant during observation period and introducing a time-dependent variable n to consider the ﬁtting error, equation (1) can be arranged as a set of linear equation systems as follows:
 
 

  v0   vðt1 Þ   1
cosr1 t1
sinr1 t1    cosrM t1
sinrM t1     nðt1 Þ 
 

 

  x1      nðt2 Þ  vðt2 Þ   1
cosr1 t2
sinr1 t2    cosrM t2
sinrM t2   y1      
3 1    : (2)    5         nðt vðt
cosr1 tN21 sinr1 tN21    cosrM tN21 sinrM tN21
   
N21 Þ   1
N21 Þ  
  xM      nðt Þ  vðt Þ   1
cosr1 tN
sinr1 tN   cosrM tN
sinrM tN
N
N
y 
M

Where t is sampling time and N is data number. This linear equation system can be written in the simple
form asZHU ET AL.

C 2014. American Geophysical Union. All Rights Reserved.
V

8042

Journal of Geophysical Research: Oceans

10.1002/2014JC009855

y5Ex1n;(3)where y is the vector for observed velocity, E is the parameter matrix, n is the vector for velocity residues
and errors that cannot be explained by the tidal constituent and mean (time independent) terms, and x isthe vector containing v0 and the harmonic constants of tidal constituents. For a classical tidal harmonic
analysis, the solutions for x are obtained using a least square ﬁtting method by minimizing the residue of
jy2Exj. Two well-known classical tidal harmonic analysis tools used with MATLAB is ‘‘T_Tide’’ [Pawlowicz
et al., 2002] and ‘‘UTide’’ [Codiga, 2011]. As shown later, for the case when sampling time is limited and theleast squares matrix become ill-conditioned, the minimizing of the residue of jy2Exj without any special
consideration may induce unexpected large errors in the solutions for the harmonic function.
Equation (3) is actually a typical linear inverse system. Solving this overdetermined system is a mathematicalproblem. Instead of a least square ﬁtting method that minimizes n, the expected solution of equation (3)can be also obtained by the ‘‘tapered least squares’’ method [Wunsch, 1996], which does not necessarily
give a minimum n but gives a reasonable solution for equation (3). The objective function for the tapered
least squares method is given byJ5ðy2Ex ÞT ðy2Ex Þ1 a2 x T x ;(4)where a2 is the weighting factor and is called a trade-off parameter. By varying the size of a2 , we can gain
some inﬂuence over the norm of the residuals relative to that of x.
If the operator matrix E has extremely small singular values, equation (3) becomes ill-posed and the solution
can become unstable. Let Dx be the error of x corresponding to the noise n. By introducing kDxk=kxk and
knk=kyk, the effect of noise is expressed by:
kDxk
knk
 CondðE Þ
;
kxk
kyk

CondðE Þ5

k1
1
kK

(5)

where Cond(E) is the condition number of matrix E, k1 is the maximum singular value, kK is the minimum
one singular value. For the tapered least squares method, the condition number Cond(E0 ) is given by replacing 1=kK using kK =ðk2k 1a2 Þ [Yamaguchi, 2005].
 
Cond E ’ 5k1 3

kK
k1
 :
ðkK 2 1a2 Þ kK

(6)

When a2 increases, the Cond(E0 ) decreases and the solution becomes insensitive to noises. However, if a2
increases too much, the solution will be excessively smoothed. So we have to ﬁnd an optimal a2 that canadjust the signal to noise ratio and give a reasonable physical solution.
There are a variety of methods to determine the optimal weighting factor a2. In the present study, we use
the L-curve method [Hansen and O’Leary, 1993] to optimally determine the a2 (Figures 3a and 3b). This factor provides a constraint to the vector norm of the solution x~ . The expected solution is determined by minimizing the above objective function J as follows:21 T
x~ 5 E T E1 a2 I
E y;

(7)

where I is the unit matrix.
Using the above inverse technique to estimate the harmonic constants of tidal constituents is different
from the classical tidal harmonic analysis method, and we will refer it as the ‘‘modiﬁed tidal harmonic analysis method’’ in this study.
Figures 2d and 2e show the shipboard ADCP sampling could not cover entire phase range of the K1 and P1
tide, suggesting the impossibility to extract them from those data using a normal tidal harmonic analysis
method unless inference is used. Therefore, we use the amplitude ratio of the O1 and K1 amplitude and theirphase difference from tide gauge records to infer K1 from O1 tidal constituent [Godin, 1972]. The same
method is also applied to the inference of P1 from O1 tidal constituent. Consequently, the equation (3) inthis study contains 33 equations (N533, i.e., repeat observations at each grid location) and 9 unknowns (v0and O1, M2, S2, MSf tidal constituents).

ZHU ET AL.

C 2014. American Geophysical Union. All Rights Reserved.
V

8043

Journal of Geophysical Research: Oceans
6

0.8

5

15
Condition number

1

k(α)

||y−Ex||

4
0.6

3

0.4
2
0.2
0

10.1002/2014JC009855

10

5

1

(a)
0

0.2
||x||

0

0.4

(b)
0

1

2

1

2

α

3

4

5

3

4

5

0

(c)
0

1

2

α

3

4

5

0.5
0.5
Total tidal energy

0.4

||x||

0.3
0.2
0.1
0

0.4
0.3
0.2
0.1

(d)
0

1

2

α

3

4

5

0

(e)
0

α

Figure 3. (a) L-curve, (b) curvature diagrams, (c) condition number versus a, (d) norm of solution versus a, (e) total tidal energy (sum ofsquares of the amplitudes) versus a. Red ‘‘solid circle’’ indicates the maximum curvature point (optimal a). The blue broken lines in Figures
3d and 3e indicate the true values.

To conﬁrm how well we can extract harmonic constants of K1 tidal constituent from the shipboard ADCP
data set and the possibility of decomposition between M2 and S2, we carry out the following simulation.
First, we obtained a 364.4 day long hourly tidal level time series for 1997 at Haikou tidal gauge station near
our observation site from Hawaii sea level center (these is not record for 2013). The harmonic analysis using
conventional method with a nodal tide correction [Foreman et al., 2009] from this data give us the harmonic
constant for 59 constituents (SSA, MSM, MM, MSf, Mf, ALP1, 2Q1, SIG1, Q1, RHO1, O1, TAU1,BET1, NO1, CHI1,
P1,K1, PHI1, THE1, J1, SO1, OO1, UPS1, OQ2, EPS2,2N2, MU2, N2, NU2, M2, MKS2, LDA2, L2, S2, K2, MSN2, ETA2, MO3,
M3, SO3, MK3, SK3, MN4, M4, SN4, MS4, MK4, S4, SK4, 2MK5, 2SK5, 2MN6, M6, 2MS6, 2MK6, 2SM6, MSK6,3MK7, M8)
and the mean value (Z0), in which the O1, K1, M2, S2 are the largest four tidal constituents among them. The
tide amplitude ratio of P1 to O1 is 0.24, and that of K2 to M2 is 0.38, respectively. We synthesized an hourly
tidal level time series for 2013 using amplitudes and phases of the most important 22 constituents (MSf, O1,
K1, M2, S2, P1, K2, 2Q1, Q1, RHO1, NO1, J1, 2N2, MU2, N2, L2, K2, MO3, SO3, MK3, MK4, M6) and Z0. Then we performconventional harmonics analyses with different inferences, using T_Tide, UTide Matlab software and our
modiﬁed tidal harmonics analyses method, respectively. The results are shown in the Table 1. It must be
noted that the inference of K1 constituent using O1 constituent (their magnitude ratio is 0.80, and theirphase difference is 54.80) is kept in all these cases. The T_Tide can only handle one inference at a time, andcannot do the inference using more than two constituents simultaneously, while the UTide can. In order tomake T_Tide to be able to infer two constituents simultaneously, we change the original equations by using
the amplitude ratio and their phase difference from tide gauge records to infer K1&P1 from O1 and K2&S2
from M2 tidal constituent and form a set of new equations. Then we use the T_Tide code to solve the new
equations.
Four main points can be known from Table 1: (1) If using the hourly data from 8 March to 2 April, the inference of constituents K1&P1 from O1 and that of K2 from S2 improve the results of O1, and the result does notdepend analysis methods; (2) If using hourly daytime data from 8 March to 2 April, the results of MSf, O1, K1,
M2, S2, Z0 are still good with inference of constituents K1&P1 from O1 and that of K2 from S2. This conclusionZHU ET AL.

C 2014. American Geophysical Union. All Rights Reserved.
V8044

Journal of Geophysical Research: Oceans

10.1002/2014JC009855

Table 1. Mean Value (Z0) and Amplitudes of Five Tidal Constituents (MSf, O1, K1, M2, S2) at Haikou and the Condition Number (Cond. N) in Each Casea,bCase

Z0 (cm)MSF(cm)O1(cm)

K1(cm)

M2(cm)

S2(cm)

Cond. N

True
Hourly, full year, Z0, MSf, O1, M2, S2 included; K1 inferred from O1 (T_Tide)
Hourly, full year, Z0, MSf, O1, M2, S2 included; K1 inferred from O1 (UTide)
Hourly, full year, Z0, MSf, O1, M2, S2 included; K1 inferred from O1 (modify)
Hourly, 8 Mar. to 2 Apr., Z0, MSf, O1, M2, S2 included; K1 inferred from O1 (T_Tide)
Hourly, 8 Mar. to 2 Apr., Z0, MSf, O1, M2, S2 included; K1 inferred from O1 (UTide)
Hourly, 8 Mar. to 2 Apr., Z0, MSf, O1, M2, S2 included; K1 inferred from O1 (modify)
Hourly, 8 Mar. to 2 Apr., Z0, MSf, O1, M2, S2 included; K1& P1 inferred from O1; K2 inferred from S2 (T_Tide)
Hourly, 8 Mar. to 2 Apr., Z0, MSf, O1, M2, S2 included; K1& P1 inferred from O1; K2 inferred from S2 (UTide)
Hourly, 8 Mar. to 2 Apr., Z0, MSf, O1, M2, S2 included; K1& P1 inferred from O1; K2 inferred from S2 (modify)
Hourly, Daytime, 8 Mar. to 2 Apr., Z0, MSf, O1, M2, S2 included; K1 inferred from O1 (T_Tide)
Hourly, Daytime, 8 Mar. to 2 Apr., Z0, MSf, O1, M2, S2 included; K1 inferred from O1 (UTide)
Hourly, Daytime, 8 Mar. to 2 Apr., Z0, MSf, O1, M2, S2 included; K1 inferred from O1 (modify)
Hourly, Daytime, 8 Mar. to 2 Apr., Z0, MSf, O1, M2, S2 included; K1& P1 inferred from
O1; K2 inferred from S2 (T_Tide)
Hourly, Daytime, 8 Mar. to 2 Apr., Z0, MSf, O1, M2, S2 included; K1& P1 inferred from
O1; K2 inferred from S2 (UTide)
Hourly, Daytime, 8 Mar. to 2 Apr., Z0, MSf, O1, M2, S2 included; K1& P1 inferred from
O1; K2 inferred from S2 (modify)
ADCP_time, Z0, MSf, O1, M2, S2 included; K1 inferred from O1 (T_Tide)
ADCP_time, Z0, MSf, O1, M2, S2 included; K1 inferred from O1 (UTide)
ADCP_time, Z0, MSf, O1, M2, S2 included; K1 inferred from O1 (modify)
ADCP_time, Z0, MSf, O1, M2, S2 included; K1& P1 inferred from O1; K2 inferred from S2 (T_Tide)
ADCP_time, Z0, MSf, O1, M2, S2 included; K1& P1 inferred from O1; K2 inferred from S2 (UTide)
ADCP_time, Z0, MSf, O1, M2, S2 included; K1& P1 inferred from O1; K2 inferred from S2 (modify)
ADCP_time, Z0, MSf, O1, M2 included; K1& P1 inferred from O1; K2 &S2 inferred from M2 (T_Tide)
ADCP_time, Z0, MSf, O1, M2 included; K1& P1 inferred from O1; K2 &S2 inferred from M2 (UTide)
ADCP_time, Z0, MSf, O1, M2 included; K1& P1 inferred from O1; K2 &S2 inferred from M2 (modify)

0.03
0.03
20.03
0.03
0.18
0.03
0.18
0.15
0.02
0.15
8.32
8.88
8.33
0.81

1.25
1.26
1.49
1.26
1.21
1.39
1.21
1.21
1.17
1.21
1.54
1.83
1.55
1.77

53.30
53.19
53.39
53.19
61.08
61.07
61.08
54.03
54.24
54.03
54.47
53.03
54.44
52.82

42.68
42.77
42.93
42.77
49.11
49.11
49.11
43.45
43.62
43.45
43.80
42.64
43.77
42.47

24.77
24.86
25.06
24.86
24.91
25.00
24.91
24.83
24.85
24.83
28.09
27.49
28.08
25.92

13.42
13.44
13.02
13.44
9.94
9.80
9.94
13.67
13.58
13.67
10.24
9.73
10.23
14.35

1.52
1.43
1.43
1.43
1.43
1.43
1.43
2.07
2.07
2.07
5.30
5.30
5.29
5.71

0.62

1.91

53.09

42.69

25.69

14.31

5.71

0.82

1.78

52.79

42.45

25.91

14.33

5.70

1.89
5.68
23.81
11.44
11.33
2.07
3.88
3.88
3.70

12.85
14.63
3.77
11.75
11.73
5.69
6.24
6.38
6.12

63.04
60.68
48.41
60.08
59.55
48.90
49.30
49.30
49.20

50.69
48.79
38.93
48.31
47.88
37.72
39.64
39.64
39.56

35.25
38.43
26.80
35.97
35.93
27.89
29.49
29.43
29.37

10.80
16.33
13.93
19.90
19.75
16.40
15.33
15.30
15.27

14.62
14.62
8.03
17.18
17.18
9.89
8.40
8.40
8.26

a
‘‘True’’ means the original harmonic constants used to synthesize the hourly data for 2013. Hourly, full year, Hourly, 8 Mar. to 2 Apr., Daytime, 8 Mar. to 2 Apr., and ADCP_time
mean Hourly data for 2013, Hourly data from 8 Mar. to 2 Apr., 2013, Daytime (07:00–18:00) data from 8 Mar. to 2 Apr., 2013 and resample data using the shipboard ADCP time, respectively. T_Tide and UTide are the classical tidal harmonic analysis tools by Pawlowicz et al. [2002] and Codiga [2011], respectively. modify means the modiﬁed tidal harmonic analysis
method given by present study.
b
Note: Bold indicates the value with a difference from its true value by 10 cm.

also does not depend on analysis method; (3) If using the data resampled using the ADCP time and using
different inferences in harmonic analysis, the results from the modiﬁed method are all good (the difference
from its true value is smaller than 5 cm). The results from conventional methods are not good (the difference from its true value is larger than 10 cm for some constituents) when using the inferences of constituents K1 from O1, or K1&P1 from O1 and K2 from S2. Using the same data with inference of constituents P1&K1
from O1 and constituents K2&S2 from M2, the results from conventional methods are remarkably improved.
The largest correlation coefﬁcients [Codiga, 2011] calculated by using UTide ranged between 0.22 for S2/M2
and 0.10 for K2/M2, and were smaller than that for K2/S2 (0.98). This may be the reason why the results were
better when S2&K2 were both inferred from M2. (4) The condition numbers from T_Tide and UTide arealways identical due to the same inferences used in the analysis. The condition number from the modiﬁed
method is generally smaller than those from the other two methods. Therefore, we include the constituentsof Z0, MSf, O1, M2, S2, and infer the constituent of K1 from O1 in the modiﬁed tidal harmonics analyses
method, include the constituents of Z0, MSf, O1, M2 and infer the constituents of K1&P1 from O1 and K2&S2
from M2 in the classical tidal harmonics analyses method for the further analyses.
To illustrate the role played by a2 in the modiﬁed tidal harmonics analyses method, we plot the condition
number, norm of solution, and total tidal energy versus a. The condition number, norm of solution, and thetotal tidal energy all decrease with increasing of a (Figures 3c–3e). The solution is close to the true value(Figure 3d), indicating that the total tidal energy loss is small (Figure 3e).
The present ADCP data have both observational error (such as the signals unresolved by the tides) anddata-processing error (such as data projection to a common section). To demonstrate the performance ofthe modiﬁed tidal harmonic analysis method in dealing with such errors, we executed an additional simulation using the tidal level data measured in the QS during the shipboard ADCP surveys. The hourly tidal levelanomaly data were ﬁrst passed through the standard harmonic analysis procedures using ‘‘T_TIDE’’ software
to obtain the solutions for the amplitude and phase of ﬁve major tidal constituents (K1, O1, M2, S2, and MSf).

ZHU ET AL.

C 2014. American Geophysical Union. All Rights Reserved.
V

8045

Journal of Geophysical Research: Oceans

10.1002/2014JC009855

Mean

0.4

Difference(m)

0.2
0
0.8
0.6
0.4
0.2
0
0.8
0.6
0.4
0.2
0
0.8
0.6
0.4
0.2
0
0.8
0.6
0.4
0.2
0
0.8
0.6
0.4
0.2
0
0.6
0.4
0.2
0
−0.2
−0.4
−0.6

O1
K1
M2
S2
MSf
Z0
0

0.04

0.08

0.12

0.16

0.2

0.24

0.28

0.32

0.36

0.4

0.44

0.48

0.52

0.56

0.6

Random error(m)Figure 4. Root mean square difference of ﬁve tidal amplitudes (O1, K1, M2, S2, and MSf) and Z0 obtained from hourly data and 33 undersampled data as a function of the given errors. The top plot shows their mean of the ﬁve tidal constituents. The red bars indicate theresults of the classical tidal harmonic analysis method and the black bars indicate the results of the modiﬁed tidal harmonic analysis
method. The center position of each bar indicates the mean of RMSD from 30 computations, and the vertical bar indicates 1 standard
deviation.

These solutions are considered as standard tidal parameters for this tidal-level anomaly data. Next, we
resampled the 33 sea level anomaly data using the time of the shipboard ADCP surveys. Using such
resampled data (hereafter called undersampled data), we carried out two parallel experiments: one using
the classical tidal harmonic analysis method (Exp. 1) and the other using the modiﬁed tidal harmonic analysis method (Exp. 2). In the Exp. 1, the constituents of Z0, MSf, O1, and M2 are directly solved while the inference of K1&P1 from O1 and that of K2&S2 from M2 are used. In the Exp. 2, the constituents of Z0, MSf, O1, M2,
and S2 are directly solved while the constituent of K1 inferred from O1. By comparing the solution from two
experiments to the standard tidal parameters calculated from the hourly data, we can evaluate the performance of the two methods dealing with the errors in undersampled data.
To consider the inﬂuence of errors in the undersampled data, we artiﬁcially added a random error to eachof the 33 tidal level anomaly data points and let the root mean square of the 33 random errors be equalled
to a given value (Rerr, the abscissa of Figure 4):vﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
ui533
uX
Rerr5terðiÞ2 =33(8)i51

where er(i) (i51,33) is the 33 random errors we added in the undersampled data.
For a given Rerr, we repeated this procedure 30 times. This repetition produced 30 sets of random errors inwhich the random error for individual data points may be different but the root mean square (Rerr) of 33
random errors for all the data is the same. By adding such random errors to the data, we executed Exp. 1
and Exp. 2, repeating each experiment 30 times.

ZHU ET AL.

C 2014. American Geophysical Union. All Rights Reserved.
V

8046

Journal of Geophysical Research: Oceans

10.1002/2014JC009855

For each calculation for the undersampled data, we obtained 10 harmonic constants and then calculated
the root mean square difference (RMSD) between the 10 harmonic constants and the standard tidal parameters calculated from hourly data. After calculating a given Rerr 30 times, we have 30 RMSDs and then caneasily calculate their mean and standard deviation. By gradually increasing Rerr, we can see the response of
the mean and the standard deviation of the RMSD from each method to increasing of Rerr (Figure 4).
When the random error of Rerr is 0.1 m, both methods give a nearly 0.05 m RMSD for MSf, O1, K1, M2 and S2,and Z0 constituents and their mean, which is approximately 4% of the tidal level range during the spring
tidal period (Figure 4). When the Rerr increases, the mean and standard deviation of the RMSDs for Exp. 1
and Exp. 2 all increase slowly, and the RMSDs for Exp. 1 are slightly larger than those for Exp. 2 (top plot of
Figure 4). Consequently, it can be concluded that the modiﬁed tidal harmonic analysis method is able towork as well as the classical tidal harmonic analysis method in controlling the contamination to the solution
from measurement errors contained in undersampled data although the modiﬁed tidal harmonic analysis
method uses less tidal constituents for inference than the classical tidal harmonic analysis method.

4. Results
4.1. Shipboard ADCP Data
The vector plots of depth-averaged current velocity obtained by shipboard ADCP from 33 cruises are shownin Figure 5. The ADCP data show only a small change in current direction in each section among the 33
cruises. The mean and standard deviation of currents from the 33 cruises is 16.1 and 8.2 (rotating anticlockwise from the east). The small change in current direction implies that it was reasonable to projectvelocity vectors from different sections onto a common section across the QS.
Vertical sections of the along-strait velocity (u) from the 33 cruises are shown in Figure 6a. Current velocities
were positive (eastward) on the 3rd, 13th, 16th, 18th, 21th–22nd and 28th–29th during the ﬂood tideperiod; and were negative (westward) for 1st–2nd, 4th–12th, 14th–15th, 17th, 20th, 23th–24th, 26th–27th,and 30th–33th during the ebb tide period. The maximum eastward velocity reached 126.1 cm s21 (16th)
and the maximum westward reached 2125.5 cm s21 (12th). During the shipboard ADCP survey from thesouthern side to the northern side of the QS (the blue arrows in each plot show ship direction), the velocities changed directions due to the changes of tidal phases (19th, 25th, and 31th). Signiﬁcant vertical shears
were noted in the deep region around two grooves in the QS.
The cross-strait velocity (v) (Figure 6b) is weaker than the along-strait velocity (Figure 6a). The velocity
ranged from 220 cm s21 to 20 cm s21, and exceeded 30 cm s21 only in a very small area in the surface
layer (29th–30th, 31th, and 33th). The velocity data (Figure 6) indicate that the tidal current is dominant in
the along-strait direction.
4.2. Tidal Currents
Depth-averaged tidal ellipses of O1, K1, M2, S2, and MSf tidal constituents are shown in Figure 7a together
with the residual current estimated by the modiﬁed tidal harmonic analysis method using depth-averaged
shipboard ADCP velocity data along the common section. The diurnal tidal current O1 and K1 tidal constituent are dominant among the ﬁve tidal constituents, and the semidiurnal tidal current S2 is the smallest
among the ﬁve tidal constituents. The amplitude of O1, K1, M2, S2, and MSf tidal constituents averaged along
the common section has a ratio of 1:0.79:0.42:0.27:0.29. Because of the use of the inference of K1 from O1
tidal constituent, the horizontal distribution of O1 and K1 tidal constituents shows the same structure along
the section (Figure 7a). The ratio of major axis of the tidal current ellipses M2 and S2 is 1:0.64, which is close
to that of their tidal level amplitude of 1:0.54. The major axis direction of the tidal ellipses of the ﬁve constituents is generally along the strait, showing a mean angle 6 its standard deviation (STD) of 11.965.7 ,11.965.7 , 13.167.0 , 23.668.5 , and 18.3610.3 (rotating anticlockwise from the east) for the O1, K1, M2,
S2, and MSf tidal constituents, respectively. The residual current is larger in northern part than in southernpart, with a mean westward velocity 6 STD of 26.062.1 cm s21 and a mean direction (from the east) 6 STD
of 203.068.0 .
The vertical distribution of the O1 and K1 tidal constituents (Figure 8a) show a uniform structure over thecommon section; the lengths and orientations of the major axis of tidal ellipses change within a small range,
presenting a mean amplitude and direction 6 STD of 55.4 6 6.1 cm s21 and 13.7 6 5.9 for the O1 andZHU ET AL.

C 2014. American Geophysical Union. All Rights Reserved.
V

8047

Distance(km)Journal of Geophysical Research: Oceans
35

Distance(km)Distance(km)Distance(km)100cm/s

100cm/s

100cm/s

1

2

3

4

5

100cm/s

100cm/s

100cm/s

100cm/s

100cm/s

6

7

8

9

10

100cm/s

100cm/s

100cm/s

100cm/s

100cm/s

11

12

13

14

15

100cm/s

100cm/s

100cm/s

100cm/s

100cm/s

16

17

18

19

20

100cm/s

100cm/s

100cm/s

100cm/s

100cm/s

21

22

23

24

25

100cm/s

100cm/s

100cm/s

100cm/s

100cm/s

26

27

28

29

30

100cm/s

100cm/s

100cm/s

20
15
10

0
35
30
25
20
15
10

0
35
30
25
20
15
10
5
0
35
30
25
20
15
10
5

Distance(km)100cm/s

25

5

0
35
30
25
20
15
10
5

Distance(km)100cm/s

30

5

0
35
30
25
20
15
10
5

Distance(km)10.1002/2014JC009855

0
35

0

5

10 15 20 25 30 35 40 0

Distance(km)30

5

10 15 20 25 30 35 40

Distance(km)

25
20
15
10
5
0

31
5

32
10 15 20 25 30 35 40 0

533
10 15 20 25 30 35 40 0

5

10 15 20 25 30 35 40

Figure 5. Vector plots of the depth-averaged current measured by the shipboard ADCP. The blue vectors show the original velocity while the red vectors show their projected current onto
the common section. The blue arrow indicates the direction of the ship movement. The bold number on the bottom-left corner of each plot is the cruise number. The tidal phase is shownwith a red dot on the sea level anomaly plot in the top right corner of each plot. The parallel gray lines indicate the axis in which the original velocity is projected onto the common section.

43.1 6 4.8 cm s21 and 13.7 6 5.9 for the K1 tidal constituents, respectively. Even near the bottom, the mean
amplitude of the O1 and K1 tidal constituents is still approximately 40 cm s21. The ratios of the minor axis
length to the major axis length (ellipticity) for the O1 and K1 tidal constituents are 0.026, indicating that thediurnal tidal current is dominated by the along-strait component.

ZHU ET AL.

C 2014. American Geophysical Union. All Rights Reserved.
V

8048

Journal of Geophysical Research: Oceans

10.1002/2014JC009855

Depth(m)

(a)
20
40
60

Depth(m)80

Depth(m)Depth(m)

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

60

20
40
60
80
20
40
60
80

Depth(m)

2

40

80

20
40
60
80

Depth(m)

1

20

20
40
60
80

29

Depth(m)

0
20

2

30
4

6

8

10 12 14 16

0

2

4

Distance(km)68

10 12 14 16 18

Distance(km)40−140 −120 −100 −80 −60 −40 −20

60
80

31
0

2

32
4

6

8

10 12 14 16

Distance(km)

0

2

6

8

10 12 14 16

Distance(km)

0

2

20

40

60

80

100 120 140

u(cm/s)

33
4

0

4

6

8

10 12 14 16 18

Distance(km)Figure 6. Vertical sections of (a) the along-strait velocity and (b) the cross-strait velocity, respectively, for the common section (looking toward the east). The positive and negative values
in Figure 6a indicate the eastward and westward currents, respectively; while those in Figure 6b indicate the northward and southward currents, respectively. The contour interval is
10 cm s21. Thick contours show zero velocity. The black arrows show the directions of the ship movement. The tidal phase is shown as a red dot on the sea level anomaly plot in the bottom right corner of each plot. The cruise number is shown on the bottom left corner of each plot.

The M2 tidal constituent is stronger than the S2 tidal constituent (Figure 8a). The sectional mean of major
axis length is 23.5 cm s21 (minor axis length 2.1 cm s21) for the M2 and 15.1 cm s21 (minor axis length
1.4 cm s21) for the S2 tidal constituents, respectively. The sectional mean of ellipticity is 0.09 for the M2 and0.10 for the S2 tidal currents; both are larger than those of the diurnal tidal constituents. The vertical distributions of M2 and S2 tidal currents also show a uniform structure over the common section, being similar tothe diurnal tidal currents. The sectional mean direction of the major axis 6 its STD is 13.1 6 7.0 for the M2
and 23.6 6 8.5 for the S2 tidal ellipses, respectively; both are close to the strait direction.
The MSf tidal current is strong in the deep central part of the QS (Figure 8a). The amplitude of the MSf tidal
current is large in the surface layers and decreases towards the bottom. The mean amplitude of the MSf
tidal current is 19.4 cm s21 and 15.1 cm s21 for the layers above and below 40 m, respectively. The majoraxes of MSf tidal ellipses are generally along the direction of QS in the center, but shift east-west in the shallow regions.
The along-strait component of residual current is all westward, large in the surface layer in the northernstrait and decreasing with depth (Figure 8a). The maximum residual current is 217.4 cm s21 at the surface
in the northern QS. The sectional mean residual current is 26.0 cm s21. The velocity is about 25.0 cm s21
at 60 m depth. The cross-strait component of residual current is much smaller than the along-strait component (ﬁgure not shown).

ZHU ET AL.

C 2014. American Geophysical Union. All Rights Reserved.
V

8049

Journal of Geophysical Research: Oceans

10.1002/2014JC009855

Depth(m)

(b)
20
40
60

Depth(m)80

Depth(m)Depth(m)

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

60

20
40
60
80
20
40
60
80

Depth(m)

2

40

80

20
40
60
80

Depth(m)

1

20

20
40
60
80

29

Depth(m)

0

2

20

30
4

6

8

10 12 14 16

02

Distance(km)

4

6

8

10 12 14 16 18

Distance(km)

40
−50

60
80

31
0

2

32
4

6

8

10 12 14 16

Distance(km)

0

2

−40

4

6

8

10 12 14 16

Distance(km)

0

2

−30

−20

−10

0

10

20

30

40

50

v(cm/s)

33
4

6

8

10 12 14 16 18

Distance(km)Figure 6. Continued.

The volume transport through the section calculated from the residual current is 20.056 Sv toward thewest. The ADCP measurements obtained velocity data only in the region deeper than 10 m. The volumetransport between two sides of the common section and the bank of QS were therefore not measured. The
residual current for such regions were obtained by extrapolation, which produces an additional volumetransport of 20.009 Sv. Therefore, the total volume transport through the QS due to the residual current is
20.065 Sv during our survey period.
For comparison, we also show the depth-averaged tidal ellipses (Figure 7b) and vertical distribution (Figure 8b)of the tidal ellipses for the O1, K1, M2, S2, and MSf tidal constituents and the residual current obtained by the classical tidal harmonic analysis method (v0, MSf, O1, M2 included; K1 & P1 inferred from O1; K2 & S2 inferred from M2).
Overall, the tidal ellipses from the classical tidal harmonic analysis are not so different from those calculated by
the modiﬁed tidal harmonic analysis method for both the depth-averaged horizontal distribution and the vertical distribution, expect for the northern part of QS, where the tidal ellipses are larger than those estimated by
the modiﬁed tidal harmonic analysis method and the residual currents are eastward as shown in the horizontal(Figure 7b) and the vertical distribution (Figure 8b), which have been never reported in the previous studies.

5. Discussion
5.1. Comparisons With Other Studies
Because the observations and data processing may generate some errors that can affect the harmonic constants of the tidal currents, we compared our results with other independent measurements to evaluate
their accuracy.

ZHU ET AL.

C 2014. American Geophysical Union. All Rights Reserved.
V

8050

Journal of Geophysical Research: Oceans(a) 35O1

Distance(km)30K1100cms−1

25

10.1002/2014JC009855

M2
100cms−1

100cms−1

20
15
10
5
0
35

S2

Distance(km)

30

−1

−1

100cms25

Residual Current

MSf
100cms

20
−1

20cms

15
10
5
0

5

10 15 20 25 30 35 40 0

5

10 15 20 25 30 35 40 0

5

10 15 20 25 30 35 40

Distance(km)(b) 35O1

Distance(km)

30

K1
−1

−1

100cms25

M2−1

100cms100cms

20
15
10
5
0
35

S2

Distance(km)

30

−1

−1

100cms25

Residual Current

MSf
100cms

20
−1

20cms

15
10
5
0

5

10 15 20 25 30 35 40 0

5

10 15 20 25 30 35 40 0

5

10 15 20 25 30 35 40

Distance(km)Figure 7. Tidal current ellipses of ﬁve tidal constituents and residual current estimated using depth-averaged shipboard ADCP velocity
data by (a) the modiﬁed tidal harmonic analysis method and (b) the classical tidal harmonic analysis method.

Shi et al. [2002] collected 37 year current observation data around the QS to examine the tidal current. Theirdata included position, depth, length of major axis and minor axis, orientation, and phase of the tidal current ellipses at four sites close to our shipboard ADCP track. However, we could not ﬁnd information on thetime of data collection. We also did not ﬁnd the values for the residual currents and therefore cannot makea quantitative comparison for the residual currents.
Comparisons between the diurnal and semidiurnal tidal currents estimated from our shipboard ADCP datausing the modiﬁed tidal harmonic analysis method (black ellipses) and those from Shi at al. [2002] (red ellipses) are shown in Figure 9 and Table 2. The distances between sites , , , and from Shi et al. [2002]
and our shipboard ADCP track (the common section) are 3.8, 3.1, 3.7, and 1.5 km, respectively (Figure 9a).
The distributions of the depth-averaged diurnal and semidiurnal tidal ellipses in the top 40 m are similar toZHU ET AL.

C 2014. American Geophysical Union. All Rights Reserved.
V

8051

Journal of Geophysical Research: Oceans

10.1002/2014JC009855

Depth(m)

(a) 0
20

O1

40
60

N

80

K1
N

−1

100cmsE−1100cms
E

20

M2

40
60

N

80

S2
N

−1

100cms
E

E

N

80

−5
−5

MSf

60

−15

−5

−10

20

−5

−5

Depth(m)040Residual
Current−1

100cms
E

0

2

4

6

8

10

12

14

−1

100cms−5Depth(m)

0

16

18

0

2

4

6

Distance(km)

8

10

12

14

16

18

Distance(km)(b)
Depth(m)

0
20

O1

40
60

N

80

K1
N

−1

100cmsE−1100cms
E

20

M2

40
60

N

80

S2
N

−1

100cms
E

−5
−1

2

4

6

100cms

8

10

12

14

16

18

0

2

4

6

8

−5

Residual
Current−5

Distance(km)

0

0

−1

E

0

−1

−5
−5

−15

−10

5

−5

0

80

−1

−5

N

−15
−15
−5

MSf

60

−10 −10

−5

20

−15

Depth(m)040−1100cmsE−10Depth(m)

0

10

12

14

16

18

Distance(km)Figure 8. Vertical structures of tidal current ellipses of the ﬁve major tidal constituents and residual current normal to the section estimated by (a) the modiﬁed tidal harmonic analysis
method and (b) the classical tidal harmonic analysis method. The thick lines in the N-E coordinate indicate the strait direction. The upward direction is north.

those of Shi et al. [2002] (Figures 9a and 9c). The vertical variations of diurnal tidal ellipses from the two setsof observations are in agreement with each other (Figure 9b). The root mean square difference (RMSD) of
the major axis, minor axis, and orientation of the tidal current ellipses at the four sites for diurnal tidal current is 13.5 cm s21, 4.8 cm s21, and 6.0 , respectively (Table 2). The vertical variations of semidiurnal tidal
ellipses from the two sets of observations are in good agreement with each other (Figure 9d). The RMSD of
the major axis, minor axis, and orientation of the tidal current ellipses for semidiurnal current is 4.2 cm s21,
3.6 cm s21, and 4.5 , respectively (Table 2).
We also compared the results from the classical harmonic analysis method (v0, MSf, O1, M2 included; K1&P1
inferred from O1; K2&S2 inferred from M2) with the data from Shi et al. [2002] (Figures 9e–9h, Table 2). The
distributions of the tidal ellipses estimated by the classical tidal harmonic analysis are roughly in agreementwith that reported by Shi et al. [2002], not only in the depth-averaged horizontal view (Figures 9e and 9g)
but also in the vertical structures (Figures 9f and 9h). The RMSD given by the classical harmonic analysis
method for the major axis, minor axis, and orientation of tidal current ellipses at the four sites are 16.1 cm s21,

ZHU ET AL.

C 2014. American Geophysical Union. All Rights Reserved.
V

8052

20
40

4

−1

50cms

(b)

Distance(km)

0
35

1

M2

30

Distance(km)

2

3

4
0

100cms

1

20

2

20

3

15
10

40
4

50cms−1

(c)

(d)

0
35

1

(K1 +O1)/2

30
25

1

80

2

3

4
0
20

2

20

3

15
10

40
4

50cms−1

(e)

(f)

35

2

3

4
0

100cms−1

1

20

2

20

3

15
10

60
80

1

M2

30
25

60

100cms−10Distance(km)

80

−1

25

0

60

40
4

(g)

Depth(m)(a)

Depth(m)

0

3

10

5

4

2

15

5

3

100cms

1

20

5

2

−1

25

5

1

(K1 +O1)/230

Depth(m)35

10.1002/2014JC009855

50cms−1(h)

5 10 15 20 25 30 35 40

Depth(m)Distance(km)Journal of Geophysical Research: Oceans

60
80

Distance(km)Figure 9. Comparisons between the diurnal ((K11O1)/2) and semidiurnal (M2) tidal current constituents estimated from the shipboard
ADCP data (black ellipses) and those from Shi at al. [2002] (red ellipses). (a) Horizontal view of depth-averaged diurnal tidal ellipses in the
top 40 m and (b) vertical distribution of diurnal tidal ellipses, (c) and (d) are the same as (a) and (b) but for M2; (e), (f), (g), and (h) are thesame as (a), (b), (c), and (d) but for the results of the classical tidal harmonic analysis method.

4.8 cm s21, 5.8 for diurnal tidal current and 4.2 cm s21, 3.4 cm s21, 5.7 for semidiurnal tidal current, respectively, being almost the same as the corresponding RMSD given by the modiﬁed tidal harmonic analysis.
It is difﬁcult for us to perform a nodal tide correction for the comparison of our results with those in Shi
et al. [2002] because Shi et al. [2002] did not give the observation time of each data and did not mention
how they treat the problem of nodal tide. Suppose the current data at these four sites were observed in thesame year, the RMSD of the major axis between our results and those of Shi et al. [2002] varies with observation year with a range between 12.0 (1997) and19.5 (1969) cm s21 for diurnal tidal currents, between 3.6(1969) and 4.3 (1997) cm s21 for semidiurnal tidal currents, respectively.
The mean phase 6 its STD estimated at the four sites using the modiﬁed tidal harmonic analysis method is163.467.6 for the diurnal tidal currents and 78.2610.6 for the M2 tidal ellipses, respectively; both are close
to those from the classical harmonic analysis method (187.165.5 for the diurnal tidal currents andZHU ET AL.

C 2014. American Geophysical Union. All Rights Reserved.
V

8053

Journal of Geophysical Research: Oceans

10.1002/2014JC009855

Table 2. Comparison of Parameters of Tidal Ellipses
Modified Tidal Harmonic Analysis
MethodShi et al. [2002]
Tide
Constituent(K11O1)/2

M2

Position

Classic Harmonic Analysis Method

Depth (m)Umax
(cm/s)

Umin(cm/s)

horient(Deg.)

hphase(Deg.)

Umax(cm/s)

Umin(cm/s)

horient(Deg.)

hphase(Deg.)

Umax(cm/s)

Umin(cm/s)

horient(Deg.)

hphase(Deg.)

3
10
20
30
3
40
3
10
20
30
40
3
10
20
3
10
20
30
3
40
3
10
20
30
40
3
10
20

57.56
57.89
51.62
46.43
58.13
39.91
78.54
65.53
79.26
74.95
37.87
31.72
45.83
49.62
28.16
29.56
25.52
25.38
32.06
24.09
23.97
27.55
25.54
23.57
17.67
24.30
27.05
25.10

4.11
2.36
4.08
3.12
4.22
2.86
6.38
3.56
4.94
1.88
0.80
4.03
2.52
1.46
3.32
2.27
0.37
1.49
2.58
2.75
3.11
0.31
4.48
1.65
0.59
1.54
1.17
3.38

23.93
10.28
4.85
8.60
10.90
9.83
7.16
16.19
16.76
14.24
15.87
2.43
7.84
21.93
14.74
22.28
20.32
11.25
11.46
12.38
12.22
18.60
14.94
13.82
14.39
1.88
2.92
3.88

125.13114.02112.45111.29
92.61
90.98
86.45
64.73
45.96
47.38
65.21
76.36
44.91
90.09
52.59
42.05
67.53
45.74
101.21
78.01
66.53
23.64
63.38
62.82
78.57
56.27
45.46
6.19

64.02
60.30
53.11
50.42
51.18
42.80
53.68
52.48
51.52
49.49
47.20
47.80
44.58
43.77
23.98
22.66
21.25
20.44
23.56
19.70
26.00
23.86
23.99
23.11
22.82
25.35
24.17
23.45

2.48
0.14
1.03
2.14
6.60
1.21
1.01
0.52
0.15
0.96
1.44
1.06
0.38
2.31
2.29
0.85
0.57
0.18
1.48
1.02
2.23
2.27
2.24
2.41
2.59
5.08
2.79
1.29

13.90
12.12
12.04
11.20
12.70
11.02
23.10
20.80
19.00
18.50
19.12
0.81
5.45
5.95
16.41
13.13
10.94
10.54
6.42
8.41
19.23
17.63
17.52
18.35
18.34
0.79
2.45
3.82

166.10172.42174.00176.23168.64165.12
159.50159.92160.74161.70162.42151.20153.32155.62
54.54
72.37
74.82
75.24
105.12
82.31
82.85
83.08
81.38
80.00
79.28
73.68
75.07
75.14

84.63
80.00
60.03
61.90
70.91
57.02
66.22
63.03
60.42
57.22
54.91
52.42
55.12
47.83
31.36
36.12
26.51
28.10
26.44
26.24
19.91
18.83
18.93
18.48
18.40
19.70
23.39
21.41

0.06
4.18
3.55
6.63
7.10
4.31
1.63
1.46
1.13
0.45
0.14
5.56
1.13
2.79
0.82
2.75
1.34
4.44
0.96
3.46
1.28
0.14
0.33
0.52
1.08
7.38
5.00
1.15

13.50
10.65
12.00
10.81
12.72
9.86
19.50
18.03
17.30
17.63
18.60
1.59
5.50
7.24
14.56
10.17
10.78
11.76
5.75
10.75
19.56
18.91
16.93
16.15
15.64
13.42
5.51
1.31

179.30191.40198.82
197.31186.92189.02185.63184.23184.82
186.35186.50180.51183.82184.8793.58
112.61129.30126.26
119.4
82.31
124.63123.94124.21125.60123.72121.53117.84114.60117.1613.4.5 for the M2 tidal ellipses). However, the phases for both the diurnal and semidiurnal tidal currents are different from those given by Shi et al. [2002] (83.4626.9 for the diurnal tidal currents and56.4623.7 for the M2 tidal ellipses). Since Shi et al. [2002] did not present the details on their treatment ofphase from data in different year, it is difﬁcult for us to make further analysis to identify the causes for thedifference of phase.
The residual currents estimated by the classical tidal harmonic analysis method are mainly westward inside
the strait expect for the northern part (Figure 7b). The eastward residual currents are against to the reports
in previous studies [Shi et al., 2002; Yang et al., 2003; Bao et al., 2005; Chen et al., 2006; Yang et al., 2006;
Chen et al., 2007; Yan et al., 2008].
The synthetic analysis of 37 years of current data by Shi et al. [2002] showed a residual current velocity of
approximately 20–40 cm s21 in the upper 10 m and approximately 10–15 cm s21 near-bottom in winterand spring. Our results show a maximum residual current velocity of 17.4 cm s21 at the surface and a mean
residual current velocity of 9.2 cm s21 in upper 10 m, both of which are slightly smaller than that given by
Shi et al. [2002]. Our results also show that the residual current velocity near the bottom reaches approximately 0–5 cm s21 in the shallow region, which is agreement with the data of Shi et al. [2002], and isapproximately 5 cm s21 in the deep region, which is slightly smaller than that reported by Shi et al. [2002].
Using a three-dimensional numerical model, Chen et al. [2009] showed that the residual currents in the QS
were stronger in the northern area than in the southern area and turned into southward near the southernarea of QS; the current magnitude is about 5–30 cm s21 in the surface layer, 3–20 cm s21 in the middle
layer, and 1–10 cm s21 near the bottom layer. Our analysis also shows that the residual currents are stronger in the northern area than in the southern area and also turn into southward near the southern side of
the section. In addition, our residual current magnitudes are within the range of their results.
Shi et al. [2002] estimated the westward residual volume transport to be approximately 0.2–0.4 Sv in winterand spring using a mean current speed of 20–40 cm s21 at the surface and of 10–15 cm s21 at the bottom,
and an average width of 29.5 km and an average depth of 46 m. Using nonsynchronous current mooring

ZHU ET AL.

C 2014. American Geophysical Union. All Rights Reserved.
V

8054

Journal of Geophysical Research: Oceans

10.1002/2014JC009855

data collected at several stations in different years, Chen et al. [2007] and Yan et al. [2008] also reported awinter westward residual volume transport of 0.055 Sv and 0.020–0.048 Sv, respectively. Our estimation,
based on the synchronous current velocity data and real bathymetry, gives a westward residual volumetransport of 0.065 Sv in spring, which is within the range of these previous reports.
Since the shipboard ADCP sampling could not cover the entire period of K1 tidal constituent, we cannot
directly extract harmonic constants of K1 tidal constituent from our data. We have to use O1 tidal constituent
to interfere K1 tidal constituent using the amplitude ratio of K1 and O1 tidal constituents and their phase difference in sea level record. It is likely that tides and tidal currents of K1 and O1 have almost same spatial distributions in the real ocean except for in a region near an amphidromic point [Guo and Yanagi, 1998] that is
far from the QS [Fang et al., 1999]. Therefore, the inference used in our analysis is likely reasonable, which isalso supported by the good agreement of our estimated diurnal (semidiurnal) tidal ellipses with thosereported by Shi et al. [2002].
5.2. Uncertainty Estimation
Using the 33 repeat shipboard ADCP data sets with the modiﬁed tidal harmonic analysis method, we estimated the tidal currents, residual currents, and associated volume transport through the QS. Here, we introduce the uncertainty in these estimations.
Using equation (7) with the tapered least squares method, we obtained the solution x~ , which is an estimation of the true value x. The uncertainty, P, which is the dispersion of x~ about the true value (sometimes
also called error covariance), can be expressed in the following form [Wunsch, 1996]:nn
21 T o
21 T 21P5 I2E E T E1 a2 I
E
E g :
(9)
R3 I2E E T E1 a2 I
The vector P contains the uncertainty of the lengths of the major axis and minor axis of tidal ellipses for theﬁve tidal constituents and the residual current velocity (Figure 10). The uncertainty for the depth-averaged
tidal currents is larger in the north than in the south (Figure 10a). The uncertainty of the major axis length isalso larger than the uncertainty of the minor axis length. The sectional mean of uncertainty of the major
axis length and its STD is 11.363.7 cm s21, 8.862.9 cm s21, 7.662.6 cm s21, 2.761.8 cm s21, and5.561.7 cm s21 (and of the minor axis length is 0.460.4 cm s21, 0.360.3 cm s21, 0.360.3 cm s21,0.260.2 cm s21, and 0.460.4 cm s21) for the O1, K1, M2, S2, and MSf tidal current constituents, respectively.
The uncertainty of major axis length is approximately 20%, 20%, 32%, 18%, and 33% of the major axis
length of the O1, K1, M2, S2, and MSf tidal current constituents, respectively.
The uncertainty of the residual current velocity is somewhat larger in the south than in the north (Figure 10a)
and also larger in the upper 10 m layer than in the bottom layer (Figure 10b). The maximum uncertainty of
the residual current velocity is 7.3 cm s21, appearing at the surface in the northern part of the QS. The sectional mean of uncertainty for the residual currents is 5.1 cm s21, which leads to an uncertainty of 0.046 Sv forthe volume transport; i.e., approximately 71% of the volume transport (0.065 Sv) due to the residual currents.
5.3. Dynamics for the Residual Currents in the QS
Using the harmonic constants of the ﬁve tidal current constituents, we can reconstruct the tidal currents(ut(x,z,t), vt(x,z,t)) at any position (x, z) within the section at any time (t). The relative contribution of the tidal
currents to the mean kinetic energy (Rmke) and to the mean eddy kinetic energy (Rmeke) can be estimated by
the following two equations [Takikawa et al., 2003]:
Rmke ðx; zÞ5

Rmeke ðx; zÞ5u2t ðx; z; tÞ1vt2 ðx; z; t Þ
u2 ðx; z; tÞ1v 2 ðx; z; t Þ
u2t ðx; z; tÞ1vt2 ðx; z; t Þ
u0 2 ðx; z; tÞ1v 0 2 ðx; z; tÞ

;

(10)

;

(11)

where u(x,z,t) and v(x,z,t) are the observed current velocities; u’(x,z,t) and v’(x,z,t) are the deviations from theresidual currents; and the top bar denotes the temporal average over the observation period.
The distributions of Rmke and Rmeke along the common section, together with the mean along-strait andacross-strait velocities are shown in Figure 11. The contribution of the tidal currents to the total Rmke and theZHU ET AL.

C 2014. American Geophysical Union. All Rights Reserved.
V

8055

Journal of Geophysical Research: Oceans

10.1002/2014JC009855

Depth(m)

(a) 0
20

O1

40
60

N

80

K1
N

100cms−1
E100cms−1E

Depth(m)

0
20

M2

40
60

N

S2
N

−1

100cms

80

E

0

−15

−5

60

N

−5

−5

MSf

40

−5

−5

−10

20

−5

Depth(m)−1

100cmsE

Residual
Current−1

100cms

80

E

0

2

4

6

8

10

12

14

16

18

0

2

4

6

Distance(km)

8

10

12

14

16

18

Distance(km)(b)
Depth(m)

0
20

O1

40
60

N

80

K1
N

−1

100cmsE

100cms−1
E

20

M2

40
60

N

80

S2
N

−1

100cms
E

−5
−1

2

4

6

8

10

12

14

16

18

0

2

4

6

8

−5

Residual
Current−5

Distance(km)

0

0

−1

100cms
E

0

−1

−5

−5

−15

−10

5

−5

0

80

−1

−5

N

−15
−15

−5

MSf

60

−10 −10

−5

20

−15

Depth(m)040−1100cmsE−10Depth(m)

0

10

12

14

16

18

Distance(km)Figure 10. (a) Uncertainty of ﬁve tidal constituents and residual current estimated by the modiﬁed tidal harmonic analysis method using depth-averaged shipboard ADCP data. (b) Vertical structures of uncertainty for the ﬁve major tidal constituents and residual current normal to the section estimated from the shipboard ADCP data. The thin and thick lines indicate
the uncertainty in the major and minor axis direction, respectively.

Rmeke is larger in the central part of the QS than on the southern and northern sides. The sectional means of
Rmke and Rmeke are 0.77 and 0.86, respectively, suggesting that the tidal currents are generally more dominant
than mean current activities and much more dominant than eddy activities. Since only ﬁve tidal current constituents (K1, O1, M2, S2, and MSf) are used to reconstruct ut(x,z,t) and vt(x,z,t), the contributions from other tidal
constituents currents in a more realistic case are likely able to increase them as well as Rmke and Rmeke.
The formation mechanism of the residual current in the QS has been proposed as the tide-induced residual
current by Shi et al. [2002], based on a simulation for tidal currents only. Alternatively, based on simulations
including tidal forcing and winds, buoyant ﬂux, and boundary currents, Gao et al. [2013] proposed thatwind-induced sea level difference between the two entrances of the QS played an important role in driving
the current through the strait. In this study, we cannot examine the role of wind-induced sea level difference along the strait direction. However, based on the ADCP measurements, we can estimate the nonlinear

ZHU ET AL.

C 2014. American Geophysical Union. All Rights Reserved.
V

8056

Journal of Geophysical Research: Oceans
−2

−10

4
−8

−2

−6

−4

−2

Depth(m)

4

−22

−18

40

4

−20

6

−2

−2

−28

20

−2

−30

−2

−22

−34
−32

−32

−4

0

10.1002/2014JC009855

2

−20

60
−18

0.92
0.
94

0.9

2

0.6

6
0.8

0.84
0.87

Depth(m)

94

0

0.96

0.

0.84

0.87

0.93

(c)

80

0.86

4

60

8

0.88

0.8

0.9

0.81

0.8

40

4

92

0.84

0.8

0.

5

0.68

0.92

0.7

81

87

20

0.

0.

0.78

0.72

0.81

0.69

0.88

0

(b)
0.8

(a)

80

(d)
2

4

6

8
10
12
Distance(km)

14

16

18 0

2

4

6

8
10
12
Distance(km)

14

16

18

Figure 11. Vertical distribution of the mean current velocity for (a) u and (b) v and contribution of tidal currents to (c) the mean kinetic energy and to (d) the total eddy kinetic energy.

effects of tidal currents that are the driving forcing of tide-induced residual currents [Shi et al., 2002; Ca’ceres
et al., 2003].
According to Ca’ceres et al. [2003], the momentum balances for the residual currents are:uð@v=@xÞ 1 v ð›v=›y Þ 1 wð@v=@zÞ 1 fu
5 21=qð@p=@yÞ 1 Az ›2 v=›z2 1 Ah ð@ 2 v=@x 2 Þ 1 Ah ›2 v=›y 2 ;(12)uð@u=@xÞ 1 v ð›u=›y Þ 1 wð@u=@zÞ 2 fv

5 21=qð@p=@xÞ 1 Az ›2 u=›z2 1 Ah ð@ 2 u=@x 2 Þ 1 Ah ›2 u=›y 2 :(13)Equation (12) is for the across-strait direction, whereas equation (13) is for the along-strait direction. The ﬁrst
three terms on the left hand side of the two equations denote the nonlinear effects of tidal current; the forthterm on the left hand side of both equations is the Coriolis term; the ﬁrst term on the right hand side of both
equations is the pressure gradient; the second term on the right hand side of both equations is vertical mixing;
the last two terms on the right hand side of the two equations are horizontal mixing. The bars over the top
denote the temporal average over one tidal cycle. Using our data, we can calculate the terms shown in boldfont. In these calculations, we set Coriolis parameter f as 5.0 3 1025 s21, the vertical eddy viscosity coefﬁcient(Az) as 0.01 m2 s21, and the horizontal eddy viscosity coefﬁcient (Ah) as 90 m2 s21. The last two parameters followed those used by Ca’ceres et al. [2003]. We consider the Coriolis term (TCor) as a reference for the other terms
and compare the advective term (Tadv) of vð@v=@yÞ ðvð@u=@yÞÞ, vertical frictional term (Tver) of Az ð@ 2 v=@z2 Þ ðAz
ð@ 2 u=@z2 ÞÞ and horizontal frictional term (Thor) of Ah ð@ 2 v=@y 2 Þ ðAh ð@ 2 u=@y 2 ÞÞ to the Coriolis term.
In the across-strait direction (Figure 12), the TCor (Figure 12b), which resembled the distribution of along-strait
mean ﬂow (Figure 11a), is of the same order as the Tadv (Figure 12c), Tver (Figure 12d), and Thor (Figure 12e).
The sectional mean for the absolute values of the TCor, Tadv, Tver, and Thor are 3.09, 10.07, 8.40, 8.71 (31026
ms22), respectively. The approximate ratio of these terms is TCor:Tadv:Tver:Thor 5 1.00:3.26:2.72:2.82. Assuming
0.01 m as the sea level difference across the QS [Shi et al., 2002], the pressure gradient due to it is approximately 5 3 1026 ms22. Therefore, the momentum balance between the sea level difference across the strait
and the Coriolis forcing is not exactly held in the QS. The nonlinear effect of tidal currents (Tadv) becomes thelargest among all these terms. The frictional forcing due to the shear of residual currents and the large viscosity coefﬁcients in both horizontal and vertical directions are also important. Therefore, although there is someZHU ET AL.

C 2014. American Geophysical Union. All Rights Reserved.
V

8057

Journal of Geophysical Research: Oceans

uncertainty in the viscosity coefﬁcientsused in our estimation, the mechanism of tide-induced residual current
given by Shi et al. [2002] is supported
by our calculation.

0.
4

0.
2
0.4

0

20

0.2

Depth(m)

0
0

40
60
80
0

10.1002/2014JC009855(a)

In the along-strait direction (Figure13), the Coriolis term (Figure 13b) is20
smaller than the Tadv (Figure 13c), Tver
40
(Figure 13d), and Thor (Figure 13e) by60
one order of magnitude. The sectional means for the absolute value of
80 (b)
the TCor, Tadv, Tver, and Thor are 1.24,
010−10
10.93, 13.94, 11.81 (31026 ms22),
20
respectively. The approximate ratio of
40
these terms is TCor:Tadv:Tver:Thor 560
1.00:8.81:11.24:9.52. Apparently, the80 (c)
Coriolis term is negligible in the0
momentum balance in the along21
3
20
strait direction. The other three terms
98
all show positive values and the
15
40
1
6
12
advective term on the left hand side60
12of equation (13) cannot be balanced
80 (d)
with the sum of vertical and horizon0
3
3
tal frictional terms on the right hand12
12
20
side of equation (13). Therefore, a40
9negative (westward) along-strait pres615
6sure gradient is necessary for the bal60
ance of momentum in the along80 (e)strait direction. In other words,
0
2
4
6
8
10
12
14
16
18
although the tidal rectiﬁcation is very
Distance(km)important to the formation of the
westward residual current in the QS
Figure 12. Contour plots for the across-strait direction of (a) the mean divergenceas suggested by Shi et al. [2002], theof the across-strait ﬂow (@v=@y, 1024 3 s21), (b) Coriolis term (fu, 1026 3 m s22), (c)advective term (vð@v=@yÞ , 1026 3 m s22), (d) vertical frictional term (Az ð@ 2 v=@z 2 Þ,
sea level difference between two
1026 3 m s22), and (e) horizontal frictional term (Ah ð@ 2 v=@y 2 Þ, 1026 3 m s22).
entrances of the QS is also necessary
for the westward residual current in
the QS as suggested by Gao et al. [2013]. Because the sea level difference between two entrances of
the QS changes with the seasonal variation in wind, a summer ﬁeld study is required to conﬁrm ourhypothesis.
−6

−8

−2

−2

Depth(m)

−4

0

−2

−1

0

20

−20

Depth(m)

10

−10

−40

10

−2

0

15

Depth(m)

6

6

15

3

Depth(m)

15

12

6

15

6

12

6. Summary
Using 33 repeat shipboard ADCP data obtained in the spring of 2013 and a modiﬁed tidal harmonic analysis
method, we estimated the harmonic constants of ﬁve major tidal currents (K1, O1, M2, S2, and MSf), the residual current and the associated volume transport through a section across the QS.
We successfully modiﬁed the classic harmonic analysis method using the inverse technique. Through a simple simulation using sea level data and comparisons with previous observations, we conﬁrmed that themodiﬁed tidal harmonic analysis method works well as the classic harmonic analysis method for tidal currents and works slightly better than the classic harmonic analysis method for the residual currents. In addition, the modiﬁed tidal harmonic analysis method is able to use less tidal constituents for inference that is
usually necessary for undersampled data. To our knowledge, our study is the ﬁrst application of the modiﬁed tidal harmonic analysis method using the inverse technique to solve an overdetermined linear equation
system to obtain the solutions for harmonic constants of tidal constituents. Baschek and Send [2001] used atidal inverse model to estimate the transport through the Strait of Gibraltar, but their tidal inverse modelZHU ET AL.

C 2014. American Geophysical Union. All Rights Reserved.
V

8058

Journal of Geophysical Research: Oceans

was different from our modiﬁed
tidal harmonic analysis method
because they assumed a fourthorder polynomial plus an exponential function in their model
functions, whereas we did not.

0.4

(a)
−2

−1

−1

40
−2

60

−10

−10

30

10

−10

10

10

40

10

10

20

−20

−1

(b)
30

60

40

16

20

8
8

40

24

(c)

80
0

16
32

60

16

(d)
40

5
5

10

40

20

20

80
0

25

20

40

10

Depth(m)−3−1Depth(m)

−4

0

1

20

80
0

Depth(m)

0.4

60
80
0

20

25

Depth(m)

−0.4

0

40

−10

Depth(m)

0
20

60

45

(e)

80
0

2

4

6

8
10
12
Distance(km)

14

16

18

Figure 13. Contour plots for the along-strait direction of (a) the mean divergence of thealong-strait ﬂow (@u=@y, 1024 3 s21), (b) Coriolis term (fv, 1026 3 m s22), (c) advective
term (vð@u=@yÞ Þ, 1026 3 m s22), (d) vertical frictional term (Az ð@ 2 u=@z 2 Þ, 1026 3 m s22),and (e) horizontal frictional term (Ah ð@ 2 u=@y 2 Þ, 1026 3 m s22).

Acknowledgments
This study is supported by the National
Basic Research Program of China(2011CB403503, 2011CB409803), theNational Natural Science Foundation ofChina (41276095, 41476020, 41276028
and 41321004), the Scientiﬁc Research
Fund of SIO under grants JT1207 andJT1402, and the project of Global
Change and Air-Sea interaction (GASI03-01-01–02). The data used in thispaper are available on request (contact
Dr. Zhu X.-H. at xhzhu@sio.org.cn).
Comments from two anonymous
reviewers were helpful in improving theoriginal manuscript.

ZHU ET AL.

10.1002/2014JC009855

With the modiﬁed tidal harmonic
analysis method, we obtained
the following results. The major
axis of tidal ellipses of ﬁve constituents is generally in the direction along the strait. The diurnal
tidal current is dominant among
ﬁve tidal constituents, and theamplitude of the diurnal tidal
current is about two times larger
than that of the semidiurnal tidal
currents. The ratio of the amplitude of O1, K1, M2, S2, and MSf
tidal currents averaged along thesection is 1.00:0.79:0.42:0.27:0.29.
The residual current is westward,
decreasing from the surface layer
to the bottom layer. The westward volume transport of theresidual current with its uncertainty through the QS is estimated to be 0.065 60.046 Sv.
This value is smaller than thevalue of 0.2–0.4 Sv reported by
Shi et al. [2002] and larger thanthe 0.055 Sv reported by Chenet al. [2007] and the 0.020–
0.048 Sv reported by Yan et al.
[2008].

Through dynamic analyses, we compared the tidal current energy with the mean kinetic energy and mean
eddy kinetic energy and showed that the tidal current activities are more dominant than mean current andeddy activities. This result is similar to that obtained in Tsushima Straits [Takikawa et al., 2003]. We also calculated the magnitude of the advective, frictional terms, and the Coriolis forcing to determine the nonlinear
tidal rectiﬁcation using the momentum equation. The ratio of TCor:Tadv:Tver:Thor in the across-strait andalong-strait direction is 1.00:3.26:2.72:2.82 and 1.00:8.81:11.24:9.52, respectively. Our results indicate that
not only the tidal rectiﬁcation but also the sea level difference between two entrances of the QS is important to the formation of westward residual current in the QS.

References
Bao, X.W., Y. J. Hou, C.S. Chen, F. Chen, and M.C. Shi (2005), Analysis of characteristics and mechanism of current system on the west coastof Guangdong of China in summer, Acta Oceanol. Sin., 24(4), 1–9.
Baschek B., and U. Send (2001), Transport estimates in the Strait of Gibraltar with a tidal inverse model, J. Geophys. Res., 106(C12), 31,033–
31,044.
Ca’ceres, M., A. Valle-Levinson, and L. Atkinson (2003), Observations of cross-channel structure of ﬂow in an energetic tidal channel, J. Geophys. Res., 108(C4), 3114, doi:10.1029/2001JC000968.
Cao D. M., and G. H. Fang (1990), A numerical model of the tides and tidal currents in Beibu Gulf [in Chinese with English abstract], Oceanol.
Limnol. Sin., 21(2), 105–113.
Chen B., J. H. Yan, D. R. Wang, and M. C. Shi (2007), The transport volume of water through the Qiongzhou Strait in the winter season [inChinese with English abstract], Periodic. Ocean Univ. China, 37(3), 357–364.

C 2014. American Geophysical Union. All Rights Reserved.
V

8059

Journal of Geophysical Research: Oceans

10.1002/2014JC009855

Chen C. L., P. L. Li, M. C. Shi, J. C. Zuo, M. X. Chen, and H. P. Sun (2009), Numerical study of the tides and residual currents in the Qiongzhou
Strait, Chin. J. Oceanol. Limnol., 27(4), 931–942.
Chen, D. S., B. Chen, J. H. Yan, and H. F. Xu (2006), The seasonal variation characteristics of residual currents in the Qiongzhou Strait [in Chinese with English abstract], Trans. Oceanol. Limnol., 2, 12–17.
Codiga, D. L. (2011), Uniﬁed tidal analysis and prediction using the UTide Matlab functions, Tech. Rep. 2011-01, 59 pp., Grad. Sch. of Oceanogr., Univ. of R. I., Narragansett. [Available at http://www.po.gso.uri.edu/~codiga/utide/utide.htm.]
Fang, G., Y. K. Kwok, K. Yu, and Y. Zhu (1999), Numerical simulation of principal tidal constituents in the South China Sea, Gulf of Tonkin
and Gulf of Thailand, Cont. Shelf. Res., 19, 845–869.
Foreman, M. G. G., J. Y. Cherniawsky, and V. A. Ballantyne (2009), Versatile harmonic tidal analysis: Improvements and applications, J. Atmos.
Oceanic Technol., 26, 806–817, doi:10.1175/2008JTECHO615.1.
Gao, J., H. Xue, F. Chai, and M. Shi (2013), Modeling the circulation in the Gulf of Tonkin, South China Sea, Ocean Dyn., 63, 979–993, doi:
10.1007/s10236-013-0636-y.
Godin, G. (1972), The Analysis of Tides, 264 pp., Univ. of Toronto Press, Toronto, Ontario, Canada.
Gordon, R. L. (1989), Acoustic measurement of river discharge, J. Hydraul. Eng., 115(7), 925–936.
Guo, X., and Yanagi, T. (1998), Three dimensional structure of tidal current in the East China Sea and the Yellow Sea, J. Oceanogr., 54(6),651–668.
Hansen, P. C., and D. P. O’Leary (1993), The use of the L-cure in the regularization of discrete ill-posed problems, SIAM J. Sci. Comput., 14,1487–1503.
Kaneko, A., S. Mizuno, W. Koterayama, and R. L. Gordon (1992), Cross-stream velocity structures and their variation of Kuroshio around
Japan, Deep Sea Res., Part A, 39, 1583–1594.
Ke, P. H. (1983), A preliminary analysis of currents and water exchanges in the Qiongzhou Strait [in Chinese with English abstract], J. Trop.
Oceanogr., 2(1), 42–46.
Pawlowicz R., B. Beardsley, and S. Lentz (2002), Classical tidal harmonic analysis including error estimates in MATLAB using T_TIDE, Comput.
Geosci., 28, 929–937.
Shi, M. C., C. H. Chen, F. Huang, and A. L. Ye (1998), Characteristics of tidal current and residual current in the Qiongzhou Straits in period
between end of winter and beginning of spring [in Chinese with English abstract], Acta Oceanol. Sin., 20, 1–10.
Shi, M. C., C. S. Chen, Q. C. Xu, H. Lin, G. Liu, H. Wang, F. Wang, and J. Yan (2002), The role of Qiongzhou Strait in the seasonal variation of
the South China Sea circulation, J. Phys. Oceanogr., 32(1), 103–121.
Takikawa, T., J.-H. Yoon, and K.-D. Cho (2003), Tidal current in the Tsushima Straits estimated from ADCP data by ferryboat, J. Oceanogr., 59,
37–47.
Teledyne RD Instruments Inc. (2013), A Teledyne RD Instruments Datasheet, pp. 2. Poway, Calif. [Available at www.rdinstruments.com.]Wang D., B. Hong, J. Gan, and H. Xu (2010), Numerical investigation on propulsion of the counter-wind current in the northern South China
Sea in winter, Deep Sea Res., Part I, 57, 1206–1221.
Wu, D. X., Y. Wang, X. P. Lin, and J. Y. Yang (2008), On the mechanism of the cyclonic circulation in the Gulf of Tonkin in the summer. J. Geophys. Res., 113, C09029, doi:10.1029/2007JC004208.
Wunsch, C., (1996), The Ocean Circulation Inverse Problem, 296 pp., Cambridge Univ. Press, N. Y.
Wyrtki, K. (1961), Physical oceanography of the southeast Asian waters, NAGA Rep. 2, Scientiﬁc Results of Marine Investigations of the South
China Sea and the Gulf of Thailand, Scripps Inst. of Oceanogr., La Jolla, Calif.
Yamaguchi, K. (2005), Study on advanced system design and data analysis of the coastal acoustic tomography, 160 pp., doctor dissertation,
Higashihiroshima, Japan.
Yan C., B. Chen S. Yang, and J. Yan (2008), The Transportation volume of water through the Qiongzhou Strait in winter season, Trans. Oceanol. Limnol., 1, 1–9.
Yang S. Y., X. W. Bao, C. S. Chen, and F. Chen (2003), Analysis on characteristics and mechanism of current system in west coast of Guangdong Province in the summer [in Chinese with English abstract], Acta Oceanol. Sin., 25, 1–8.
Yang S. Y., B. Chen, and P. L. Li (2006), A study of the characteristics of water transport from the South China Sea into Beibu bay via the
Qiongzhou Strait in summer in terms of temperature and salinity data [in Chinese with English abstract], Trans. Oceanol. Limnol., 1, 1–7.
Yu, M. G., and J. F. Liu (1993), The system and pattern of the South China Sea circulation [in Chinese], Mar. Predict, 10(2), 13–17.
Zu, T. T., J.P. Gan, and S. Y. Erofeeva (2008), Numerical study of the tide and tidal dynamics in the South China Sea, Deep Sea Res., Part I, 55,
137–154.

ZHU ET AL.

C 2014. American Geophysical Union. All Rights Reserved.
V8060