1 Introduction

The atmosphere is constantly in motion and the computation of the displacement of an air particle is important in both trajectory analysis and the designation of numerical models. The Lagrangian approach was developed to trace the motion of air flow by following a fluid parcel as it moves with the flow (Bowman et al., 2013). The fluid parcel can be either air or seawater (van Sebille et al., 2014; Qin et al., 2015). The four-dimensional trajectory of the fluid parcel represents the temporal and spatial motion of the fluid. Passive tracers that are carried by and/or dispersed with the fluid parcels have also received the attention of researchers. These tracers may be radioactive pollution (Stohl et al., 2012; Fei et al., 2014), volcanic ash (Webster et al., 2012), CO₂ emissions (Trusilova et al., 2010), dust (Guo et al., 2013; Notaro et al., 2013), vapor (Zhang and Li, 2014; Li and Dolman, 2016) or ozone in the atmosphere, or tracers (Smith et al., 2016) and drifters (Shchekinova et al., 2016) in the ocean.

The meridional–zonal resultant displacement method (MRDM) is the simplest computation model used in trajectory studies. Bowman et al. (2013) reviewed six trajectory models. Of these models, only the Stochastic Time-Inverted Lagrangian Transport Model (STILT; Lin et al., 2003) is designed for the regional domain and the other models—including the Flexible Particle Dispersion Model (FLEXPART; Stohl et al., 2005), the Hybrid Single Particle Lagrangian Integrated Trajectory Model (HYSPLIT; Draxler and Hess, 1997, 1998; Draxler, 1999), the Lagrangian Analysis Tool Trajectory Model (LAGRANTO; Wernli and Davies, 1997), the Next-Generation Atmospheric Dispersion Model (NAME; Jones et al., 2007), and the Three-Dimensional Trajectory Model (TRAJ3D; Bowman, 1993; Bowman and Carrie, 2002)—have been developed for both the regional and global domains. The MRDM model in spherical coordinates is a convenient selection with which to investigate the trajectories of fluid particles in middle and lower latitudes. In such a coordinate, the latitude circles have a large curvature and the longitude lines are convergent at high latitudes.

It is essential to compute the mean displacement of fluid particles accurately when their motion trajectories are traced. The MRDM model is a traditional method of computing the arrival position of an air particle on a sphere. An air particle without meridional motion will move in a curvilinear manner rather than along a straight line if its displacement at high latitudes is examined by the MRDM. This motion is completely different from Newtonian inertial motion, in which a particle physically moves along a straight line. Similar problems arise when the mean displacement is solved under curved grid coordinates. To solve this polar singularity problem, a plane tangent to the sphere at the pole measured from the map projection is used (Harris and Kahl, 1994) and an exact latitude spacing within 0.25° of the pole has been proposed (Draxler and Hess, 1997).

The sources of error in trajectory models have been investigated by a number of researchers. Stohl (1998) reviewed the Lagrangian trajectory model in terms of the calculation method, precision, and application. The trajectory error varies between examples and a position error of 20% of the travel distance is considered to be typical, limiting the general applicability of the calculated trajectories; errors of the order of 30% are typical of the forecast trajectories. In general, the accuracy of the trajectory computation had been steadily increasing. Bowman et al. (2013) suggested that the input data for trajectory model will affect the accuracy of the model. One of the requirements is for an increased output frequency to match the continuing increase in spatial resolution. This is because the accuracy is highly dependent on the time interval: changing the wind sampling interval from 3 to 1 h reduces the error by about one order of magnitude (Bowman et al., 2013, their Fig. 1). Because the truncation errors from an approximate scheme to solve the particle’s new position are generally small (Bowman et al., 2013), the dominant errors in the trajectory calculations typically arise from errors in the wind velocities themselves or from the limited spatial and temporal resolution of the gridded wind fields (Stohl et al., 2001).

Costa and Blanke (2004) assessed the trajectory errors in the North Atlantic Ocean current and found that the errors increased roughly linearly with the travel time. As a general rule, three-dimensional trajectories limited the trajectory errors and neglecting the vertical velocity component yielded results with unacceptably large errors of about 40% or more of the travel distance. In this case, the vertical shear of the horizontal velocity was the dominant cause of deviations in the horizontal trajectory. Riddle et al. (2006) reviewed many studies and reported errors between 2% and 30% of the total trajectory length for trajectories between several hours and several days in duration. During the International Consortium for Atmospheric Research on Transport and Transformation 2004 campaign, the mean trajectory errors for the balloon trajectory over 2 and 12 h were about 26% and 34% of the flight distance for ECMWF-based and GFS-based trajectories, respectively. An important feature of using balloons is that they allow horizontal trajectory errors to be isolated from the more complex problem of accurately determining the vertical trajectory motion because they do not follow the vertical motion of air. Thus, the two trajectory error percentages given here can be regarded as errors relative to the horizontal trajectory.

The displacement of air particles is also important in numerical semi-Lagrangian schemes for the computation of the departure position of the fluid particle (e.g., Bates et al., 1990; Ritchie and Beaudoin, 1994; Melvin et al., 2010). Many researchers have used rotation matrices to compute the departure position (e.g., Temperton et al., 2001; Staniforth et al., 2010; Wood et al., 2010). As reviewed by Wood et al. (2010), most of these methods are applied in scenarios assuming a shallow atmosphere, although Staniforth et al. (2010) and Wood et al. (2010) achieved a uniform departure point scheme for the assumption of a deep atmosphere.

A general routine to compute the departure point has been designed that uses iterative methods and two- or three-time level schemes (e.g., Temperton et al., 2001; Staniforth et al., 2010; Wood et al., 2010). The procedure is similar to each iterative step and the estimated departure point and mean velocity of displacement tend to be more accurate with each step of the iteration. The accuracy of the iteration is thought to depend on the accuracy of each step within the iterative circle. If only one iterative step is considered, then the circling procedure is simplified to a one-step problem to isolate the impact of other sources—for example, in the MRDM and great circle displacement computing method (GCDM), the initial wind might be considered as a one-step wind during the iteration in which the wind gives more accurate trajectory as iterative step amount increases. The “error” in the MRDM seen by the GCDM is the difference between them, which is not the real error. We evaluated the error inhomogeneity of the MRDM with the GCDM in non-iterative examples and then compared the performance of the MRDM and GCDM using the trajectories of Temperton et al. (2001) as the “real” reference trajectories.

This paper is arranged as follows. In Section 2, the GCDM is deduced through geometric relations and the evaluation elements are introduced. Section 3 presents the non-iterative case, wherein the mean displacement inhomogeneity of the MRDM is accessed through the GCDM with the NCEP/NCAR long-term mean monthly wind fields (Kalnay et al., 1996) and idealized wind fields. In Section 4, the iterative versions of the MRDM and GCDM are compared based on the “real” trajectories from a three-time level trajectory scheme (see Temperton et al., 2001, their appendix). The discussion and conclusions are presented in Section 5.

2 Deduction of the GCDM and evaluation elements 2.1 Deduction of GCDM

Under the approximation of a shallow atmosphere, assuming that an air particle moves above the spherical surface at a one-step speed of V ( = u 2 + v 2 + w 2 ) and that its departure position A is defined as ( λ 1 , φ 1 , h 1 ) in spherical coordinates, then its arrival position B ( λ 2 , φ 2 , h 2 ) after a lapse time Δ t can be solved by the MRDM via the following equation

λ 2 = λ 1 + u Δ t R cos φ 1 , φ 2 = φ 1 + v Δ t R , h 2 = h 1 + w Δ t ,

(1)

where λ, φ, and h are the longitude, latitude, and height above mean surface sea level and u, v, and w are the zonal, meridional, and vertical components of air velocity, respectively, and R, the distance from A to the center of the earth, is the summation of the radius of the earth and h₁.

The GCDM assumes that the air particle moves horizontally along the great circle that is tangential to its horizontal wind velocity vectors at the departure position and that the arc length between the departure and arrival positions is equal to the mean displacement defined by the product of the horizontal wind speed of the departure position and lapse time. The coordinates of A are transformed from spherical to rectangular coordinates with the origin point at the center of the earth, the z-axis along the axis of the earth through the North Pole, and the x-axis along the radius with λ = 0, φ = 0 pointing away from the earth. The position A in the rectangular coordinate ( x A , y A , z A ) is

x A = R cos φ 1 cos λ 1 , y A = R cos φ 1 sin λ 1 , z A = R sin φ 1 .

(2)

Each spherical velocity component causes a displacement along the three directions of the rectangular coordinate. d x | u denotes the displacement in the x direction caused by the horizontal wind component u, and similarly for v and w. The symbol x can be changed to y or z to represent the displacement in the y or z direction. The x-displacement can be written as:

d x = d x | u + d x | v + d x | w .

(3)

We first considered the contribution of the zonal wind to the rectangular displacement. Figure 1 is a sketch map of the latitude circle plane in spherical coordinates. The plane is parallel to the x–y plane in rectangular coordinates. The zonal wind A F → = u i can be decomposed as the vectors A L → and L F → , i.e., A F → = A L → + L F → , where A L → = u cos λ 1 n and L F → = − u sin λ 1 m are the y and x displacements caused by the wind u, respectively. In these equations and hereafter, ( i , j , k ) denotes the unit vector along the zonal, meridional, and vertical directions in spherical coordinates, and ( m , n , l ) denotes the unit vector along the ( x , y , z ) direction in rectangular coordinates. The latitude circle plane is perpendicular to l, thus

d x | u = − u sin λ 1 Δ t ,

(4)

d y | u = u cos λ 1 Δ t ,

(5)

d z | u = 0.

(6)

The contributions of the meridional wind were investigated in a similar manner. Figure 2 is a sketch map of a longitude great circle. The meridional wind A D → = v j is decomposed as A Q → and Q D → , i.e., A D → = A Q → + Q D → , where A Q → is the projection of A D → onto l and Q D → is the projection of A D → onto O M → . The unit vector of O M → is cos λ₁ m + sin λ₁ n. Thus the displacement contributions of wind v are:

d x | v = − v sin φ 1 cos λ 1 Δ t ,

(7)

d y | v = − v sin φ 1 sin λ 1 Δ t ,

(8)

d z | v = v cos φ 1 Δ t .

(9)

The contribution of the vertical velocity w can also be demonstrated in Fig. 2. A T → represents the velocity in the w direction, i.e., w k = A T → , A T → = A P → + P T → , where A P → is the projection of A T → onto l and P T → is the projection of A T → onto O M → . Therefore, the displacements caused by w are:

d x | w = w cos φ 1 cos λ 1 Δ t ,

(10)

d y | w = w cos φ 1 sin λ 1 Δ t ,

(11)

d z | w = w sin φ 1 Δ t .

(12)

According to Eq. (3) and Eqs. (4)–(12), the displacements in the rectangular coordinate are:

d x = ( − u sin λ 1 − v sin φ 1 cos λ 1 + w cos φ 1 cos λ 1 ) Δ t ,

(13)

d y = ( u cos λ 1 − v sin φ 1 sin λ 1 + w cos φ 1 sin λ 1 ) Δ t ,

(14)

d z = ( v cos φ 1 + w sin φ 1 ) Δ t .

(15)

Equations (13)–(15) are consistent with the unit vector forms of Ritchie and Beaudoin (1994). These equations satisfy the relationship for displacement, velocity, and lapse time as follows:

d x 2 + d y 2 + d z 2 = u 2 + v 2 + w 2 Δ t .

(16)

If the air particle is traveling along a straight line, then its arrival position B in the rectangular coordinate can be expressed as:

x = x A + d x , y = y A + d y , z = z A + d z .

(17)

The rectangular coordinate has been temporarily introduced to calculate the arrival position. Position B is then further transformed back to the spherical coordinate using the inverse relation of Eq. (2). If the radius of B is:

R * = x 2 + y 2 + z 2 ,

(18)

then the arrival position of B ( λ 2 , φ 2 , h 2 ) in spherical coordinates can be denoted as:

h 2 = R * − R ,

(19)

λ 2 = arctan ( y x ) ,

(20)

ϕ 2 = arcsin ( z R * ) .

(21)

Equations (19)–(21) are the analytical mathematical solution of the arrival position of the air particle after a straight line displacement. The vertical height of the air particle, as shown in Eqs. (18) and (19), is determined by the vertical speed (w) and the horizontal speeds (u, v). The vertical displacement caused by the horizontal movement is large. It is estimated that a horizontal movement at the equator of 2.5° causes about 6 km of vertical displacement. If this were true, then the air particles at high levels in the troposphere would all penetrate the tropopause. This scenario is never observed and is inconsistent with the horizontal movement characteristics of the atmosphere. In the traditional MRDM, Eq. (1) shows that the horizontal motion of the air particle has the following features: when the air particle is anywhere on the earth apart from the polar points with u = 0 or along the equator with v = 0, it moves along a great circle with its distance of movement equal to the length of the arc instead of the tangent line of the great circle. The displacement along the great circle is more reasonable than the displacement along a tangent. Therefore, the displacement is corrected to the great circle displacement. As a first step, following the treatment of the MRDM in Eq. (1), the vertical height is determined by the vertical speed only:

h 2 = h 1 + w Δ t .

(22)

The horizontal displacement is then corrected from the tangent line length to the arc length and the arrival position in the GCDM is corrected. Figure 3 is a sketch map of the great circle tangential to the horizontal velocity vector. The great circle is OAG; A is the departure point, and B is the arrival point, A B ⊥ A O ; E is on the straight line AB, OE intersects the great circle at point C; the arc A C ⌢ and line AB have equivalent length and OD is the radius of great circle, O D ∥ A B . The arc length is:

A C ⌢ = u 2 + v 2 Δ t .

(23)

The following relation is shown in Fig. 3 (the key triangle is shaded gray), i.e., $ AE/R = \tan \left( {\angle AOC} \right) = $ tan ( A C ⌢ / R ) results in A E = R tan ( u 2 + v 2 Δ t / R ) and α = u 2 + v 2 Δ t / R . The displacement in each direction of the rectangular coordinate [Eqs. (13)–(15)] is corrected to meet the arc length displacement definition, d x E = ( A E / A B ) d x :

d x E = tan α α d x , d y E = tan α α d y , d z E = tan α α d z .

(24)

The coefficient of Eq. (24) has a limitation, lim α → 0 tan α / α = 1. This relation shows that when the horizontal displacement is much less than the radius of the earth, then the displacements under the definitions of the arc length and tangent line length are equivalent.

Thus, the GCDM is deduced. The straight line displacements in the rectangular coordinate are computed with Eqs. (13)–(15). The displacements are then corrected to the arc length definition using Eq. (24). The arrival position in rectangular coordinates is successively computed by adding the corrected displacements to the departure position through Eq. (17). The coordinates of the arrival position are then transformed back to spherical coordinates, the horizontal positions of which are calculated by using Eqs. (18), (20), and (21) and the vertical height is obtained from Eq. (22).

2.2 Evaluation elements of the MRDM

The contribution of vertical winds to horizontal displacement is negligible in the one-step computation of trajectories. Under the approximation of a shallow atmosphere, the movement of an air particle from A to Bmay be separated into two steps: first, horizontal movement from points A to M, then vertical movement from points M to B. The ratio of the horizontal displacement caused by the vertical speed of the particle equals the ratio between the vertical displacement and the distance to the center of the earth with a ceiling value of 1/300. However, if each step of the computation is the same, then it is unclear why the three-dimensional wind field causes such a large change in the trajectory. The reason is that the vertical wind transports the air particle to the position of the new horizontal wind, i.e., the vertical shear of the horizontal wind is an essential factor in the trajectory (Costa and Blanke, 2004). Thus, assuming that the air particle moves horizontally along the earth’s surface, i.e., w = 0, h₁ = 0, then the departure position is A ( λ 1 , φ 1 , 0 ) and the velocity vector is ( u , v , 0 ) . After a lapse time of Δ t , the arrival positions computed through the MRDM and GCDM are B ( λ 2 , φ 2 , 0 ) and C ( λ 3 , φ 3 , 0 ) , respectively. The corresponding arc length of the MRDM and GCDM is:

d m = R cos − 1 [ sin φ 1 sin φ m + cos φ 1 cos φ m cos ( λ 1 − λ m ) ] ,

(25)

where m = 2 and 3, respectively. Correspondingly, the magnitude of the displacement error is defined as:

D m = | d m − u 2 + v 2 Δ t | .

(26)

Using Eq. (2), three positions of the air particle are transformed to rectangular coordinates, i.e., A ( x 1 , y 1 , z 1 ) , B ( x 2 , y 2 , z 2 ) , and C ( x 3 , y 3 , z 3 ) . Therefore, A B → = ( x 2 − x 1 , y 2 − y 1 , z 2 − z 1 ) and A C → = (x₃–x₁, y₃–y₁, z₃–z₁) are the displacement vectors in the rectangular coordinate computed by the MRDM and GCDM. Based on the dot product rule of vectors, i.e., A B → ⋅ A C → = | A B → | ⋅ | A C → | cos θ , the difference in angle between the vectors is defined as:

θ = cos − 1 [ ( x 2 − x 1 ) ( x 3 − x 1 ) + ( y 2 − y 1 ) ( y 3 − y 1 ) + ( z 2 − z 1 ) ( z 3 − z 1 ) ( x 2 − x 1 ) 2 + ( y 2 − y 1 ) 2 + ( z 2 − z 1 ) 2 ( x 3 − x 1 ) 2 + ( y 3 − y 1 ) 2 + ( z 3 − z 1 ) 2 ] .

(27)

In this equation, if | A B → | ⋅ | A C → | = 0 , we set θ = 0 .

The geometrical image of the displacement errors and difference in angle are shown in Fig. 4, where N is the North Pole, x–y are lines of longitude and the concentric circles are circles of latitude. The departure point for the air particle is A. The horizontal velocity vector is A E → = A F → + A D → , i.e., V → = u i + v j and ADEF is a rectangle. After Δ t 1 > 0 , the arrival point is G, with A B ⌢ = u Δ t 1 , A C ⌢ = v Δ t 1 and ACGB is a curvilinear quadrilateral. After Δ t 2 ( Δ t 2 > Δ t 1 ), the arrival point is T, with A H ⌢ = u Δ t 2 , A P ⌢ = v Δ t 2 , and the curvilinear quadrilateral is AHTP A P T H . In the extreme case, after Δ t 3 ( Δ t 3 > Δ t 2 ), the air particle arrives at the North Pole N, with A M ⌢ = u Δ t 3 , A N ⌢ = v Δ t 3 and the curvilinear quadrilateral changes to the curvilinear triangle A N M instead and the displacement vector is A N → . In the MRDM, the air particle moves along the diagonal line of the curvilinear quadrilateral enclosed by two pairs of longitude- and latitude-circle borders. The lengths of the latitude circle borders are usually unequal. When and only when Δ t → 0 , , the curvilinear quadrilateral infinitely tends to a rectangle and the displacement and direction will be equal to those computed by the GCDM. Otherwise, the displacement is not equal to the distance defined by Eq. (23) in most cases, except for the air particle with u = 0 anywhere in the world or with v = 0 along the equator. As shown in Fig. 4, for an air particle at a given departure position, the larger the displacement, then the larger the MRDM displacement error and the angle difference between the MRDM and GCDM displacement vectors. The extreme case is that the arrival position is the North Pole and the curvilinear quadrilateral then becomes a curvilinear triangle.

To give a notional result, three parameters were set as constant in the following computations: w = 0 m s^–1, h₁ = 0 m, and R = 6370 km at each grid point.

3 Evaluation of MRDM: Non-iterative case

In Section 2, Eq. (16) shows that when the air particle moves along the tangent line of a great circle, its physical progress satisfies the velocity–time–displacement relation. Equation (24) indicates that when the ratio between the displacement of the air particle and the radius of the earth is small, then the GCDM generally satisfies the relation. Because the mean horizontal displacement between the departure and arrival positions is assumed to be equal to the real displacement defined by the velocity multiplied by the lapse time in the deduction of the GCDM, the displacement computed through the GCDM is accurate in theory and its displacement error is only a computational error. However, the displacement computed through the MRDM has both analytical and computational errors (Fig. 4). In the following discussion, the MRDM is evaluated in terms of the distance error [Eq. (26)] and the angle difference [Eq. (27)] using long-term monthly mean wind field and ideal wind field datasets.

3.1 Long-term mean monthly wind field

The NCEP/NCAR long-term mean monthly wind field is derived from data during 1981–2010 (Kalnay et al., 1996). The dataset includes all 12 months of the year and the horizontal resolution is 2.5° × 2.5° with a global coverage. It contains 17 standard pressure levels in the vertical direction. January is selected to represent the boreal winter and July to represent the boreal summer. The 850-, 500-, and 200-hPa levels are selected to represent the low, middle, and top levels of the troposphere, respectively. The polar region is excluded and the region within 80°S–80°N is selected. The general circulation of the atmosphere in the boreal winter shows that there are strong westerly winds over East Asia at 500 hPa (Fig. 5), with a maximum wind speed > 36 m s ^–1 (about 130 km h^–1); the wind speed at 200 hPa is higher and there are more vortexes at 850 hPa.

In the region of interest, each grid point is treated as the departure position of the air particle, which moves at the wind speed at this location. For the one-step displacement is small compared with the midlatitude wavelength on the climate mean chart (Fig. 5) and the detailed structure is only slightly distorted under scale consideration. The mean displacements and their errors and the angle difference between the MRDM and GCDM are calculated and set at the departure grid. Thus, the distribution pattern of these variables can be plotted. The displacement lapse times are set to 1, 3, 6, and 12 h. The mean distance errors of the MRDM, the mean distance errors of the GCDM, and the angle differences between the displacements computed by the MRDM and GCDM for the 500-hPa wind in January are shown in Figs. 6–8. To give general information, the data are smoothed once by nine points before being plotted. When the lapse time is 1 h, the displacement errors computed by the MRDM and GCDM have similar distribution patterns (Figs. 6a, 7a). Because the displacement error of the GCDM (Fig. 7a) is only a computational error and Figs. 6a and 7a are similar, it might be conjectured that the 1-h displacement errors of the MRDM are dominated by computational errors. The displacement error of the GCDM has a decreasing trend with increasing lapse times (Figs. 7a–d), whereas that of the MRDM has an increasing trend (Figs. 6a–d). The GCDM displacement error has a decreasing trend because the GCDM displacement is corrected by Eq. (24), in which the correction coefficient is related to the magnitude of displacement. The mathematically accurate form is tan α α . For the computation procedure, both the numerator and denominator inevitably have computational errors of tan α + ε 1 α + ε 2 , where ε 1 and ε 2 are the computational errors. When the displacement is small, although mathematically lim α → 0 tan α / α = 1 , both the numerator and denominator are small and tend to zero, which means that effect of the computational errors on the coefficient is large. When the displacement is large, the computational errors are less important in magnitude than the numerator and denominator, which makes the coefficient more accurate and thus the computational error of the magnitude of displacement decreases with the lapse time. The region with a large angle difference is located in the middle and high latitudes and its magnitude increases as the lapse time increases (Figs. 8a–d).

To include the whole layer and seasonal information, the confinement of the vertical wind field is extended to include the low, middle, and upper troposphere, and the temporal coverage is expanded to cover both the boreal winter and summer. The mean and maximum regional displacement errors and the maximum angle differences are shown in Table 1. The maximum displacement errors of both methods are small with the 1-h lapse time and the GCDM is slightly better than the MRDM for the mean displacement error. The displacement error of the MRDM increases with the lapse time. To determine the acceptable upper limit of the displacement error, Scheele and Siegmund (2001) applied the Royal Netherlands Meteorological Institute TRAJKS trajectory model (Scheele et al., 1996) to compute the three-dimensional displacement during the 12-h interval of the ECMWF ensemble predictions using an iterative scheme. A horizontal distance between neighboring iterative steps of < 300 m is one of the preconditions that stops the iteration. Therefore, 300 m was regarded as the maximum acceptable displacement error in this study. Table 1 shows that the regional mean displacement error of the MRDM is > 300 m for lapse times at 200, 500, and 850 hPa of 3/3, 6/6, and 12 h/6 h in the boreal winter (January)/summer (July), respectively. This result indicates that the MRDM should be applied with caution when the lapse time is > 3 h. For the MRDM, the 12-h regional maximum displacement errors at 200 hPa are about 70 and 50 km in January and July, respectively, with corresponding regional mean errors of about 8 and 5 km. The maximum angle difference is up to about 10° in July with a lapse time of 12 h. However, the regional mean (maximum) displacement error of the GCDM is < 100 m (< 3 km) for all time lapses.

3.2 Ideal wind field

To test the dependence of the trajectory error on the latitude and the definition range of the wind direction, an ideal wind field is designed for the one-step trajectory problem. There are two centers of large MRDM displacement errors in front of and behind the large North American trough (Figs. 6b–d). As discussed in Section 2.2, the displacement error of the MRDM is determined not only by the magnitude, but also by the wind direction at the departure position. Moreover, the latitude bound is limited to 80°S–80°N. The wind speed in this region is consistently 10 m s^–1 and the defined range of the wind direction is 0°–360°. Thus, the ideal wind field is constructed in the space of latitude and wind direction. The 0°–80°N and 0°–360° wind direction region is selected for the symmetrical pattern of the Northern and Southern Hemispheres (Fig. 9). The lapse times are still 1, 3, 6, and 12 h. The displacement errors of the MRDM and GCDM are investigated, together with the angle difference. As the displacement is the product of the wind speed and lapse time, the idealized error distribution figures are similar to a look-up table using the relation between the displacement vector and the latitude.

The displacement error of the GCDM also decreases with the lapse time. The displacement error of the MRDM (Fig. 10) and the angle difference (Fig. 11) increase with the lapse time (Figs. 10 and 11 are smoothed as in Figs. 6–8). The displacement error of the MRDM has a four-wave distribution pattern in the defined range of wind direction, with maximum values at 45°, 135°, 225°, and 315° and sharp wave peaks (Fig. 10). The angle difference has a two-wave feature, the pure meridional displacements have zero angle differences, and there are two flat peaks between them in two sections of wind direction (0°–180° and 180°–360°) (Fig. 11). Table 2 gives the variation of the displacement errors and angle differences with the lapse time. The characteristics are similar to those in Table 1. As the ideal wind field is used, the error is more orderly distributed. For example, the maximum angle differences increase with the lapse time. For the given region (here 0°–80°N), the regional maximum angle difference is larger for larger displacements. The differences in angle increase nonlinearly with latitude (Fig. 11).

3.3 Applicability of the GCDM at the North Pole

The applicability of GCDM at the North Pole is investigated. The data used are the long-term monthly mean NCEP/NCAR wind field at the North Pole. The wind field data consist of the wind components at 144 grid points overlapping at the North Pole. The wind speed is almost constant and the wind components vary with longitude as a sine or cosine function (Fig. 12a). The distribution of the displacement error in the space of longitude and lapse time is shown in Fig. 12b. The distribution pattern of the displacement error and its statistical information is similar to the examples discussed in the preceding sections. A close study of the data shows that the mean wind speeds at 850 hPa in January are not strictly equivalent at the North Pole, with a difference of about 0.03 m s^–1 between the maximum and minimum values. The difference will result in 108, 324, 648, and 1296 m of displacement with lapse times of 1, 3, 6, and 12 h, respectively; the corresponding maximum displacement errors of the GCDM are 127.0, 64.3, 37.9, and 17.8 m, respectively. They are of the same magnitude with a time lapse of 1 h, but the errors in the GCDM are much smaller when the time interval is > 3 h. This result suggests that the GCDM is applicable in the polar region.

4 Evaluation of MRDM: The iterative case

In the preceding section, the trajectory is only dependent on the “initial” wind and the GCDM acts as the “real” trajectory when assessing the angle difference with the MRDM. In fact, there is no exact arrival position of the air particle when only the initial wind vector of the departure position is given. For example, a Northern Hemisphere air particle with a wind speed of u = 10 m s^–1 and v = 0.0 m s^–1 will move to the east with some northward deviation if it is to the south of a small low pressure center, but some southward shift will occur if it is to the north of a small high. That is to say, all the trajectories are dependent on the wind field. To solve this dependence, iterative methods are introduced to compute the displacement more accurately: not only the initial wind vectors, but also the iterative wind vectors at the other (departure) position are used. Because of the different physical conditions, a one-step trajectory cannot be evaluated by an iterative method. This section compares the MRDM and GCDM using iteration.

4.1 Introduction of test designation

Iterative procedure. The iterative procedure can significantly improve the computational accuracy of the trajectory (Stohl, 1998). The general form of iteration is:

X 1 ( t 1 ) ≈ X ( t 0 ) + ( Δ t ) X ˙ ( t 0 ) , X 2 ( t 1 ) ≈ X ( t 0 ) + 1 2 ( Δ t ) [ X ˙ ( t 0 ) + X ˙ 1 ( t 1 ) ] , ⋮ X i ( t 1 ) ≈ X ( t 0 ) + 1 2 ( Δ t ) [ X ˙ ( t 0 ) + X ˙ i − 1 ( t 1 ) ] ,

(28)

where X is the position vector, X ˙ is the wind vector, t is time, the subscripts 0 and 1 denote two time levels and Δ t is the lapse time between the two time levels. X ˙ i ( t 1 ) is taken at position X i ( t 1 ) and the superscript i denotes the ith iterative step of the corresponding vector.

Avoidance of interpolation errors. In general, the interpolation errors include both spatial and temporal interpolation errors (Bowman et al., 2013) and both are major sources of error in the computation of trajectories. The temporal variation of the analytical wind field of cases 1–3 of Nair and Lauritzen (2010) are introduced to avoid interpolation errors:

Case 1:

{ u ( λ , φ , t ) = k sin 2 ( λ / 2 ) sin ( 2 φ ) cos ( π t / T ) , v ( λ , φ , t ) = k 2 sin ( λ ) sin ( φ ) cos ( π t / T ) .

(29)

Case 2:

{ u ( λ , φ , t ) = k sin 2 ( λ ) sin ( 2 φ ) cos ( π t / T ) , v ( λ , φ , t ) = k sin ( 2 λ ) cos ( φ ) cos ( π t / T ) .

(30)

And case 3:

{ u ( λ , φ , t ) = − k sin 2 ( λ / 2 ) sin ( 2 φ ) cos 2 ( φ ) cos ( π t / T ) , v ( λ , φ , t ) = k 2 sin ( λ ) cos 3 ( φ ) cos ( π t / T ) .

(31)

The analytical wind fields of Nair and Lauritzen (2010) are theoretical and the parameters are re-appointed according to the atmospheric conditions. The flow parameter k is roughly the scale of wind component globally and is taken as 30 m s^–1. Parameter T is the timescale of wind revision and is taken as 72 h (3 days).

Reference method. The method of Temperton et al. (2001) is used as the reference method. The three-time level trajectory scheme of Temperton et al. (2001) had been used as a standard method in the ECWMF operational model. Their algorithm is so accurate and robust that its modifications are still used by world-leading weather forecast centers. Thus, the trajectories computed by using the method of Temperton et al. (2001) are treated as the “real” trajectories.

Test designation. Because the MRDM is not suitable for high latitudes, the test region for the departure position is set as the spherical band between the ±80° latitudes. The departure positions are distributed evenly with a resolution of 1° × 1°. A one-step trajectory is tested from t = 0 with lapse times of 1, 3, 6, and 12 h, respectively. To avoid the computational singularity of Eq. (24), the correction coefficient is set to 1 when the displacement is <1 mm (corresponding angle α approximately 10^–10). The assessment factors are same as the preceding errors in magnitude and angle, but the “real” trajectory is that of Temperton et al. (2001).

4.2 MRDM and GCDM error distributions

The magnitude error and angle error of the MRDM and GCDM are computed with respect to the results of Temperton et al. (2001). Cases 1 and 2 are examples of non-divergent deformation flow and case 1 presents a large gyre centered at 180° of the equator view in a latitude–longitude map projection [Eq. (29) and Fig. 13a]. Case 2 presents two gyres centered at 90° and 270° from the equator, respectively [Eq. (30) and Fig. 14a]. Case 3 is an example of divergent deformation flow, a saddle-form system in meteorology with a singular point at 180° from the equator [Eq. (31) and Fig. 15a].

In each case of the ideal wind field, the distribution patterns of both errors are unchanged with the lapse time, although their magnitudes increase. Considering the 6-h lapse time for wind field in case 1, two centers of angle error of the MRDM (Fig. 13b) are located at the most northern and southern boundaries where the winds are at their maximum along the boundary latitudes. Because of computational errors, there may be singular values of the angle error at positions where the winds tend to zero. The singular angle error is characterized by one or more extremely large angle errors (≤ 180°) surrounded by angle error regions with smaller values; the angle errors change abruptly spatially. Figure 13b shows a singular angle error of 180° at 0° N, 180°. The singular values of angle difference are related to computational errors and are not discussed further. The magnitude error of the MRDM has two close maximum centers at 60° N and 60° S due to the shape of the sphere and the two centers shift in the upwind direction. The GCDM errors have a similar distribution pattern to the MRDM (Fig. 13b), but only half the magnitude (Fig. 13c). The error distribution of case 2 is very similar to that of case 1, except that two gyres are introduced (Figs. 14b, c). The magnitudes of the change in wind direction along the southern and northern boundaries are doubled (one gyre in Fig. 13a to two gyres in Fig. 14a) and the upwind shift in the centers of the angle and magnitude errors are more obvious than those in case 1. The angle error distribution pattern in both cases 1 and 2 is similar to the ideal test in Fig. 11; the maximum values are located poleward of the gyre centers and the difference in the distance of the northeast/southwest wind shown in Fig. 10 is not obvious. For case 3, the wind speed is small at high latitudes (Fig. 15a) and the maximum angle error and magnitude error are much less than those of cases 1 and 2. The error centers are located between 30° and 60° in the Southern/Northern Hemisphere and the error centers slightly shift in the upwind direction (Figs. 15b, c). All cases 1–3 indicate that the errors of both methods are related to latitudinal location, wind speed, and wind direction.

The general information about the MRDM and GCDM errors in Table 3 shows that the magnitude error and angle error are related to the wind pattern. In case 2, both methods have the largest error, and in case 3, both methods have the least error. Unlike the prior non-iterative cases, the magnitude error of the GCDM, together with all the other errors, increases with the lapse time. For a lapse time of 1 h, the global mean distance errors of GCDM are 282.4, 380.2, and 56.8 m for cases 1, 2, and 3. The last mean distance error is a little larger than the preceding 300 m critical value. In general, the GCDM result is close to the trajectory of Temperton et al. (2001) at a timescale of 1 h. If a shorter time lapse is given, the trajectory of the GCDM is even closer to the trajectory of Temperton et al. (2001) than the MRDM. The magnitude and angle errors of the MRDM are about twice as large as the corresponding errors of the GCDM.

5 Conclusions and discussion

We investigated the one-step computation problem of trajectories on a sphere with and without iteration. For the non-iteration cases, the inhomogeneity error of the computation of displacement using the MRDM was shown to depend on the wind speed and direction and the latitude of the departure position. The error increases with wind speed and absolute latitude. The magnitude of the error and the direction difference has four- and two-wave patterns in the definition range of the wind direction. In terms of the general circulation of the atmosphere, the wind field has a pattern in which westerly wind bands alternate with easterly wind bands, especially at high levels. When the displacement is large, the displacement computed by the GCDM has a component toward the equator in pure westerly/easterly wind bands, but the MRDM ensures that the air particle moves along the circle of latitude in this situation. Thus, the general circulation of the atmosphere favors the MRDM through its circulation pattern. When the displacement is close to zero, the GCDM and MRDM give the same results under this extreme condition (Figs. 4, 6a, 7a). At the lower levels of the atmosphere, winds may occur with arbitrary directions, and the GCDM gives better results.

The clearest difference between the MRDM and GCDM is the dependence of the displacement error of the MRDM on the wind direction and the latitude location of the air particle. Under some special situations, the physical concepts of the GCDM and MRDM are similar and the displacements computed by the two methods are equivalent. The comparison between the GCDM and MRDM is thus a comparison between the physical figure of the MRDM under special and common situations. The comparison explains why the displacement error of the MRDM depends on the wind direction. The GCDM improves the MRDM by extending its physical figure from examples with a few special wind directions and locations to all wind directions and the global domain. The GCDM has a consistent physical figure and is a suitable reference method for assessing the inhomogeneity of the MRDM.

As the displacement computed by the MRDM is different in both magnitude error and angle differences from the GCDM, it is natural to ask under what conditions the MRDM can be used to compute the displacement. The MRDM can be used when the displacement is small and the region of displacement is far from the pole. The Courant number is a criterion for the legality of displacement computation (Kuo et al., 1985; Innocentini, 1999). The precondition is a Courant number that is much smaller than unity. In spherical coordinates, the meridional grid space decreases rapidly away the equator and the precondition is hard to satisfy at high latitudes. Thus, the MRDM has a polar singularity problem and the GCDM is an alternative method of computing the mean displacement in polar regions.

For the iterative cases, temporal variations in the analytical ideal wind fields are introduced to avoid interpolation errors. The magnitude errors and angle errors of the MRDM and GCDM with respect to the trajectories of Temperton et al. (2001) are compared. The results show that the GCDM is closer to the trajectories of Temperton et al. (2001) than the MRDM, their corresponding errors have a similar distribution pattern and the mean GCDM errors are half those of the MRDM. The lapse time is usually set to a critical value at which an air particle cannot move across more than one grid space. The chosen stream factor is 30 m s^–1 and the 1-h displacement is about 108 km, which is on the scale of the ideal wind grid spaces near the equator. The error distribution with a lapse time of 6 h is given in the iterative section because the other three-time lapse examples have similar error distribution patterns.

The results of this study indicate that the curvilinear coordinate itself might be a source of error in trajectory models when the displacement is computed through a velocity resultant method (e.g., the MRDM). The error is not negligible in the region where at least one of the coordinates has a large curvature (e.g., the polar region of spherical coordinates) and/or the displacement is large. The GCDM is an alternative displacement computation method with less error than the MRDM and can be used to solve the trajectory at high latitudes.

Acknowledgments. The authors thank the anonymous reviewers for their thoughtful critiques and suggestions, which substantially improved the paper. The NCEP/NCAR reanalysis data are provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, at www.esrl.noaa.gov/psd/.