HTML
--> --> -->As a new method of atmospheric remote sensing, the GPS radio occultation (RO) technique has the advantages of high precision, high vertical resolution, global quasi-uniform coverage, all-weather conditions, and long-term stability, among others, and has received extensive attention. In particular, the occultation project named the Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC), jointly implemented by Taiwan and the United States, was able to obtain nearly 2000 profiles a day during the five years of its mission (Liou et al., 2007). These occultation profiles have important applications in the study of global climate change and numerical weather forecasting.
In view of the high vertical resolution characteristics of occultation data, many scholars use them to determine ABL heights. (von Engeln et al., 2005) proposed the determination of ABL height by using the truncation height in full spectrum inversion. (Sokolovskiy et al., 2006) estimated ABL depths (i.e., ABL top heights) by determining the maximum lapse (breakpoint) of the refractivity vertical gradient. Subsequently, in a follow-up study (Sokolovskiy et al., 2007), they proposed that the ABL depth can also be determined by the bend angle profile. (Basha and Ratnam, 2009) used the high vertical resolution radiosonde data from the tropical station Gandaki to determine the ABL height by calculating the minimum value of the refractivity vertical gradient, and compared their results with those determined using COSMIC RO data, showing good consistency. (Guo et al., 2011) applied the breakpoint method proposed by (Sokolovskiy et al., 2006) to the refractivity data from COSMIC RO, and obtained the global distribution and seasonal variation of the marine ABL height. (Ao et al., 2012) defined the minimum vertical gradient of the profiles of refractivity and vapor pressure as the ABL top, and introduced the relative minimum gradient (or sharpness parameter) to characterize the quality of the ABL height. (Chan and Wood, 2013) improved the breakpoint method proposed by (Sokolovskiy et al., 2006) through fulfilling more constraints to ensure the quality of the results, and presented the seasonal cycle characteristics of the global ABL height. (Liao et al., 2015) and (Liu et al., 2016) also used a similar method to obtain and evaluate the seasonal and diurnal variations of the global marine ABL height.
The work mentioned above mainly used refractivity profile data, with bending angle (BA) profile data rarely employed. However, the BA, as an intermediate product in the occultation data processing chain, has some unique advantages. For example, the BA is a rawer product than refractivity, and thus contains less noise in data processing. Furthermore, because the BA is calculated directly from the observational data (and not by the Abel transform), there is no limitation on atmospheric ducting. In addition, since refractivity is obtained through integrating the BA according to the Abel inversion formula, refractivity profiles are smoother than BA profiles, which means that BA profiles contain more information than refractivity profiles (Fig. 1). This was also pointed out by (Rieder and Kirchengast, 2001), who showed in their Fig. 1 that the BA profile exhibited the highest signal-to-noise ratio. Therefore, BA profile data are more suitable for detecting ABL heights.
Figure1. COSMIC RO refractivity and BA profiles (atmPrf_C001.2008. 275.03.49.G31_2013.3520_nc).On the other hand, most of the methods mentioned above for determining ABL height are reduced to calculating the vertical derivatives of refractivity or BA profiles and carried out usually by the finite difference method. This can be addressed in the framework of inverse problems. From the viewpoint of the theory of inverse problems, calculating derivatives of functions using limited and discrete observational data (often known as numerical differentiation), belongs to inverse problems. Plus, numerical differentiation is generally ill-posed as an inverse problem, which is characterized mainly by numerical instability (Tikhonov and Arsenin, 1977). Hence, if conventional methods, such as the finite difference method, are applied to solving numerical differentiation problems, noise in the observational data will be amplified, resulting in derivatives far from true values that are sometimes even completely useless (Hanke and Scherzer, 2001; Ramm and Smirnova, 2001). To overcome the ill-posedness (mainly numerical instability) in numerical differentiation, some regularization strategies, such as Tikhonov regularization, have been introduced, and some stable numerical differentiation methods, such as the cubic spline interpolation method, mollification method, and variational regularization method, have been developed to obtain stable approximate derivatives (Ramm and Smirnova, 2001; Cheng et al., 2003).
In the above context, the present paper aims at determining ABL heights from BA profiles. Unlike in previous work, the present paper adopts another approach, i.e., combining the finite difference algorithm with the Tikhonov regularization technique, to obtain the vertical derivatives of BA profiles and then determine the ABL heights.
This paper is structured as follows: Section 2 introduces the numerical differential method and data used in this paper. Section 3 is the validation of the method by comparing with the ABL heights from COSMIC products, and results obtained by the finite difference method upon refractivity. Finally, a summary and discussion are presented in section 4.
2.1. Methods
The ABL top is a transition layer from the ABL to the free atmosphere, which is usually accompanied by temperature inversion and/or a sharp decrease in water vapor (Sokolovskiy et al., 2007, Fig. 1). According to the Smith-Weintraub equation (Smith and Weintraub, 1953), \begin{equation} \label{eq1} N=77.6\frac{P}{T}+3.73\times 10^5\frac{e}{T^2} , \ \ (1)\end{equation} where T is the temperature (K), P is the total atmospheric pressure (hPa), and e is the water vapor pressure (hPa), the refractivity will drop significantly at the ABL top and, correspondingly, there will also be a steep reduction in the BA (α). However, (Sokolovskiy et al., 2007) did not directly calculate the minimum vertical gradient (negative) of the BA; instead, they worked out its maximum lapse (positive) in a height interval ? z in the form ?α/(?α/? z) (i.e., max{?α/(?α/? z)}) to determine the ABL height, as small-scale structures in BA profiles may produce spurious gradient minima. The method of (Sokolovskiy et al., 2007) essentially calculates the vertical derivative of the BA via the finite difference method, but with the BA smoothed. Because the Tikhonov regularization technique has the function of smoothing, we apply it directly to the BA profile data, obtaining the vertical gradient (first derivative) of the BA. The height of the minimum vertical gradient is determined as the ABL height. The formulation of this problem is given as follows:Suppose α(z) to be the BA profile, z b≤ z≤ z t, where z b and z t represent the bottom and top of the α(z) profile, respectively. Grid points {zi,i=1,?,m} divide [z b,z t] into m-1 parts (z b=z1<z2<?<zm=z t) with an equal interval of h=(z t-z b)/(m-1). The values of α(z) at zi, α(zi), are known by the COSMIC RO data, and the next goal is to obtain the first derivative of α(z) at zi, α'(zi).
For convenience, some notation is introduced as follows: $\alpha'(z)\equiv\varphi(z)$, $\alpha(z_i)\equiv\alpha_i$, $\varphi(z_i)\equiv\varphi_i$. By the Newton-Leibniz formula, there is \begin{equation} \label{eq2} \alpha_{i+1}=\alpha_{i-1}+\int_{z_{i-1}}^{z_{i+1}}{\varphi(z)dz},\quad i=2,\cdots,m-1 \ \ (2)\end{equation}and the integral of the right-hand side of Eq. (3) can be approximately obtained by the Simpson formula, \begin{equation} \label{eq3} \alpha_{i+1}-\alpha_{i-1}\approx \frac{h}{3}(\varphi_{i-1}+4\varphi_i+\varphi_{i+1}),\quad i=2,\cdots,m-1 \ \ (3)\end{equation} which can be written in matrix form as \begin{equation} \label{eq4} {AX}={B} ,\ \ (4)\end{equation} where \begin{eqnarray*} {A}&\!\!=\!\!&\left( \begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c} 1 & 4 & 1 & & \\[-0.5mm] & \ddots & \ddots & \ddots & \\[-0.5mm] & & 1 & 4 & 1 \end{array} \right){X}=(\varphi_1,\varphi_2,\cdots,\varphi_{m-1},\varphi_m)^{\rm T}\\ {B}&\!\!=\!\!&\bigg(\,\!\dfrac{3}{h}(\alpha_3\!-\!\alpha_1),\!\dfrac{3}{h}(\alpha_4\!-\!\alpha_2),\!{\cdots}\,\!,\!\dfrac{3}{h}(\alpha_{m-1}\!-\!\alpha_{m-3}),\,\! \dfrac{3}{h}(\alpha_m\!-\!\alpha_{m-2})\!\bigg)^{\!\rm T}. \end{eqnarray*} Now, the problem is down to solving these linear equations.
The linear equations [Eq. (5)] can also be obtained from the finite difference method, and so solving them involves simply calculating the vertical derivative of a profile of the BA by the finite difference method. Since it is an inverse problem, solving these linear equations is ill-posed, and especially numerical unstable for observational bending angle. To overcome the ill-posedness, Eq. (5) is transformed into solving the minimization problem of the Tikhonov functional with a regularization functional included, \begin{equation} \label{eq5} \min J,\quad J=\|{AX}-{B}\|^2+\gamma\|{LX}\|^2 \ \ (5) \end{equation} where γ||LX||2 is the regularization functional and γ>0 is the regularization parameter that should be given in advance. The matrix $L=\left(\begin{array} -1 & 1 & \quad & \quad & \quad \\ \quad & -1 & 1 \quad & \\ \quad & \quad & \ddots & \ddots \\ \quad & \quad & \quad & -1 & 1 \end{array} \right)$ is the first derivative operator that can control the smoothness of the solution x via choosing the correct value of the regularization parameter γ. It is simple to find that, as γ increases, ||LX||2 will decrease and then X will be smoother.
Equation (6) is reduced to solving the linear equations (well posed) as \begin{equation} \label{eq6} ({A}^{\rm T}{A}+\gamma {L}^{\rm T}{L}){X}={A}^{\rm T}{B} \ \ (6)\end{equation} and its solution, \begin{equation} \label{eq7} {x}_0=({A}^{\rm T}{A}+\gamma {L}^{\rm T}{L})^{-1}{A}^{\rm T}{B} \ \ (7)\end{equation} is an optimal stable approximate solution of Eq. (5), where the regularization parameter γ is determined by the L-curve method (Hansen, 1992), which balances the two items in Eq. (6). Also, γ ranges from tens to hundreds according to the smoothness of the first derivative of the BA. The position of the minimum value of the BA vertical gradient is defined to be the ABL height.
To ensure the reliability of estimated ABL heights for refractivity data, (Ao et al., 2012) introduced the sharpness parameter, and (Chan and Wood, 2013) introduced the relative distinctness of minima, which is defined as the ratio of the global minimum in a refractivity gradient to the mean of all local minima in the same refractivity gradient. The present paper also defines a sharpness parameter to measure the reliability of estimated ABL heights for BA data, which is given as \begin{equation} \label{eq8} \lambda=\dfrac{\alpha'_{{\rm global}}}{\mathop{{\rm aver}}\limits_{i\le 5}(\alpha'_{\min-i})} \ \ (8)\end{equation} where α' global is the global minimum (negative) of a BA vertical gradient, ${\mathop{{\rm aver}}\limits_{i\le 5}(\alpha'_{\min-i})}$ is the average of the first five minimum values (negative) of the BA vertical gradient (with the global minimum included; and if less than five, it is calculated according to the actual number of minima of the BA vertical gradient). The average in Eq. (9) takes only the first five minor minima of the BA vertical gradient, as the minima in different occultation profiles may differ greatly in number. It is obvious that Λ≥ 1 and, generally, as Λ increases, the ABL top will be more significant and the ABL height more reliable correspondingly.
2
2.2. Data
The occultation data used in this paper are from CDAAC (the COSMIC Data Analysis and Archive Center; http://www.cosmic.ucar.edu/), including four months of global BA and refractivity profiles (cosmic2013 atmPrf files) in January, April, July and October of 2008. The corresponding ECMWF high-resolution gridded analysis data (echPrf files) are also used for the validation of the method. Since the ABL height is usually below 6 km, we only deal with profiles below 6 km, and require the difference between the minimum height of a profile and the topographic height at the occultation point to be less than 500 m. The topographic information is from the GTOPO30 global digital elevation model (https://lta.cr.usgs.gov/GTOPO30). Accordingly, the number of eligible profiles is 11 263, 12 489, 12 066 and 11 748, in the four months respectively (47 566 in total). To facilitate data processing and calculation, we interpolate the profiles to equidistant grid points with an interval of h=10 m.2
3.1. Comparison with ABL heights from COSMIC products
3.1.1. Case of BAFirst, we compare the ABL height Z BA_r with Z BA_c under the constraint Λ>1.75. From Figs. 2a-d, the average bias between Z BA_r and Z BA_c in January, April, July and October is -0.04 to -0.02 km (Z BA_r is slightly lower than Z BA_c), and the correlation coefficient is about 0.98. Though the two are in good consistency, there still exist some cases with large differences. As shown in Figs. 2e-f, with reference to the red point in Fig. 2a, both an inversion and a decrease in vapor exist at the height of 0.72 km, and the Z BA_r seems more acceptable. This example represents the majority of cases with large differences in the upper-left area of the line y=x. That is, Z BA_c sometimes fails to detect ABL heights near the endpoint of the profile. Essentially, the fixed 0.3 km sliding window for sharp ABL tops, aiming at avoiding small-scale structures (Sokolovskiy et al., 2007), might neglect any thin ABL tops of less than 0.3 km. From Figs. 2g and h, with reference to the cyan point in the lower-right area of Fig. 2a, this represents another kind of situation where there are multiple sharp peaks in the BA gradient (Fig. 2g). From Fig. 2h, an obvious ABL top cannot be found according to the temperature and vapor pressure profiles and different methods can give out different results.
Figure2. (a-d) Scatter diagrams of Z BA_r and Z BA_c for Λ≥ 1.75 in January, April, July and October. The brown line is y=x. (e-f) COSMIC RO BA profile (atmPrf_C006.2008.027.05.43.G13_2013.3520_nc, with reference to the red point in (a), and the temperature and vapor pressure profiles from its corresponding echPrf file (echPrf_C006.2008.027.05.43. G13_2013.3520_nc). The BA gradient is solved by the method in the present paper (the following are the same). (g-i) COSMIC RO BA and refractivity profiles from the atmPrf file (atmPrf_C005.2008.015. 04.27.G02_2013.3520_nc, with reference to the cyan point in Figs. 3a, 5a and 6a) and profiles from its corresponding echPrf file. The refractivity gradient is solved by the finite difference method (the following are the same).When the constraint is weakened——say, the sharpness parameter Λ≥ 1.5 (Fig. 3a)——the correlation coefficient is still over 0.95, and there are more eligible occultation profiles (46.2% of profiles in January can be used in the comparison) without a larger bias between the two results. When the constraint is strengthened——say, the sharpness parameterΛ≥ 2 (Fig. 3b)——the correlation coefficient is up to 0.99, and there is less bias, indicating the best consistency.
Figure3. Scatter diagrams of Z BA_r and Z BA_c in January for (a) Λ>1.5 and (b) Λ>2. The brown line is y=x.3.1.2. Case of refractivity
Next, we compare the ABL height Z BA_r with Z Re_c for Λ>1.75. From Figs. 4a and b, the average bias between Z BA_r and Z Re_c is -0.260 to -0.220 km, and the correlation coefficient is higher than 0.85. That is, Z BA_r is less than Z Re_c systematically, but there is still a high correlation between Z BA_r and Z Re_c. In addition to the points of good correlation, there are also some points with large deviations.
Taking the red point in Fig. 4a as an example, ZBA_r=0.45 might be more reasonable according to the temperature and vapor pressure profiles from Fig. 4f. The same line, y=x, divides the diagrams into two areas, as in section 3.1.1. The example in Figs. 5e-g can also represent the cases with large differences in the upper-left area. As the sliding window mentioned above, (Guo et al., 2011) also use it to avoid small-scale irregularities in refractivity, resulting in the negligence of the ABL tops near the surface. For another reason, the refractivity profiles retrieved from the BA profiles by the Abel inverse transform will be smoothed during the integrating process, i.e., the gradient patterns of the BA profiles will not be retained, which yields the differences. As for the lower-right area in Fig. 4a, the magenta point and the cyan point illustrate the large difference (4.12 km in Fig. 2h) and the relatively smaller difference (1.67 km in Fig. 4i), respectively. These cases are usually excluded because the refractivity gradients do not meet criterion b in (Chan and Wood, 2013). It is also hard to determine these non-sharp ABL tops according to the temperature or the vapor pressure profiles.
Figure4. (a-d) Scatter diagrams of ZBA_r and ZRe_c for Λ≥ 1.75 in January, April, July and October. Panels (e-g), with reference to the red point in (a), are from atmPrf & echPrf_C004.2008.013.05.48.G27_2013.3520_nc. Panels (h-j), with reference to the magenta point in (a), are from atmPrf & echPrf_C003.2008.018.18.56.G11_2013.3520_nc. The cyan point in (a) denotes the case in Fig. 2 (g)-(i). The brown line is y=x.2
3.2. Comparison with ABL heights from refractivity by the finite difference method
In research on determining ABL heights by the gradient method, the finite difference method is widely utilized (Basha and Ratnam, 2009; Seidel et al., 2010; Ao et al., 2012; Chan and Wood, 2013). Therefore, the ABL height obtained in this paper is also compared with that obtained by applying the finite difference method to the refractivity profile data.From Fig. 5 (for Λ≥ 1.75), the average bias between Z BA_r and Z Re_fd is 0.071-0.082 km, and the correlation coefficient is higher than 0.95. That is, Z BA_r is also slightly higher than Z Re_fd, and there is still a high correlation between Z BA_r and Z Re_fd. It seems that the finite difference method upon refractivity profiles is able to detect lower ABL heights compared with the two kinds of COSMIC products. However, cases of large deviation still exist. As shown in Figs. 5e-g, with reference to the red point in Fig. 5a, this case is located at an inland area where there is less moisture content (60.0°N, 66.5°E). The temperature inversion at about 0.87 km may play a more important part in Eq. (2) than the vapor pressure at that height, and the Z BA_r might be more reliable to some extent. Thus, apart from detecting lower ABL heights near the bottom of the profile, compared with methods upon refractivity profiles, this new method can also determine the thin ABLs, as some gradient patterns in the BA profiles are smoothed during the Abel inverse transform. On the other hand, the cases when Z BA_r is significantly higher than Z Re_fd appear to be ambiguous according to the ECMWF comparison. Taking the cyan point and the magenta point in Fig. 5a for example, though it is difficult to recognize an ABL top by the refractivity (no refractivity gradient is less than -50), by the temperature (no inversion) or by the vapor pressure (no sharp decrease), an ABL top can be provided for reference in the BA profile. Thus, the detected ABL tops by the refractivity profiles in the lower right of the line y=x are usually less credible, for they do not meet some criteria [relative distinctness from (Ao et al., 2012) or sharpness constraint from (Chan and Wood, 2013)] and the ABL tops provided by the method in this paper might be a better choice.
Figure5. (a-d) Scatter diagrams of ZBA_r and ZRe_c for Λ≥ 1.75 in January, April, July and October. Panels (e-g), with reference to the red point in (a), are from atmPrf and echPrf_C005.2008.021.14.17.G03_2013.3520_nc. Panels (h-j), with reference to the magenta point in (a), are from atmPrf and echPrf_C005.2008.016.04.02.G02_2013.3520_nc. The cyan point in (a) denotes the case in Figs. 2g-i. The brown line is y=x.2
3.3. Results analysis
The above compares the ABL heights determined by the numerical differentiation method combined with the regularization technique on BA profiles Z Re_fd with the results obtained by other methods and data (refractivity profiles). The method in this paper will not be influenced by seasonal factors and are stable in the four months of comparison. Under the sharpness parameter constraint Λ≥ 1.75, the result using the finite difference method ZRe_fd obtains the minimum value of the ABL height, while COSMIC provides the highest ABL height by using refractivity data based on the breakpoint method ZRe_ c. Our developed method upon BA profiles gives medium ABL heights.From Figs. 2, 4 and 5, ZBA_r and ZBA_c have the best consistency, and ZRe_ c has significantly higher deviations than the other two methods when compared with ZBA_r, which is related to the definition, for another reason, that the ZRe_ c traces the top of the interfacial layers (Guo et al., 2011). If the constraint of the sharpness parameter is weakened (Λ≥ 1.5) or strengthened (Λ≥ 2), the above conclusion still holds. What is different is that, when the sharpness parameter is weakened, the numbers of eligible occultation profiles will significantly increase, the three kinds of bias (section 3.1 and 3.2) will increase, and the correlation coefficients will decrease; conversely, the numbers will markedly decrease, the biases will be reduced, and the correlation coefficients will further increase when the sharpness parameter is enhanced. Further analysis will be made in the next section on the effect of the sharpness parameter on the bias between the results obtained by different methods.
2
3.4. Relationship between the sharpness parameter and ABL height uncertainty
In order to study the effects of the sharpness parameter on estimating the ABL height, we use the COSMIC RO BA data and refractivity profile in January 2008 (similar in April, July and October) to show the relationship between the deviations of ABL height by different methods and the sharpness parameter Λ calculated in the present paper (Fig. 6).
Figure6. (a-c) Relationship between the sharpness parameter and the bias of ZBA_r from ZRe_c, ZBA_c and ZRe_fd in January respectively, respectively, where the red solid line represents the mean bias when Λ≥ i and the red dashed lines point out the corresponding 1-σ region.The data points in Fig. 6a distribute the most symmetrically, while the points in Figs. 6b and c are not as symmetrical as those in Fig. 6a, indicating that ZBA_r and ZBA_c are the best conformed, and the consistency of ZBA_r and ZRe_ c, ZBA_r and ZRe_fd, is relatively poor. The bias between ZBA_r and ZBA_c is generally zero from Λ=1 to Λ=4, while ZBA_r and ZBA_fd have positive deviation. There is negative deviation almost all the time between ZBA_r and ZRe_c, due to the definition of ZRe_c. The three kinds of standard deviation basically decrease with the increase in the sharpness parameter, and when Λ≥ 1.75 relatively little uncertainty exists (Figs. 6a and c). Generally speaking, there is small uncertainty when using different methods to estimate the ABL height as the sharpness parameter increases.
As the BA has obvious advantages over refractivity data, the numerical differentiation method combined with regularization techniques is applied to the COSMIC RO BA to estimate the ABL height. To verify the validity of this method, our result, ZBA_r, is compared with the ABL heights of ZRe_ c (based on the breakpoint method upon the refractivity) and ZBA_c (based on the maximum lapse of the BA) from COSMIC, and ZRe_fd (based on the finite difference method on refractivity). The comparison shows that, on average, the ZRe_fd ABL height is the lowest, while that from COSMIC (ZRe_ c) is the highest. The developed ABL height, ZBA_r, is a medium height. Of the cases with large differences, a lower ZBA_r, i.e., in the upper-left area, can often detect thin ABL tops that are sometimes neglected by methods based on the sliding window to avoid noise, or by the Abel integral process from refractivity to BA profiles. A higher ZBA_r can usually provide an alternative ABL height for the non-sharp cases that are usually excluded by other methods.
To give a metric for the reliability of the ABL height, the sharpness parameter is introduced, as in previous work. Statistical analysis shows that, for large sharpness parameters (sharp ABLs), the biases of the ABL heights by different methods is not too large——namely, there is small uncertainty in using different methods to determine the ABL height; whereas, for small sharpness parameters (non-sharp ABLs), the results of different methods are quite different——that is, the uncertainty in the ABL height determined by different methods is quite large.
Since the true value of the ABL height cannot be known, strictly speaking, we cannot judge which method is better or the best. However, the BA is rawer than the refractivity in the processing chain of COSMIC RO data, and therefore has less noise than the refractivity and contains more information. In terms of theory and numerical experiments, the regularization techniques have the function of smoothing, so, it is reliable to apply the numerical differentiation method combined with regularization techniques on BA profiles to determine the ABL height.
