Thursday, December 5, 2013

Cost Functional


The proper specification of a cost functional, which is typically defined as a weighted sum of lack of fits to observations and other physical constraints, is essential for obtaining a satisfactory solution of an objective analysis problem (Legler and Navon 1991).

The data that I am typically dealing with in my research is derived from scanning Doppler radars, e.g., measurements from X-band (95-GHz) and C-band (5.6-GHz) radars. These radars scan the atmosphere in the so-called plan position indicator (PPI) mode, which is best described by the spherical coordinate system as each gate location can be represented by the coordinate-triplet \((r, \phi, \theta)\). This data is usually then interpolated onto a more natural and intuitive 3-D Cartesian analysis domain. From here, numerous geophysical parameters can be estimated, e.g., the 3-D wind field \((\mathbf{u}, \mathbf{v}, \mathbf{w})\). And this is where a cost functional is often needed.

The Cost Functional

A cost functional \(J\) is typically defined as the sum of independent constraints, one of these relying on observations from one or more instruments, and more often than not some sort of constraint on the smoothness of the analysis. All of these individual constraints represent so-called lack of fits between observations or a priori assumptions. We arrive at our best estimate of the truth when the global minimum of \(J\) has been found. So for example, if we're attempting to estimate the 3-D wind field, we're looking for the simultaneous satisfaction of

\frac{\partial J}{\partial \mathbf{u}} = 0, \; \frac{\partial J}{\partial \mathbf{v}} = 0, \; \frac{\partial J}{\partial \mathbf{w}} = 0 .

As a result, we are expected to compute the gradient of the cost functional at one point or another, which, for some constraints is easily defined analytically, but for others (e.g., smoothness) not so much. So let's focus on the smoothness constraint.

The Smoothness Constraint

The smoothness constraint can be defined multiple different ways, e.g., using first-order or second-order partial derivatives. For the sake of this post, let's focus on the smoothness constraint defined in Potvin et al. (2012), in particular only the first term defined in equation (5). This term looks something like the following

J_s = \frac{1}{2} \sum_{i,j,k} \lambda_s
             \left[\left(\frac{\Delta^2 u}{\Delta x^2}\right)^2 +
                    \left(\frac{\Delta^2 u}{\Delta y^2}\right)^2 +
                    \left(\frac{\Delta^2 v}{\Delta x^2}\right)^2 +
                    \left(\frac{\Delta^2 v}{\Delta y^2}\right)^2    

The summation is over all the Cartesian grid points in the analysis domain, and \(\lambda_s\) represents the weight of this term on the entire analysis.

For the sake of simplicity we'll assume \(u\) is only a function of \(x\), i.e. \(u_i = u(x_i)\). Therefore the second order partial derivative can be approximated by a finite difference,

\left( \frac{\partial^2 u}{\partial x^2} \right)_i^2 = \left( \frac{u_{i+1} - 2 u_i + u_{i-1}}{\Delta x^2} + \mathcal{O}( \Delta x^2) \right)^2

Equation \(\ref{Eq:2}\) becomes an approximation when the \(\mathcal{O}(\Delta x^2)\) term is dropped. This is referred to as a centered difference scheme with second-order accuracy. Note how a scheme like this can only be used for interior grid points since it requires information from left and right grid points. Therefore forward or backward difference schemes must be employed at the grid boundaries.

Now, since \(J_s\) is the summation of all these partial derivatives over all the Cartesian grid points, it is actually a function of several terms involving \(u_i\). To illustrate this, consider the following,

\( \left( \frac{\partial^2 u}{\partial x^2} \right)_{i-2}^2 \approx \left( \frac{u_{i-1} - 2 u_{i-2} + u_{i-3}}{\Delta x^2} \right)^2 \)

\( \left( \frac{\partial^2 u}{\partial x^2} \right)_{i-1}^2 \approx \left( \frac{u_{i} - 2 u_{i-1} + u_{i-2}}{\Delta x^2} \right)^2 \)

\( \left( \frac{\partial^2 u}{\partial x^2} \right)_{i}^2 \approx \left( \frac{u_{i+1} - 2 u_{i} + u_{i-1}}{\Delta x^2} \right)^2 \)

\( \left( \frac{\partial^2 u}{\partial x^2} \right)_{i+1}^2 \approx \left( \frac{u_{i+2} - 2 u_{i+1} + u_{i}}{\Delta x^2} \right)^2 \)

\( \left( \frac{\partial^2 u}{\partial x^2} \right)_{i+2}^2 \approx \left( \frac{u_{i+3} - 2 u_{i+2} + u_{i+1}}{\Delta x^2} \right)^2 \)

Notice how \(u_i\) manifests itself in three of the above equations. However, from the finite difference scheme we chose to approximate the partial derivative, we see that the influence of the \(u_i\) term vanishes when we move past the \(i - 2\) and \(i + 2\) derivatives. Therefore, in general for a second-order-accurate centered difference scheme, the partial derivative of \(J_s\) with respect to \(u\) at the \(i\)-th grid point is,

\left( \frac{\partial J_s}{\partial u} \right)_i^2 = \frac{\lambda_s}{\Delta x^2} \left[ \left(\frac{\partial^2 u}{\partial x^2}\right)_{i-1} - 2 \left(\frac{\partial^2 u}{\partial x^2}\right)_i + \left(\frac{\partial^2 u}{\partial x^2}\right)_{i+1} \right]

As stated before, equation \(\ref{Eq:3}\) was derived using the second-order-accurate centered difference scheme we chose at the beginning. The chain rule has been used here to calculate the derivatives. Fundamentally, equation \(\ref{Eq:3}\) is highly dependent on the finite difference scheme we choose to approximate the partial derivatives. As a result we need to keep track of which finite difference schemes are being utilized. For example, for grid points at the boundaries of the grid, say at \(i = 1\), we now have to use a forward difference scheme because we only have information of \(u\) in the forward direction. Notice that now we would have,

\( \left( \frac{\partial^2 u}{\partial x^2} \right)_{i=1}^2 \approx \left( \frac{u_{1} - 2 u_{2} + u_{3}}{\Delta x^2} \right)^2 \)

\( \left( \frac{\partial^2 u}{\partial x^2} \right)_{i=2}^2 \approx \left( \frac{u_{3} - 2 u_{2} + u_{1}}{\Delta x^2} \right)^2 \)

\( \left( \frac{\partial^2 u}{\partial x^2} \right)_{i=3}^2 \approx \left( \frac{u_{4} - 2 u_{3} + u_{2}}{\Delta x^2} \right)^2 \)

where a forward difference scheme of \(\mathcal{O}(\Delta x)\) is used at \(i = 1\) and a centered difference scheme of \(\mathcal{O}(\Delta x^2)\) has been used for all subsequent derivatives. Notice how \(u_1\) manifests itself in two of the above equations, where its influence vanishes once we get past the second grid point \((i \gt 2)\). So now we would have a slightly different equation than that of equation \(\ref{Eq:3}\) to describe the derivative of \(J_s\) with respect to \(u\) at the boundary \(i = 1\),

\left( \frac{\partial J_s}{\partial u} \right)_{i=1}^2 = \frac{\lambda_s}{\Delta x^2} \left[ \left(\frac{\partial^2 u}{\partial x^2}\right)_{i=1} + \left(\frac{\partial^2 u}{\partial x^2}\right)_{i=2} \right]


Once again, please keep in mind that equations \(\ref{Eq:3}\) and \(\ref{Eq:4}\) have been produced using the particular finite difference schemes that we chose to approximate the partial derivatives. The take home message is that if a constraint cannot be expressed analytically, like the smoothness constraint which has to be approximated by finite differences, then in turn the derivative of that constraint is highly dependent on those finite difference schemes.


Legler, D. M. and I. M. Navon, 1991: VARIATM - A Fortran program for objective analysis of pseudostress wind fields using large-scale conjugate-gradient minimization. Computers & Geosciences, 17, pp. 1-21.

Friday, August 2, 2013

Steiner Algorithm for Echo Classification


I have been using the Steiner et. al (1995) algorithm for echo classification for some time now. It has its drawbacks, like any algorithm. For example, there are many free parameters that are region dependent, like the automatic convection parameter and the peakedness parameter (\(\Delta Z\)). These parameters, as one might suspect, would have entirely different values for convection near Darwin, Australia than convection near Oklahoma City, Oklahoma.

This algorithm was designed and tuned for radar data collected in Darwin, Australia during February 1988. Therefore, if I want to apply this algorithm to data collected at ARM's Southern Great Plains (SGP) site near Lamont, Oklahoma, I need to first check the results using the default parameters suggested by Steiner et al. (1995), and then tune these parameters accordingly if need be.

It should be noted that this algorithm is applied to gridded reflectivity data at a height chosen to be below the melting layer. In the case of Steiner et al. (1995), they used grids 240 x 240 km in \(x\) and \(y\), with 2 km resolution in each dimension. Here I will be using 100 x 100 km grids in \(x\) and \(y\), with 500 m resolution in each dimension.


The algorithm is based on three criteria: intensity, peakedness, and surrounding area. The intensity criterion is a simple threshold, where any grid point with a reflectivity of at least 40 dBZ is automatically convective. The peakedness criterion (\(\Delta Z\)) checks the difference between the grid point and the average reflectivity taken over the surrounding background (\(Z_{bg}\)). If this difference is greater than a certain value, then the grid point is labelled convective. Finally, any grid points within an intensity-dependent radius (\(R_{sa}\)) from those already labelled as convective by either the intensity or the peakedness criteria are also labelled convective.

The default parameters defined in Steiner et al. (1995) are as follows:

  1. Automatically convective: 40 dBZ.
  2. Background radius: 11 km.
  3.  \(\Delta Z\) = \(\begin{cases} 10, & \mbox{if } Z_{bg} \lt 0 \newline 10 - \frac{Z_{bg}^2}{180}, & \mbox{if } 0 \le Z_{bg} \lt 42.43 \newline 0, & \mbox{if } Z_{bg} \ge 42.43 \end{cases}\)
  4. There are three different relations for \(R_{sa}\), the small, medium, and large relations, but each has 1 km and 5 km as the smallest and largest radii, respectively.


So let's first check how the algorithm, with its default parameters, classifies a squall line event during May 20th, 2011. Near 1033 UTC the squall line, with an axis slightly offset from a north-south direction, was located approximately in the middle of the analysis domain, as it moved through in a west-east direction. There was a leading anvil region to the east and a trailing stratiform region to the west. The results are shown in Fig. 1.

Fig. 1. C-SAPR reflectivity and corresponding echo classification at 1033 UTC on May 20th, 2011. The working level (or separation altitude) was a constant 1.5 km. The medium \(R_{sa}\) relation was used for this classification.

The algorithm does a decent job at classifying the convective squall line feature, however it does miss its leading edge to the east. The trailing stratiform region to the west is another story though. The algorithm is severely affected by brightband contamination here, and a large swath of this region is labelled convective.

One obvious reason for some of the poor classification in the trailing stratiform region is due to the intensity criterion, which was set to 40 dBZ. This value is too weak for convection observed in Oklahoma, which consistently produces reflectivities of 55+ dBZ. Therefore, my first change will be to up this value to 45 dBZ. Another, albeit less obvious reason for the poor classification in the stratiform region is due to the peakedness criterion. Examining the peakedness relation, which varies depending on the intensity of the background reflectivity \(Z_{bg}\), I see that it does not require a sharp enough gradient between the grid point and \(Z_{bg}\), and as a result, the brightband contamination, which has a high intensity but at the same time is spatially smooth, gets labelled convective. Therefore, I will use a new relation for \(\Delta Z\):

$$ \Delta Z = \begin{cases} 14, & \mbox{if } Z_{bg} \lt 0 \newline 14 - \frac{Z_{bg}^2}{180}, & \mbox{if } 0 \le Z_{bg} \lt 42.43 \newline 4, & \mbox{if } Z_{bg} \ge 42.43 \end{cases} $$

Basically I have shifted the \(\Delta Z\) curve up by 4 dBZ, or, in other words, a sharper gradient between the grid point and the background reflectivity is required. The results of these two changes are shown in Fig. 2.

Fig. 2. Same as Fig. 1 but tuning the intensity criterion and the peakedness criterion.


Steiner, M., and R. A. Houze Jr., and S. E. Yuter, 1995: Climatological Characterization of Three-Dimensional Storm Structure from Operational Radar and Rain Gauge Data. J. Appl. Meteor.34, 1978-2007

Saturday, July 27, 2013

Biggerstaff Algorithm for Echo Classification


I recently took a look at the Biggerstaff and Listemaa (2000) algorithm for echo classification, which is supposed to be an improved scheme from the Steiner et al. (1995) algorithm. To examine this algorithm, I use reflectivity observations from a scanning ARM C-band radar (C-SAPR), and apply the algorithm to this data. All data used in this post was recorded during the Midlatitude Continental Convective Clouds Experiment, or MC3E for short. Furthermore, I will be using gridded (or interpolated) C-SAPR data, on a Cartesian grid 100 x 100 x 17 km in \(x\), \(y\), and \(z\), with a resolution of 500 m for each dimension.


This modified algorithm focuses on two estimates: the magnitude of horizontal reflectivity gradient \(\left(\left|\left|\nabla_h Z_e\right|\right|\right)\) at the working level, and the vertical lapse rate of reflectivity \(\left(\nabla_z Z_e\right)\). The vertical lapse rate of reflectivity is defined as the decrease in reflectivity with increasing height in a 3 km layer directly above the height of maximum reflectivity.

These two estimates are put through a thresholding, as described below:

  • \(\left|\left|\nabla_h Z_e\right|\right| \ge\) 3 dB km-1 \(\rightarrow\) grid point possibly convective
  • \(\left|\nabla_z Z_e\right| \le\) 3.5 dB km-1 \(\rightarrow\) grid point possibly convective

I say possibly in the above criteria because there are some other criteria involved in this algorithm (e.g. brightband fraction and a 2-D window filter). These other criteria are ignored in this post.

Since the vertical lapse rate is computed for a 3 km layer above the height of maximum reflectivity, I used the following equation for each column,

$$\nabla_z Z_e \approx \frac{Z_e|_{\text{max}} - Z_e|_{+3 \; \text{km}}}{3 \; \text{km}}$$.

Note with the above equation that we should expect a positive vertical reflectivity lapse rate to be computed, without taking the absolute value.

Finally, note that finite differences used in the gradient calculation were that of 1st and 2nd-order accurate centered, forward, and backward differences.


To test this algorithm, I want to look at the results of both the magnitude of horizontal reflectivity gradient and the vertical lapse rate of reflectivity, and see if these estimated values reflect the expected values described in Biggerstaff and Listemaa (2000). In order to do this, I picked an observation time during the passage of a squall line which had well-defined anvil and stratiform regions.

Fig. 1. C-SAPR reflectivity at the working level (1.5 km) during 1040 UTC on May 20th 2011. Reflectivity values below 0 dBZ were masked.

The cross section of reflectivity shown in Fig. 1 shows a north-south squall line located in the middle of the analysis domain, which was moving west-east. There is an anvil region ahead of the squall line, as well as a trailing stratiform region. Since Fig. 1. is a cross section at 1.5 km altitude, the anvil region is largely not seen.

First I will show the results of the lapse rate calculation.

Fig. 2. Vertical reflectivity lapse rate within the 3 km layer above the height of maximum reflectivity for the same time as shown in Fig. 1.

Notice that the upper limit of the color bar in Fig. 2 is 3.5 dB km-1. From the criteria shown above, we expect all areas shown in red to possibly be stratiform, and the rest possibly convective. The minimum and maximum values computed were 0.02 and 30.3 dB km-1, respectively. The missing data seen around 97.5 W and 36.8 N is due to the cone of silence of the C-SAPR. Also, since the lapse rate was computed for each column, the anvil region has filled in because there is data above the working level.

Now let's look at the results of the gradient calculation.

Fig. 3. Magnitude of horizontal reflectivity gradient at the working level for the same time as shown in Fig. 1.

Here I have set the upper limit of the color bar to highlight areas of possible convection. Any areas shaded in red are possibly convection according to Biggerstaff and Listemaa (2000). The minimum and maximum values computed were 0.001 and 14.1 dB km-1, respectively. Since this calculation is at the working level, the anvil region is once again not available.


A first look at the results shows that most areas will not be well classified using these two criteria alone. By any standard the current implementation of the lapse rate estimation does a very poor job at highlighting both convective and stratiform grid points. The gradient criterion does not appear to fair much better, and is severely affected by poor attenuation correction (the ray-like features in Fig. 3). The gradient criterion does do a good job at classifying the leading edge of the squall line as convective, since the precipitation mass in this region would likely be a direct result of convective updrafts.

Further analysis and discussion with colleagues is in order here!


Biggerstaff, M. I., and S. A. Listemaa, 2000: An Improved Scheme for Convective/Stratiform Echo Classification Using Radar Reflectivity. J. Appl. Meteor., 39, 2129-2150

Steiner, M., and R. A. Houze Jr., and S. E. Yuter, 1995: Climatological Characterization of Three-Dimensional Storm Structure from Operational Radar and Rain Gauge Data. J. Appl. Meteor., 34, 1978-2007