Up-scaling short-term process-level understanding to longer timescales using a covariance-based approach

General experience in hydrologic modelling suggests that the parameterisation of a model changes over different time and space scales. As a result, hydrologists often re-parameterise their models whenever different temporal or spatial resolutions are required. Here, we investigate theoretical aspects of this issue in a search for the cause(s) of the need for re-parameterisations. Based on Taylor series expansion, we present a mathematical approach for temporal up-scaling that involves covariance-based corrections. We apply the theory using a unique database of half-hourly pan evaporation measurements (comprising 237 days) and examine how the model parameters change when integrating from half-hour to daily and then monthly integration periods. We show that the model parameters change over different integration periods because of changes in the covariance between the model variables. In our model system, we find that the covariance-based correction is highly variable from day to day but settles down to a reasonably constant value over periods longer than about 15 days. The 15 days timescale is likely to be specific to our model system, nonetheless the underlying principle that there is a characteristic timescale for the covariance-based scaling correction of a particular hydrologic process might be general. If that proved true it would open up the possibility of systematically searching for characteristic integration periods for the key covariance-based scaling terms in other key hydrologic processes. That would in turn enable the development of more generalised hydrologic closure scheme(s).


Introduction
Thirty years ago, hydrologists had a reasonable empirical knowledge of the typical rates of many basic processes like rainfall, evaporation, infiltration, etc.One of the central challenges at that time was to use that knowledge to say something useful at the larger temporal and spatial scales of typical hydrological interest, e.g. annual runoff from an entire catchment or river basin, flood peak estimation, etc. Thirty years on, the term "up-scaling" invokes images of networks of bores, weirs and flux towers that are all tied together in a spatial database using satellite imagery to describe the land cover.These technological tools have proved immensely useful in the scaling of site-specific measurements.
In contrast to ongoing technological developments, much less attention has been given to the theoretical side of the up-scaling problem (Blöschl and Sivapalan, 1995;Viney and Sivapalan, 2004;McDonnell et al., 2007).The theoretical side of the problem is dominated by a key task -the need to correctly calculate the space-time averages (e.g.Milly and Eagleson, 1987;Rastetter et al., 1992;Parlange et al., 1993;Famiglietti and Wood, 1995;Hu andIslam, 1997a, b, 1998;Choudhury, 1999;Koster and Suarez, 1999;Roderick, 1999;Hansen and Jones, 2000;Vogel and Sankarasubramanian, 2003;Albertson and Montaldo, 2003;Rodriguez-Iturbe et al., 2006).That problem is not unique to hydrology -many other related disciplines are also grappling with the same basic problem.For example, climate scientists use the term, "sub-grid variability", to summarise two key concepts involving the calculation of space-time averages.The first key Published by Copernicus Publications on behalf of the European Geosciences Union.
W. H. Lim and M. L. Roderick: Scaling theory concept relates to the representation of various intensive state variables (e.g.pressure, temperature, specific humidity, etc.) in climate models that typically have very large grid cells, and it is known a priori that the intensive state variables vary spatially within a grid cell.The second concept is the calculation of fluxes, and thereby the changes over time, when many of the underlying processes are known to be non-linear.This latter issue is of fundamental importance, because ignoring the non-linearity results in biased predictions (Larson et al., 2001).Of course such bias can be eliminated using "tunable" model parameters but this has the immediate problem that one is never sure if the same tuning would apply in the future.
Different approaches to handling this bias can be envisaged.The traditional approach is to use larger computers to support smaller grid cells and shorter integration periods (e.g.minutes instead of hours, or hours instead of days, etc.) in numerical models.This is a brute-force approach, but at best, it can only reduce the bias, because to account for it completely using a numerical approach would require both infinitesimally small grid cells and time periods which is not possible due to practical computing constraints.For an infinitesimally resolved model to be correct, we also need infinitesimally resolved data for the initial and boundary conditions, and for the parameters in the highly resolved space.Such knowledge would not typically exist.An alternative approach is to account for the "scale bias" by directly estimating it, but to do that requires a clear theoretical understanding of what the "scale bias" actually is.For example, hydrologic modellers are well-aware that model parameters can change with the temporal resolution of rainfallrunoff models (e.g.Littlewood and Croke, 2008;Wang et al., 2009;Kavetski et al., 2011).Woods and Sivapalan (1999) examined space-time variability during storm events and showed that the scaling problem could be formulated using the covariance.Further extension of their work (Viglione et al., 2010a) enabled assessment of the dependence of the catchment flood response to the space-time interactions between rainfall, runoff generation and routing mechanisms followed by its application to characterise the type of flood events based on different precipitation patterns (see details in Viglione et al., 2010b).Although these works focused on single event flood response, a more general idea of explicitly incorporating covariance terms into the theory may assist (i) the description of the scale bias and (ii) development of a more generalised approach in hydrologic modelling.
By comparison, atmospheric scientists have long used a formal closure framework (see Chapter 19 in Shuttleworth (2012) for an accessible and relevant discussion of closure schemes used to model turbulence) to guide their research.In that framework, the covariance terms are explicit and are therefore open to further experimental investigations.However, it is difficult to investigate the feasibility of developing closure schemes for hydrology using a system as complex as a catchment.Hence, there are numerous advantages in using a simpler system to gain theoretical insights into the up-scaling issue.
In this paper we use a unique half-hourly database of high quality pan evaporation measurements (Lim et al., 2012) to examine scale bias in model predictions as a function of the model integration period.To do that we examine how the parameters of a previously formulated and tested evaporation model change when integrating from half-hourly to daily and monthly time periods.We then use those results to develop a covariance-based approach to up-scale evaporation.

Statement of the problem: evaporation hypothesis testing
Dalton's Law (Dalton, 1802) for evaporation from a wet surface can be expressed as where E [m s −1 ] is the evaporation rate of liquid water in traditional hydrologic units of depth per unit time, e s (T s ) (Pa) is the vapour pressure at the evaporating surface and e a (T a ) (Pa) is the air vapour pressure at the same height that air temperature is measured.
The scaling of Eq. ( 1) over a given period (0 to τ ) can be written as where is the aerodynamic function whose numerical value depends on the integration period τ .Here f (τ ) v can be understood as an effective parameter (McNaughton, 1994).It is effective in the sense that it will give the correct estimate of E (τ ) because it is calculated from observations using E (τ ) /(e s (T s ) − e a (T a )).This procedure means that if the time period of the model was changed f (τ ) v might also change.
The correct scaling procedure is to begin with the full equation, where f v [m s −1 Pa −1 ] is the aerodynamic function that is independent of the integration period.By inspection we see that the correct procedure is to calculate the mean of the product (Eq.3) and not the product of the means (Eq.2).To investigate the differences between Eqs.

Taylor series expansion
Based on Taylor series expansion, the theory of integrating the product of two variables (e.g.z = ab) over a given time period (0 to τ ) can be expressed as (see Appendix A for formal derivation) where σ ab is the covariance between a and b, r [a,b] (range: −1.0 to 1.0) is the correlation between a and b, σ a is the population standard deviation of a and σ b is the population standard deviation of b.Note that this result is independent of the distribution of the variables.If r [a,b] → 0 then σ ab → 0 and therefore ab → a b.In other words, when the variables are uncorrelated, the mean of the product is equal to the product of the means.
In the more general circumstance, the variables are correlated.We can simplify Eq. ( 4) by incorporating the coefficient of variation for both a and b (i.e.C a ≡ σ a a and where χ [a,b] is a correction factor arising from the covariance between a and b.Application of Eq. ( 5) for a product involving more than two variables is demonstrated in Sect.3.

Model system
Following our previous study (Lim et al., 2012), we formulate pan evaporation E pan [m s −1 ] as where M w [kg mol −1 ] is the molecular mass of water, R [J mol −1 K −1 ] is the ideal gas constant, ρ w [kg m −3 ] is the density of liquid water, D v [m 2 s −1 ] is the diffusion coefficient for water vapour in air, T a [K] is the air temperature and z is the boundary layer thickness (see Appendix C for details).In our original research, we parameterised Eq. ( 6) using half-hourly data.
Following Eq. ( 3), the scaling of Eq. ( 6) over longer term periods with constants M w and R, and near-constant ρ w removed outside the integration, results in Following Eq. ( 5), we rearrange Eq. ( 7) as a product of the means, , that is multiplied by a collection of additive terms that are associated with covariance-based correction factors.In general, for a product involving five variables, there will be ten covariance-based correction factors.In this instance, there is one less because e s (T s ) and e a (T a ) are related by a sum (and the mean of a sum equals the sum of the means) and not a product.Using these nine covariance-based correction factors, we have (see Appendix B for detailed derivation) This is the general equation for temporal up-scaling to be examined in later sections.

Data and methods
The data were collected from an experimental US Class A pan (with bird guard) located at the Bureau of Meteorology (BoM) field station at Canberra Airport in Australia over the 2007-2010 period (see Sect. 3 in Lim et al., 2012, for full details).We present a brief summary below.
The evaporation rate was calculated from water height measurements made using a magnetostrictive linear displacement transducer (MLDT) (MagneRule Plus, MRU-4001-015, Schawitz Sensors, Hampton, VA, USA) with a spherical float.The float was installed in a stilling well connected to one side of the pan.The water surface temperature was measured using an infrared thermometer (Model: M50-1C-06-L, Mikron Instrument Co. Inc., Oakland, NJ, USA).Wind speed, air temperature, air vapour pressure and atmospheric pressure were measured based on standard installations.Atmospheric pressure was measured using a Vaisala Pressure Transmitter (Model: PTB101B).Wind speed was measured using a cup anemometer at 2 m a.g.l.(above ground level) (u 2 ) and air temperature and air vapour pressure were measured with a Vaisala Humitter (Type 50Y, Vaisala, Helsinki, Finland).All data were collected at 5 min intervals and integrated to half-hourly intervals.
We use high quality (half-hourly) data collected over 237 days for evaluating the magnitude of the covariancebased corrections per Eq. ( 8).The vapour pressure at the evaporating surface was calculated assuming the air immediately adjacent to the surface was saturated at water surface temperature (Monteith and Unsworth, 2008, Eq. 2.27).The boundary layer thickness was estimated using theory developed in the previous paper (see formulation in Appendix C).The half-hourly data are aggregated into days, and then into months.Occasional missing data made it difficult to obtain a complete month of high quality pan evaporation measurements, thus we consider months with a minimum of 16 days as valid months.The resulting 11 months included all seasons (2 months in spring, 3 months in summer, 3 months in autumn and 3 months in winter).

Evaluating the scaling corrections
We examine the use of Eq. ( 8) in up-scaling from half-hourly to daily or monthly time periods in Sects.5.1 and 5.2, respectively.From there, we approximate Eq. ( 8) into a simpler form and show how aerodynamic functions arise from scaling corrections for different integration periods (Sect.5.3).We assess alternative approaches to relate the scaling corrections (Sect.5.4).Finally, we quantify the scaling corrections for practical application (Sect.5.5).

Scaling from half-hourly to daily
The magnitude of the nine covariance-based correction factors (Eq.8) when scaling from half-hourly to a daily basis are shown in Fig. 1.We found most (seven out of nine) of the covariance-based correction factors (i.e.χ) to be approximately zero (Fig. 1a-g).The main reason for the near zero covariance-based correction factors in seven instances is that the coefficient of variation of two variables (D v and 1 T a ) is close to zero (results not shown).In contrast, the covariancebased correction factors for the remaining two terms involving the inverse of the boundary layer thickness ( 1 z ) and the two vapour pressures (e s (T s ), e a (T a )) were relatively large (Fig. 1h, i).Physically, that makes intuitive sense because we observed a strong diurnal variation where the wind speed tended to increase each afternoon thereby increasing 1 z (Lim et al., 2012) and resulting in a correlation between these variables.
The contribution of each covariance-based correction term (Eq.8) to the overall evaporation rate from the pan is shown in Fig. 2. The results show that the product of the means (i.e.E * pan ) accounted for about 80 % of the daily integrated evaporation rate (Fig. 2a).Only those corrections involving the inverse of the boundary layer thickness ( 1 z ) and the two vapour pressures (e s (T s ), e a (T a )) made any significant difference to the integration (Fig. 2i, j).With those correction terms the final results showed excellent agreement with observations (Fig. 2k; R 2 = 0.97, n = 237, RMSE = 0.50 mm d −1 ).

Scaling from half-hourly to monthly
We repeated the above analysis but this time we integrated from half-hourly to a monthly period.The results were virtually identical with the earlier analysis based on integration to a daily time period (Fig. 3).Again, only those covariancebased correction factors involving the inverse of the boundary layer thickness ( 1 z ) with the two vapour pressures (e s (T s ), e a (T a )) were of practical importance (Fig. 3h, i).If all the covariance-based correction terms were ignored, the bias would result in an underestimate of monthly pan evaporation of around 17 % (Fig. 4a).Including these correction terms gave excellent agreement with the integrated observations (Fig. 4k; R 2 = 0.99, n = 11, RMSE = 0.32 mm d −1 ).

Comparing half-hourly, daily and monthly aerodynamic functions
The previous results have demonstrated that in our application, most of the covariance-based correction factors make little practical difference.Retaining the two important correction factors χ that relate the inverse of the boundary layer thickness ( 1 z ) (which increases with the wind speed) with the two vapour pressures (e s (T s ), e a (T a )), we can rewrite Eq. ( 8) as This approximation allows us to express the aerodynamic function f where z .We use that expression to calculate the numerical values of f * v at daily (n = 237) and monthly (n = 11) integration periods and compare those with the original half-hourly results reported by Lim et al. (2012).As anticipated, the long-term (daily, monthly) aerodynamic functions are generally (but not always) larger than the aerodynamic function computed using half-hourly data because of the previously noted correlations between (1 z ) with e s (T s ) and e a (T a ) (Fig. 5).

Attempt to relate scaling corrections using alternative approaches
If short-term data were available it would be straight forward to numerically estimate the two key covariance-based correction factors (i.e.χ ).Without doing that, one could ask whether these correction factors are themselves correlated?If so, one could make further simplifications and possibly avoid the need for detailed calculations.We tested that proposition and found no relationship at either daily or monthly integration periods (Fig. 6).An alternative approach is to seek a physically justified relationship between these correction factors and a key environmental variable that is routinely measured.By defining a new variable, h * ≡ e a (T a ) e s (T s ) , we can rewrite Eq. ( 9) as Fig. 2. Estimated E pan (Eq.8; prime symbol (') denotes partial results) versus observed E pan for 237 days (1 : 1 line shown): and Eq.(10) as We note that e s (T s ) is not normally observed in standard operational practice.Accordingly, we compare the resulting scaling correction with the observed relative humidity h ≡ e a (T a ) e s (T a ) that is routinely measured over both daily and monthly integration periods (Fig. 7).Over daily integration periods, the resulting scaling corrections varied from −10 to 40 % with a mean of ∼ 13 %.Over monthly integration periods, the resulting scaling corrections vary from 0 to 20 % with a mean of ∼ 11 %.The overall (but weak) relation to emerge is that the resulting scaling correction approaches zero when the relative humidity approaches saturation (100 %).At the other extreme, when the relative humidity approaches 30 %, the resulting scaling correction approaches 25 %.However, this relationship is not accurate enough for practical applications.
In summary, the magnitude of the scaling correction relative to the product of the means (i.e.E * pan , f * v ) remain substantial and there does not appear to be a simple way of accurately estimating that as a function of a readily measured environmental variable (e.g.relative humidity).

Diagnostic estimates of the scaling corrections
Our results (Figs. 2, 4) showed that there was little difference in the scaling correction at either daily or monthly integration periods.That implies that the underlying physical reason for the covariance correction is the diurnal co-variation of the wind speed, the vapour pressure of the evaporating surface and the air vapour pressure.In that sense, the accuracy of scaling correction depends on the ability to properly capture the diurnal cycle of these key variables.
Our results also suggest that the mean diurnal cycle of longer integration period might be sufficient to estimate a "representative" total correction factor (i.e. 1 + χ All ; where χ All is the sum of all scaling corrections) that is multiplied with the product of the means (i.e.E * pan , f * v ).To test this idea, we first assess whether the correlation between the key variables (above) remains in the mean diurnal cycle of each season; and identify the number of days needed to calculate the mean diurnal cycle and hence the total correction factor (using summer as an example) in Sect.5.5.1.Using the sampling theory, we then assess whether a mean diurnal cycle obtained from a limited number of days would adequately estimate the total correction factor for each season in Sect.5.5.2.

Mean diurnal cycle
We first classified all daily data from the 2007-2010 period (i.e.237 days) into four seasons and used those days to calculate a mean diurnal cycle for each season (Fig. 8).The results confirm that the correlation between wind speed and the two measures of vapour pressure holds over all seasons.
Using summer data (64 days available) as an example, we assess the accuracy of the total correction factor as a function of the number of days used to estimate the mean diurnal cycle.In Fig. 9, the x axis represents the number of days used (left to right: 1, 2, . . ., 64 days) to obtain a single mean diurnal cycle for each of the three key variables.For instance, when "1 day" is used to calculate the total correction factor, we use all 64 days and therefore have 64 total correction factors (1 + χ All , top panel of Fig. 9).When "2 days" are used, we combine 2 consecutive diurnal cycles into one single mean diurnal cycle and therefore have 32 total correction factors.This aggregating process is continued until we have a single diurnal cycle computed using data from 64 days.Considering the "1 day" case first, the total correction factor varies from 1.02 to 1.37 (mean = 1.16, sample standard deviation = 0.095).Note that the accuracy of the "1 day" case is the highest possible because each of the 64 days has their own accurate total correction factor.At the other extreme, using all 64 days to compute a single diurnal cycle led to a total correction factor of 1.16 that is similar to the mean of the 64 individual daily estimates as would be expected.There is a gradient between these extremes with the root mean square error (RMSE) increasing from 0.6 ("1 day" case) to 0.9 mm d −1 at a threshold of ∼ 16 days (Fig. 9).The practical implications are as follows: (i) while the total correction factor varies widely from day to day, (ii) the total correction factor for a longer period will approach the mean of the individual daily corrections, (iii) and in this example, as few as 15 days are needed to accurately estimate the total correction factor for the entire season.These results are generalised in the following section using sampling theory.

Using sampling theory
Referring to the previous section, the standard error of the mean (SEM) is the sample standard deviation σ s divided by the square root of the number of observations.In the previous summer example we have SEM = 0.095/ √ n.Hence if we have a total correction factor computed using, say 15 days, the SEM = 0.024.Given that the mean total correction factor is ∼ 1.16 this represents an error of around 2 %.Clearly, if  σ s is known then the number of daily observations needed to achieve a given level of accuracy can be estimated.This turns out to be very convenient in our example because the σ s was more-or-less the same in each season (spring: 0.087; summer: 0.095; autumn: 0.086; winter: 0.091).These results are v versus u 2 (wind speed at 2 m a.g.l.) for the experimental pan for different integration periods.readily generalised.Using summer as an example (Fig. 10), the curves of SEM for the total correction factor (i.e.±1 SEM, ±2 SEM, ±3 SEM) narrow sharply when n increases beyond 10 days and are consistent with the pattern shown in Fig. 9 (top panel).By approximating σ s ∼ 0.10 for all the seasons, we estimate the SEM for the total correction factor using the diurnal cycle averaged over 15 days as ±0.05 at a .The results imply that ∼ 15 days of high frequency data would estimate the total correction factor for a season within ±4 % (at the 95 % confidence level).
To test the practical applicability of the above approach we randomly selected 15 days in each season and used that data to calculate a single mean diurnal cycle.With that we estimated the total correction factor and applied that to all days in each season.The results for each season are shown in Fig. 11 and confirm that one can make an accurate seasonal correction using a single mean diurnal cycle based on only 15 days.

Discussion
Starting from known process-level physics, we used a Taylor series expansion to make explicit the covariance terms when integrating from shorter to longer integration periods.We find that short-term (half-hourly) measurements from an evaporation pan can be accurately up-scaled to longer (days, months) integration periods by directly accounting for the covariance using covariance-based corrections per Eq. ( 8).The underlying physical reason for the emergence of persistent covariance at our study site was a general tendency for wind speed to increase in the mid-afternoon when the air temperature and pan evaporation rate are also generally higher (Lim et al., 2012(Lim et al., , 2013)).The details of that relation are no doubt specific to our site but the principle that the temporal (spatial) scaling depends on the temporal (spatial) covariance between model variables is general (e.g.Parlange et al., 1993;Famiglietti and Wood, 1995;Albertson and Montaldo, 2003;Vogel and Sankarasubramanian, 2003).
Our first attempt to make a prognostic estimate of the scaling correction using readily available measurements of relative humidity was physically justified but not accurate enough for practical applications (Fig. 7).The remaining approach is to estimate the covariance-based scaling correction using the mean diurnal cycle over some specific period.The number of days needed to make a correction of a given accuracy depends directly on the day-to-day variability in the numerical value of the scaling correction.In our example, we found that a scaling correction based on the mean diurnal cycle calculated using averages over 15 days could be used to estimate the total correction factor to within 4 % (at the 95 % confidence level) (Fig. 9).Note that errors in the range 10-40 % can be anticipated if the covariance is not explicitly considered (Fig. 8).
The basic idea presented here is similar to the use of covariance terms when examining space-time variability in flood event hydrology (Woods and Sivapalan, 1999;Viglione et al., 2010a, b).Our results extend that by showing that one can use the covariance-based approach in a prognostic way to estimate a scaling correction.Our study was based on a highly simplified system -an evaporation pan.Extending these results to a more realistic setting, e.g. a catchment, remains a grand challenge for the future.In a catchment one also has to deal with numerous thresholds, e.g.runoff occurs when soil saturates, and actual evapotranspiration is less than the potential rate because of water supply restrictions or wind shelter effect associated with natural or artificial barriers (Hipsey and Sivapalan, 2003;Hipsey et al., 2004).Despite that real world complexity, the use of a simple model system has demonstrated that a general covariance-based framework could be developed for up-scaling.

Summary and conclusions
The model system used here, i.e. a half-hourly database of high quality pan evaporation measurements, is simpler than the problems typically faced in hydrology, e.g.catchmentscale water balances, etc.By using a simple system, we showed that the model parameters change for different integration periods because the covariance between model variables changed over those periods.Exact mathematical expressions, based on a Taylor series expansion were developed and used to quantify how the change in covariance propagates into a change in the numerical value of model parameters.
Using our pan evaporation database we found that not all covariance-based correction terms actually matter.In our example, there were nine covariance-based correction terms, yet only two of those made any numerical difference to the  results.The key physical factor in both was the inverse of the boundary layer thickness ( 1z ) (which increases with the wind speed).The two important covariance-based correction terms arose from the covariance between (i) 1 z and the vapour pressure at the evaporating surface e s (T s ) and (ii) 1 z and the air vapour pressure e a (T a ).With those two correction terms we showed that at this site, the numerical value of the aerodynamic function was generally (but not always) larger at both daily and monthly integration periods compared to the original half-hourly data (Fig. 5).That arose because of a strong diurnal cycle in the pan evaporation data where the wind speed usually peaks in the mid-afternoon each day (Lim et al., 2012).
We found that the resulting scaling correction in the pan evaporation application could be readily understood as a function of the relative humidity (Fig. 7).However, that relation was not sufficiently accurate for routine practical applications.We also found that a single covariancebased correction could be implemented for averaging periods longer than about 15 days in our pan evaporation example (Figs. 9,10).This finding points to the possibility that characteristic "covariance timescales" might exist for some other key hydrologic processes.More investigations will be needed to confirm or refute this idea, but if it is found to be general then this might prove a useful way to start developing closure schemes along the lines of those routinely used in the atmospheric sciences (Shuttleworth, 2012).

Appendix A Taylor series expansion for a product of two variables
Let z i = a i b i (for i = 1, . . ., n) be a product of two variables.We now examine the mean z = a b by applying a Taylor series expansion in two variables (see any calculus text, e.g.Adams, 1991) about the point (a, b) to express z in terms of a and b: to the other.Following Appendix A, we let z i = a i b i c i d i (for i = 1, . . ., n) be a product of four variables.With that, the mean z = a b c d in terms of a, b, c and d becomes ≡ g ρ a ρ a L is the "characteristic" speed of air in the vertical direction.We note that g [m s −2 ] (= 9.81 m s −2 ) is the gravitational acceleration and that ρ a and ρ a are calculated using ] is the molecular mass of water.We set ρ a = 0 for conditions when the air density at the evaporating surface is greater than that at the reference height, i.e. minimal buoyancy forces; subsequently u V,C = 0.
term pan evaporation measurements from Eq. (2) in terms of the product of the means, denoted f * v and the same covariance-based correction factors, i.e.

Fig. 1 .
Fig. 1.Relation between the covariance-based correction factors (per Eq. 8) and the relevant product of the means for all nine combination for 237 days as follows: (a) D v , 1 T a ; (b) D v , 1 z ; (c) 1 T a , 1 z ; (d) D v , e s (T s ); (e) D v , e a (T a ); (f) 1 T a , e s (T s ); (g) 1 T a , e a (T a ); (h) 1 z , e s (T s ); (g) 1 T a , e a (T a ); (h) 1 z , e s (T s ); (i) 1 z , e a (T a ).
a (T a )] −e a (T a ) e s (T s ) − e a (T a ) s (T s )] e s (T s ) e s (T s ) − e a (T a ) a (T a )] −e a (T a ) e s (T s ) − e a (T a ) E * pan , (k) sum of (a) to (j) (y = 0.99 x − 0.14, R 2 = 0.97, n = 237, RMSE = 0.50 mm d −1 ).

Fig. 3 .
Fig. 3. Relation between the covariance-based correction factors (per Eq. 8) and the relevant product of the means for all nine combination for 11 months as follows: (a) D v , 1 T a ; (b) D, 1 z , (c) 1 T a , 1 z ; (d) D v , e s (T s ); (e) D v , e a (T a ); (f) 1 T a , e s (T s ); (g) 1 T a , e a (T a ); (h) 1 z , e s (T s ); (g) 1 T a , e a (T a ); (h) 1 z , e s (T s ); (i) 1 z , e a (T a ).

Fig. 6 .
Fig. 6.Relation between the two key covariance-based correction factors for different integration periods.

Fig. 9 .
Fig. 9. Accuracy of the total correction factor (1 + χ All ) versus the number of days used to calculate the mean diurnal cycle.(Data are for summer only.)

Fig. 10 .
Fig. 10.Standard error of the mean of the total correction factor (SEM) versus the number of days (n).The curves (±1 SEM, ±2 SEM, ±3 SEM) were calculated as sample standard deviation σ s (i.e.±1σ s , ±2σ s , ±3σ s ) divided by square root of n. (Analysis based on summer data only.)

+
2 b i − b (c i − c) ∂ 2 z ∂b i ∂c i + 2 b i − b d i − d ∂ 2 z ∂b i ∂d i + 2 (c i − c) d i − d ∂ 2 z ∂c i ∂d i + higher order derivatives, (B2)where all derivatives are evaluated at(a, b, c, d).Here i | (a,b,c,d) = b c, ∂ 2 z ∂b i ∂c i | (a,b,c,d) = a d, ∂ 2 z ∂b i ∂d i | (a,b,c,d) = a c, ∂ 2 z ∂c i ∂d i | (a,b,c,d)= a b and that first and all higher order derivatives are zero.Hence Eq. (B2) becomesz = a b c d = a b c d + σ ab c d + σ ac b d + σ ad b c + σ bc a d + σ bd a c + σ cd a b = a b c d 1 + χ [a,b] + χ [a,c] + χ [a,d] + χ [b,c] + χ [b,d] + χ [c,d] (B3)where χ is a correction factor arising from the covariance between two variables (see simple example in Eq. 5).Similarly we have a b c e = a b c e [1 + χ [a,b] + χ [a,c] + χ [a,e] + χ [b,c] + χ [b,e] + χ [c,e] ].From there, we express Eq.(B1) asE pan = M w R ρ w a b c d − e 1 + χ [a,b] + χ [a,c] + χ [b,c] D v ,e s (T s )] e s (T s ) e s (T s ) − e a (T a ) + χ [D v ,e a (T a )] −e a (T a ) e s (T s ) − e a (T a ) + χ [ 1 Ta ,e s (T s )] e s (T s ) e s (T s ) − e a (T a ) + χ [ 1 Ta ,e a (T a )] −e a (T a ) e s (T s ) − e a (T a ) + χ [ 1 z ,e s (T s )] e s (T s ) e s (T s ) − e a (T a ) + χ [ 1 z ,e a (T a )] −e a (T a ) e s (T s ) − e a (T a ), s (T s ) − e a (T a )).This is Eq.(8) in the main text.
− e a (T a )) M a + e a (T a ) M w T a − (P a − e s (T s )) M a + e s (T s ) M w − e a (T a )) M a + e a (T a ) M w T a + (P a − e s (T s )) M a + e s (T s ) M w T s (C3) respectively, where R [J mol −1 K −1 ] is the ideal gas constant, T s [K] is the water surface temperature, T a [K] is the air temperature, P a [Pa] is the atmospheric pressure, e a (T a ) [Pa] is the air vapour pressure, e s (T s ) [Pa] is the vapour pressure at the evaporating surface, M a [kg mol −1 ] is the molecular mass of air and M w [kg mol −1