Residual kriging analysis of airborne measurements: application to the mapping of atmospheric boundary-layer thermal structures in a mountain valley

Residual kriging (RK) is applied to data from measurements performed with an instrumented motorglider in an Alpine Valley, where a coupled lake and valley wind develops. Results provide an insight into a typical mean vertical structure of the valley boundary layer, displaying a rather shallow convective mixed layer, surmounted by a deep stable layer. Mapping with RK allows to detect spatial temperature anomalies associated with the local development of diurnal winds over each explored valley cross section. Local scale features, such as cross-valley thermal asymmetries, amenable to surface inhomogeneities and their effects on heat ﬂuxes, are also identiﬁed. Copyright 2013 Royal Meteorological


Introduction
Instrumented light airplanes are ideal tools to investigate fine-scale structures of the atmospheric boundarylayer (ABL), especially in mountainous regions, where a variety of daily-periodic winds typically develop under fair weather conditions (Zardi and Whiteman, 2013). Indeed several projects, investigating valley wind phenomena, included research flights in their intensive field observations. In particular, airborne measurements performed under the projects DISKUS (Hennemuth, 1985(Hennemuth, , 1986, MAP-Riviera (Weigel and Rotach, 2004;Rotach and Zardi, 2007) and ALP-NAP (de Franceschi and Zardi, 2009;Gohm et al ., 2009;Harnisch et al ., 2009) provided an unprecedented basis to obtain high-resolution pictures of threedimensional (3D) fields of temperature, wind, turbulent fluxes and chemical species throughout a valley volume. A crucial step in the analysis of measurements collected along the flight trajectories is represented by the choice of an appropriate method to obtain a suitable remapping of data on regular grids. In fact, this is essential to extract readable views of the 3D complex structures explored by the airplane. To this purpose various interpolation techniques were used in the literature: Egger (1983) adopted an inverse distance weighted interpolation, Hennemuth (1985) an exponentially-weighted inverse squared distance (exp-ISD) interpolation, while de Wekker (2002) and Weigel and Rotach (2004) used a natural neighbour interpolation with Delaunay triangulation.
In this study we apply the geostatistical technique called residual kriging (RK) (Ahmed and de Marsily, 1987), which relies on intrinsic stochastic properties of the dataset. Kriging techniques are widely applied in meteorology, but mainly for the interpolation of surface variables from sparse data, such as rainfall, air temperature, radiation and soil moisture (see Hudson and Wackernagel, 1994 for an example).
Here, we show that kriging is also a suitable tool for the integrated mapping of airborne and surface measurements, allowing for the retrieval of 3D finescale features of the observed meteorological fields.

Target area and field measurements
The present analysis concentrates on the lower Sarca river Valley and the adjacent Lakes Valley in the Alps (Figure 1). These two valleys form a corridorconnecting the area north of Lake Garda to the Adige river Valley -where a coupled lake and valley breeze typically occurs on sunny days, the so-called Ora del Garda. This wind arises over the northern shorelines of Lake Garda, and propagates northward up-valley, until it breaks into the adjacent Adige Valley, near the city of Trento, through an elevated saddle.
Typical patterns of this wind were previously investigated on the basis of surface measurements only (de Franceschi et al ., 2002). For this reason a research project was later organized to explore the 3D thermal structures associated with the Ora del Garda by means of airborne measurements. Accordingly, a flight was performed on 24 September 1998, under weak synoptic forcing and cloudless conditions, securing the full development of local wind systems. The flight trajectory was designed so as to follow ascending and descending spiralling paths over selected valley cross sections ( Figure 1).
Details on the measurement platform and its operation are reported in . Notice that all the sensors were slow-response with respect to the typical timescales of turbulent fluctuations. Hence they naturally filtered out high-frequency temporal and spatial fluctuations of temperature and pressure, and recorded mean values only. Figure 2 shows potential temperature profiles obtained from plotting airborne data versus height (i.e. neglecting horizontal coordinates) for each valley section, along with the respective flight timings. Flying each spiral required less than 30 min, which is half the typical timescale of ABL structure daytime evolution (Stull, 1988): this allows to consider the whole atmosphere under a quasi-steady state during each spiral. Grey lines are pseudo-soundings obtained from moving-window (width: 230 m) vertical averages of airborne data, extended below the sampled region under the assumption of a constant potential temperature in the lowest ABL, the so called diurnal mixed layer (ML). Also reported in Figure 2 are observations recorded around the overflight time at surface stations nearest to the area flown over by each spiral.

Materials and methods
The pseudo-soundings represent the dominant vertical structures of the valley ABL. Small-scale anomalies, reflecting local patterns associated with specific flow features, driven by complex topography and land cover heterogeneities, are superimposed as a stochastic component. These thermal anomalies can be retrieved and visualized by mapping the potential temperature field on regular grids by means of a kriging technique.
Kriging is a class of geostatistical methods devised to predict the best estimate of a regionalized variable Z at any unexplored point x 0 . The estimate consists in a weighted average of observed values at nearby points, where the weights are determined from the spatial covariance of observations, minimizing the estimate variance under an unbiasedness condition (Matheron, 1971;Journel and Huijbregts, 1978;Isaaks and Srivastava, 1989;Goovaerts, 1997). Among the various kriging options, the most appropriate for this purpose is RK (Ahmed and de Marsily, 1987;Odeh et al ., 1994;Odeh et al ., 1995;Goovaerts, 1999). RK assumes that the target field Z can be decomposed into a deterministic drift µ and stochastic residuals δ, as follows: The former is usually estimated from a best-fit approximation of the original data to some structure, presumed a priori from the physics of the observed phenomena. In this case the pseudo-soundings ( Figure 2) naturally provide the required drift: this Figure 2. Potential temperature measurements (black dots) at each valley cross section compared with routine surface observations (averaged on overflight times indicated in the panels) from nearby stations (grey bullets; labels refer to Figure 1). Grey lines are the pseudo-soundings representative of the mean vertical structure (drift) used for the evaluation of residuals in the kriging procedure. The local valley-floor height for each cross-section is indicated by the grey shading. All the heights are above mean sea level (MSL).
implies that the deterministic component depends mainly on the vertical coordinate, as typical in the ABL. Residuals instead are separately estimated by means of a standard kriging procedure. In the present case the residuals, resulting from subtraction of the above drift from the original data, displayed local second order stationarity, thus allowing for the application of an ordinary kriging algorithm (Journel and Huijbregts, 1978).
The latter assumes the following expression for the residual estimate δ (x 0 ) based on observations at n nearby points x α : Kriging techniques imply that the weights, λ α , be evaluated by means of the so called semivariogram function γ (h), which estimates the spatial dependence of the observations covariance in terms of lagaveraged pairwise dissimilarity: where h is the separation lag vector between the positions of observation pairs, N (h) is the number of pairs separated by h and Z (x α ) the observed value at x α (Goovaerts, 1997). In practice, semivariograms may be empirically estimated by a variety of estimator functions. Here the robust estimator proposed by Cressie and Hawkins (1980) was found to give the best results in terms of semivariogram stability. In fact empirical semivariograms are experimental realizations of the intrinsic covariance structure of the observed field, which is assumedly representable by a mathematical function to be determined. To this purpose, various theoretical models may be used to fit empirical semivariograms, and then used to compute the weights λ α (see Goovaerts, 1997 for further details). For the present case we adopt a spherical model: Here the nugget n (i.e. semivariogram value at the origin) accounts for variability at distances smaller than the shortest sampling interval or lower than the measurement resolution, the sill s is the asymptotic value for increasing lags, and the practical range r is the lag distance at which the semivariogram practically reaches the sill value (Journel and Huijbregts, 1978). Beyond this distance, correlation between data is negligible. Notice, however, that spatial covariance may depend not only on the modulus, but also on the direction of the lag h between pairs. In the present case, the two principal directions are along the vertical and the cross-valley coordinates, respectively. Therefore, horizontal and vertical directional semivariograms of residuals were evaluated separately. As an example these are shown in Figure 3 for spiral B/w. As expected, vertical and cross-valley characteristic ranges are different. Similar results are found on all the explored cross sections, with ratios between the two ranges varying between 3.5 and 6.0. The two ranges also provide the appropriate scaling factors for normalizing the space coordinates to get an isotropic field. As the flight patterns were performed closely around vertical cross-valley planes, the (weak) variability in the third direction, i.e. along-valley, could not be derived from the data. Instead the same range was assumed for this direction as for the cross-valley one. On the basis of the above reasoning, omnidirectional semivariograms were then computed and finally used for RK mapping.
From the above description it can be argued that kriging offers a series of advantages compared with other interpolation methods. First, it is an exact interpolator, i.e. it returns the observed values at measurement points. Second, kriging assigns the interpolation weights according to the characteristic covariance of the target field, not on the sole basis of the geometric configuration of observation points. Third, it allows to account for the intrinsic anisotropy of the target field, by using the appropriate scaling factors for the principal space directions. Moreover the estimate variance represents an objective measure of confidence in the  estimate at each prediction point. A preliminary comparison between RK and ISD method applied to the analysis of airborne measurements  confirmed the better performance of the former.
For the sake of completeness, it must be mentioned that RK is computationally more expensive than the other mentioned methods.

Results and discussion
Some common features of the ABL structures in the explored valleys can be outlined from the pseudosoundings shown in Figure 2. In particular, at all the cross sections the vertical thermal structure displays a rather shallow ML, surmounted by a deep stable layer (lapse rate ∼5 K km −1 ). Such a structure has been observed in many other valleys (Weigel and Rotach, 2004), and is associated with daytime subsidence of potentially warmer air from the free atmosphere. Indeed the latter compensates for the air removed at the valley floor by thermally driven upslope winds occurring at the valley sidewalls. Subsidence also inhibits the convective ML growth, and stabilizes the upper part of the profile Kimura, 1995, 1997;Rampanelli et al ., 2004;Weigel and Rotach, 2004;Rotach and Zardi, 2007;Serafin and Zardi, 2010a, 2010b.
The upper stable layer for section A displays an overall stability (∼5 K km −1 ) slightly larger than the other ones (∼4 K km −1 ): this is consistent with the stronger downward velocity associated with subsidence in narrower valley sections (B and C). Despite the relatively early stage of the convective boundarylayer, the incipient development of a warmer quasimixed layer at the crest level can be envisaged. This is very likely produced by the reconnection to the valley centre of the airflow generated by upslope flows at the sidewall crests, as outlined by Kimura (1995, 1997) and Serafin and Zardi (2010a, 2010b. Potential temperature fields were interpolated by means of RK on regular grids (50 × 50 × 50 m 3 spacing) around the vertical planes shown in Figure 1. For each cross section the basis for RK application was provided by residuals of airborne data, as explained above, and nearby surface observations, where available. Examples of the interpolated fields are discussed below. The mapping areas include only grid points where the interpolation can be reliably trusted (i.e. nodes whose distance from observation points is smaller than the characteristic semivariogram range) and are further restricted to the interior of each valley basin. For all the represented cross sections, the estimate standard deviation (i.e. the square root of the estimate variance; not shown) was everywhere smaller than 0.35 K, tending to null values when approaching trajectory points, as expected. Data from surface stations complement the airborne dataset for the lowest layers not covered by the flights and display a very good matching. Figure 4(a) shows the RK-interpolated potential temperature field on the vertical cross section A. The lowest layers display a uniform temperature distribution up to 590 m above ground level (AGL) throughout the cross-valley section, consistently with the occurrence of a ML on the wide (∼5 km) valley floor. A minor East-West asymmetry is detected in the intermediate layer below ∼1650 m AGL, reasonably associated with up-slope flows stronger on the East facing sunlit sidewall at the time of the flight. Figure 4(b) shows the cross sections including both the area over Cavedine Lake on the western side (spiral B/w), and the elevated Cavedine Valley on the eastern half of the valley section (spiral B/e). Interpolated fields clearly highlight an asymmetric thermal structure: on the western side the sunlit steep sidewall in bare rock is likely to be warmer than the valley center; here Cavedine Lake seems to favour colder conditions, as suggested by time series of surface measurements at from the nearby station labeled CAV in Figure 1 (not shown). In section B/e, again a 400-m deep ML occurs above the elevated valley floor, and at upper levels a minor asymmetry, detectable from depressed isentropes, is promoted by downward subsidence at the valley center. Figure 4(c) shows the potential temperature field at spiral C. The thermal structure displays remarkable cross-valley asymmetries, clearly reversing around 1300 m MSL. This picture reflects a different heating on the two sidewalls: indeed, the North-West facing sidewalls are never directly sunlit during the day, whereas the opposite side of the valley gets direct solar radiation since early morning. Such a contrast is expected to produce a cross-valley circulation (Hennemuth, 1986) with an up-slope flow occurring on the South-East facing sidewall, and a down-slope flow on the opposite side. The two circulations are part of a unique cross-valley cell, which tilts isentropes downward towards the sunlit side. Such an asymmetric pattern is confirmed by the complementary view offered by horizontal sections at different levels ( Figure 5). Possible interference from upper winds may be excluded, on the basis of synoptic observations clearly showing a high-pressure situation and a light south-westerly wind over the area (not shown). Instead the local eastward bending of the valley axis may enhance the purely thermally-driven cross-valley flow. Indeed local curvatures force alongvalley winds to develop secondary circulations over the whole valley cross sections. In the MAP Riviera project (Weigel and Rotach, 2004) such an effect was invoked to explain the anomalous down-slope flow on a sunny sidewall. Different from that case, where the dynamically driven secondary flow was acting to contrast the thermally-driven slope-flow, here it would be in favor of the cross-valley circulation produced by the thermal asymmetry.

Conclusions
ABL structures in an Alpine Valley were investigated by RK interpolation of airborne and surface measurements at selected valley cross sections.
The results cast new light on the thermal structure associated with the development of the Ora del Garda wind. In particular, pseudo-soundings from airborne data allowed to identify a well-defined common vertical structure for the valley ABL, including a shallow ML surmounted by a deep stable layer. The RK-reconstruction of 3D fields revealed further local features, not detectable from vertical profiles only, associated with complex terrain and surface inhomogeneities. In particular, the effects of the differential heating of the valley sidewalls, and/or of topography irregularities, such as the bending of the valley axis, were outlined. Even small-scale effects determined by a small lake (Cavedine Lake) on the thermal structure were detected.
RK appears to be a promising technique for the investigation of 3D structures, not only of potential temperature, but also of other atmospheric variables, such as concentrations of water vapor, air pollutants and other species, as well as wind velocity and turbulent quantities. These analyses could usefully complement results from surface measurements for a more complete characterization of typical features concerning turbulent fluxes de Franceschi et al ., 2009), pollutant transport (de Franceschi and, atmospheric convection and associated precipitation (Bertò et al ., 2004) and urban climate (Giovannini et al ., 2011(Giovannini et al ., , 2013 in valley environments. Compared with other methods available in the literature, kriging presents a series of advantages: its implementation implicitly accounts for the characteristic covariance and for the intrinsic anisotropy of the target field, and it also naturally provides an estimate of the interpolation error. Furthermore RK allows to merge in a combined analysis both surface and airborne observations. Therefore, it is a suitable technique for integration of data from different surface measurements (e.g. at different locations: valley floor and sidewall slopes), as well as from remote sensing observations. Finally, RK-interpolated fields can also provide a suitable basis to validate high-resolution numerical simulations, in view of a more comprehensive investigation of processes involved in the development of local winds and related ABL processes.