Toward Deducing the 3-D Thermodynamic Structure from Only Wind Fields-A Validation Study using Simulated Dataset

This paper presents a newly designed thermodynamic retrieval algo­ rithm whereby one can deduce the potential temperature and pressure gra­ dient fields in a three dimensional space using only 3-D wind measurements. The latter can be observed by Doppler radar in real cases. In order to achieve this purpose, all three equations of motion are implemented in a single cost function. Then, given the detailed wind information, and through the pro­ cess of variational minimization, these equations are solved simultaneously in a least square sense. The products of the approximated solutions are the three dimensional potential temperature, and the pressure gradients along any direction. Some preliminary studies using model-generated data sets show that this is a feasible tool for the retrievals, and satisfactory results can be obtained when the performance of this method against various types of input data errors is investigated. Compared to the traditional method proposed by Gal-Chen, the advantage of this scheme is that the vertical structure of the thermodynamic variables can be determined without re­ quiring a point of independent observation of the pressure and tempera­ ture for each layer. In addition, this method can easily be applied to an area with topography. It is believed that this algorithm can be of particularly useful in many Doppler weather radar applications. (


INTRODUCTION
It is known that although Doppler radar can provide information about wind with a fine spatial and temporal resolution, the structure of thermodynamic variables, such as temperature and pressure, cannot be directly detected.The pioneer work by Gal-Chen ( 1978, hereafter GC) demonstrated that given the detailed Doppler radar derived 3-D wind observations, and their time derivatives through consecutive volume scans, one could extract the information of the pressure and temperature perturbations.This powerful method turned out to be particularly popular, and was later widely adopted by many researchers to deal with various weather phe nomena (e.g., Hane et al. 1981;Gal-Chen and Kropfli 1984;Hane and Ray 1985;Lin et al. 1986;Parsons et al. 1987, among others).Despite the encouraging results, it should be empha sized that GC only retrieves the deviations (from the horizontal mean) of the pressure and temperature perturbations (with respect to a base state) for each horizontal plane .That is, only 8' -< 8 ' > and p' -< p'> are solved uniquely at each horizontal layer, where<> stands for a spatial average over that layer.In order to specify the values of< 8 ' >and< p'> , and deter mine the pressure and temperature fluctuations 8 ' and p', at least one point of independent field measurement of the temperature and pressure for each altitude must be available.Al though a sounding may provide such information, unfortunately, considering the scenario of daily operations, difficulties persist.First, one has to synchronize the sounding operation with the radar observations.But it is known that the former is usually conducted only twice a day.Secondly, the trace of the sounding may not always remain within the region of the radar volume scans.The consequence of this drawback is that, whenever an interpretation of the vertical structure of the thermodynamic fields is required, one must proceed with caution.
A new approach was proposed by Roux (1985, 19 88), whereby the pressure and tempera ture perturbations can be solved uniquely up to a single constant.To deduce this unknown constant requires only one independent pressure and temperature observation at a single point somewhere in the retrieval domain.Roux and Sun (1990) modified the method by Roux (1988) by taking into account a simplified form of the thermodynamic equation everywhere in the domain, so that the temperature gradient was provided both horizontally and vertically.Addi tional improvements in Roux's retrieval technique can be found in Roux et al. (1993).On the other hand, with the rapid development of the adjoint method, Sun and Crook (1996) demon strated the advantages of using the so-called 4D-Var technique over the traditional GC scheme for retrieving the thermodynamic field.In contrast to the previous techniques, an important feature shown by the 40-Var formulation is that the temperature perturbations can be inferred without the auxiliary of any extra, independent measurements .
In this research, a different approach is suggested.The basic methodology involves imple menting the horizontal and vertical equations of motion into a single cost function.Through the method of variational analysis, a set of optimally-determined potential temperature fluc tuations and pressure perturbation gradients, that minimize this cost function, can be deduced simultaneously.
This paper is organized as follows.In the next section the proposed method is discussed in detail.Section 3 gives the data sets used for verification.Comparisons and discussions of the retrieval results are presented in Section 4 , followed by the conclusions in Section 5.

BASIC FORMULATION OF THE RETRIEVAL SCHEME
In order to fully investigate the performance and properties of this method, the validation is carried out using artificial data sets generated by a numerical model.To simplify the problem, the air is assumed to be dry.The basic momentum equations may be written as follows: The Lagrangian time derivative is defined as:

D
a a a a -=-+u-+v-+w-Dt at

ax (}y
Ck ' ( 1) (2) (3) in which Pis the pressure, P0=100 kPa, R is the gas constant, and CP is the specific heat capacity at a constant pressure.The subscript "O" represents a horizontally homogeneous ba sic state, from which the non-hydrostatic perturbations are expressed by variables with a single prime.It can be seen that F, G and H are functions of the winds, as well as the basic state temperature 80• The three-dimensional air motion can be estimated through either the mul tiple-Doppler radar synthesis procedure (e.g., Armijo 1969;Ray et al. 1980), or single Dop pler wind retrievals (e.g., Sun et al. 1991;Qiu and Xu 1992;Laroche and Zawadzki 1994;Shapiro et al. 1995;Liou 1999, and others), while 80 can be determined by an environmental sounding.
A cost function can be formulated based on (1)-( 3), as shown in the following: where J = � J Jf 1, ( a/�2 + a);z + a3�2 + tX4fi2 + asPs2 + a6�2 ) dn. , ( 6 )   dQ. = dxdydzdt , P. =(an ' _ F ) ((J n' 8' ) In (6), the first two terms ( Pi.2, �2) represent the differences between the retrieved hori zontal pressure gradients and the calculated F and G (using the Doppler derived 3-D wind observations in real applications), respectively.P3 2 denotes the difference in the retrieved vertical pressure gradient, the buoyancy, and the quantity H.Although this coefficient is specified empirically, the general patterns of the retrieved fields are not too sensitive to the value of a4• It is known that although the total potential tempera ture is increasing with height, the vertical gradient of the potential temperature perturbations can be either upward or downward.In other words, () 8' I dz< O is physically true.However, our experimental results show that this constraint can effectively eliminate the "spurious" thermal instability, while the region of "true" thermal instability remains intact.Furthermore, since P.i 2 is added as a "weak" constraint, therefore, when adequate weighting is given, it would not make ()()'/Ck< O disappear completely.Although interpreted differently, this constraint was also included in the retrieval technique proposed by Roux and Sun (1990), as well as in the improved version suggested by Roux et al. (1993).Sun and Crook (1996) also adopted this constraint for the same purpose in their 4D-Var retrieval formulations.
Finally, constraints P/-� 2 represent low-pass filters to suppress noise.Based on Testud and Chong (1983), the weighting coefficients, a5 and a6, are estimated by: --k-: µ <Xs -a 6 -N -N, (9) where k0 denotes the cutoff wave number, and N is the total number of grid points in the retrieval domain.If 2 � is chosen as the cutoff wavelength, where � is the horizontal grid resolution, then the corresponding value ofµ equals (�/ Jr)4 -0.01 � 4 .
From the first three cost function constraints, it is clear that the potential temperature appears in its undifferentiated form, but the pressure is always expressed by first order derivatives.Consequently, when ( 6) is minimized, the results will provide the structure of the buoyancy field itself in a three-dimensional space, as well as the pressure fluctuation gradients, not only for each horizontal plane, but also in the vertical direction.In other words, the poten tial temperature field itself is determined, as is the pressure field up to a "volume wide" constant.
No additional information is needed.We believe this to be an important advance over the traditional GC scheme.However, when it is desirable to compare the retrieved pressure per turbations against the true solutions, one should remove the average over the entire domain from each field, respectively.
In order to minimize (6), which is a cost function with very large number of dimensionalities, a quasi-Newton conjugate-gradient algorithm, named VA15AD by Liu and Nocedal (1989), is employed.This method, however, requires the user supplied knowledge of the cost function gradients with respect to Jr ' , ()' at each point (i.e., ()] / d!r' and ()] / ()()' ).This can be de rived by attributing the variation of the cost function, represented by 8J, to the variations of Jr', e', denoted by oJr ' and o e' respectively, throughout the interior region as well as over the boundaries.That is: The readers can refer to Liou (1999) for the details of th .e derivations.For example, for the interior points, the gradients are expressed by: ( Using an initial guess (usually zero) for the unknown variables n ' and 8 ', and the gradi ents for J computed by (1 1)-( 12), subroutine VA15AD determines a search direction whereby it finds a new set of n' and e ' that reduce the value of the cost function J.This new solution becomes another initial guess, and the searching process continues until the reduction of the cost function converges to zero.Unlike the traditional GC algorithm in which the spatial de rivatives of n ' is specified along the boundaries using wind observations in order to solve a Poisson equation for n ' , one advantage of this method is that the boundary conditions for n ' and e ' are not required in the minimization procedure.
All existing thermodynamic retrieval techniques start by using real Doppler radar derived wind data, in which observational errors are already embedded, to estimate those terms on the left hand side of momentum equations shown by ( 1)-(3) (i.e., the local time rate of change, advection, Coriolis force, and sub-grid parameterizations).The estimations are usually con ducted in finite difference form, and it is known that this procedure would inevitably amply the errors.The author would like to point out that the minimization algorithms developed by Roux (1 985,1 988), and Roux and Sun (1990) require additional spatial derivatives of the mo mentum equations, and this procedure is likely to further amplify the impacts of the observa tional errors.By contrast, in this research, the momentum equations are directly applied in their original forms without any further derivatives.Thus, the problem mentioned above can be avoided.

INPUT DATA SETS AND VERIFICATION
In this paper the data sets utilized for validation come from the numerical simulation of a group of thermal bubbles.The model domain contains a total of 55 X 55 X 55 grid points.The horizontal and vertical resolutions are specified to be 1.0 km and 0.15 km, respectively.Thus, the simulation is performed within a three-dimensional domain with a volume of approxi mately 55 X 55 x 8 km3• The background atmosphere is assumed to be neutral, with a basic state potential temperature of 300°k.The temperature perturbations are given at the layer of Z = 225m.The formulation of each thermal bubble is expressed by: where e: denotes the maximum magnitude of the bubble, ( x e , Y e ) are the coordinates of the bubble's center, and ( x,, y ,) represent the half-widths of the perturbation in x and y directions respectively.In this study they both are given the value of 10 km.A group of five bubbles is given to the initial temperature perturbation field, with their centers located at (28km, 28km), (14km, 28km), (42km, 28km), (28km,14km) and (28km, 42km) respectively.The value of e: is specified by 3°k, except for the bubble at (28km, 28km), where a stronger e: equals to 6°k is provided.The model is numerically integrated forward in time for 40 minutes, using leap-frog scheme with a time step M of 6 seconds.No perturbations are imposed to the initial pressure field.The dimensions of the domain, as well as the size and magnitude of the initial conditions are carefully specified so that the fields of the pressure and temperature fluctua tions remain well inside the domain throughout the entire simulation.Figures 1 and 2 show the simulated pressure and potential temperature perturbations on a vertical x-z cross section passing through the center of the domain after 35 minutes of simulation time, while their horizontal distributions at the height of 1 km are illustrated in Figs. 3 and 4. Note that in Figs. 1 and 3, the volume wide average has been removed.These four fields are selected as true solutions for verifying the retrieval results.
In additional to the visual comparisons, the accuracy of the retrieval results is investigated by evaluating several quantitative indices.If "A" represents either the potential temperature or the pressure, and the subscripts "t" or "r" denote the "true" or "retrieved' quantities respectively, we then obtain the following definitions: Spatial correlation coefficient (SCC) : Here, A can be an average taken either over a three-dimensional volume, or a two-dimen sional plane.N is the total number of grid points in the volume/plane.
Reliability parameter for retrievals ( E,) : fJJ (a:-F) + -a:-G + ---a: -a: + g 7/t-g � dxdydz The SCC score measures the "similarity" between two fields, and reaches 1 when they are identical.The purpose of calculating this index is that in many diagnostic studies of weather phenomena using Doppler radar data, it suffices to construct the "pattern" of the thermody namic variables.In other words, one is interested in the relative positions of the high/low pressure, as well as the warm/cold centers, rather than the exact numbers at each grid point.As for the reliability parameter E,, it estimates the relative errors of the retrieval, and is analo gous to the "momentum-checldng" parameter by Gal-Chen and Kropfli (1984).They gave a threshold of 0.5, above which the retrieval is believed to be virtually meaningless.Hane and Ray (1985) also suggested that an E, �0.25 would be considered to be a successful retrieval.

4.RESULTS
Most of the retrieval experiments were conducted within a smaller domain of 45 X 45 X 33 grid points in order to avoid the influence from the numerical model's boundary conditions.Only preliminarily investigations are executed in this study.Basically, these experiments are conducted to explore the feasibility of the proposed formulation, as well as its sensitivity to two commonly encountered degraded input data in radar meteorology, namely, the discrete nature of the Doppler radar volume scans, and the wind data with embedded observational errors.We will also focus on the vertical structure of the retrieved thermodynamic variables.

Reference Run
The first experiment represents a reference run in which the tendency terms in the input data sets F, G and H (i.e., du I at, dv I at, aw I at) are calculated accurately by finite differencing the model outputs of the velocity fields at T = 35 min-6 s and T = 35 min+ 6 s.Note that 6s is the length of the leap-frog time step in the model integration.For the advection, Coriolis f orce, and turbulence parameterizations, we use the model results at T = 35 min.The F, G and Hare then substituted into the retrieval algorithm without contamination.Note that even with perfect input data, one should not expect an exact solution from the proposed scheme.This is because the strategy as outlined in Section 2 only iteratively seeks an approximated solution, which will minimize the differences between the momentum equations and the input data.By contrast, the traditional methods, such as GC, directly solve f or the exact thermody namic variables from the model equations.With this in mind, Figs.5-8 display the retrieved pressure and potential temperature perturbations along the vertical and horizontal planes respectively.Their degree of agreement to the true solutions (Figs.1-4) is very encouraging.The author would like to emphasize that not only the horizontal distributions of the pressure and temperature are correct, the vertical structures can also be determined with certainty.
To make a comparison quantitatively, the indices introduced in ( 14)-( 15) are calculated globally over the entire retrieval domain.Table 1 shows the resulting statistics.The value of E, (=0.149) is far below 0.25, indicating a successful retrieval.Of particular satisfaction is the good sec scores (> 0.9) for all retrieved variables, which imply a close similarity between the retrieved and the true fields over a three-dimensional space.The SCC parameters are also estimated for the selected vertical and horizontal layers as displayed in Figs.5-8  For example, Figure 9 displays the resulting distributions of potential temperature perturba tion slice by slice along the y coordinate, and one immediately recognizes that the high scores assure that the vertical structures of thermodynamic perturbations can be determined without ambiguity.

Sensitivity to Time Intervals Between Radar Volume Scans
As discussed in GC, talcing the time derivatives over a long time interval is equivalent to simulating errors due to non-simultaneous radar observations.A criteria for an allowable radar scan time lit is suggested in Gal-Chen (1979) by: Ulit where U is a typical velocity, and L stands for the spatial resolution.In Experiment 2, the sensitivity of this retrieval scheme to the time intervals between each wind observation is explored.As mentioned earlier in the previous section, the tendency terms (i.e., (Ju I at, dv I at, aw I at) in the input data sets F, G, Hare prepared by finite differencing the model outputs of the velocity fields at T = 35 min-lit and T = 35 min+ lit, where /j.f is the length of the leap frog time step.For the advection, Coriolis force, and turbulence parameterizations, they are evaluated using the model-produced wind outputs at T = 35 min.In this experiment, we pro long the time intervals between two radar data sets gradually from 12s (Experiment 1), to 180s (Experiment 2), and to 300s (Experiment 3).Note that some research Doppler radar (e.g., Doppler on Wheel) is capable of conducting 60s to 180s fast scans, while 300s is approxi mately the time interval between WSR-880 volume scans under operational mode.
Table 1 lists the results of the quantitative comparisons.It is obvious that the proposed method works very well as /).. t increases from 12s to 300s.The magnitude of E, varies from 0.149 to 0.200, but still under the threshold value of 0.25.It may also be noted that in this experiment the SCC scores for both pressure and temperature remain close or above 0.90 for different time intervals, which implies, at least qualitatively, that the retrieved thermodynamic perturbation fields are robust for this kind of degraded data_ The vertical and horizontal struc ture of the retrieved pressure and temperature in Experiment 3 are illustrated in Figs.10-13, along with the calculated SCC scores for that particular slice of the domain.The agreements with their counterparts shown in Figs.1-4 are very satisfactory.Therefore, Experiments 2 and 3 state that this retrieval algorithm does have a certain tolerance to the errors associated with the discreteness of the radar measurements.In real case applications, one should realize that the importance of the time derivative term is case dependent, but the author believes that an intensive scanning strategy would be beneficial in terms of improving the quality of the retrievals.
In Experiment 4 the time derivatives are completely dropped out from the calculations of F, G, and H to test the validity of local steady-state assumption.From Table 1 it can be seen that the quality of the retrievals deteriorates severely, as is shown by the value of E, which exceeds the threshold, and becomes 1.049.Such a large number indicates that this retrieval is virtually a failure.Significant differences are also found in the plots of the retrieved fields.For example, Figure 14 denotes the potential temperature perturbations on the same x-z cross section as in Fig. 2, but the SCC score is found to be as low as 0.476.One should realize that the result found in this experiment is in agreement with that from GC, in which the commonly used local steady-state hypothesis was actually disqualified.In additional, GC also suggested that the meaning of steady-state should be interpreted statistically.In other words, it is the statistical properties, which can be obtained by some sort of spatial, temporal or ensemble average, that may appear time-independent.It is also known that the tendency terms Jul dt, ()v I dt, aw I dt comprise the effects of advection and evolution.For the former, it is possible to account for its influence by applying the concept of a moving frame of reference, as dis cussed by Gal-Chen (1982).However, the only way to mitigate errors due to evolution is to conduct fast radar volume scans so that the tendency terms can be estimated accurately.

Sensitivity to Contaminated Datasets
Experiment 5 is designed to study the impact of using imperfect data sets, which can be obtained by superimposing random errors of significant magnitudes onto F, G and H.In order to reduce the noise, filtering can be applied either on the input data (a priori) or on the result ing thermodynamic fields (a posteriori).The latter was advocated by GC when the artificial model-produced data were used for investigation, since he argued that filtering the input data would remove not only the noise, but useful signals as well.Therefore, he suggested applying the filtering only to the results, not the input.However, when GC was applied to a real bound ary layer case study (PHOENIX project 1978) in Gal-Chen and Kropfli (1984), a somewhat different conclusion was addressed that the smoothing should be applied not only before, but also after the retrieval was completed.This is called "double smoothing" in their paper.In this research, we employ this conclusion.
To prepare the input data sets, let A represent either the F, or G, or H field at each grid point; it is then superimposed by a random number as shown in the following: A(i,j, k) =A(i,j, k) + t: X(2xrdn -1) XA(i,j, k), (17) where "rdn" denotes a random number generator, whose range is from 0 to 1, and t: specifies the maximum magnitude of the perturbations to imitate the observational errors.Then, the perturbed field is smoothed by a three-dimensional filter which has the following form: A(i,j,k);:::; LL :Liwt(i+l,j+m,k+n)xA(i+l,j+m,k+n), 1=1 m=-1 n=-1 where wt is the weighting coefficient assigned to the points surrounding the center (i,j, k).For each grid point, a total of 27 points are involved in smoothing the variable.Let S = I / I + I m I +In I, then when S = 0 (the center point), wt= 1/ 8; when S = 1 (6 points), wt== 1/ 16; when S = 2 (12 points), wt= 1 1 32; S = 3 (8 points), wt== 1/ 64.In order to study the performance of this filter, the low-pass smoothing constraints represented by P5 and P6 in (6) are turned off in this experiment.
Let the maximum magnitude of the random error ( t:) be 25% of the value at each point, the statistics in Table 1 demonstrates that SCC scores remain high (0.895-0.9 60), while the value of E ( =0.188) indicates that useful information has been extracted through the retrievals.

Application to Topography
In a mountainous area such as the island of Taiwan, understanding how a weather system (e.g., typhoon, Mei-Yu front, squall line, etc.) evolves and interacts with the topography is always an important issue.To accomplish this goal, it turns out to be a necessary task to conduct thermodynamic retrievals for a non-flat region.Instead of reformulating the govern ing e,quations on a terrain-following coordinate system, which is usually the case for doing numerical model simulations, the method proposed in this research can easily adopt the topog raphy by setting the gradients of the cost function (J) with respect to the unknowns ( n' and 8 ') to zeros at these points where the altitudes are below the topography.Experiment 6 dem onstrates an example in which the retrieval scheme is utilized in a domain with an idealized bell shaped .mountain.The peak ofthe mountain reaches 900m, and the half-widths is 10 km.

+22 +22
Figures 19-20 depict the retrieved vertical pressure and potential temperature distributions for a x-z cross-section.It can be seen that the positions of the high/low pressure centers, and the structures of the triple thermal bubbles, are all successfully recovered above the terrain.

CONCLUSIONS AND FUTURE WORK
A different thermodynamic retrieval scheme is developed in this research to extract the pressure and potential temperature field information from 3-D wind observations.The latter, in real cases, can be derived by Doppler radar.Based on the aforementioned results, the fol lowing is a list of the conclusions: • Numerical experiments have shown that the algorithm outlined in this research is a feasible tool to perform thermodynamic retrievals with sufficient accuracy.
• The proposed method can tolerate the errors introduced by using finite differences to ap proximate the local temporal derivatives, and extract meaningful results from seriously con taminated input data sets, provided an appropriate filter is applied.However, completely neglecting the time rate of change of the wind fields leads to unsuccessful retrievals.
• The retrieval scheme can easily be applied to a region with topography.
• Only 3-D wind data are needed for the retrievals, and the products of the proposed scheme are the three-dimensional potential temperature fluctuations, and the pressure perturbation gradients, in any direction.Therefore, the interpretation of the thermodynamic structures in the vertical direction should no longer cause confusiion.We believe that this is one advan tage over the traditional scheme of GC, and should be of particular usefulness in many meteorological applications.
The current data sets focus on simulated data sets, and only preliminary tests are conducted.Future work, some already underway, should include more investigations using Doppler radar observations obtained from real field experiments, and extend the method to incorporate the thermodynamic equation.It is believed that this equation can provide additional information to constrain the gradients of the potential temperature in a three dimensional space.

[Fig. 2 .
Fig. I.The model-generated "true" pressure field on a vertical cross section crossing the center of the domain.The volume wide mean has been removed.
Fig. 9.The 2-D SCC scores of potential temperature perturbations ( e') on each x-z cross section along the y coordinate for Experiment 1.