Sensitivity of large-aperture scintillometer measurements of area-average heat fluxes to uncertainties in topographic

Introduction Conclusions References

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 (3 where z is the height above the ground, ζ ≡ z/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 Wyngaard et al. (1971) and in Andreas (1989).These values may not be appropriate for all field sites.As can be surmised from Eqs.
(3) and (2), 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 and homogeneous 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 of 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 Figures

Back Close
Full 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.
In using the Monin-Obukhov similarity hypothesis we are assuming that the flow is stationary and that the terrain is homogeneous.While the topography is not flat, we will assume that it is nearly so and that the surface features are horizontally homogeneous.Heterogeneous terrain implies spatial gradients in fluxes; in this case many authors make the assumption that the scintillometer beam is above the blending height where gradients in fluxes are negligible (e.g., Wieringa, 1986;Mason, 1987;Claussen, 1990Claussen, , 1995;;Meijninger et al., 2002;Hartogensis et al., 2003;Lu et al., 2009).Sensitivity studies have so far been restricted to single values of beam height as in Andreas (1989) and in Appendix A of Hartogensis et al. (2003).Hartogensis et al. (2003) anticipated the quantification of sensitivity in H S to variable topography for a largeaperture scintillometer strategy with independent u measurements.
We are thus considering a large-aperture scintillometer strategy with independent u measurements as in Andreas (1989) and in 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 Eqs. ( 3) and ( 2) at any point along the beam path where z is taken to be the local height of the beam above the underlying terrain, 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 Introduction

Conclusions References
Tables Figures

Back Close
Full 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 normalized distance along the beam path.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 the effective beam height 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 normalized Introduction

Conclusions References
Tables Figures

Back Close
Full path length u, G(u) is the optical path weighting function, g is gravitational acceleration, and κ is the von Kármán constant.
The propagation of uncertainty from source 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 good estimate of the uncertainty in the derived variables that results from small errors in source measurements is given by where derived variable f is a function of source measurement variables x 1 , x 2 , . . ., x N with respective systematic error σ x s 1 , σ x s 2 , . . ., σ x s N and with respective independent Gaussian distributed uncertainties with standard deviations σ x r 1 , σ x r 2 , . . ., σ x r N as seen in Taylor (1997).The numerical indices indicate different independent variables, such as temperature T , pressure P , or z, for example.Computational error on f due to the inaccurate solution of the theoretical equations is represented by σ f c .The first and last terms in Eq. ( 9) represent an offset from the true solution (inaccuracy), whereas the Introduction

Conclusions References
Tables Figures

Back Close
Full It is convenient to re-write Eq. ( 9) as where the sensitivity functions S f ,x are defined as 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 relative error on the individual source measurement variable x.
The problem of resolving the uncertainty on the derived variables is a matter of identifying the magnitude and character of the source measurement uncertainties, and then solving for the partial derivative terms in Eqs. ( 9) and ( 11).
We seek a solution to the sensitivity function of sensible heat flux as a function of topography S H S ,z .The sensitivity function S H S ,z is a function only of ζ in the flat and homogeneous terrain case as seen in Andreas (1989).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 a function of both ζ and u.We are thus aiming to expand S z seen in Fig. 4 of Andreas (1989) (our S H S ,z in Fig. 8) from one dimension to two.This extra dimension will come from the fact that some derived variables such as z eff are functions of an integral over continuous variables z(u) and G(u), where we consider for generality that z(u) has a continuous Introduction

Conclusions References
Tables Figures

Back Close
Full uncertainty σ z (u).To illustrate this we re-write for example Eq. ( 6) as where subscript i indicates that u = (i /N).The propagation of error defined by Eq. ( 9) involves derivatives of the dependent variables as functions of all the independent variables, where each z i is independent.For such a partial derivative term we have where δ i k is the Kronecker delta, and 1 ≤ k ≤ N.This anticipates the definition of a new differential operator that is the same as a normal partial derivative in all steps including chain rule, product rule, etc., except upon the application of the Leibniz rule while differentiating the primary variable, whereupon we multiply the integrand by a Dirac delta function.We can name this operator the "Dirac-Leibniz" derivative, and we denote it here as with respect to x.Using continuous notation and again differentiating Eq. ( 12), we can write Note that ; it can also be derived that .
We thus define Under stable conditions the set of Eqs. ( 1)-( 6) is coupled in l through ζ ; we begin de-coupling them by combining Eqs. ( 4) and ( 5) to obtain Since ζ > 0, the unsolved sign is positive.With the substitution we re-arrange Eq. ( 16) to obtain where z eff in the stable case is determined by a-priori known functions z(u) and G(u) through Eq. ( 6).The value of Λ, including C 2 T , is directly determined from the source measurements.The solution of Eq. ( 18) follows by re-writing it as a fourth degree algebraic equation in ζ 2/3 : or more practically, it can be solved through fixed point recursion on the function

Conclusions References
Tables Figures

Back Close
Full where we must consider the positive root.Note that since Eq. ( 19) 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 source 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. ( 18) as a starting point.The tree diagram for any set of source measurements in stable conditions is seen in Fig. 1.The source 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.For example, we have We will need some derivatives that we are not able to directly retrieve from explicit definitions.By implicitly differentiating Eq. ( 18) under the guidance of the tree diagram seen in Fig. 1, we derive The Dirac-Leibniz derivative term δz eff δz (u) for stable conditions has been evaluated in Eq. ( 14).

Conclusions References
Tables Figures

Back Close
Full , this leads to We substitute Eq. ( 24) into Eq.( 8) 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. ( 21).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. ( 24): We can solve for δζ δz (u) by implicit differentiation of Eq. ( 25).In finding where, from Eqs. ( 25) and ( 28), we have From Eq. ( 28), we have that such that, from Eqs. ( 25), (28), and (29), we derive

Conclusions References
Tables Figures

Back Close
Full All the information we need to solve for S H S ,z (u) is now resolved.Introduction

Conclusions References
Tables Figures

Back Close
Full

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

Imnavait creek basin field campaign
As an example we use topographic data from the Imnavait Creek Basin field site 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 is 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 ζ in unstable conditions through Eq. ( 25).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 ζ for ζ < 0.
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. ( 33) and ( 34) lead to the sensitivity function seen in Fig. 7.Note that S H S ,z (u) is a function of only u and ζ , since, for any one beam height transect z(u), Λ is mapped bijectively with respect to ζ through Eq. ( 25), as seen in Fig. 6.Introduction

Conclusions References
Tables Figures

Back Close
Full Note that if we consider a constant ratio of z(u) , the term in, for example Eq. ( 10), can be re-written as The term in square brackets on the right of Eq. ( 35) 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.

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 the Dirac-Leibniz derivative.The solution to S H S ,z (u) is given in Eqs. ( 15), ( 33) and ( 34).
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 in unstable conditions, and it is Introduction

Conclusions References
Tables Figures

Back Close
Full both concentrated in the center of the path and in areas where the underlying topography approaches the beam height.This finding intuitively makes sense for two reasons: firstly scintillometers are more sensitive to C 2 T at the center of their beam path, and secondly, C 2 T decreases nonlinearly in height above the surface and it 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 concentrations 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 on 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 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 greatly decrease our uncertainty by taking high precision topographic measurements in these areas, especially for Arctic beam paths, which must be low due to thin boundary layers.Note also that S H S ,z (u) is not analogous to a footprint or to a path weighting function; the scintillometer is still sensitive to C 2 T along the whole path.S H S ,z (u) should not be interpreted beyond being a measure of how uncertainties in topographic measurements propagate through to the derived sensible heat flux.
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) in stable conditions where z eff is de-coupled from ζ , and nearly identical results (depending on the path) in 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.Introduction

Conclusions References
Tables Figures

Back Close
Full 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 (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 results between these two studies is not due to the differences between single and double wavelength strategies.The Obukhov length (denoted by Hartogensis et al., 2003) is a function of z LAS through 2003 Eqs. ( 5) and ( 6).The addition of chain rule terms to reflect the dependence of l on z in Hartogensis et al. (2003) Eq. (A2) resolves differences between Hartogensis et al. (2003) Fig. A1 and Andreas (1989) Fig. 4; the flat terrain sensitivity function for ζ < 0 is which is given correctly in Andreas (1989).Eqs. ( 9), ( 11), (33), and (34) 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 source measurements.Thus, for any set of source 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 on the two sets.
In the application of Eq. ( 9), 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, which 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. ( 25), which is in fixed point form.The solution to Eq. ( 25) 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) Introduction

Conclusions References
Tables Figures

Back Close
Full recursively using Eq. ( 20), where F (ζ ) demonstrates convergence properties that are similar to those of F (ζ ) in Eq. ( 25).It was found 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 including pathaveraged 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;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,1.3) in Andreas, 1992) to substitute into Eqs.( 18) and ( 25), 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 replaces u as a source measurement; u is derived through information from the full set of source measurements.Note that if u is derived through source 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 in 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 is coming 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 instrument footprint (e.g., Meijninger et al., 2002;Liu et al., 2011).Additionally, new theory may be developed for heterogeneous terrain involving complex distributions of water availability and roughness length such as the terrain in Imnavait Creek Basin.Introduction

Conclusions References
Tables Figures

Back Close
Full We will need some derivatives that we are not able to directly retrieve from explicit definitions.
By implicitly differentiating Eq. ( 18  We substitute Eq. ( 24) into Eq.( 8) to obtain This single equation is in the single unknown ζ since z(u), G(u) and Λ are known; it is also in ) and (4).
scintillometers measure 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 Figures Back Close Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | central square root term represents the breadth of uncertainty due to random error (imprecision).The source measurement variables being considered here are P , T , C 2 n , beam wavelength λ, z(u), and u .
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | δ δx (as opposed to ∂ ∂x ) when differentiating Figures Back Close Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 3 Solution of the sensitivity function S H S ,z (u) 3.1 Stable conditions (ζ > 0) Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 1 .
Fig. 1.Variable inter-dependency tree diagram for the stable case (ζ > 0).The source measurement variables are at the end of each branch; all other variables are derived.

Fig. 1 .
Fig. 1.Variable inter-dependency tree diagram for the stable case (ζ > 0).The source measurement variables are at the end of each branch; all other variables are derived.

Fig. 2 .
Fig. 2. Variable inter-dependency tree diagram for the unstable case (ζ < 0).The source measurement variables are at the end of each branch; all other variables are derived.

Fig. 5 .
Fig. 5. Graphical visualization of the fixed point solution of Eq. (25).The recursively defined series [ F (ζ guess ), F ( F (ζ guess )), F ( F ( F (ζ guess ))), . . .] converges monotonically for any ζ guess < 0. A typical value of Λ = 3/4 is used representing unstable conditions in the atmospheric surface layer.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.