An Application of the Immersed Boundary Method for Recovering the ...

MAY 2012

LIOU ET AL.

1603

An Application of the Immersed Boundary Method for Recovering the Three-Dimensional Wind Fields over Complex Terrain Using Multiple-Doppler Radar Data YU-CHIENG LIOU AND SHAO-FAN CHANG Department of Atmospheric Sciences, National Central University, Jhongli City, Taiwan

JUANZHEN SUN National Center for Atmospheric Research, Boulder, Colorado (Manuscript received 27 June 2011, in final form 17 October 2011) ABSTRACT This study develops an extension of a variational-based multiple-Doppler radar synthesis method to construct the three-dimensional wind field over complex topography. The immersed boundary method (IBM) is implemented to take into account the influence imposed by a nonflat surface. The IBM has the merit of providing realistic topographic forcing without the need to change the Cartesian grid configuration into a terrain-following coordinate system. Both Dirichlet and Neumann boundary conditions for the wind fields can be incorporated. The wind fields above the terrain are obtained by variationally adjusting the solutions to satisfy a series of weak constraints, which include the multiple-radar radial velocity observations, anelastic continuity equation, vertical vorticity equation, background wind, and spatial smoothness terms. Experiments using model-simulated data reveal that the flow structures over complex orography can be successfully retrieved using radial velocity measurements from multiple Doppler radars. The primary advantages of the original synthesis method are still maintained, that is, the winds along and near the radar baseline are well retrieved, and the resulting three-dimensional flow fields can be used directly for vorticity budget diagnosis. If compared with the traditional wind synthesis algorithm, this method is able to merge data from different sources, and utilize data from any number of radars. This provides more flexibility in designing various scanning strategies, so that the atmosphere may be probed more efficiently using a multiple-radar network. This method is also tested using the radar data collected during the Southwest Monsoon Experiment (SoWMEX), which was conducted in Taiwan from May to June 2008 with reasonable results being obtained.

1. Introduction Although a Doppler radar is capable of providing wind and reflectivity information inside a precipitation system with high spatial and temporal resolution, the so-called Doppler radial velocity, however, is merely the projected component of the complete wind field along the radar beam. Several methods of using multiple-radar measurements to reconstruct a full three-dimensional wind field have been developed since the late 1960s (Armijo 1969; Miller and Strauch 1974; Ray et al. 1975; Doviak

Corresponding author address: Dr. Yu-Chieng Liou, Department of Atmospheric Sciences, National Central University, 320, Jhongli City, Taiwan. E-mail: [email protected] DOI: 10.1175/MWR-D-11-00151.1 Ó 2012 American Meteorological Society

et al. 1976; Brandes 1977; Ray et al. 1978). Given a prescribed boundary condition, these methods basically solve a unique solution of the wind field from a set of closed equations, which usually consists of the geometric relations between the Cartesian wind (u, y, w) and the radial wind components detected by each radar, the equation of mass continuity, and an estimation of the particle terminal fall speed. The variational technique has also proven to be an effective means for adjusting errors introduced by inadequate specification of the top or bottom boundary conditions for the vertical velocity, or for enforcing the satisfaction of the continuity equation (O’Brien 1970; Ray et al. 1978, 1980; Chong and Testud 1983; Ziegler et al. 1983). The variational approach has started to play an even more important role in various studies such as those by Scialom and Lemaitre

1604

MONTHLY WEATHER REVIEW

(1990), Protat and Zawadzki (1999), Shapiro and Mewes (1999), Gao et al. (1999), and Gao et al. (2004). In these methods, a cost function is usually formulated first. Within this cost function, the major ingredients for performing multiple-Doppler wind synthesis, such as the radial velocity observations, the mass conservation equation, or some forms of low-pass filters, are employed as either weak or strong constraints. By minimizing this cost function, one obtains an optimally adjusted wind field. However, all the multiple-Doppler synthesis methods mentioned above are conducted over a flat surface. This limits their practical application in regions with topography. Georgis et al. (2000) developed a new technique for using airborne or ground-based Doppler radars to observe precipitation systems over areas with complex orography. Chong and Cosma (2000) improved the formulation of a method called the multiple-Doppler synthesis and continuity adjustment technique (MUSCAT; Bousquet and Chong 1998), which can be applied to resolve the flow fields over both flat and complex terrain. Chong et al. (2000) combined the methods of Georgis et al. (2000) and Chong and Cosma (2000). They reported results from a real-time and automated multiple-Doppler synthesis algorithm utilized data acquired during the Mesoscale Alpine Programme. Bousquet et al. (2007, 2008a,b) described the successful implementation of the MUSCAT method to establish an operational multipleDoppler wind retrieval algorithm with the French radar network. Drechsel et al. (2009) showed the threedimensional wind field retrieved by MUSCAT with dual-Doppler lidar data observed during the 2006 Terraininduced Rotor Experiment (TREX). Recently, the first author of this current work designed a variational-based multiple-Doppler three-dimensional wind synthesis method (Liou and Chang 2009, hereafter LC09). A simplified vertical vorticity equation, due to its advantages in providing additional constraints for wind analysis (Protat and Zawadzki 2000; Mewes and Shapiro 2002; Shapiro et al. 2009; Alaoui et al. 2010), is also incorporated in LC09. Thus, the resulting wind fields produced by LC09 can be readily applied to conduct vorticity budget analysis or derive the thermodynamic structure of the weather system. In addition, similar to Chong and Bousquet (2001), the method of LC09 is also able to reconstruct the wind field along and near the radar baseline. However, in Chong and Bousquet (2001), the successful recovery of the baseline winds by their method is attributed to the smoothing term, which involves a low-pass filter with second-order derivatives of the winds. On the other hand, the experiments in LC09 show that the continuity equation and the vorticity equation play even important roles in retrieving the wind field along the baseline.

VOLUME 140

This study presents an extension of LC09, for the purpose of developing a multiple-Doppler radar synthesis method to retrieve the three-dimensional flow fields over complex terrain, while still retaining the original advantages reported in LC09. To take into account the influence imposed by a nonflat surface, the immersed boundary method (IBM) is adopted. The IBM can provide realistic topographic forcing, without the necessity to change the Cartesian grid configuration into a terrainfollowing coordinate system, therefore is particularly suitable for our purpose. In the next section, the formulation of the multipleradar synthesis method is introduced. In sections 3–4, a series of (observing system simulation experiment) OSSE-type tests are conducted to explore the performance of the proposed method. In section 5, this method is applied to analyze a heavy rainfall event using radar data collected during the Southwest Monsoon Experiment (SoWMEX), which was conducted from May to June 2008 in Taiwan. In section 6, some conclusions and plans for future work are discussed.

2. Methodology a. Variational-based wind synthesis The cost function formulated in this study basically follows the LC09 design, with the exception of some modifications for selecting the top and bottom boundary conditions. The first constraint represents the geometric relation relating the retrieved Cartesian wind components and the radial velocity directly measured by the individual radar, expressed as (Vr )i,t 5

(y 2 Piy ) (x 2 Pix ) ut 2 yt ri ri 2

ri 5

(z 2 Piz ) (wt 1 WT ,t ), ri

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (x 2 Pix )2 1 (y 2 Piy )2 1 (z 2 Piz )2 ,

(1a)

(1b)

where the subscripts i (from 1 to N) and t (from 1 to 2) are the indices for each radar and each time level, respectively. In (1), (Vr)i,t is the radial winds observed by the ith radar at time t; (ut, y t, wt) denote the three-dimensional wind at the location (x, y, z) at time t; (Pix , Piy , Piz ) depict the coordinates of the ith radar; ri stands for the distance from each grid point to the ith radar; and WT,t stands for the terminal falling speed of the raindrop, which can be estimated using the radar reflectivity and an empirical equation (e.g., see Shapiro et al. 1995). A coefficient is

MAY 2012

assigned to this constraint, and is given a value of either 1 or 0 to mark whether a specific grid point is covered by the ith radar or not. The second and the third constraints are the anelastic continuity equation and a simplified vertical vorticity equation, with the mixing and baroclinic terms neglected, as follows: ›(r0 u) ›(r0 y) ›(r0 w) 1 1 5 0, ›x ›y ›z

(2)

›j ›j ›j ›j ›u ›y 2 (j 1 f ) 52 u 1y 1w 1 ›t ›x ›y ›z ›x ›y ›w ›y ›w ›u 2 , (3a) 2 ›x ›z ›y ›z ›y ›u 2 , ›x ›y

J5

(3b)

å å [w(i, j, k 5 NZ)]2, i

(4b)

j

where i, j, k are the grid indices along the x, y, and z directions, respectively; and NZ is the index for the highest model level. The treatment of the bottom boundary condition for w over terrain will be discussed in the next section. For the horizontal wind fields (u, y), we consider a scenario wherein the radar data availability would decrease after reaching a certain altitude. If at a given layer the amount of data falls below a certain threshold (e.g., ,50%), then information from the background fields would be utilized. Here two options are provided. If the background field is generated by a sounding, one can let the horizontally averaged u and y approach a single value. That is

j5

1605

LIOU ET AL.

NZ

J5

å f[u(k) 2 usounding(k)]2 1 [y(k) 2 ysounding(k)]2g,

k5N1

In (2), r0 is the air density, which varies only with height. In (3a), j is the vertical vorticity, while f represents the Coriolis parameter. The first, second, and third terms on the right-hand side of (3a) denote the advection, divergence, and tilting terms of the vorticity equation, respectively. Finally, the overbar means a temporal average over the two time levels. The background term and a Laplacian smoothing filter, which involves second-order derivatives of the wind fields have the same forms as those in LC09, thus they are not shown here. The background term is implemented to fill in the data-void regions, and serves as the initial guess for the retrieved variables. It can be prepared using in situ observations (e.g., radiosonde), outputs from a mesoscale numerical model, reanalysis data, or a combination of these. The method for determination of the top and bottom boundary conditions for the winds, however, is different from that of LC09. Instead of enforcing the horizontally averaged w approaching zero, as suggested in LC09, other alternatives are provided. In this study, in the regions with flat terrain, we extend the first level of the domain to the ground or a small altitude, so that the vertical velocity at this level can be reasonably set to zero. A constraint for this purpose is formulated in (4a), and applied only to areas with flat surface in the analysis domain. A similar top boundary condition expressed by (4b) can also be specified, if the upper limit of the analysis domain reaches a sufficient altitude (e.g., close to tropopause): J5

å å [w(i, j, k 5 1)]2 , i

j

(4a)

(5) where N1 denotes the starting layer to apply this constraint and the overbar stands for a horizontal average. On the other hand, if the background field carries a threedimensional spatial structure, then the constraint can be rewritten as NZ

J5

å å å [(u 2 ubkgd )2 1 (y 2 ybkgd )2 ].

k5N1

i

(6)

j

Constraint (6) simply states that (u, y), from layers N1 to NZ, are adjusted to approach their background values (ubkgd, y bkgd) at each grid point. It should be pointed out that all the constraints are incorporated into a single cost function in the least squares forms. Applying the three-dimensional variational data assimilation (3DVAR) analysis, the expressions for the gradients of the cost function with respect to six unknown control variables (i.e., three wind components at two time levels) can be derived. Starting from an initial guess (usually zero) for these variables, and using the information of the gradients to determine the descent direction, a new set of us and ys, which can reduce the cost function to a smaller magnitude is solved. This procedure is repeated until the cost function converges to a minimum, which implies that in the least squares sense all the constraints are satisfied simultaneously. As discussed in LC09, the weighting coefficients for each constraint are determined through a trial-and-error process. A series of tests are conducted to ensure that the retrieved fields are not sensitive to this set of weighting coefficients.

1606


FIG. 1. A two-dimensional schematic example of the IBM configuration. The thick solid line indicates the immersed boundary, the solid circles stand for the so-called ghost-cell points (G), while the open circles represent their image points (I). The ghost-cell point and its corresponding image point have the same distance to the immersed boundary. The lines connecting G and I intersect the boundary perpendicularly at point B. The crosses stand for grid points inside the terrain, and the open squares mark the points in the flow region.

b. The immersed boundary method The IBM was first introduced by Peskin (1972) for the study of cardiac mechanics. The IBM is capable of simulating fluid patterns over a complex geometry, but the entire computation can be carried out on a Cartesian grid, without the need to use a terrain-following coordinate. This is accomplished by modifying the nodes near the immersed boundary to simulate the force exerted by the surface on the flow. Numerous ways of performing the modifications and refinements have been proposed. The manuscript by Mittal and Iaccarino (2005) provided an overview of the variants of the IBM. Recently, Lundquist et al. (2010) implemented the IBM into the Weather Research and Forecasting Model (WRF). The coupled IBM–WRF is able to simulate flows over highly irregular surface, such as an urban skyline with nearvertical slopes. This would be, however, difficult for the WRF to simulate using its original s coordinates. Since the algorithm reported in LC09 can only be applied over

VOLUME 140

a flat surface, and the purpose of this study is to extend the LC09 method for synthesizing a three-dimensional wind field on complex terrain, thus the IBM method is chosen and implemented, because of its simplicity and the advantages discussed above. The IBM algorithm used in this study is similar to Tseng and Ferziger (2003), but with a different treatment of the near-boundary nodes. A two-dimensional example is shown in Fig. 1 to illustrate how the grid points are grouped into two categories: the flow region and the ghost cell (Tseng and Ferziger 2003). The terrain elevation is specified as a function of x and y only, and has the same resolution as the analysis domain. At each column, the first point below the immersed boundary is identified as the ghost cell. In Fig. 1, G stands for the ghost cell, and B is the location where the boundary condition is to be satisfied. To prevent numerical instability, the values at nearby points in the flow region are first interpolated to the image point of the ghost cell crossing the boundary, labeled I in Fig. 1. Each ghost-cell point and its corresponding image point have the same distance to the immersed boundary. The straight line connecting G and I intersects the boundary perpendicularly at point B. In this study, the wind components are assumed to change linearly in a threedimensional space, and the interpolation is performed using multiple linear regression fitting. For Dirichlet-type boundary condition, an overdetermined system can be formulated as in the following: 8 u0 5 uI 1 ux X0 1 uy Y0 1 uz Z0 > > > > > > u1 5 uI 1 ux X1 1 uy Y1 1 uz Z1 > > > < u2 5 uI 1 ux X2 1 uy Y2 1 uz Z2 (7) , > > > . > .. > > > > > :u 5 u 1 u X 1 u Y 1 u Z N

I

x

N

y

N

z

N

where u0 stands for the boundary condition of u, which needs to be satisfied at point B over the immersed surface (see Fig. 1); ui with i 5 1 ; N are the values of u at the nearest N points in the flow region surrounding the image point uI; and (ux, uy, uz) denote the spatial derivatives along the x, y, and z directions, respectively. The distances from the image point to point B and the other N points are represented by (X0, Y0, Z0) and (Xi, Yi, Zi, i 5 1 ; N), respectively. Applying the multiple linear regression method, one obtains a set of equations in matrix form: AX 5 B, where

(8)

MAY 2012

1607

LIOU ET AL.

2

N

å Xi

6 (N 1 1) 6 6 6 6 N 6 6 6 X 6 i50 i 6 A56 6 N 6 6 Y 6 6 i50 i 6 6 6 N 4 Zi

The unknown parameters in matrix X can be solved by

å Zi

7 7 i50 7 7 7 N 7 7 Xi Zi 7 7 i50 7 7 7 N 7 Yi Zi 7 7 7 i50 7 7 N 7 5 Zi Zi

i50

N

3

N

å Yi

i50

å

N

X 5 A21 B.

å Xi Xi å Xi Yi å

å

i50

i50

N

N

i50

i50

N

N

8 ›u > > ^ ^ ^ > > ›n 5 0 uI 1 ux nx 1 uy ny 1 uz nz > > > > > u1 5 uI 1 ux X1 1 uy Y1 1 uz Z1 > > > < u2 5 uI 1 ux X2 1 uy Y2 1 uz Z2 , > > > > > .. > > . > > > > > :u 5 u 1 u X 1 u Y 1 u Z

å Xi Zi å Yi Zi å

i50

2

i50

i50

3

N

å

6 u 7 7 6 6 i50 i 7 7 6 7 6 N 2 3 uI 7 6 6 ui Xi 7 6u 7 7 6 6 x7 7 6 i50 6 7. 6 7 X 5 6 7 B 56 7 N u 4 y5 7 6 7 6 6 ui Yi 7 uz 7 6 i50 7 6 7 6 N 7 6 5 4 ui Zi

å

N

å

2

N

2

N

å Xi

N

N

z

N

n^x n^x 1

å Xi Xi

i51

n^x n^y 1

n^z n^x 1

å Yi Xi

n^y n^y 1

N

å Yi Yi

i51 N

å Zi Xi

i51

å X i Yi

i51 N

i51

n^z n^y 1

å Zi Yi

i51

7 7 7 7 7 N 7 Xi Zi 7 n^x n^z 1 7 7 i51 7, 7 N 7 7 Yi Zi 7 n^y n^z 1 7 i51 7 7 N 7 5 Zi Zi n^z n^z 1 i51

N

N

å

å Zi

i51

N

n^y n^x 1

å

å

7 6 ui 7 6 i51 7 6 7 6 N 7 6 7 6 ›u 7 6 ^ n 1 u X 6 ›n i i 7 i51 7 6 B 56 7. N 7 6 7 6 ›u 6 n^ 1 ui Yi 7 7 6 ›n y i51 7 6 7 6 N 7 6 5 4 ›u ui Zi n^z 1 ›n i51

å

i51

uG 5 uI 2

å

(13)

å å

The unknown matrix X can now be solved by (10). The value of u at the ghost cell (uG) is computed by

(12)

å

3

N

y

3

N

å Yi

i51

å

i51

x

where (^ nx , n^y , n^z ) are the components of the unit vector perpendicular to the immersed boundary and ›u/›n denotes the normal derivative of u crossing the boundary. The expression of (8) can still be used, but with the matrices A and B replaced by the following:

i50

å

I

(11)

(9)

å

6 N 6 6 6 6 N 6 6 X 6 6 i51 i A56 6 N 6 6 6 Y 6 i51 i 6 6 6 N 4 Zi

(10)

Then, the value of u at the ghost cell (uG) is obtained by uG 5 2u0 2 uI. Similarly, for a Neumann-type boundary condition, the overdetermined system is written as

å Yi Xi å Yi Yi å

å

i50

N

›u L, ›n

(14)

with L representing the distance between the ghost-cell and its image point. In (10), for solving the four unknowns (uI, ux, uy, uz), N is given a number greater than four to avoid the ill conditioning of the matrix. In the following experiments, the eight nearest points are selected. Finally, in this study, the Neumann boundary condition is applied to u and y, and the Dirichlet boundary condition is used for w. That is ›u ›y ›h ›h 5 5 0 and w 5 u 1 y , ›n ›n ›x ›y

(15)

1608


VOLUME 140

where h is the terrain height. The right part of (15) is also known as the terrain-induced vertical velocity.

3. Indices for verification, model configuration, and simulated results for the observing system simulation experiment To give a quantitative assessment of the accuracy of the retrieved results, the root-mean-square error (rmse), relative rms error (r_rmse), and spatial correlation coefficient (SCC) between the retrieved variable and its true counterpart are computed. The indices used by Gao et al. (1999) are employed and extended. For the horizontal wind (VH), these parameters are calculated by combining the u and y components at both time levels together:

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u2 u 2 2 u t i51 [(ur 2 ut )i 1 (y r 2 y t )i ] , (16) rmse(VH ) 5 43M vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u2 u 2 2 u u i51 [(ur 2 ut )i 1 (y r 2 y t )i ] u , (17) r_rmse(VH ) 5 u 2 u t [(ut )2i 1 (y t )2i ]

åå

åå

åå

i51

2

å å(ur 2 aur )i (ut 2 aut)i 1 (yr 2ayr )i (yt 2 ayt )i

SCC(VH) 5 i51

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Ar 3 At

, (18)

where 2

Ar,t 5

å å[(ur,t 2 aur,t )2i 1 (yr,t 2 ayr,t )2i ], 2

aur,t 5

(19)

i51

2

å å(ur,t )i

i51

2M

and ayr,t 5

å å(yr,t )i

i51

2M

.

(20)

In (16)–(20), the subscripts ‘‘r’’ or ‘‘t’’ denote the retrieved or true values, and the subscripts i (51, 2) represent time levels 1 or 2. The summation can be either over a two-dimensional plane or a threedimension space, and M is the total number of grid points used in the computation. The same calculation is also applied for the vertical wind, but the summation only considers w.

FIG. 2. The model-generated true horizontal wind field (vectors) at Z 5 750 m. Three radar stations—RCCG, RCKT, and RCGI—are illustrated by the square, circle, and triangle, respectively. The thin solid line is the coast line of Taiwan. The terrain heights are represented by the shading, with an interval of 500 m.

All the datasets used for the OSSE-type numerical experiments are generated from the WRF model. The model domain comprises 101 3 101 grid points in the horizontal and 31 s-layers in the vertical. The horizontal grid size and vertical resolution are 2.0 km and approximately 0.5 km, respectively. The model domain covers southern Taiwan, with realistic terrain height available at each grid point. The wind data at two time levels, with an interval of 3 min, are extracted from the model simulation then projected onto radar sites. The results are utilized as the simulated radial wind observations from two volume scans measured by each radar. No attempt is made to represent observational errors. It should also be pointed out that the synthesized results could depend on the data sampling frequency. The problem regarding the correction of the errors induced by advection and nonsimultaneous radar observations, as discussed in Gal-Chen (1982), Chong et al. (1983), Caillault and Lemaitre (1999), and Matejka (2002), is not addressed in the current study. Figure 2 shows the horizontal wind vectors (averaged over two time levels) on a horizontal plane close to the surface (Z 5 0.75 km), while Figs. 3a–c display the u, y, and w fields along a vertical cross section (Y 5 174 km), where the maximum vertical velocity is identified. The locations of the three radars (RCCG, RCKT, and RCGI) are also marked in Fig. 2. They are actual existing radar stations currently operated by the Central Weather Bureau and Air Force of Taiwan in this area. It can be

MAY 2012

LIOU ET AL.

1609

seen from Fig. 2 that the prevailing wind is coming from the southwest. The flows are deflected when approaching the topography of Taiwan. Figure 3a illustrates that the u component of the wind has a maximum speed about 16.5 m s21 at the top of the mountain. In Fig. 3b, the y component of wind reveals a maximum magnitude (.15.0 m s21) to the west of the mountain. After examining the flow structures along other east–west cross sections, this maximum turns out to be an indication of the deflection of the flow by the topography, resulting in an enhancement of the wind speed in the north–south direction. The w field, depicted in Fig. 3c, exhibits a clear wave pattern above the complex terrain. The strongest updraft, which occurs along the upslope of the mountain, can reach 3.5 m s21. In the next section a series of OSSE experiments is conducted. The primary purposes are to validate the implementation of the IBM method, to examine the performance of the newly designed synthesis algorithm when different numbers of radars and/or the terrain blockage are considered, and to investigate the role played by the background field, especially in data-void regions. Finally, a hypothetical design for a combined radar scanning strategy is proposed to optimize the usage of a multipleDoppler radar network. In all OSSE experiments, the background wind field is constructed using the sounding data employed for initializing the WRF simulation discussed above.

4. Results of the multiple-radar wind synthesis obtained using model-simulated datasets a. Retrieval with three radars, no terrain blockage (experiments NB_3R and NB_3R_w0)

FIG. 3. The model-generated true wind field over a vertical cross section at Y 5 174 km: (a) u, (b) y, and (c) w. The contour interval is 1.5 m s21 for (a) and (b), and 0.5 m s21 for (c). The solid (dashed) line denotes a positive (negative) value, which represents the westerly (easterly) wind in (a), the southerly (northerly) wind in (b), and updraft (downdraft) in (c). The topography is plotted by the dark shading.

In experiment NB_3R an idealized scenario is proposed in which all the grid points in the analysis domain, except those within the mountains, are covered by three radars, as shown in Fig. 2. The purpose is to test the performance of this new scheme after implementing the IBM technology for recovering the wind field when a set of flawless multiple-Doppler radial wind data is provided everywhere throughout the domain. Figures 4 and 5 depict the retrieved horizontal wind vectors and u, y, and w fields along the same cross sections as in Figs. 2–3. It can be seen that the retrieved three-dimensional wind fields are highly consistent with their ‘‘true’’ counterparts. The major features, such as the deflection of the horizontal flows by the topography, the locations of the maximum u and y, as well as the wave structures exhibited by w, are all successfully retrieved. However, the magnitudes of the vertical velocity tend to be underestimated, with a maximum value reaching only about

1610


VOLUME 140

FIG. 4. As in Fig. 2, but for the retrieved horizontal wind field (vectors) at Z 5 750 m from experiment NB_3R.

2.3 m s21. A quantitative comparison is shown in Table 1 and indicates that the rms errors of the horizontal and vertical winds are below 0.25 m s21. The spatial correlations between the true and retrieved VH and w can reach as high as 0.987 and 0.869, respectively. To further validate the synthesized results, the variables involving first-order derivatives of the winds are examined. These fields are known to be highly sensitive to the accuracy of the retrieved flow fields. Figures 6 and 7 demonstrate two examples in which the spatial distributions of the simulated and retrieved horizontal divergence (›u/›x 1 ›y/›y) and vertical vorticity (›y/›x 2 ›u/›y) along an east–west direction at a selected level (Z 5 2.25 km, Y 5 174 km) are compared. It can be seen that there are variations with significant amplitudes near the mountains (depicted by the broken sector of the curves), and that these features are successfully reproduced in the retrieved results. On the west side of the mountains, there is an identifiable convergence zone near X ; 102 km, associated with positive vorticity. The updraft (not shown) in this region can reach 1.5 m s21. Another characteristic of our method is that the retrieved three-dimensional winds will approximately satisfy the vertical vorticity equation, as expressed by (3a), in a least squares sense. Consequently, the vorticity budget diagnosis becomes rather straightforward. Figure 8 shows the retrieved vorticity tendency (›j/›t), or the left-handside term of (3a), at the same level as in Figs. 6 and 7. The right-hand side of (3a), namely, the advection term, divergence term, and tilting term, is illustrated in Fig. 9. It can be found that the summation of these three terms, indicated by the thick solid line in Fig. 9, is nearly

FIG. 5. As in Fig. 3, but for the retrieved wind field over a vertical cross section at Y 5 174 km from experiment NB_3R: (a) u, (b) y, and (c) w. The contour interval is 0.3 m s21 for (c).

identical to that shown in Fig. 8. It should be pointed out if the traditional wind synthesis method is applied at two time levels for conducting the same vorticity budget diagnosis, there is no constraint to enforce the satisfaction

MAY 2012

1611

LIOU ET AL.

TABLE 1. The rms error, r_rms error, and SCC for 11 OSSE experiments. The prefix B or NB indicate whether the terrain blockage is activated or not, respectively. The number in front of the letter R represents the number of radars used for computation. NO or ONLY stands for the exclusion or inclusion of a specific radar, and NOBG means that the background field is not utilized. The numbers inside the parentheses in NB_2R_NOKT, NB_2R_NOGI, and NB_2R_NOCG denote the values of RAD [defined by (21) in the text] computed for RCKT, RCGI, and RCCG, respectively. rms (m s21)

r_rms

SCC

Expt

VH

w

VH

w

VH

w

NB_3R NB_3R_w0 NB_2R_NOKT (64%) NB_2R_NOGI (68%) NB_2R_NOCG (72%) NB_3R_NOBG NB_2R_NOKT_NOBG NB_1R_ONLYCG NB_1R_ONLYCG_NOBG B_3R B_3R_CMBN

0.249 0.256 0.423 0.527 0.546 0.0544 0.451 1.060 6.228 0.710 0.991

0.243 0.295 0.257 0.262 0.305 0.177 0.270 0.316 0.829 0.296 0.320

0.0272 0.0279 0.0462 0.0575 0.0597 5.94 3 1023 0.0492 0.116 0.680 0.0776 0.108

0.585 0.712 0.618 0.631 0.734 0.427 0.650 0.761 1.997 0.714 0.771

0.987 0.986 0.984 0.982 0.981 0.989 0.984 0.960 0.439 0.974 0.962

0.869 0.726 0.854 0.841 0.765 0.903 0.785 0.714 0.219 0.747 0.714

of the vorticity equation. Thus, the right-hand and lefthand sides of (3a) may not always balance each other and will produce a so-called residual term. In some cases, the magnitude of this residual term turns out to be comparable to that of the other terms, which makes it difficult to provide a reasonable interpretation of the changes in vorticity. Another experiment (NB_3R_w0) is designed in which the w is set to zero along the ground. In other words, we pretend the surface is flat, even though it is not. A comparison of the retrieved wind fields from experiments NB_3R and NB_3R_w0 reveals that the major

FIG. 6. The spatial variation of horizontal divergence (›u/›x 1 ›y/›y) along an east–west direction at Z 5 2.25 km, Y 5 174 km. The solid line is the true solution from the model simulation. The dashed line shows the retrieved result from experiment NB_3R. The magnitude of the divergence has been multiplied by 103. The discontinuity of the curves is caused by the terrain blockage.

differences (not shown) taking place at layers close to the terrain. The overall accuracy of the retrieved threedimensional wind, as expected, turns out to be lower than that from experiment NB_3R, particularly for w. The statistics are also included in Table 1. For example, it can be seen that the relative rms error for w increases significantly from 0.585 in NB_3R to 0.712 in experiment NB_3R_w0. This experiment demonstrates the usefulness of imposing the new boundary condition for w. As a supplemental test, the retrieval experiment is repeated by withholding the zero-normal-gradient condition for u, y as shown in (15). It is found that in the resulting three-dimensional wind fields (not shown) all the major features exhibited in the true solutions are still well

FIG. 7. As in Fig. 6, but for the spatial variation of vertical vorticity (›y/›x 2 ›u/›y) from the model simulation and experiment NB_3R.

1612


VOLUME 140

FIG. 8. The spatial variation of the temporal derivative of the vertical vorticity along an east–west direction at Z 5 2.25 km, Y 5 174 km from experiment NB_3R. The magnitude of the curve has been multiplied by 106. This is also the left-hand-side term from (3a) in the text.

FIG. 9. As in Fig. 8, but for the right-hand side of (3a) in the text. The thin solid line is the advection term, the long dashed line denotes the divergence term, the dotted line represents the tilting term, and the thick solid line is the summation of all three terms.

recovered. However, a quantitative assessment reveals that the rms errors of the retrieved horizontal and vertical velocities in this experiment increase from those of NB_3R by about 14% and 17%, respectively. This result demonstrates the value of using the immersed boundary method and a proper treatment of the bottom boundary conditions. The aforementioned results suggest that the IBM has been successfully implemented into the wind synthesis scheme and functions well. The statistics computed above are used as reference for the following experiments.

speed. This index is certainly a function of the magnitudes and directions of the flow field, and the geometric location of the radar. The formula used for this purpose is

b. Retrieval with two radars, terrain blockage not included (experiments NB_2R_NOKT, NB_2R_NOGI, and NB_2R_NOCG) This experiment consists of three tests (NB_2R_NOKT, NB_2R_NOGI, and NB_2R_NOCG). In each test, one of the three radars illustrated in Fig. 2 is excluded, meaning that only two radars are used to conduct the wind synthesis. By examining the retrieved three-dimensional winds (not shown), it is found that similar to the finding by LC09, the flow structures along or near the baseline formed by any two radars, can still be well recovered. This is a major advantage of the variational approach over the traditional method. Nevertheless, Table 1 reveals that the accuracy of the synthesized wind fields depends on the different combinations of the radars, as shown by the rms errors, as well as the SCC values of the horizontal and vertical winds. To explain this discrepancy, we estimate the ratio (RAD) between wind information detectable from each radar (i.e., the radial wind), and the total wind

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (Vr )2 ffi. RAD 5 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (u2 1 y 2 1 w2 )

å

å

(21)

The summation is carried out over all grid points within the flow region. Thus, if a grid point is blocked by the mountains, and cannot be reached by the radar, as in the next experiment, then the radial wind at this location is given a value of zero, meaning no wind information is available at this point. The computed RADs for the RCKT, RCGI, and RCCG radars are 64%, 68%, and 72%, respectively. The difference in the ratio of RAD implies that for this particular flow field and the relative locations of the three radars, RCCG (RCKT) is able to detect the most (least) amount of wind information. Thus, if the wind synthesis is performed without RCCG (RCKT) (i.e., experiments NB_2R_NOCG and NB_2R_NOKT), the computation would be subjected to stronger (weaker) impact, and is expected to produce the worst (best) retrievals. This is indeed confirmed by the statistics shown in Table 1, as NB_2R_NOCG (NB_2R_NOKT) has the lowest (highest) SCC and largest (smallest) rms errors. However, it should be pointed out that although NB_2R_NOCG produces the worst retrieved results, the major features described above, such as the enhancement of the north–south winds and the wave structure in w, are still well reproduced. This is particularly true for the horizontal winds, as

MAY 2012

LIOU ET AL.

revealed by their high SCC values (0.981). The quality of the retrieved vertical winds is also comparable to that reported by other studies (e.g., Gao et al. 1999).

c. The role of the background wind field (experiments NB_3R_NOBG, NB_2R_NOKT_NOBG, NB_1R_ONLYCG, and NB_1R_ONLYCG_NOBG) In this section, the role played by the background field is investigated using four experiments. Experiment NB_3R_NOBG is the same as NB_3R, but the information from the background wind field is not included. In other words, the constraint associated with the background term in the cost function is turned off, and the initial guess for all retrieved variables are set to zero. Similarly, in the second test, NB_2R_NOKT_ NOBG is conducted by repeating NB_2R_NOKT, but without the background term. Note that judging from the radial winds that are detectable by each radar (see the ratio RAD in section 4b), this test uses the observations from RCCG and RCGI that can provide more information to the retrievals. In the third and fourth tests, the radial wind data from only RCCG (which provides the most wind information) are utilized, and NB_1R_ONLYCG and NB_1R_ONLYCG_NOBG stand for the experiments with and without the background term, respectively. It is found that NB_3R_NOBG takes a longer time to converge, but the statistics shown in Table 1 indicate an even closer agreement with the true solutions than in NB_3R. In contrast, when two radars are involved, better results are produced from NB_2R_NOKT than those from NB_2R_NOKT_NOBG. For experiment NB_1R_ONLYCG, Table 1 shows that the results are further degraded, but still remain at an acceptable level. An examination of the retrieved fields (not shown) reveals that the principal features embedded in the flow region can still be recovered, even though data from only single Doppler radar are employed. Finally, if the background information is removed (NB_1R_ONLYCG_NOBG), the accuracy of the retrieved fields drops very dramatically. It is known that although the background term represents a description of the flow’s general pattern, it contains errors, since the detailed small scale features are usually not resolved. When the grid points are covered by a sufficient number of radars (e.g., three), the retrieval algorithm receives enough information to accurately reconstruct the wind fields. Under this condition, the inclusion of the background field turns out to worsen (although not significantly) instead of improving the retrievals. There is a trade-off between the accuracy and the computational efficiency. When a grid point is

1613

covered by two radars, having the background term does help to improve the accuracy. However, when only a single Doppler radar is available, our study shows that the extra information provided by the background term will become extremely crucial for the success of the retrievals. An interesting finding from NB_1R_ONLYCG is that when the general structure of the flow field can be provided by a background wind field, our method has the potential to recover the flow field using observations from only single Doppler radar. This is also mentioned in Gao et al. (1999).

d. Retrieval with three radars, terrain blockage is considered (experiment B_3R) In Fig. 2, the radars illustrated in the analysis domain are actually three existing operational radar stations. Among them, the RCKT and RCCG are S-band Doppler radars both operated by the Central Weather Bureau, and the RCGI is a C-band Doppler radar run by the Air Force of Taiwan. In this experiment, we let all radars conduct volume scans with discrete elevation angles. For RCKT and RCCG, nine elevation angles (0.58, 1.48, 2.48, 3.48, 4.38, 6.08, 9.98, 14.68, and 19.58) are employed. For the RCGI, the scan comprises 13 elevation angles (0.48, 1.08, 1.78, 2.68, 3.68, 4.88, 6.08, 7.58, 9.58, 12.58, 16.58, 24.08, and 38.08). All the radars are assumed to have a 18 beamwidth. The maximum detectable range for RCKT and RCCG is 230 km, and for RCGI it is 120 km. In addition, the blockage of the radar beams by the topography is also taken into account, which means that grid points behind the mountains become data-void regions. The design of this experiment represents a more realistic scenario, in which the radar coverage is much less than in the previous experiments. A quantitative estimation shows that the ratio of the grid points covered by RCKT, RCCG, and RCGI is 66%, 58%, and 51% of the entire domain, respectively. Furthermore, the values of RAD, as defined by (21), are 54%, 57%, and 50% for RCKT, RCCG, and RCGI, respectively. Even with the incomplete radar data coverage, the retrieved flow structures are found to be in good agreement with their true counterparts. Figure 10 presents one example, in which one of the key features exhibited in the true solution, namely, the wave motion in the w field, is still well retrieved. The statistics for the retrieved three-dimensional wind field are revealed in Table 1. Although somewhat degraded from the previous experiments, they remain rather satisfactory. For example, the relative rms error of the horizontal winds is only about 7.76%. In addition, the SCC between the retrieved and true w reaches above 0.747. It was suggested from the results discussed in sections 4b,c that the accuracy of the retrievals is dependent on

1614


FIG. 10. As in Fig. 5c, but for the retrieved vertical velocity field over a vertical cross section at Y 5 174 km from experiment B_3R. The contour interval is 0.3 m s21.

coverage by a sufficient number of radars. This argument can be further elaborated by showing the retrievals along the vertical direction. We compute the ratio between the number of points covered by at least two radars and the total number of grids at each horizontal plane. The relative rms error and SCC used for verification are also displayed as a function of height. Figure 11 indicates that when the ratio of points with multiple ($2 radars) radar coverage is one (in experiment NB_3R), which represents an idealized case with data from all three radars available everywhere, the relative rms errors and SCC of the horizontal wind remains below 5% and close to 1.0, respectively, as revealed in Fig. 12. In contrast, in experiment B_3R, the ratio exceeds 0.9 only at approximately Z 5 2–4 km, then decreases to about 0.2 at 12–14 km. It is noticed that the accuracy of the retrieval has peak values near 4 km, where the radar coverage is the best, then starts to deteriorate with height. A larger relative rms error and smaller SCC are found when the altitudes of less radar coverage are reached. The accuracy of the retrieved w along the vertical direction behaves similarly, thus the results are not shown.

e. A hypothetical design for a combined radar scanning strategy (B_3R_CMBN) In this section, a hypothetical combined radar scanning strategy is proposed. The basic concept is to let each radar focus on probing different areas of the domain. An example is demonstrated in which we let two of the radars (RCKT and RCGI) scan the lower to middle portion of the domain from approximately 0.5 to 6.5 km, while the third radar (RCCG) is allowed to observe the middle to upper layers from about 5.5 to 11.5 km. The

VOLUME 140

FIG. 11. The ratio between the number of grid points covered by at least two radars and the total number of points at each height. The solid line is the ratio from experiment NB_3R and the dashed line is for experiment B_3R.

purpose of this design is to form multiple-Doppler ($2) coverage in the lower to middle atmosphere, so that the weather systems existing in this region can be well resolved. However, from the middle to upper layers, the retrievals reply on the information provided by the single Doppler observed radial velocity and/or the background fields. Apparently, when more radar stations are available, one could certainly expand the area covered by multiple-Doppler radars. Figure 13 illustrates the retrieved vertical velocity, superimposed by gray scales representing the coverage by different numbers of radars. It is seen that even though the data used for the retrievals in each region come from different sources, the major wave pattern that characterizes the w field can be smoothly retrieved over the entire domain. The rms errors and SCC values listed in Table 1 also indicate the consistency of the retrieved three-dimensional winds with their modelgenerated counterparts. The authors believe that because of the capability of our method to merge observational data of different natures, it should offer greater flexibility to design various scanning strategies when a multipleDoppler radar network is used.

5. Testing the multiple-radar synthesis method using field experiment datasets The method used in LC09 only deals with a flat bottom boundary. Its lowest layer has to be higher than the mountains. Thus, it can only recover the wind fields above the mountain peaks, even though the radar data might be available along the mountain slopes. The results from section 4 indicate that the newly designed method in this

MAY 2012

LIOU ET AL.

1615

FIG. 13. The retrieved vertical velocity field from experiment B_3R_CMBN. The blank, lightly, moderately, and heavily hatched regions represent the radar coverage by zero, one, two, and three radars, respectively. The topography is denoted by the dark shading.

portrays the relative positions of these three radars, the range of the analysis domain, as well as the topography of southwestern Taiwan. The analysis domain covers an area of 328 km 3 546 km 3 12.25 km. The

FIG. 12. The vertical distribution of the (a) relative rms error and (b) SCC of the retrieved horizontal wind. The solid and dashed lines are for experiments NB_3R and B_3R, respectively.

research has the potential to retrieve the wind field everywhere in the flow region, including the mountain slopes. In this section, the performance of the proposed multiple-radar synthesis method equipped with the immersed boundary method is investigated using the radar data collected during SoWMEX. This field experiment was conducted from May to June 2008 in the southwestern part of Taiwan. The scientific purpose was to understand the mechanisms that trigger the heavy precipitation in the South China Sea and Taiwan during the Asian summer monsoon season, so as to improve the accuracy of the quantitative precipitation forecasts in this area. The case selected in this study is from the intensive observing period 8 (IOP8) on 14 June 2008. Observational data from three S-band radars are utilized: S-POL (operated by the National Center for Atmospheric Research); and RCKT and RCCG (operated by the Central Weather Bureau of Taiwan). Figure 14

FIG. 14. The retrieved horizontal wind field (vectors) at Z 5 1.0 km, superimposed on Taiwan’s topography (shading at 1000-m intervals), and the locations of the data collected from different sources. The RCKT, S-POL, and RCCG are denoted by a solid circle, triangle, and square, respectively. The open circle indicates the radiosonde station. The cross inside a circle represents the surface station. The soundings released at Lyu-Dao (LD), PingTung (PT), and Yong-Kang (YK) are used for verification.

1616


FIG. 15. The u component of the wind at the (a) LD, (b) PT, (c) YK stations. The solid and dashed lines are from retrievals and sounding observations, respectively. The solid triangle, circle, square, and diamond represent that this particular layer is covered by zero, one, two, and three radars, respectively. The area above Z 5 5 km at LD station is mainly covered by only one radar. In contrast, the PT and YK sites are covered by two or three radars.

horizontal and vertical resolutions are 2.0 and 0.5 km, respectively. Also shown in Fig. 14 are the locations of nine surface stations and six radiosonde stations. The data collected by soundings released at Lyu-Dao (LD), Ping-Ting (PT), and Yong-Kang (YK) are used for verifications. Wind data measured by the other three soundings and the surface stations are used to prepare the background wind field. They are first interpolated to the model levels, then to the model grid points using a two-dimensional Barnes scheme with a 200-km radius

VOLUME 140

FIG. 16. As in Fig. 15, but for the y component of the wind.

of influence. Thus, the background wind field in this real case contains certain horizontal structure, which does not exist when only one sounding is used for setting up the background field. The multiple-Doppler retrieved horizontal wind at Z 5 1.0 km is superimposed on Fig. 14. It is shown that the prevailing wind moves from the oceans to the southwestern part of Taiwan. Tai et al. (2011) has pointed out that the warm and moist air brought by this southwesterly flow produces heavy precipitation in southern Taiwan (;200 mm day21). The flow is deflected northward at the western flank of the mountain. To validate the wind fields, Figs. 15 and 16 illustrates a comparison of the retrieved u and y components of wind against the

MAY 2012

LIOU ET AL.

1617

FIG. 17. The retrieved vertical velocity field (contour interval 0.3 m s21) for IOP 8 of SoWMEX at Y 5 276 km. The blank, lightly, moderately, and heavily hatched regions represent the grid points covered by zero, one, two, and three radars, respectively. The terrain is denoted by the dark shading.

measurements taken by the LD, PT, and YK radiosondes. The comparison is made only for the lower to middle layer of the atmosphere (i.e., ,7.0 km), since the soundings may drift too far away from the site. Also graphically displayed is the number of radars covering each layer. It can be seen that due to the terrain blockage and geographic location, the LD site is only observed by one radar at the middle altitude. In contrast, the areas above PT and YK stations are better covered by at least two radars. Consequently, the retrieved horizontal winds at PT and YK are expected to have higher accuracy. This can be verified by computing the differences between the retrieved and the observed horizontal winds available at LD, PT, and YK, for which the rms errors are 2.87, 2.69, and 1.98 m s21, respectively. It should be pointed out that although the wind profile at LD is obtained by relying on single-Doppler data, the error is still comparable to that in other relevant studies, such as Sun et al. (2010) or Tai et al. (2011). Finally, Fig. 17 displays the vertical velocity field along a east– west-oriented vertical cross section at Y 5 276 km, where the retrieved strongest downdraft in the domain is identified. Lack of direct measurement of the vertical wind makes verification difficult. However, the coverage by a different number of radars is exhibited for reference. Based on the results from previous experiments, it is believed that the retrieved flow field in those regions with a sufficient number ($2) of radar coverage should be rather reliable. Updrafts are found on the western (upstream) side of the mountains. Unfortunately, the eastern (downstream) side of the mountains is not reachable by

any radar. Nevertheless, the downward motion near the lee side of the mountains appears to be reasonable.

6. Conclusions and future work In this study we have implemented the immersed boundary method to a variational-based multipleDoppler radar wind synthesis algorithm. The performance of the newly designed method is investigated using model-simulated datasets, and multiple-Doppler radar observations collected during IOP8 of the 2008 SoWMEX field project. In the latter, the verifications are conducted against sounding measurements. The experimental results show that the wind fields over a complex topography can be well retrieved, while the advantages of the original method reported in LC09, are still maintained. These advantages include its capability to recover the flow fields along and near the radar baseline, the direct application of the retrieved wind fields to vorticity budget diagnosis, and more selections in terms of specifying the top and bottom boundary conditions for the horizontal and vertical winds. Furthermore, this method is able to merge data available from different sources, and utilize data from any number of radars. This provides us with more flexibility in designing various scanning strategies, so that one may probe the atmosphere more efficiently using a multipleradar network. In future work, the accuracy of the retrieved vertical velocity needs to be carefully validated using instruments, such as a vertically pointing radar. This

1618


method should have the potential to produce a complete three-dimensional wind field using data from all radars in Taiwan. Acknowledgments. The authors are very grateful to the participants of the SoWMEX project. Special thanks go to the crew members who operated the S-POL, RCCG, and RCKT. This research is supported by the National Science Council of Taiwan under Grants NSC100-2625-M-008-002 and NSC100-2111-M-008-005.

REFERENCES Alaoui, A. H., L. White, and A. Shapiro, 2010: Estimation of vertical wind velocity using Doppler radar data with vorticity constraints. Appl. Math. Comput., 215, 4119–4131. Armijo, L., 1969: A theory for the determination of wind and precipitation velocities with Doppler radars. J. Atmos. Sci., 26, 570–573. Bousquet, O., and M. Chong, 1998: A multiple-Doppler synthesis and continuity adjustment technique (MUSCAT) to recover wind components from Doppler radar measurements. J. Atmos. Oceanic Technol., 15, 343–359. ——, P. Tabary, and J. Parent du Chatelet, 2007: On the use of operational synthesized multiple-Doppler wind fields. Geophys. Res. Lett., 34, L22813, doi:10.1029/2007GL030464. ——, T. Montmerle, and P. Tabary, 2008a: Using operationally synthesized multiple-Doppler winds for high resolution horizontal wind forecast verification. Geophys. Res. Lett., 35, L10803, doi:10.1029/2008GL033975. ——, P. Tabary, and J. Parent du Chatelet, 2008b: Operational multiple-Doppler wind retrieval inferred from long-range radial velocity measurements. J. Appl. Meteor., 47, 2929– 2945. Brandes, E., 1977: Flow in severe thunderstorms observed by dualDoppler radar. Mon. Wea. Rev., 105, 113–120. Caillault, K., and Y. Lemaitre, 1999: Retrieval of three-dimensional wind fields corrected for the time-induced advection problem. J. Atmos. Oceanic Technol., 16, 708–722. Chong, M., and J. Testud, 1983: Three-dimensional wind field analysis from dual-Doppler radar data. Part III: The boundary condition: An optimum determination based on a variational concept. J. Climate Appl. Meteor., 22, 1227–1241. ——, and S. Cosma, 2000: A formulation of the continuity equation of MUSCAT for either flat or complex terrain. J. Atmos. Oceanic Technol., 17, 1556–1564. ——, and O. Bousquet, 2001: On the application of MUSCAT to a ground-based dual-Doppler radar system. Meteor. Atmos. Phys., 78, 133–139. ——, J. Testud, and F. Roux, 1983: Three-dimensional wind field analysis from dual-Doppler radar data. Part II: Minimizing the error due to temporal variation. J. Climate Appl. Meteor., 22, 1216–1226. ——, and Coauthors, 2000: Real-time wind synthesis from Doppler radar observations during the Mesoscale Alpine Programme. Bull. Amer. Meteor. Soc., 81, 2953–2962. Doviak, R. J., P. Ray, R. G. Strauch, and L. J. Miller, 1976: Error estimation in wind fields derived from dual-Doppler radar measurement. J. Appl. Meteor., 15, 868–878.

VOLUME 140

Drechsel, S., M. Chong, G. J. Mayr, M. Weissmann, R. Calhoun, and A. Dornbrack, 2009: Three-dimensional wind retrieval: Application of MUSCAT to dual-Doppler lidar. J. Atmos. Oceanic Technol., 26, 635–646. Gal-Chen, T., 1982: Errors in fixed and moving frame of references: Applications for conventional and Doppler radar analysis. J. Atmos. Sci., 39, 2279–2300. Gao, J., M. Xue, A. Shapiro, and K. K. Droegemeier, 1999: A variational method for the analysis of three-dimensional wind fields from two Doppler radars. Mon. Wea. Rev., 127, 2128– 2142. ——, ——, K. Brewster, and K. K. Droegemeier, 2004: A threedimensional variational data analysis method with recursive filter for Doppler radars. J. Atmos. Oceanic Technol., 21, 457– 469. Georgis, J. F., F. Roux, and P. H. Hildebrand, 2000: Observation of precipitation systems over complex orography with meteorological Doppler radars: A feasibility study. Meteor. Atmos. Phys., 72, 185–202. Liou, Y.-C., and Y.-J. Chang, 2009: A variational multiple-Doppler radar three-dimensional wind synthesis method and its impacts on thermodynamic retrieval. Mon. Wea. Rev., 137, 3992– 4010. Lundquist, K. A., F. K. Chow, and J. K. Lundquist, 2010: An immersed boundary method for the Weather Research and Forecasting Model. Mon. Wea. Rev., 138, 796–817. Matejka, T., 2002: Estimating the most steady frame of reference from Doppler radar data. J. Atmos. Oceanic Technol., 19, 1035– 1048. Mewes, J. J., and A. Shapiro, 2002: Use of the vorticity equation in dual-Doppler analysis of the vertical velocity field. J. Atmos. Oceanic Technol., 19, 543–567. Miller, L. J., and R. G. Strauch, 1974: A dual-Doppler radar method for the determination of wind velocities within precipitating weather systems. Remote Sens. Environ., 3, 219–235. Mittal, R., and G. Iaccarino, 2005: Immersed boundary methods. Annu. Rev. Fluid Mech., 37, 239–261. O’Brien, J. J., 1970: Alternative solutions to the classical vertical velocity problem. J. Appl. Meteor., 9, 197–203. Peskin, C. S., 1972: Flow patterns around heart valves: A numerical method. J. Comput. Phys., 10, 252–271. Protat, A., and I. Zawadzki, 1999: A variational method for realtime retrieval of three-dimensional wind field from multipleDoppler bistatic radar network data. J. Atmos. Oceanic Technol., 16, 432–449. ——, and ——, 2000: Optimization of dynamic retrievals from a multiple-Doppler radar network. J. Atmos. Oceanic Technol., 17, 753–760. Ray, P. S., R. J. Doviak, G. B. Walker, D. Sirmans, J. Carter, and B. Bumgarner, 1975: Dual-Doppler observation of a tornadic storm. J. Appl. Meteor., 14, 1521–1530. ——, K. K. Wagner, K. W. Johnson, J. J. Stephens, W. C. Bumgarner, and E. A. Mueller, 1978: Triple-Doppler observations of a convective storm. J. Appl. Meteor., 17, 1201–1212. ——, C. L. Ziegler, W. Bumgarner, and R. J. Serafin, 1980: Singleand multiple-Doppler radar observations of tornadic storms. Mon. Wea. Rev., 108, 1607–1625. Scialom, G., and Y. Lemaitre, 1990: A new analysis for the retrieval of three-dimensional mesoscale wind fields from multiple Doppler radar. J. Atmos. Oceanic Technol., 7, 640–665. Shapiro, A., and J. J. Mewes, 1999: New formulations of dualDoppler wind analysis. J. Atmos. Oceanic Technol., 16, 782– 792.

MAY 2012

LIOU ET AL.

——, S. Ellis, and J. Shaw, 1995: Single-Doppler velocity retrievals with Phoenix II data: Clear air and microburst wind retrievals in the planetary boundary layer. J. Atmos. Sci., 52, 1265–1287. ——, C. K. Potvin, and J. Gao, 2009: Use of a vertical vorticity equation in variational dual-Doppler wind analysis. J. Atmos. Oceanic Technol., 26, 2089–2106. Sun, J., M. Chen, and Y. Wang, 2010: A frequent-updating analysis system based on radar, surface, and mesocale model data for the Beijing 2008 Forecast Demonstration Project. Wea. Forecasting, 25, 1715–1735.

1619

Tai, S.-L., Y.-C. Liou, J. Sun, S.-F. Chang, and M. C. Kuo, 2011: Precipitation forecasting using Doppler radar data, a cloud model with adjoint, and the Weather Research and Forecasting model: Real case studies during SoWMEX in Taiwan. Wea. Forecasting, 26, 975–992. Tseng, Y., and J. Ferziger, 2003: A ghost-cell immersed boundary method for flow in complex geometry. J. Comput. Phys., 192, 593–623. Ziegler, C. L., P. S. Ray, and N. C. Knight, 1983: Hail growth in an Oklahoma multicell storm. J. Atmos. Sci., 40, 1768–1791.