Functional derivatives applied to error propagation of uncertainties in topography to large-aperture scintillometer-derived heat fluxes

Scintillometer measurements allow for estimations of the refractive index structure parameter $C_n^2$ over large areas in the atmospheric surface layer. Turbulent fluxes of heat and momentum are inferred through coupled sets of equations derived from the Monin--Obukhov similarity hypothesis. One-dimensional sensitivity functions have been produced that relate the sensitivity of heat fluxes to uncertainties in single values of beam height over flat terrain. However, real field sites include variable topography. We develop here, using functional derivatives, the first analysis of the sensitivity of scintillometer-derived sensible heat fluxes to uncertainties in spatially distributed topographic measurements. Sensitivity is shown to be concentrated in areas near the center of the beam path and where the underlying topography is closest to the beam height. Relative uncertainty contributions to the sensible heat flux from uncertainties in topography can reach 20\,\% of the heat flux in some cases. Uncertainty may be greatly reduced by focusing accurate topographic measurements in these specific areas. A~new two-dimensional variable terrain sensitivity function is developed for quantitative error analysis. This function is compared with the previous one-dimensional sensitivity function for the same measurement strategy over flat terrain. Additionally, a~new method of solution to the set of coupled equations is produced that eliminates computational error.

Abstract. Scintillometer measurements allow for estimations of the refractive index structure parameter C 2 n over large areas in the atmospheric surface layer. Turbulent fluxes of heat and momentum are inferred through coupled sets of equations derived from the Monin-Obukhov similarity hypothesis. One-dimensional sensitivity functions have been produced that relate the sensitivity of heat fluxes to uncertainties in single values of beam height over flat terrain. However, real field sites include variable topography. We develop here, using functional derivatives, the first analysis of the sensitivity of scintillometer-derived sensible heat fluxes to uncertainties in spatially distributed topographic measurements. Sensitivity is shown to be concentrated in areas near the center of the beam path and where the underlying topography is closest to the beam height. Relative uncertainty contributions to the sensible heat flux from uncertainties in topography can reach 20 % of the heat flux in some cases. Uncertainty may be greatly reduced by focusing accurate topographic measurements in these specific areas. A new two-dimensional variable terrain sensitivity function is developed for quantitative error analysis. This function is compared with the previous one-dimensional sensitivity function for the same measurement strategy over flat terrain. Additionally, a new method of solution to the set of coupled equations is produced that eliminates computational error.

Introduction
Large-aperture scintillometers infer the index of refraction structure parameter C 2 n over large areas of terrain in the atmospheric surface layer. The structure parameter for temperature C 2 T is resolved, and this information solves for the sensible heat flux H S through the application of equations derived from the Monin-Obukhov similarity hypothesis (Hartogensis et al., 2003;Moene, 2003). The sensible heat flux in the atmospheric surface layer is given by where ρ is the density of air, c p is the heat capacity at constant pressure, u is the friction velocity, and T is the temperature scale (e.g., Monin and Obukhov, 1954;Obukhov, 1971;Sorbjan, 1989;Foken, 2006). The temperature scale T is resolved by where z eff is the effective beam height above the ground, ζ ≡ z eff / l, where l is the Obukhov length (e.g., Sorbjan, 1989), and a, b and c are empirical parameters. The values of the empirical parameters are taken to be a = 4.9, b = 6.1, and c = 2.2, as seen in Andreas (1989) after an adjustment from the original values seen in Wyngaard et al. (1971).
These values may not be appropriate for all field sites. We will assume that C 2 T is resolved by neglecting the influence of humidity fluctuations, although this does not influence our results.
As can be imagined from Eqs. (2) and (3), it is important to know the height z at which C 2 T is being sampled; this corresponds to the scintillometer beam height. The beam height usually varies along the beam path. Even if turbulence is being sampled above an extremely flat field, uncertainty in z will still be present. Previous studies such as Andreas (1989) and Hartogensis et al. (2003) have quantified the sensitivity of H S to uncertainties in z over flat terrain. It is the goal of this study to extend the theoretical uncertainty analysis of Andreas (1989) and Hartogensis et al. (2003) to take into account variable terrain along the path. The value of this is in the ability to evaluate uncertainty estimates for scintillometer measurements over variable terrain, as well as to study the theoretical effect that the underlying terrain has on this uncertainty.
The studies of Andreas (1989) and Hartogensis et al. (2003) assume an independently measured friction velocity u . With large-aperture scintillometers, u may be inferred through the Businger-Dyer relation of wind stress, which is coupled to the Monin-Obukhov equations (e.g., Hartogensis et al., 2003;Solignac et al., 2009). Alternatively, with displaced-beam scintillometers, path-averaged values of the inner-scale length of turbulence l o can be measured (in addition to C 2 n ), which are related to the turbulent dissipation rate , which is in turn related through coupled Monin-Obukhov equations to u (e.g., Andreas, 1992). As a first step towards a variable terrain sensitivity analysis for large-aperture scintillometers, we will assume independent u measurements such that the Businger-Dyer equation will not be considered. Additionally, in order to take into account thick vegetation, the displacement distance d is often introduced. We will not consider this for the purposes of this study.
We are thus considering a large-aperture scintillometer strategy with independent u measurements as in Andreas (1989) and Appendix A of Hartogensis et al. (2003), and we consider the line integral effective beam height formulation from Hartogensis et al. (2003) and Kleissl et al. (2008). The effective height formulation is also discussed in Evans and De Bruin (2011) and in Geli et al. (2012). The assumptions behind this line integral approach are that the profile of C 2 T above the ground satisfies the Monin-Obukhov profile at any point along the beam path, and also that H S is constant vertically and horizontally within the surface layer region sampled by the beam. In this case, two coupled effects must be taken into account. Firstly, the scintillometer is most sensitive to fluctuations in the index of refraction towards the center of its beam. This is due to the optical configuration of the scintillometer system; a unit-less optical path weighting function takes this into account (e.g., Ochs and Wang, 1974;Hartogensis et al., 2003). The second effect is that, in areas where the topography approaches the beam, the C 2 T being sampled is theoretically more intense than in areas where the terrain dips farther below the beam.
In Sect. 2 of this paper, we define the sensitivity function S H S ,z (u) for the sensible heat flux H S as a function of variable topography z(u), where u is the relative path position along the beam. In Sect. 3, we solve for S H S ,z (u) for any general given z(u). In Sect. 4 we visualize the results by applying the resulting sensitivity function to the topography of a real field site in the North Slope of Alaska. We then apply the resulting sensitivity function to examples of synthetic beam paths. In Sect. 5 we discuss our results, and we conclude in Sect. 6.

Definition of the sensitivity function S H S ,z (u)
Under stable conditions (ζ > 0), the set of equations to consider consists of Eqs. (1) and (3), as well as where z eff is derived in Kleissl et al. (2008) based on the theory from Hartogensis et al. (2003), z(u) is the height of the beam along the relative path position u, T is the temperature, G(u) is the optical path weighting function, g is gravitational acceleration, and κ is the von Kármán constant (0.4). For unstable conditions (ζ < 0), Eqs. (1), (2) and (4) are considered, but Eq. (5) is replaced by where z eff is derived in Hartogensis et al. (2003). The propagation of uncertainty from measurements such as z(u) to derived variables such as H S will be evaluated in the context of the inherent assumptions behind the theoretical equations. A standard approximation (e.g., Taylor, 1997) to the uncertainty in estimating the derived variable f = f (µ), µ = (µ 1 , µ 2 , . . . , µ N ), byf = f (x), a function of measurement variables x = (x 1 , x 2 , . . . , x N ), is The numerical indices indicate different independent (measurement) variables, such as T , P , C 2 n , u , and beam wavelengths λ and z. It is convenient to re-write Eq. (7) as where the sensitivity functions Sensitivity functions such as these are developed in Andreas (1989) and Andreas (1992). They are each a measure of the portion of relative error in a derived variable f resulting from a relative error in the individual measurement variable x i . The problem of resolving the uncertainty in the derived variables is a matter of identifying the magnitude and character of the measurement uncertainties, and then solving for the partial derivative terms in Eqs. (7) and (9).
Here we seek a solution to the sensitivity function of sensible heat flux as a function of topography S H S ,z . In the flat terrain case, the sensitivity function S H S ,z has a single component, corresponding to the single measurement variable z (Andreas, 1989). In our situation, however, we may imagine that since z(u) is distributed over one dimension instead of a single value of z, S H S ,z will be composed of a spectrum of components: We are thus aiming to expand the sensitivity function denoted "S z " in Fig. 4 of Andreas (1989) Fig. 8) from one dimension to infinitely many, owing to the fact that some derived variables such as z eff are functions of an integral over continuous variables z(u) and G(u) (we consider for generality that z(u) has a continuous uncertainty σ (u) 2 ). In other words, Being dependent on a continuum of measurement variables, the sensitivity function S H S ,z (u) here requires the calculation of a so-called functional derivative, δz eff /δz(u) (e.g., Courant, 1953;Greiner and Reinhardt, 1996). Functional derivatives have a long history of application to statistical error analysis (e.g., Fernholz, 1983;Beutner, 2010, and many references therein).
For our purposes, a heuristic derivation of δz eff /δz(u) results from an interpretation of the integral in z eff as the limit of Riemann sums. That is, where subscript i indicates that u = (i/N ). Upon discretizing the input variables, we have Letting k = arg min k |z(u) − z k | and taking the limit N → ∞, the desired functional derivative is given by We thus define as the sensitivity function of sensible heat flux H S to uncertainties in variable topography z(u). It is our goal to evaluate Eq. (14).
3 Solution of the sensitivity function S H S ,z (u)

Stable conditions (ζ > 0)
Under stable conditions, the set of Eqs.
(1), (3), (4) and (5) is coupled in l through ζ ; we begin de-coupling them by combining Eqs. (3) and (4) to obtain Since ζ > 0, the unsolved sign is positive. With the substitution we re-arrange Eq. (15) to obtain where z eff in the stable case is determined by a priori known functions z(u) and G(u) through Eq. (5). The value ofˆ , including C 2 T , is directly determined from the measurements. The solution of Eq. (17) follows by re-writing it as a fourthdegree algebraic equation in ζ 2/3 : or more practically, it can be solved through fixed-point recursion on the function where we must consider the positive root. Note that since Eq. (18) is fourth degree, Galois theory states that it has an explicit solution form (e.g., Edwards, 1984). It is thus possible in theory to write H S = h(z(u), C 2 n , P , T , λ, u ), where h is an explicit function of the measurements; however, it would be quite an unwieldy equation.
We do not need an explicit solution in order to study the sensitivity; we can use the chain rule and implicit differentiation as in Gruber and Fochesatto (2013). We establish the variable inter-dependency using Eq. (17) as a starting point. The tree diagram for any set of measurements under stable conditions is seen in Fig. 1. The measurements are at the ends of each branch, and all other variables are dependent.
The required global partial derivatives are now defined through the variable definitions, the above equations, and the tree diagram. We have We will need one derivative that we are not able to retrieve directly from explicit definitions. By implicitly differentiating Eq. (17) under the guidance of the tree diagram seen in Fig. 1, we derive The functional derivative term δz eff δz(u) for stable conditions has been evaluated in Eq. (13).

Unstable conditions (ζ < 0)
Under unstable conditions, the set of Eqs. (1), (2), (4) and (6) is coupled in l through ζ ; note that z eff is coupled to ζ in the unstable case. We combine Eqs. (2) and (4) to obtain Since ζ < 0, the sign is negative. With the substitution , this leads to We substitute Eq. (23) into Eq. (6) to obtain This single equation is in the single unknown ζ , since z(u), G(u) and˘ are known; it is also in the fixed-point form ζ = F (ζ ). The tree diagram for the unstable case is seen in Fig. 2. Evaluation of global partial derivatives proceeds analogously to the stable case as in Eq. (20). Now we have To pursue the solution of S H S ,z (u), we will need to solve for ∂z eff ∂ζ by the differentiation of Eq. (23): We can solve for δζ δz(u) by implicit differentiation of Eq. (24). In finding δζ δz(u) , it is useful to define where, from Eqs. (24) and (27), we have From Eq. (27), we have such that, by implicitly differentiating Eq. (28) and then substituting, we derive All the information we need to solve for S H S ,z (u) is now resolved.

Full expression for the sensitivity function S H S ,z (u)
Since we are considering an independent u measurement, we have S T ,z (u) = S H S ,z (u) = z(u) T δT δz(u) . We obtain

Application of the results for the sensitivity function
S H S ,z (u)

Imnavait Creek basin field campaign
As an example, we use topographic data from the Imnavait Creek basin field site (UTM 5N 650220.5 East, 7615761.5 North), where there is a campaign to determine large-scale turbulent fluxes in the Alaskan tundra; it is seen in Figs. 3a and 4. We assume for simplicity that vegetation patterns, water availability, and other changes across the basin that could affect the flow in the atmospheric surface layer do not represent a significant source of surface heterogeneity. The elevation data seen in Fig. 3a are from a 5 m resolution digital elevation map (DEM), which has a roughly 0.5 m standard deviation in a histogram of the difference between the DEM elevations and 50 randomly distributed GPS ground truth points, as seen in Fig. 3b. Note that the systematic offset between the DEM and the GPS ground truth measurements does not contribute to systematic error in z(u). Note also that some of this spread in data may be due to an active permafrost layer. For this field site, we can solve for ζ under unstable conditions through Eq. (24). As can be seen in Fig. 5, we arrive at the solution for ζ with the recursively defined series [F (ζ guess ),F (F (ζ guess )),F (F (F (ζ guess ))), . . . ] that is guaranteed to converge monotonically for any ζ guess < 0.
A plot of ζ as a function of˘ for this field site is seen in Fig. 6. Note that the relationship between ζ and˘ is bijective; any value of˘ is uniquely associated with a value of ζ .
Considering the field case study of the Imnavait Creek basin, where the height of the beam over the terrain z(u) and the standard path weighting function G(u) are seen in Figs. 3a and 4, Eqs. (31) and (32) lead to the sensitivity function seen in Fig. 7. Note that S H S ,z (u) is a function of u and ζ only, since, for any one beam height transect z(u),˘ is  mapped bijectively with respect to ζ through Eq. (24), as seen in Fig. 6. Note that if we consider a constant ratio of σ z (u) z(u) , systematic error propagation can be re-written as The term in square brackets on the right of Eq. (33) is plotted in Fig. 8.

Synthetic scintillometer beam paths
It is interesting to examine the sensitivity function over synthetic paths that are representative of commonly used paths in scintillometry. Two synthetic paths can be seen in Fig. 9. They include a slant path as well as a quadratic path representing a beam over a valley.
The sensitivity function S T ,z (u) = S H S ,z (u) for synthetic path 1 (the quadratic path) seen in Fig. 9 is seen in Fig. 10. For synthetic path 2 (the slant path), the sensitivity function is seen in Fig. 11.   Atmos. Meas. Tech., 7, 2361-2371, 2014 www.atmos-meas-tech.net/7/2361/2014/   ) and (24) produced with a monotonically converging series as explained in the text and as visualized in Fig. 5. The Imnavait Creek basin terrain and beam path are used for z(u), along with the standard path weighting function G(u) as seen in Figs. 3a and 4. The mapping between ζ and˘ and between ζ and is bijective. Note that the solution of ζ for˘ = 1/4 corresponds to the intersection ofF with ζ in Fig. 5.

Discussion
A sensitivity function mapping the propagation of uncertainty from z(u) to H S has been produced for a large-aperture scintillometer strategy incorporating independent u measurements, and the line integral footprint approach to variable topography developed in Hartogensis et al. (2003) and Kleissl et al. (2008). This was accomplished by mapping out the variable inter-dependency as illustrated in the tree diagrams in Figs. 1 and 2, and by applying functional derivatives. The solution to S H S ,z (u) is given in Eqs. (14), (31) and (32).
As seen in Figs. 3a, 4, and 7, our results for S T ,z (u) = S H S ,z (u) show that sensitivity to uncertainties in topographic heights is generally higher under unstable conditions, and it is both concentrated in the center of the path and in areas where the underlying topography approaches the beam height. This finding intuitively makes sense, since scintillometers are more sensitive to C 2 T at the center of their beam path, and C 2 T decreases nonlinearly in height above the surface and strengthens with greater instability. For the Imnavait Creek basin path, the value of S H S ,z (u) increases to 3 at small dips in the beam height beyond the halfway point of the path, as seen in Fig.7. Note that the asymmetry along u of S H S ,z (u) corresponds to the asymmetry of the path, which is mostly at a higher (> 6 m) height in the first half, and at a lower height (≈ 4 m) in the second half, as seen in Fig. 4. Also note that the local maxima in S H S ,z (u) occur at roughly u ≈ 60 % and u ≈ 65 %; these correspond directly to topographic protuberances seen in Figs. 3a and 4. Note that the total error in H S is contributed from the whole range of u along S H S ,z (u), so even though we may have values of up to 3 in the sensitivity given by 1 0 S T ,z (u)du, and the flat terrain sensitivity function S z derived in Andreas (1989) (for ζ > 0, the functions are identical). For stable conditions (ζ > 0), S T ,z (u) is given by Eq. (31). For unstable conditions (ζ < 0), S T ,z (u) is given by Eq. (32), where values for ζ as a function of˘ are obtained through a numerical solution of Eq. (24), which may be visualized in Fig. 6. The Imnavait Creek basin terrain and beam path are used for z(u), along with the standard path weighting function G(u) as seen in Figs. 3a and 4. functions, our error bars may still be reasonable. The average value of S H S ,z (u) along u is never higher than 1, as seen in Fig. 8. Knowledge of where the concentration in sensitivity is allows us to decrease our uncertainty greatly by taking high-accuracy topographic measurements in these areas, especially for Arctic beam paths, which must be low due to thin boundary layers.
For example, if the random error in z(u) in the Imnavait Creek basin were 0.5 m, the relative error resulting in H S due to uncertainty in z(u) alone would be just 2 % under slightly unstable conditions where˘ = 1/4 and ζ ≈ −3.75, whereas if we reduce the uncertainty in z(u) to 0.1 m, the relative error in H S due to uncertainty in z(u) would be just 0.3 %, so with a reasonable number of survey points (100), the error can be quite small. However, if we look at Fig. 3b, we see that there is significant systematic error, perhaps due to shifting permafrost. If we have a perfectly even systematic error across the whole map, then this error is not propagated. However, if we have even a small amount of systematic error such as 0.5 m distributed around the center of the beam path near the local maxima in sensitivity, we can easily achieve 10 % to 20 % relative error in H S . In comparison to other variables, the values for S H S ,u are similar in magnitude to S H S ,z under unstable conditions, smaller under neutral conditions, and larger under stable conditions (Andreas, 1989). Under unstable conditions, error from u may therefore be similar in magnitude to error from z(u); however, for path-averaged u scintillometer strategies, this is not an issue. For C 2 n , the sensitivity functions are usually smaller, but in isolated regions they are larger (Andreas, 1989).
The average value of S H S ,z (u) over the beam path reduces to identical results to the flat terrain sensitivity function S z from Andreas (1989) (which would be denoted S T ,z here) under stable conditions where z eff is de-coupled from ζ , and nearly identical results (depending on the path) under unstable conditions where z eff is coupled to ζ , as seen in Fig. 8. It is unknown as to whether the addition of equations for path-averaged u measurements such as the Businger-Dyer relation seen in Hartogensis et al. (2003) and Solignac et al. (2009), or displaced-beam scintillometer strategies as seen in Andreas (1992), would change these results significantly.
We note that the study of Hartogensis et al. (2003) evaluated a function similar to S H S ,z for flat terrain with an independent u measurement (the 2003 Eq. 7 is ignored); however, at ζ ≈ 0 they found a sensitivity of 1/2 instead of 1/3 as found in Andreas (1989). The difference in the results between these two studies is not due to the differences between single-and double-wavelength strategies. The Obukhov length (denoted by L MO in Hartogensis et al., 2003) is a function of z LAS through the 2003 Eqs. (5) and (6). The addition of chain rule terms to reflect the dependence of l on z in Hartogensis et al.'s (2003) Eq. (A2) resolves differences between Hartogensis et al.'s (2003) Fig. A1 and Andreas et al.'s (1989) Fig. 4; the flat-terrain sensitivity function for ζ < 0 is which is given correctly in Andreas (1989). Equations (7), (9), (31), and (32) may be implemented into computer code for routine analysis of data. It is worth noting that the sign of ζ is an a priori unknown from the measurements. Thus, for any set of measurements, we should calculate the set of all derived variables and their respective uncertainties assuming both stable and unstable conditions, and if uncertainties in the range of ζ overlap with ζ = 0 for either stability regime, we should then consider the combined range of errors in the two sets.

Atmos
In the application of Eq. (7), we must recognize computational error σ f c . Previous studies have incorporated a cyclically iterative algorithm that may not converge, as seen in Andreas (2012), or that may converge to an incorrect solution, as illustrated in the section on coupled nonlinear equations in Press et al. (1992). We have developed techniques to eliminate this error. For unstable cases (ζ < 0), the solution of ζ follows from Eq. (24), which is in fixed-point form. The solution to Eq. (24) is guaranteed to converge monotonically with the recursively defined series [F (ζ guess ),F (F (ζ guess )),F (F (F (ζ guess ))), . . . ] as seen in Traub (1964) and in Agarwal et al. (2001), and as demonstrated in Fig. 5. We may solve for the stable case (ζ > 0) recursively using Eq. (19), whereF (ζ ) demonstrates convergence properties that are similar to those ofF (ζ ) in Eq. (24). It was found to be practical to make ζ guess = ±1.
Future expansions of the results presented here should focus on including multiple wavelength strategies to evaluate the latent heat flux and H S , as well as on including path-averaged u measurements using l o and C 2 n scintillometer strategies as in Andreas (1992) or using a point measurement of wind speed and the roughness length via the Businger-Dyer relation (e.g., Panofsky and Dutton, 1984; Figure 11. Sensitivity function S H S ,z (u) = S T ,z (u). For stable conditions (ζ > 0), S T ,z (u) is given in Eq. (31). For unstable conditions (ζ < 0), S T ,z (u) is given by Eq. (32), where values for ζ as a function of˘ are obtained through a numerical solution of Eq. (24), which may be visualized with Fig. 6. Synthetic beam path 2 (the slant path) is used for z(u), along with the standard path weighting function G(u) as seen in Figs. 9 and 4. Solignac et al., 2009). Modification of the analysis for including path-averaged u measurements involves the addition of one or two more equations (e.g., Eq. 8 in Solignac et al., 2009, or Eqs. 1.2 and 1.3) in Andreas, 1992) to substitute into Eqs. (17) and (24), as well as the definition of new tree diagrams to reflect that u is now a derived variable. In these cases, either the turbulence inner-scale length l o or a point measurement of wind speed and the roughness length replaces u as a measurement; u is derived through information from the full set of measurements. Note that if u is derived through measurements including z, Eq. (1) implies that S H S ,z = S T ,z + S u ,z . It is worth investigating whether computational error can still be eliminated in these cases.
We have considered here the effective height line integral approach derived in Hartogensis et al. (2003) and Kleissl et al. (2008) to take into account variable topography. Even if we assume a constant flux surface layer, under realistic wind conditions, turbulent air is advected in from nearby topography. For example, in the Imnavait Creek basin path seen in Fig. 3a, if wind comes from the west, the turbulent air being advected into the beam path comes from a volume that is higher above the underlying topography than if wind came from the east. Sensitivity studies should be produced for two-dimensional surface integral methods that take into account the coupling of wind direction and topography on an instrument footprint (e.g., Meijninger et al., 2002;Liu et al., 2011). Additionally, a new theory may be developed for heterogeneous terrain involving complex distributions of water availability and roughness length such as the terrain in the Imnavait Creek basin.

Conclusions
Sensitivity of the sensible heat flux measured by scintillometers has been shown to be highly concentrated in areas near the center of the beam path and in areas of topographic protrusion. The general analytic sensitivity functions that have been evaluated here can be applied for error analysis over any field site as an alternative to complicated numerical methods. Uncertainty can be greatly reduced by focusing accurate topographic measurements in areas of protrusion near the center of the beam path. The magnitude of the uncertainty is such that it may be necessary to use high-precision LIght Detection And Ranging (LIDAR) topographic data as in Geli et al. (2012) for Arctic field sites in order to avoid large errors resulting from uneven permafrost changes since the last available DEM was taken. Additionally, computational error can be eliminated by following a computational procedure as outlined here.