Power-function expansion of the nondimensional complementary relationship of evaporation: the emergence of dual attractors

The polynomial form of the nondimensional complementary relationship (CR) follows from an isenthalpic process of evaporation under a constant surface available energy and unchanging wind. The exact polynomial expression results from rationally derived first and second-order boundary conditions (BC). By keeping the BCs, the polynomial can be extended into a two-parameter (a and b) power function for added flexibility. When a = b = 2 it reverts to the polynomial version. With the help of Australian FLUXNET data it is demonstrated that the power-function formulation excels among CR-based two-parameter models considered, even when a = 2 is prescribed to reduce the number of parameters to calibrate to two. The same powerfunction approach (a = 2) is then employed with a combination of different gridded monthly potential evaporation terms across Australia, while calibrating b against the multiyear simplified water-balance evaporation rate on a cell-by-cell basis. The resulting bi-modal histogram of the b values peaks first near b = 2 and then at b 1 (secondary modus), confirming earlier findings that occasionally a linear version (i.e., b = 1) of the CR yields the best estimates. It is further demonstrated that the linear form emerges when regional-scale transport of moist air is negligible toward the study area during its drying, while the more typical nonlinear CR version prevails otherwise. A thermodynamic-based explanation is yet to be found as to why the flexible power function curves (i.e., b [?] 2) converge to the polynomial one (b = 2) in such cases.


37
The complementary relationship (CR) of evaporation is a powerful tool [see the latest global 38 studies by Ma et al. (2021), Brutsaert et al. (2020)] for predicting actual land evaporation (E) 39 rates with the help of basic meteorological variables (i.e., air temperature, humidity, net surface 40 radiation and wind speed) all obtained at a single elevation above the ground. Since its original 41 formulation by Bouchet (1963), it has evolved into various versions [see Han & Tian (2020) for 42 a brief overview] based on different heuristic arguments. 43 After almost six decades of the groundbreaking study by Bouchet (1963), Szilagyi (2021)    been originally inspired by the study of Brutsaert (2015). 60 Here a brief summary of this thermodynamical approach is provided. First the nondimensional 61 linear as well as the polynomial CR equations are derived. Then the latter is expanded by a 62 power function formulation to make it more flexible. The resulting power function with two 63 additional parameters (a and b) is to be applied with daily measurements of air temperature (T), 64 pressure (p), vapor pressure deficit (VPD), net radiation (Rn), ground heat conduction (G) and Note that this work is not meant as a calibration/verification analysis of a preferred two-77 parameter approach over other existing similar (or single parameter) approaches. That is why the 78 steps required for such a study (i.e., validation with data separate from calibration, sensitivity 79 analysis of the parameters, etc.) are deliberately not repeated here, specifically because it would 80 blur the focus of the present work which is the investigation/demonstration, by the help of a 81 recently discovered power-function expansion, of how and when the horizontal advection of 82 humidity over a drying surface produces/affects the thermodynamically-derived linear and 83 nonlinear forms (Szilagyi, 2021) of the CR of evaporation and the typical environmental 84 conditions under which, one or the other, emerges.   (Stull, 2000). The wet-bulb 105 temperature, Twb, is the lowest temperature the air at the measurement height can attain by 106 evaporation, but this temperature is rarely reached during natural processes due to large-scale 107 vertical mixing of free tropospheric air into the boundary layer (Brutsaert, 1982). Instead, a wet 108 environment air temperature, TPT ≥ Twb, generally occurs. TPT however is not known during 109 drying conditions of the environment (i.e., when Ta > TPT), but it can be estimated by the wet-110 surface temperature, Tws, because in humid conditions air temperature changes mildly with 111 elevation above the ground (Laikhtman,1964;Stull, 2000;Szilagyi, 2014).

112
Tws can be estimated (Szilagyi & Jozsa, 2008)  given by where Δ denotes the slope of the saturation vapor pressure curve (hPa C°-1 ) at the measured air 123 temperature Ta, and the empirical wind function, fu (mm d -1 hPa -1 ), is traditionally specified as fu 124 = 0.26(1 + 0.54u2) (Brutsaert, 1982). Here u2 (m s -1 ) is the horizontal wind speed at 2-m above 125 the ground and can be estimated by a power function (Brutsaert, 1982) from measurements (uh) 126 at h meters above the surface as u2 = uh (2 / h) 1/7 . The e * (Ta) -ea expression in the aerodynamic 127 term of Eq. 3 is often referred to as the vapor pressure deficit (VPD).

128
With the two isenthalps anchored to the saturation vapor pressure curve, one may notice that The right-hand-side of Eq. 4 can be further expanded into The combination of Eqs. 5 and 6 yields (Szilagyi, 2021) 144 which via the corresponding evaporation terms in Fig. 1 can be written as  complementarity in the CR means that E and Ep change in opposite ways (Bouchet, 1963), best 167 seen in Eq. 10 between E and EpX. When Ep increases (i.e., the environment dries), wi decreases As an area dries, large-scale horizontal advection of more humid air from the surrounding larger 170 region may occur. This is especially true for areas lying downwind from a sea, or other large 171 body of water, or areas surrounded by mountains having much wetter conditions. This influx of 172 external humid air suppresses or completely eliminates the weak vertical humidity gradient that 173 would otherwise exist. This means that the resulting suppressed and therefore vanishing es -ea 174 term in Eq. 7 would not respond anymore to changes in es and therefore to the ensuing ea that a 175 change in es would normally generate, leaving the left-hand-side of Eq. 7 unresponsive to any 176 changes in (the transported) ea itself, thus causing dy / dX → 0 when X → 0 (Szilagyi, 2021).

177
Note that e * ws and ePT are conservative (invariant) quantities only under isenthalpic processes and

190
Eq. 11 has already been applied on a monthly basis in a calibration-free mode, employing a   211 The polynomial in Eq. 11 can be expanded by a power-function approach using the same BCs.

212
The resulting function has two parameters additional to Eq. 11, a and b. most practical applications the parameter ranges can be narrowed to 1 < a ≤ 2 and 1 < b < 10.

220
In order to reduce the number of parameters to just two (the PT-α, and b) in Eq. 12 for a 221 meaningful comparison with other existing two-parameter CR-based methods, a = 2 is 222 prescribed in this study for evaporation estimation. It makes also possible that the power-1.001 indeed has a vanishing slope at X = 0, as BC iv) requires, but it is indistinguishable from 226 the y = X line of Eq. 10 by the naked eye. For this reason, during calibration of b in the ensuing 227 analysis with a = 2 imposed, a value of b = 1 will be allowed for practicality, even though it 228 violates BC iv).  Table 4

270
A 5-day aggregation instead of a weekly one is chosen, because Morton (1983) argues that it is 271 the shortest time-interval over which any effect of passing weather systems, temporarily 272 upsetting the dynamic equilibrium between the surface and the overlying air, can be expected to 273 be substantially subdued.

274
Performance of the calibrated models is summarized in Table 3. The four (plus Eq. 11) models 275 behave similarly in terms of the root-mean-square error (RMSE), but Eq. 12 produces the best 276 results in seven out of the nine cases considered, followed by Eq. 11 (four occasions, provided 277 Eq. 12 is excluded) and KB06 (twice). In fact, Eq. 12 is always the best performing model with forests serious aridity never occurs as the majority of the points are situated at X > 0.5 (Fig. 6b), 284 with corresponding evaporation rates, EEC > 1 mm d -1 (Fig. 6a), therefore the effect of any in several points at X < 0.2 (Fig. 6d), and EEC < 0.5 mm d -1 (Fig. 6c). Any horizontal moisture 289 transport to the grass sites somewhat drier than their environs will leave the eddy-covariance 290 measurements largely unaffected in the beginning of drying when vertical gradients of the vapor 291 pressure over the grass are still significant, but nonetheless, will depress the value of Ep (which is 292 KB06 HT12 * GX21 Parameters: α, d pressure curve at high temperatures), and thus boost the wetness index, wi, within X, which then 294 moves the measurement points to the right horizontally in Fig. 6d, away from the 1:1 line for 0.2 295 < X < 0.45. The measurement points however will follow the diminishing slope of Eq. 12 at 296 extreme low X (< 0.2) values (as seen in Fig. 6d) and get closer to the 1:1 line again when large-297 scale horizontal moisture advection itself weakens as arid conditions probably spread spatially.
298 Figure 6. Regression plots of the modeled (Emod) 30-day evaporation rates against eddy-covariance measurements 300 (EEC) at two forested (a) and two grass (c) sites of FLUXNET in Australia (see Fig. 5 for locations) together with the 301 least-squares-fitted straight lines. Graphical representation of the calibrated (see Table 3 Ew expectation. While such is the case mostly for the grass sites (Fig. 6d), it is not so for the 307 forested ones (Fig. 6b), due to their incorrect scaling that produces xKB and x (Table 2)   Testing the power-function approach with gridded simplified water-balance data 319 Eq. 12 is further tested across Australia for the spatial distribution of its b value, employing 0.25- is negative. This is to be expected, as the measurement points (Ewb or EEC) are fixed in the 361 nondimensional graph once the value of α is set within X. An overestimation (i.e., when the 362 curve is above a given marker point in e.g., Fig. 6d) in Eq. 12 can only be corrected by moving 363 the curve to the right which is achieved by increasing the value of b (Fig. 4), and vice versa for

378
An interesting property of the histogram is that it is bimodal, with a secondary peak near b = 1.

379
As discussed above, a unity value of b and a linear relationship between y and X (except in the 380 vicinity of X = 0 where the slope must vanish due to the BCs in Eqs. 11 and 12) result in theory 381 when in Fig. 1  Tasmania, all where annual precipitation rates are the highest and aridity the lowest (Fig. 10), 389 forcing the calibrated values of b to remain unity (Fig. 8a).

399
As seen in Fig. 9, about 95% of the histogram values are less than three. In fact, b > 3 occurs 400 predominantly along the dry southern and western seashore (Fig. 8a)

433
Table 4 also indicates that the value of b and its spatial behavior with gridded data are not Naturally, the preservation of a constant relative speed between the two isenthalps' state 489 coordinates cannot be expected to exist in a strict sense, at all times, due to unavoidable changes 490 in Qn, air pressure, and/or wind conditions during the averaging period (typically from day to 491 month), but rather in a statistical sense, as a mean behavior over the averaging period.

492
Eq. 12 may be preferrable over the existing single-parameter (and calibration-free when applied 493 with gridded data of a large domain) polynomial approach of Eq. 11, due to its built-in flexibility 494 when calibration is made possible by available measured (e.g., eddy-covariance derived) or 495 water-balance based E estimates and/or the possibility exists that a linear CR approach (i.e., 496 when a = 2 and b = 1 in Eq. 12) yields (even temporarily, during wet conditions that appear in Data availability All data used in this study are publicly available from the following sites.