Determining Atmospheric Boundary Layer Height with the Numerical Differentiation Method Using Bending Angle Data from COSMIC

Shen YAN,Jie XIANG?,2,and Huadong DU

1College of Meteorology and Oceanography,National University of Defense Technology,Nanjing 211101,China

2Key Laboratory of Mesoscale Severe Weather of Ministry of Education,Nanjing 210093,China

ABSTRACT This paper presents a new method to estimate the height of the atmospheric boundary layer(ABL)by using COSMIC radio occultation bending angle(BA)data.Using the numerical differentiation method combined with the regularization technique,the first derivative of BA profiles is retrieved,and the height at which the first derivative of BA has the global minimum is defined to be the ABL height.To reflect the reliability of estimated ABL heights,the sharpness parameter is introduced,according to the relative minimum of the BA derivative.Then,it is applied to four months of COSMIC BA data(January,April,July,and October in 2008),and the ABL heights estimated are compared with two kinds of ABL heights from COSMIC products and with the heights determined by the finite difference method upon the refractivity data.For sharp ABL tops(large sharpness parameters),there is little difference between the ABL heights determined by different methods,i.e.,the uncertainties are small;whereas,for non-sharp ABL tops(small sharpness parameters),big differences exist in the ABL heights obtained by different methods,which means large uncertainties for different methods.In addition,the new method can detect thin ABLs and provide a reference ABL height in the cases eliminated by other methods.Thus,the application of the numerical differentiation method combined with the regularization technique to COSMIC BA data is an appropriate choice and has further application value.

Key words:atmospheric boundary layer height,numerical differentiation method,COSMIC,bending angle,regularization

1.Introduction

The atmospheric boundary layer(ABL)is the closest layer to the Earth's surface.Characterized by significant diurnal variations and turbulent motions,it is of great importance to the vertical transport of heat,momentum,water vapor and chemical composition,the evolution of large-scale weather processes and the study of climatology(Garratt,1994).The ABL height is one of the important parameters of the ABL,which is usually used in the parameterization of the physical processes related to the ABL.An accurate ABL height has great influence on improving the precision of numerical weather forecasting models and general circulation models(Medeiros etal.,2005).Therefore,estimating the ABL height has been a challenging topic and one of the frontier research fields in atmospheric science(Hong et al.,1998;Xu et al.,2002;Li et al.,2006;Mao et al.,2006;Seidel et al.,2010;Dai et al.,2014).

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.

Fig.1.COSMIC RO refractivity and BA profiles(atm-.2008.275.03.49.G312013.3520nc).

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 illposedness(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.Methods and data

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),

whereTis the temperature(K),Pis the total atmospheric pressure(hPa),andeis 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 Δzin 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 verti-cal 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,zb≤z≤zt,wherezbandztrepresent the bottom and top of the α(z)profile,respectively.Grid points{zi,i=1,···,m}divide[zb,zt]intom-1 parts(zb=z1＜z2＜ ···＜zm=zt)with an equal interval ofh=(zt-zb)/(m-1).The values of α(z)atzi,α(zi),are known by the COSMIC RO data,and the next goal is to obtain the first derivative of α(z)atzi,α′(zi).

For convenience,some notation is introduced as follows:α′(z)≡ φ(z),α(zi)≡ αi,φ(zi)≡ φi.By the Newton—Leibniz formula,there is

and the integral of the right-hand side of Eq.(2)can be approximately obtained by the Simpson formula,

which can be written in matrix form as Now,the problem is down to solving these linear equations.

The linear equations[Eq.(4)]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.(4)is transformed into solving the minimization problem of the Tikhonov functional with a regularization functional included,

where γ‖LX‖2is the regularization functional and γ ＞ 0 is the regularization parameter that should be given in advance.The matrix

is the first derivative operator that can control the smoothness of the solutionxvia choosing the correct value of the regularization parameter γ.It is simple to find that,as γ increases,will decrease and thenwill be smoother.

Equation(5)is reduced to solving the linear equations(well posed)as

and its solution,

is an optimal stable approximate solution of Eq.(4),where the regularization parameter γ is determined by the L-curve method(Hansen,1992),which balances the two items in Eq.(5).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

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 ofh=10 m.

3.Validation of the method

The main purpose of this section is to verify the scheme for determining the ABL height proposed in section 2 by comparing the ABL heights obtained from the numerical differentiation method upon the BA profile datawith those obtained by other methods.Since the true values of the ABL heights are unknown,we mainly consider the case of the sharp ABLs.The information on the sharp ABLs is easier to extract,and the ABL heights obtained by different methods are in good consistency,which can be used to test the applicability of the new method.Due to the constraints of the sharpness parameter on ABLs,the numbers of eligible occultation profiles in the four months are 2924,3666,3417 and 3617,respectively.Of the profiles successfully extending to below 500 m above the surface,about 28.7%on average are used.Other ABL heights for comparison include:two kinds of ABL height from COSMIC data products[obtained by the BA based on the maximal BA lapse at 0.3 km intervals(Sokolovskiy et al.,2007),and the refractivity based on the break point with maximum dN/dzlapse(Guo et al.,2011),denoted asand]and the ABL height determined by the finite difference method upon COSMIC refractivity data(denoted as).As mentioned above,there usually occurs a temperature inversion and/or a sharp decrease in water vapor at ABL tops.Thus,in some cases with large differences,the echPrf files are used to assist with the validation.

3.1.Comparison with ABL heights from COSMIC products

3.1.1.Case of BA

Fig.2.(a—d)Scatter diagrams ofandfor λ ≥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.,with reference to the red point in(a),and the temperature and vapor pressure profiles from its corresponding echPrf file(.2008.027.05.43..3520nc).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(.2008.015.04.27..c,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).

First,we compare the ABL heightwithunder the constraint λ ＞ 1.75.From Figs.2a—d,the average bias betweenandin January,April,July and October is-0.04 to-0.02 km(is slightly lower than),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 theseems more acceptable.This example represents the majority of cases with large differences in the upper-left area of the liney=x.That is,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.

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.

3.1.2.Case of refractivity

Next,we compare the ABL heightwithfor λ＞1.75.From Figs.4a and b,the average bias betweenandis-0.260 to-0.220 km,and the correlation coefficient is higher than 0.85.That is,is less thansystematically,but there is still a high correlation betweenand.In addition to the points of good correlation,there are also some points with large deviations.

Taking there d point in Fig.4 a as an example,=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 smallscale 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.

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 betweenandis 0.071—0.082 km,and the correlation coefficient is higher than 0.95.That isis also slightly higher than,and there is still a high correlation betweenand.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.(1)than the vapor pressure at that height,and themight 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 whenis significantly higher thanappear 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 liney=xare 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.

Fig.3.Scatter diagrams of andin January for(a)λ＞1.5 and(b)λ＞2.The brown line is y=x.

Fig.4.(a—d)Scatter diagrams ofand for λ ≥ 1.75 in January,April,July and October.Panels(e—g),with reference to the red point in(a),are from atmPrf&.2008.013.05.48.Panels(h—j),with reference to the magenta point in(a),are from atmPrf&.2008.018.18.56.The cyan point in(a)denotes the case in Fig.2(g)—(i).The brown line is y=x.

3.3.Results analysis

Fig.5.(a—d)Scatter diagrams of and for λ ≥ 1.75 in January,April,July and October.Panels(e—g),with reference to the red point in(a),are from atmPrf and .2008.021.14.17..Panels(h—j),with reference to the magenta point in(a),are from atmPrf and.2008.016.04.02..The cyan point in(a)denotes the case in Figs.2g—i.The brown line is y=x.

The above compares the ABL heights determined by the numerical differentiation method combined with the regularization technique on BA profileswith 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 methodobtains the minimum value of the ABL height, while COSMIC provides the highest ABL height by using refractivity data based on the breakpoint method.Our developed method upon BA profiles gives medium ABL heights.From Figs.2,4 and 5,andhave the best consistency,andhas significantly higher deviations than the other two methods when compared with,which is related to the definition,for another reason,that thetraces 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.

Fig.6.(a—c)Relationship between the sharpness parameter and the bias of from,andin 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.

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).

sically 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.

4.Summary and discussion

The ABL is a transitional area between the free atmosphere and land or ocean areas.It is of great importance in air—sea interaction,large-scale synoptic processes,and the prediction of atmospheric pollutants.The ABL height is one of the key parameters of the ABL,and plays an indispensable role in the parameterization of its physical processes.Occultation technology,based on the global navigation satellite system,provides us with a lot of high-accuracy and high resolution atmospheric profiles,which opens up a new data source for detecting the ABL's structure.

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,,is compared with the ABL heights of(based on the breakpoint method upon the refractivity)and(based on the maximum lapse of the BA)from COSMIC,and(based on the finite difference method on refractivity).The comparison shows that,on average,theABL height is the lowest,while that from COSMICis the highest.The developed ABL height,,is a medium height.Of the cases with large differences,a lower,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 highercan 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.

Acknowledgements.This study was supported by the National Natural Science Foundation of China(Grant No.41475021).Thanks to the UCAR COSMIC project for providing occultation data,and thanks to Dr.HE for his help in downloading data.Finally,thanks to the anonymous reviewers for their valuable comments.