Methodology for determining multilayered temperature inversions

Temperature sounding of the atmospheric boundary layer (ABL) and lower troposphere exhibits multilayered temperature inversions specially in high latitudes during extreme winters. These temperature inversion layers are originated based on the combined forcing of localand large-scale synoptic meteorology. At the local scale, the thermal inversion layer forms near the surface and plays a central role in controlling the surface radiative cooling and air pollution dispersion; however, depending upon the large-scale synoptic meteorological forcing, an upper level thermal inversion can also exist topping the local ABL. In this article a numerical methodology is reported to determine thermal inversion layers present in a given temperature profile and deduce some of their thermodynamic properties. The algorithm extracts from the temperature profile the most important temperature variations defining thermal inversion layers. This is accomplished by a linear interpolation function of variable length that minimizes an error function. The algorithm functionality is demonstrated on actual radiosonde profiles to deduce the multilayered temperature inversion structure with an error fraction set independently.


Introduction
The atmospheric boundary layer (ABL) is the lowest part of the troposphere in permanent contact with the earth surface and responds to thermal and roughness surface forcing in timescales of minutes to hours (Stull, 1988).Synoptic largescale meteorological processes condition the ABL state by driving relatively rapid horizontal air mass exchange for example under a cyclonic air mass advection (Fochesatto et al., 2002) or simply by imposing, on upper tropospheric levels, an adiabatic compression through the formation of a highpressure anticyclone (Byers and Starr 1941;Bowling et al., 1968;Curry, 1983;Cassano et al., 2011).These large-scale synoptic processes constrain the temperature of the air mass in tropospheric levels without initial connection to the local surface temperature forming vertically localized positive upward thermal gradients called elevated inversion (EI) layers (Csanady, 1974;Milionis andDavies, 1992, 2008;Mayfield and Fochesatto, 2013).
Although the ABL timescale for surface response ranges from minutes to hours (Garrat and Brost, 1981;Stull, 1988), the ABL response to synoptic forcing at local scale initiated by the radiative and dynamic interaction in the presence of the EIs can vary from hours to several days depending on a number of factors related to topographic conditions, temperature and wind distribution, thermal stratification and synoptic air mass type (e.g., moist cyclone advection or dry anticyclone situations, Mayfield and Fochesatto, 2013).
In polar atmospheres, the surface energy balance during winter is dominated by surface radiative exchange (Shulski and Wendler, 2007;Huff et al., 2010).This exchange occurs between a shallow layer close to the ground capped by the surface-based inversion (SBI) top height.The SBI also plays a thermalizing role in the surface cooling process by promoting radiative equilibrium between surface long-wave emission and the air mass radiative emission contained within the SBI layer depth (Malingowski et al., 2014).The SBI is present normally under non-advective conditions over iceand snow-covered surfaces and is a common feature in continental places as well as in the Arctic Ocean, sea ice and maritime environments (Overland and Guest, 1991;Ueno et al., 2005;Eastman and Warren, 2010;Bintanja et al., 2011).The SBI formation and persistence is also strongly dependent on terrain morphology and orientation (e.g., formation of cold pools Mahrt et al., 2001;Clements et al., 2003;Whiteman and Zhong, 2008;Fochesatto et al., 2013).Of particular interest is the formation and breakup of the SBI layer (Billelo, 1966;Hartman and Wendler, 2005), its thermal regime across time as it forms (i.e., fast cooling at the formation process followed by a growing depth phase) (Malingowski et al., 2014) and its thermodynamic characteristics controlling the sensible and latent heat exchange rate (Raddatz et al., 2013a, b).Moreover, simultaneous determination and characterization of EI layers and SBI has the potential to improve assessment on local air pollution meteorology (André and Mahrt, 1982), to improve forecasting of freezing rain episodes (Robbins and Cortinas, 2002), to determine heat exchanges in unconsolidated sea-ice surfaces (Raddatz et al., 2013 a, b) and potential impact on subarctic and Arctic climate (Ueno et al., 2005;Eastman and Warren, 2010;Bintanja et al., 2011).
The most widely used instrumentation to observe thermal inversion layers is direct temperature determination by highresolution radiosondes.However the methodologies to analyze the radiosonde temperature profile (see for example Seidel et al., 2010) still do not account for complex embedded thermal structures.On the other hand some authors have classified these layers as being short-living structures and therefore discarded from their analyses (Kahl 1990;Serreze et al., 1992).However, it was recently demonstrated during the Winter Boundary-Layer Experiment (Wi-BLEx) (Fochesatto et al., 2013;Mayfield and Fochesatto, 2013;Malingowski et al., 2014) that complex thermal structures persist with some changes through the day depending on environmental conditions.
Nonetheless, trying to detail more in-depth turbulent mixing processes in the ABL, the observations demand instrumentation that can profile the atmosphere at higher spatial and temporal resolution than what a radiosonde instrument can provide.Thus active remote sensing instruments give a better description of the vertical structure and dynamics of the ABL such as the case for example of lidar (Light Detection And Ranging, Fochesatto et al., 2001a, b, 2004, 2005, 2006) and sodar (SOund Detection And Ranging; Beyrich and Weill, 1993;Fochesatto et al., 2013).Particularly interesting is the high temporal and spatial resolution at which lidars can describe the initiation of convection close to the surface as well as the dynamic of the residual layer (Fochesatto et al., 2001 a, b).Similarly, insightful description of the turbulent and dynamic structure of the ABL when shallow cold flows penetrate the stable ABL can be gained by using sodar profilers (Fochesatto et al., 2013).On the other hand, passive ground-based remote sensors i.e., microwave radiometers, based on background sky microwave radiation sensing are another option to profile the ABL at ∼ minute resolution.Methodologies to determine the ABL height using lidars are based on detecting the changes in the optical backscattering level when the laser beam crosses the interface between the air mass in the ABL and the free troposphere (Menut et al., 1999;Bianco and Wilczak, 2002;Brooks, 2003;Lam-mert and Bosenberg, 2006).However, determining the ABL height using sodar becomes a little more complicated because the observations give access to the temperature structure of turbulence C 2 T profile that changes depending upon the ABL phase, terrain and multilayer structure (Holmgren et al., 1975;Beyrich and Weill, 1993).
Radiosonde determination of thermal inversion layers to study the Arctic Ocean ABL has been done using function fitting on radiosonde in the Arctic Ocean (Serreze et al., 1992) to determine the large-scale climate signature in the SBI time series of polar regions of the Arctic and Antarctic (Seidel et al., 2010;Zhang and Seidel, 2011;Zhang et al., 2011).In short-term studies, Khal (1990) and Khal et al. (1992) implemented a methodology to explore the temperature profile and then, applying some constraints, deduce the SBI top height, using similar methodologies trying to link time series of SBI heights to large-scale climate processes in continental Alaska (Bourne et al., 2010).
In this article a mathematical procedure is described, and its numerical implementation is documented to determine the multilayer thermal structure present on a given temperature profile in the lower troposphere.Section 2 describes the mathematical procedure and numerical implementation, Sect. 3 determines the calibration factor, Sect. 4 applies the numerical procedure to routine and high-resolution radiosonde observations and Sect. 5 gives a final assessment of the implemented numerical methodology.

Methodology
A thermal inversion layer represents a region in the atmosphere defining a positive upward temperature gradient dT dz > 0, where T is the temperature and z is the height.As described in the Introduction this positive upward temperature gradient can have two origins: the first one related to surface cooling (i.e., surface net radiation loss) giving origin to the formation of the so-called SBI, and the second one originated after synoptic large-scale processes (i.e., upper level warming also known as EI layers) (Mayfield and Fochesatto, 2013).Under these conditions the temperature profile in the lower troposphere is expected to describe either a neutral thermal condition (no inversion present), isothermal layers, presence of a surface-based inversion (SBI), development of stratified layers within the SBI layer depth, shallow inversion layers above the SBI inversion top or more complex situations in which upper-level thermal inversions are also present with base and top in the outflow ABL (Khal, 1990;Mayfield and Fochesatto, 2013).Figure 1 depicts the more relevant situations.
To deduce the existing thermal inversion layers in a temperature profile, a computational routine was implemented to identify the heights and temperatures of the base and top of inversion layers.In this case the temperature profile was reduced to a mathematical expression based on a linear func- tion calculated on the basis of an iterative linear interpolation routine.The numerical routine is based on fitting and refitting a piecewise linear interpolating function with variable length to define layer by layer the temperature profile composition.An error quote controls the algorithm convergence by layer and is calculated based on the largest singular value ε by which thermal features in the temperature profile can be reproduced.The numerical method actually reduces the number of points present in the original temperature profile.Therefore after applying the numerical methodology, the temperature profile became a data structure defining thermodynamic characteristics of thermal layers found in the temperature profile.Further on, exploring the data structure we can define a thermal inversion layer as the layer in which the temperature gradient turns positive upward to negative upward.
Equation (1) indicates the mathematical formula describing the piecewise linear approximation of the temperature profile: where k represents each linear interval fitting and N is the total number of temperature measurements.The sequence of linear fitting is executed by selecting two temperature points at the time -one from the top height and the other from bottom height maintaining bottom fixed -until the linear piecewise function minimizes the error quote against the temperature profile as indicated in Eq. ( 2): Here ε represents the mathematical norm or Euclidean distance between the linear piecewise mathematical representation φ (z) and the original temperature profile T (z).This way the temperature profile is assimilated to the series of sequentially fitted linear intervals defining the temperature function.It has to be noted that after the layer (k) has been detected, the layer (k+ 1) considers its lowermost point as being the top height of layer (k) and starts over again "connecting" this point with the top-height point of the temperature profile and, from there, descending in height until the ε value reaches the preset threshold.This way the procedure continues until the entire temperature profile has been scanned and assimilated to frames defined by piecewise linear functions.
Basically the number of piecewise linear functions (k) will be defined by the ultimate number of temperature measurements composing the temperature profile (N ); thus max(k) = N − 1.However, increasing the error fraction (ε) will therefore reduce the max(k) reached.In such a case the final approximation formula is expressed in Eq. (3): where p is the total number of points from the original temperature profile that has been retained by the fitting process.
Table 1 indicates a five-layer temperature structure in the ABL up to 1 km height based on six-point temperature simulation.The profile defines a stratified SBI with top height at 250 m topped by an EI at 670 m. Figure 2 illustrates the sequences in which the algorithm operates to search and define each layer from 1 to 5. In this case we retain all points in the temperature profile.

Convergence factor
Since the ability of the algorithm to extract thermal inversion layers depends upon the value chosen for ε, then it is reasonable to determine for example what threshold value of temperature profile gradient dT /dz in • C/100 m can be retrieved when a given value of ε is preset.To deduce this relationship a numerical simulation was performed scanning values of ε from 0.1 to 20 over a temperature profile displaying a single thermal inversion which was varied in temperature strength T from 0.1 to 20 • C and vertical depth z from 1 to 500 m.For each step on ε values the temperature gradient threshold dT /dz, able to determine that specific layer, was deduced.
Figure 3 panel on the left illustrates the simulated temperature inversion and the applied temperature strength and vertical depth variability.Similarly the temperature gradient threshold as a function of the convergence error ε is shown in Fig. 3 in the panel on the right.This relationship was determined in an ideal scenario of a single temperature inversion layer.Following Fig. 3 (right panel), if a constant inversion layer depth is analyzed (e.g., z = 200 m) and starting on ε = 20, the minimum temperature strength across the given layer will be 6.5 • C/100 m while reducing ε = 5 allows detection of 2 • C/100 m; if ε = 1 is considered, then the detection limits drop to less than 1 • C/100 m.It has to be noted that ε is a convergence factor that gives to the algorithm more  or less sensitivity to detect thermal structures present in the given temperature profile.

Numerical application
In this section the numerical method is applied to the restitution of thermal layers from radiosonde data.Two cases are examined.The first one comes from the radiosonde data archive of NOAA-ESRL database for the Fairbanks International airport WMO station code 70261 retrieved from the Wyoming radiosonde database (http://weather.uwyo.edu/upperair/sounding.html).The sounding was taken on 15 January 2014 at 00:00 UTC (Fig. 4 -left panel).The second case corresponds to one of the GPS sounding cases taken during an intensive observing period in the framework of Wi-BLEx (Malingowski et al., 2014).In this case a morning sounding was chosen during 13 October 2009 at 15:00 UTC (Fig. 4 right panel) during formation process of the SBI layer.
Table 2 summarizes all retrieved layers extracted presetting ε = 0.1.This setting allowed detecting p = 24 layers over 32 temperature points below 3 km in the case of 15 January 2014 at 00:00 UTC, while for 13 October 2009 at 15:00 UTC the sounding presents 677 points below 3 km and p = 56 layers were detected.
From the analysis of detected layers one can initially search in the data structure the changes in m k .m k > 0 will be assigned to positive upward dT /dz, while negative will indicate the opposite.From the analysis of the resulting data it can be deduced that the inversion layer, defined when m k changes from positive to negative, will be assigned as SBI when the base of that inversion is located at the surface, regardless of its stratification condition, while, on the other hand, any other inversion layer with base off the ground will be called EI.In the case of 15 January 2014 at 00:00 UTC, it was found that the SBI top height is located at 165 m above the ground, while for 13 October 2009 at 15:00 UTC is was located at 44 m above the ground.On the other hand two EI layers were detected in the first case with inversion top heights at 925 and 1751 m from the surface.However in the second case a number of EI layers have been identified with top inversion heights at 289, 425, 541, 654, 1098 m, etc.The retrieved thermal inversion SBI and EI layers are summarized in Table 2.
Furthermore, the algorithm performance to extract the multilayered thermal structure was evaluated using the abovementioned examples in terms of the number of layers detected as function of the preset convergence error ε and the final temperature retrieved error (%).In this case the algorithm was executed in a loop varying ε in the interval from 0.1 to 1.0 in steps of 0.1 and from 1.0 to 15 in a 0.5 step.The final temperature error was computed based on the norm of the final temperature vector resampled in terms of the original vector height.This explains why the actual preset value ε, which applies to every step in the fitting and re-fitting process, is different from the final temperature profile error after fitting the entire temperature profile.Figure 5 depicts two curves indicating the number of thermal layers retrieved and the final overall temperature error profile as a function of the ε preset value for both sounding cases.Table 2. Results of application of the numerical routine to extract thermal layers from radiosonde database NOAA-ESRL PAFA 70261 for 15 January 2014 at 00:00 UTC and IOP WiBLEx GPS sounding on 13 October 2009 at 15:00 UTC.In this case ε = 0.1 for both profiles.The temperature profile contains 32 points below 3 km height and the program detected p = 24 layers; for 13 October 2009 the sounding contains 677 points below 3 km and p = 56 layers.In the classification column, the thermal inversion layers detected are indicated as SBI or EI.First column represents the layer number, second column is the layer base height Z b (m), third column is the layer top height Z t (m), fourth column is the temperature at Z b , fifth column is the temperature at Z t , sixth column is the temperature gradient m k ( • C/100 m), seventh column is the depth of the layer (m) and eighth column is the layer classification.In order to summarize the application of this methodology and visualize the effects of the final error over the reconstructed temperature profile, Fig. 6 illustrates the resampled temperature profile based on the case of four selected error levels for both cases.
It is clear that the lowest error (%) is associated with lowering the convergence factor ε, while on the other hand increasing the preset convergence factor ε increases the overall temperature retrieved error; therefore, the fine structure in the temperature profile vanishes.This is very important to consider.For example, for more than 10 % of error the entire ABL thermal structure of local forcing disappears indicating just the presence of the EI layers.

Summary
A simple numerical iterative routine was developed to deduce thermal inversion layers composing a given atmospheric temperature profile.The numerical routine converts the temperature profile in a reduced data structure containing base and top heights and temperatures as well as thermal gradients allowing identification of most important air mass changes accounting for surface-based inversion (SBI) and elevated inversions (EI) present in the temperature sounding.
This methodology has been applied to the study of 10 years of upper air data in the interior of Alaska to deduce the statistical significance of the presence of SBI and multilevel EIs.In addition, the application of this methodology allowed classification of EI layers based on the type of synoptic air mass by means of the dew-point temperature gradient across the inversion layer depth (Mayfield and Fochesatto, 2013).
The method was also applied to study the temporal evolution of SBIs using high-resolution GPS soundings during the formation and breakup cycle in spring and fall seasons (Malingowski et al., 2014).
The method does not introduce new temperature points; it rather reduces the number of them according to the preset convergence factor.After the data structure is obtained, the determination of temperature inversion layers follows the search of changes in the m k coefficient.
Based on this methodology temporal series of temperature inversion height can be formed analyzing the data structure output and, after applying a specific set of criteria, identify a time series of a particular nature.For example, in Mayfield and Fochesatto (2013), the interest was to perform a statistical analysis of heights of SBIs and the presence of EI-1 to −4 layers.Using this numerical routine, it was possible to classify them in terms of their dew-point temperature gradient across the inversion depth.For that, it was necessary to include the dew-point temperature in the data structure after reading and processing the radiosonde data so that the retrieved EI layer could be classified.The relationship between threshold thermal gradient dT /dz, the preset convergence factor ε and the overall final error between the resample temperature profile and the original temperature profile is not straightforward.Section 3 describes the relationship between dT /dz and ε for a theoretical profile to deduce this relationship.Section 4 applies the methodology, and then a resampled version of the temperature profile was compared with the original temperature profile to calculate an overall percent error.This is important to differentiate since ε applies internally as a convergence factor to increase fidelity in the temperature fitting process over a Euclidean norm that applies over a variable vector length step by step.On the other hand, the overall error indicating how accurate the resampled temperature profile reproduces the original profile accounts for the entire profile at once.
The application of such methodology has numerous venues.For instance, air pollution meteorology and pollutant transport processes are strongly dependent on meteorological conditions near the emission source -in particular at high latitudes where the SBI layers and stratified SBI show a persistent presence during the winter in the absence of or under a little diurnal cycle.The application of this methodology is not restricted to low-level layers exclusively since it has been demonstrated by Mayfield and Fochesatto (2013) that EI layers also play an important role in the high-latitude ABL because of the dynamic-radiative forcing that drives surface warming and pollutant recirculation.Another application is the evaluation of large-scale synoptic meteorological changes and their impact on high-latitude environments.In this case it is important to consider a methodology that can scrutinize thermal layers present in the temperature profile over ocean, maritime and continental places.In order to be able to evaluate the impact of prescribed synoptic meteorological patterns over a specific area, it is necessary to retrieve simultaneously thermal layers affected by local-scale and those dominated by large-scale atmospheric motion.Therefore evaluation of temporal series containing this multilayered information could assist in improving our understanding of how large-scale feedbacks are introduced and affect the local environment.Since the radiative interaction between SBI and EI layers is controlled by the vertical stratification, then the analysis of surface warming in high-latitude continental and oceanic environments in particular over sea-ice surfaces will be one future application of this methodology.Of course databases in which the method can be applied need to have sufficient vertical resolution to pull level details close to the surface to be able to define the multilayered structure.
The algorithm was proposed for radiosonde profiles, but it can also process temperature outputs from microwave radiometers and large-scale reanalysis databases as well as cli-mate outputs downscaled to temperature profiles.Particularly interesting is the use of this algorithm in climatological determination of the temperature inversion layers SBI and EI and their connection to large-scale synoptic meteorology.

Figure 1 .
Figure 1.Illustration of temperature profiles indicating different thermodynamics conditions from left to right: (a) no inversion present; (b) presence of an SBI; (c) presence of a stratified SBI; (d) SBI stratified in the presence of an EI and (e) SBI stratified and multiple EIs.Modified after Khal (1990).

Figure 2 .
Figure 2. Sequence of the numerical routine to identify the thermal layers indicated in Table1.Fitting function is given by the thick trace with "o" markers while the temperature profile is indicated by thin black trace with "*" markers.Top panel indicate the sequence to detect Layer 1, second panel from top is sequence to detect Layer 2, third panel from top is the detection of Layer 3, fourth panel from top is to detect Layer 4 and bottom panel is the sequence to determine Layer 5.

Table 1 .Figure 3 .
Figure 3. Convergence error simulation: the left panel illustrates the simulated temperature inversion varying Z in meters and T in degrees Celsius.Right panel is the minimum thermal gradient detected expressed in • C/100 m on a single inversion layer as function of the preset error norm ε varying from 0.1 to 20.The parameter in this curve is the layer depth represented in this case from 1 to 200 m.

Figure 4 .
Figure 4. Analysis of radiosonde temperature profile: panel on the left represents temperature observation made on 15 January 2014 at 00:00 UTC and panel on the right is a GPS radiosonde observation during IOP WiBLEx on 13 October 2009 at 15:00 UTC both at the National Weather Service (NWS) Fairbanks International Station (code 70261).Vertical axis is height (m) above ground level.

Figure 5 .
Figure 5. Retrieved number of thermal layers (o) and final temperature profile error ( ) as function of the ε convergence preset value.Panel on the left corresponds to 15 January 2014 at 00:00 UTC and panel on the right corresponds to 13 October 2009 at 15:00 UTC.

Figure 6 .
Figure 6.Resampled temperature profile after running the numerical methodology over different preset convergence factor values ε = 0.1, 1.0, 5.5 and 10.5 as function of the final retrieved error (%).Panel on the left is for 15 January 2014 at 00:00 UTC and panel on the right is for 13 October 2009 at 15:00 UTC.