###### Abstract

Scintillometer measurements of the turbulence inner-scale length and refractive index structure function allow for the retrieval of large-scale area-averaged turbulent fluxes in the atmospheric surface layer. This retrieval involves the solution of the non-linear set of equations defined by the Monin-Obukhov similarity hypothesis. A new method that uses an analytic solution to the set of equations is presented, which leads to a stable and efficient numerical method of computation that has the potential of eliminating computational error. Mathematical expressions are derived that map out the sensitivity of the turbulent flux measurements to uncertainties in source measurements such as . These sensitivity functions differ from results in the previous literature; the reasons for the differences are explored.

guess

## Section 1 Introduction

Scintillometers detect fluctuations in the intensity of a beam of light that passes through a path length of
50 m to 5000 m of near-ground turbulence in the surface layer [2008]. These fluctuations are related to
the structure function of the index of refraction , and the turbulence inner-scale length
[1961, 1988, 1994]. The index of refraction is a function of temperature and humidity; thus
can be decomposed into structure functions of temperature and humidity as , and . Scintillometer
wavelengths are selected that are each more sensitive to fluctuations in one variable (such as temperature) than others (such as
humidity), so that , and may be resolved.
For example, intensity fluctuations of visible and near-infrared beams are more sensitive to temperature fluctuations than humidity
fluctuations, while microwave beams are more sensitive to humidity fluctuations [1990].
Structure functions such as are described in
? (?), and represent the strength and spacial frequency of perturbations in variables;
thus is a measure of turbulence intensity weighted by the susceptibility of the index of refraction of the medium to changes in variables such as temperature and humidity.

The goal of this study is to solve for the sensible heat flux and the momentum flux as functions of source measurements such as and , as well as to quantify the propagation of uncertainty from source measurements to the calculated values of and . Another type of turbulent flux is the latent heat flux . The turbulent fluxes are given by

(\theequation) | |||||

(\theequation) | |||||

(\theequation) |

where and are the temperature and humidity scales, is the friction velocity,
is the density of the air, is the specific heat at constant pressure, and is the latent heat of vaporization.
Determining area-averaged turbulent fluxes involves solving for and , which are related to the path-length scale
structure-function measurements through the non-linearly coupled
Monin-Obukhov similarity equations [1989]. This procedure also involves solving for
in Eqs. 1, 1 and 1. The friction velocity
can be related either to path-length scale measurements
as with displaced-beam scintillometer strategies described in ? (?),
or to the wind profile and roughness length with large-aperture scintillometer strategies via the Businger-Dyer relation
[1984, 1989, 2002, 2003].

We consider here a displaced-beam scintillometer strategy in which path-averaged measurements of and
are obtained. Other required measurements
include temporally-averaged pressure , temperature , humidity , as well as the
height of the beam above the underlying terrain .
Thus , , , , and are referred to as the
source measurements. Each of these measurements demonstrates temporal and spacial variability as well as
measurement uncertainty. Uncertainty
propagates from the source measurements to the derived variables via the
set of equations being considered.
Uncertainties in and are described in ? (?), while uncertainties in , and depend on the particular instrument being used.
Here, we explore the use of scintillometers over flat and homogeneous terrain, thus the height of the beam is considered to be a single value with its associated uncertainty.
While and are representative of
turbulent fluctuations along the whole beam, , and are typically point measurements representative of localized areas near
their respective instruments.

Applications for scintillometers
include agricultural scientific studies such as ? (?) and ? (?), and aggregation of surface measurements to satellite-retrieval scales for weather
prediction and climate monitoring as in ? (?) and in ? (?).
The unique spacial scale of scintillometer measurements gives them the potential for a key role in bridging the gap between
ground-based instruments with footprints on the order of and model and satellite-retrieval scales on the
order of .

The scale of scintillometer measurements introduces an additional complexity in the retrieval of the turbulent fluxes.
This retrieval combines the large-scale scintillometer measured variables and with source measurements that are not necessarily
representative of the same scale. The only exception to be considered is the atmospheric pressure .
In particular, measurements of and may be representative of smaller footprints around their respective instruments.
Specifically, assuming that variables such as average temperature T represent the
entire beam path introduces a form of uncertainty.
This uncertainty is somewhat similar to a systematic error, although it may be difficult to
quantify because of its temporal variability.

Of previous scintillometer sensitivity studies, some stand out as possibly contradicting each other.
For instance, the conclusion of the error analysis in ? (?) for a and strategy was that
“The Monte Carlo analysis of the propagation of the statistical errors shows that there is only moderate sensitivity of the flux calculations to the initial errors in the measured quantities.”
The error analysis of ? (?), however, results in sensitivity functions that feature singularities.
The sensitivity functions presented there
imply that the resolution of and consequently of , and by scintillometer
and measurements
is intrinsically restricted to low precision over a certain range of environmental conditions.
While these two studies use different methods and present results over slightly different ranges in variables,
they produce sensitivity functions that for the same range differ significantly.

In Sect. 2 below, we decouple the set of equations including those of the Monin-Obukhov similarity hypothesis for and scintillometer strategies for the example of unstable surface-layer conditions to arrive at single equations in single unknowns. The variable inter-dependency is mapped out as illustrated by tree diagrams. In Sect. 3, we take advantage of the mapped out variable inter-dependency to guide us in using the chain rule to solve the global partial derivatives in sensitivity functions to investigate error propagation. We produce sensitivity functions for , and as functions of both and . In Sect. 4 we explore the ramifications of our results and compare them to previous literature, and we give conclusions in Sect. 5.

## Section 2 Measurement Strategy Case Study: Displaced-Beam Scintillometer System in Unstable Conditions

We consider here a two-wavelength system as introduced in ? (?), where one of the scintillometers measures both and as in ? (?).
With this strategy, our measurements can resolve humidity and temperature fluctuations separately since
the two scintillometers have different wavelengths and that have differing sensitivities in the index of refraction
to humidity and temperature. This technique therefore requires fewer assumptions than the
corresponding single-wavelength strategies as seen in ? (?).

The following set of equations determines , and from the source measurements, and subsequently determines the turbulent fluxes:

(\theequation) | |||||

(\theequation) | |||||

(\theequation) | |||||

(\theequation) | |||||

(\theequation) | |||||

(\theequation) |

where is the local acceleration due to gravity, is the Gamma function, is the turbulent energy dissipation rate,
is the specific gas constant, is the von Kármán constant, , where is the Obukhov length,
is the Obukhov-Corrsin constant, is the viscosity of air and is the thermal diffusivity of air
(Andreas, 1989; 1992; 2012)
and are structure functions of the refractive index for the separate wavelengths and .
Eqs. 2 and 2 determine directly from and the other source measurements.
Inherent in Eqs. 2 and 2 is the assumption that , which is validated previously [1989, 1990].

The similarity functions and are given by

(\theequation) | |||

(\theequation) |

for
which corresponds to unstable conditions. The form of the similarity functions and their parameters follow from ? (?) and ? (?); the values are
taken to be , , and [1988].

The source measurements may not determine the sign of , which is unknown a priori for every set of source measurements at any one time interval. We follow ? (?) in solving for and from Eqs. 2 and 2, making sure to note that the signs of are not yet solved by introducing unknowns and :

(\theequation) | |||||

(\theequation) |

where the roots on the left-hand side are considered to be positive. Following ? (?), these can be re-arranged to isolate and with the as yet undetermined signs:

(\theequation) | |||||

(\theequation) |

where

(\theequation) |

It is useful to include the definition of the Bowen ratio as

(\theequation) |

We can solve for as

(\theequation) |

where . It is useful to consider as well as as unit-less independent variables in our sensitivity analyses that represent certain meteorological
regimes. They represent the ratio of the
sensible to latent heat fluxes and an indicator of surface-layer stability, respectively.

Since we are considering unstable conditions, we have since , so from Eq. 2 we have

(\theequation) | |||

(\theequation) | |||

(\theequation) |

We begin decoupling the set of equations by taking Eqs. 2 and 2 and substituting into Eq. 2, then cubing the resulting equation as well as squaring Eq. 2 to arrive at

(\theequation) | |||||

(\theequation) |

where and are defined as

(\theequation) |

where

(\theequation) |

is determined directly from the source measurements. Here we note that the left-hand side is negative, and so the term in square brackets in is negative as well.
From any set of measurements we know the sign of , and we also know the values of the two terms that multiply the unknown signs.
Occasionally these relations are enough to determine all the signs; otherwise the signs remain ambiguous and they are evaluated from observations of the temperature and humidity stratification
as seen in ? (?).

Eq. 2 can be solved with a fixed-point recursive technique as illustrated in Fig. 1. The recursive function

(\theequation) |

is used.

A good estimate of the uncertainty in the derived variables that results from small errors in source measurements is given by

(\theequation) |

where the derived variable is a function of source measurement variables with respective systematic error and
with respective independent Gaussian distributed uncertainties with standard deviations as seen in ? (?).
The numerical indices indicate different independent variables, such as , , or , for example.
Computational error due to the inaccurate solution of the theoretical equations is represented by . The first and last terms in Eq. 2 represent an offset
from the true solution (inaccuracy), whereas the central square-root term represents the breadth of uncertainty due to random error (imprecision).

It is practical for the purpose of a sensitivity study to rewrite Eq. 2 as

(\theequation) |

where are unitless sensitivity functions defined by

(\theequation) |

The sensitivity functions are each a measure of the portion of the error
in the derived variable resulting from error on each individual source measurement . In addition to the error on source measurement variables, we
can also recognize
that , and have been resolved to some level of certainty by fitting field data.
We thus treat them here in the same way as source measurements.

In the application of Eqs. 2 and 2, we recognize the addition of the computational error .
In previous field and sensitivity studies [2002, 2002, 2009, 2012],
the full set of equations has been incorporated into
a cyclically iterative algorithm which cycles through the full set of equations, allowing multiple variables to change.
This numerical algorithm sometimes fails to converge, as demonstrated in ? (?).

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. 2 and 2. These derivatives are global^{1}^{1}1
Global partial derivatives are those which propagate from the dependent (derived) variable down to the independent (source measurement) variable through
the entire tree
diagram, whereas local partial derivatives propagate as if the equation being differentiated were independent of the rest of the equations in the set. An alternative to direct evaluation
of global partial derivatives via the chain rule is
a total-differential expansion (where all derivatives are local) of each equation in the set.
This approach can be used to solve for global partial derivatives by re-grouping all total-differential terms into one equation. Readers may refer to ? (?).
; that is, they take into account
all the relationships in all of the relevant equations through which the variable is derived. Without an analytic solution of the set of coupled equations
we could either solve for the partial derivatives
through a total-differential expansion of each equation individually, followed by a re-grouping of all differential terms as seen in
Andreas (1989; 1992)
or we could use numerical
error propagation techniques as in the
Monte Carlo analysis of ? (?) or as in the analysis of ? (?).

We investigate
inter-variable sensitivity analytically via Eq. 2, using Eq. 2 as a starting point. We use Eq. 2
to determine the details of the variable inter-dependency to define
our use of the chain rule.
A tree diagram representing the variable inter-dependency is broken into three parts shown in Figs. 3, 4, and 5.

Eq. 2 can be reduced to a choice of two algebraic equations

(\theequation) | |||

(\theequation) |

with the substitution

(\theequation) |

Galois theory implies that, since Eqs. \theequation and \theequation are ninth order, there is no way to write for any general values of and , where is an explicit function of the source measurements [1984]. It is thus simplest to extract by implicit differentiation of Eq. 2; the results are in given in Appendix A.

## Section 3 Results: Derivation of Sensitivity Functions

Following the solution method described above, we solve for global partial derivative terms in Eqs. 2 and 2 through use of the general chain rule guided by the variable inter-dependency tree diagrams seen in Figs. 3, 4 and 5. We will obtain sensitivity functions of the sensible heat flux and the momentum flux as functions of and . From Eqs. 1, 2 and 2 we have

(\theequation) | |||||

(\theequation) |

thus we seek solutions for , , , and .

We first obtain with guidance from the tree diagram depicted in Fig. 4:

(\theequation) |

The individual terms of Eq. \theequation are given in Appendices A and B. Combining them, we obtain

(\theequation) |

We now obtain :

(\theequation) |

The individual terms of Eq. \theequation are developed in Appendices A and C. Combining them, we obtain

(\theequation) |

We now obtain with guidance from the tree diagram depicted in Fig. 5. We have

(\theequation) |

The individual terms in Eq. \theequation are developed in Appendices A and D. Combining them, we obtain

(\theequation) |

We now obtain . We have

(\theequation) |

The individual terms in Eq. \theequation are developed in Appendices A and E. Combining them we obtain

(\theequation) |

Combining our results in Eqs. \theequation, \theequation, \theequation, and \theequation, we can obtain and from
Eqs. \theequation and \theequation; the results are seen in Fig. 6.

The absolute value of our results for given by Eqs. \theequation, \theequation and \theequation is similar to the sensitivity multiplier found in ? (?) as seen in their Fig. 10. The absolute value of our result of given by Eqs. \theequation and \theequation is also compatible with the results of ? (?) seen in their Fig. 9. However, our result for in Eq. \theequation differs from that obtained in ? (?) as seen in Fig. 7.

Similarly, our result for in Eq. \theequation differs from that obtained in ? (?) as seen in Fig. 8.

## Section 4 Discussion

The reason for the difference between our results and those of ? (?) in Figs. 7 and 8 can be seen to have arisen in Eqs. A.7 and A.10 of ? (?) . Even though there is a typographical error in Eq. A.7 in the application of the product rule (it should be

(\theequation) |

where the second term contained originally), this is not the origin of the reason since the result in Eq. A.8 follows from the modified Eq. A.7. The reason is found to be that Eqs. A.7 and A.8 are not differentiated locally with respect to Eq. 1.3 of ? (?) as they should be in a total-differential expansion. The local derivative is

(\theequation) |

keeping constant regardless of the relationship between and . The relationship between and is taken into account when we re-group the full set of locally expanded equations (which are coupled in and ). The second term on the right-hand side of Eq. \theequation and Eq. A.7 of ? (?) is thus not necessary and does not appear in Eq. \theequation. Taking into account the relationship between and via the chain rule is appropriate for direct evaluation of global derivatives, but not in individual derivatives of a total-differential expansion of the full set of equations. Eqs. A.10 and A.11 of ? (?) have the same issues of not being differentiated locally with respect to Eq. 1.3 of ? (?). The local derivative there is

(\theequation) |

A re-analysis of the ? (?) differential expansion including the local derivatives in Eqs. \theequation
and \theequation is reproduced in Appendix F; the results for and
are identical to those found here in Eqs. \theequation and \theequation.
Note that the left-hand side of Eq. \theequation contains the terms and instead of
and as in Eq. A.16 of
? (?). These differences also influence the ? (?) sensitivity functions for and .

The technique presented here for the direct evaluation of partial derivatives can be applied to evaluate sensitivity functions for other variables involved in this scintillometer strategy for both stable and unstable conditions, however we will now focus on the implications of our results on other previous studies. Another instance where we found divergence in results is in the study of ? (?) where in Eq. A2 and Fig. A1 should be the same as the results of ? (?) in Fig. 4, regardless of the differences between a single and double wavelength strategy. Note that in ? (?), for , it was found that

(\theequation) |

for a scintillometer strategy involving independent measurements, whereas a value of was found in ? (?). The issue here is not due to the differences in scintillation strategies (note that the Businger-Dyer relation is ignored in the sensitivity study of ? (?)). The issue is that Eq. A1 of ? (?) is coupled to Eqs. 5-6 of ? (?) in . In the derivation of Eq. A1, ? (?) essentially have considered to be the same as in ? (?), and they have considered similar equations that assume an independent measurement (Eq. 7 of ? (?) is ignored). Including the coupling of Eq. 7 of ? (?) (the Businger-Dyer relation) in adds complication; however if we continue to assume an independent measurement, we achieve the same results as in ? (?), viz:

(\theequation) |

A similar example is in the analysis of ? (?), when the sensitivity of to is being examined. Eq. 13 of ? (?) is not a “direct” relation of to source measurements, since is a derived variable. There is coupling to and thus we may investigate the sensitivity with

(\theequation) |

where is modified for the single scintillometer and strategy. Also in ? (?), it is stated that errors in are attenuated in deriving (here denoted ) due to the square-root dependence; however we can go a step further by realizing that Eq. 9 of ? (?) is not yet decoupled from . As follows from our analysis applied to the case considered in ? (?) (modifying Fig. 4 for a single-wavelength strategy), we obtain

(\theequation) |

Note that there may be no way to actually obtain “direct” relationships between the source
measurements and the derived variables if the implicit equation in (such as Eq. 2) is fifth order or higher.

## Section 5 Conclusions

A new method of deriving sensitivity functions for and scintillometer measurements of turbulent fluxes has been produced by mapping
out the variable inter-dependency and solving for partial derivatives with the chain rule. We have bypassed the need for an explicit solution
to the theoretical equations by including one implicit differentiation step on Eq. 2, which is
a bottleneck on the tree diagrams seen in Figs. 4 and 5.
This allows for the evaluation of
sensitivity functions that are useful not only for optimizing the measurement strategy and selecting the most ideal wavelengths,
but the closed, compact form of sensitivity functions produced using the method presented here is convenient
to incorporate into computer code for the analysis of data.
It is noteworthy that the actual functional relations change at , which corresponds to neutral conditions.
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.
If errors
on overlap with for either stability regime, we should then consider the combined range of errors.

In addition to the source measurements, the empirical parameters , and have been included in the tree diagrams.
Future study should quantify the sensitivity of derived variables to these parameters.
In considering errors on the empirical parameters or on other source measurements such as ,
a total-differential expansion such as in
Andreas (1989; 1992)
may become intractable, whereas an analysis of the type presented here remains compact.

Results obtained here have resolved some issues in the previous literature.
For example, we have confirmed the conclusion of ? (?) that and scintillometers can obtain fairly precise measurements of turbulent fluxes.
In the range of , the results derived here for and
are similar to those in ? (?);
however for the separate results differ greatly in both magnitude
and in the shape of the curves as seen in Figs. 7 and 8.
These sensitivity functions in ? (?) contain singularities near ; this effectively implies that it is impossible to resolve
in this stability regime. The sensitivity functions derived here demonstrate a small magnitude for typical values of including the range .
The sensitivities of the sensible heat flux to uncertainties in and are found in Eqs. \theequation and
\theequation and are seen in Fig. 6; they are compatible with the results of ? (?) and they imply that, with optimal wavelengths, we can arrive at
reasonably precise measurements of path-averaged turbulent fluxes and friction velocity.

An advantageous byproduct of having reduced the system of equations into a single equation in a single unknown is that the error in the actual computation of the derived variables can be
essentially eliminated, or it can be estimated. Eqs. \theequation and \theequation are polynomials; numerical methods for their accurate solution are well established.
Using fixed-point recursion, the maximum computational error can be resolved, and monotonic convergence can be guaranteed as seen in ? (?) and more recently in ? (?).

In contrast,
the classical iterative algorithm
(Andreas, 1989; 2012; Hartogensis, 2003; Solignac, 2009)
may diverge or alternate about a potential solution. At worst, techniques such as the classical algorithm
may stop at a “bottleneck” and converge to a false solution as illustrated in ? (?). In their section on
non-linear coupled equations, it is stated:

“We make an extreme, but wholly defensible, statement: there are no good, general (numerical) methods for solving systems of more than one non-linear equation. Furthermore,
it is not hard to see why (very likely),
there never will be any good, general (numerical) methods…”

In ? (?), similar one-dimensional iterative methods of numerical computation of were used to eliminate computational error,
however the fixed-point algorithm we have presented converges for any (with the correct sign).
We argue that at least some of the spread of data in Figs. 5 and 6 in ? (?) may be due to
computational uncertainty as well as the incorporation of , , and measured at the
scale of an eddy covariance system’s footprint while being forced to assume that they are representative of the beam path scale.
The scatter in these plots may not be entirely due to unreliable and measurements.

Future expansions of the sensitivity analysis presented here may focus on taking into account field sites with heterogeneous terrain and variable topography. For stationary turbulence with beams above the blending height, the line integral formulation for effective beam height given by Eq. B2 in ? (?) and Eqs. 10-12 in ? (?) could be incorporated. Two-dimensional footprint analyses involving surface integrals that take into account variable roughness length and wind direction as in ? (?) and in ? (?) may be incorporated for flat terrain that is heterogeneous enough to force the scintillometer beam to be below the blending height [1986, 1987]. Further theoretical developments may be anticipated that take into account both heterogeneity and variable topography. It is hoped that the general mathematical approach presented here can help to keep track of uncertainty for any scintillometer application, as well as to eliminate the byproducts of iteration.

## Acknowledgements

The authors thank the Geophysical Institute at the University of Alaska Fairbanks for its support, Derek Starkenburg and Peter Bieniek for assistance with editing, two anonymous reviewers, one in particular, for very helpful comments. In addition, the authors thank Flora Grabowska of the Mather library for her determination in securing funding for open access fees. GJ Fochesatto was partially supported by the Alaska Space Grant NASA-EPSCoR program award number NNX10N02A.

## Appendix

## Section A Relations between and

(\theequation) | |||||

(\theequation) | |||||

(\theequation) |

## Section B Individual terms in for unstable conditions ()

(\theequation) | |||||

(\theequation) |

## Section C Individual terms in for unstable conditions ()

(\theequation) | |||||

(\theequation) | |||||

(\theequation) |

## Section D Individual terms in for unstable conditions ()

(\theequation) | |||||

(\theequation) | |||||

(\theequation) |

## Section E Individual terms in for unstable conditions ()

(\theequation) | |||||

(\theequation) | |||||

(\theequation) |

## Section F Total differential expansion as in Andreas (1992) for unstable conditions ()

Here we reproduce the analysis of ? (?). Subscripts indicate the equation that is being differentiated locally. The coupled equations are

(\theequation) | |||

(\theequation) | |||

(\theequation) | |||

(\theequation) |

We expand Eqs. \theequation and \theequation as

(\theequation) | |||

(\theequation) |

Combining these, we obtain

(\theequation) | |||