Analyzing energy–water exchange dynamics in the Thar desert

Abstract

Regions of strong land–atmosphere coupling will be more susceptible to the hydrological impacts in the intensifying hydrological cycle. In this study, micrometeorological experiments were performed to examine the land–atmosphere coupling strength over a heat low region (Thar desert, NW India), known to influence the Indian summer monsoon (ISM). Within the vortex of Thar desert heat low, energy–water exchange and coupling behavior were studied for 4 consecutive years (2011–2014) based on sub-hourly measurements of radiative–convective flux, state parameters and sub-surface thermal profiles using lead-lag analysis between various E–W balance components. Results indicated a strong (0.11–0.35) but variable monsoon season (July–September) land–atmosphere coupling events. Coupling strength declined with time, becomes negative beyond 10-day lag. Evapotranspiration (LE) influences rainfall at the monthly time-scale (20–40 days). Highly correlated monthly rainfall and LE anomalies (r = 0.55, P < 0.001) suggested a large precipitation memory linked to the local land surface state. Sensible heating (SH) during March and April are more strongly (r = 0.6–0.7) correlated to ISM rainfall than heating during May or June (r = 0.16–0.36). Analyses show strong and weak couplings among net radiation (Rn)–vapour pressure deficit (VPD), LE–VPD and Rn–LE switching between energy-limited to water-limited conditions. Consistently, +ve and −ve residual energy [(dE) = (Rn − G) − (SH + LE)] were associated with regional wet and dry spells respectively with a lead of 10–40 days. Dew deposition (18.8–37.9 mm) was found an important component in the annual surface water balance. Strong association of variation of LE and rainfall was found during monsoon at local-scale and with regional-scale LE (MERRA 2D) but with a lag which was more prominent at local-scale than at regional-scale. Higher pre-monsoon LE at local-scale as compared to low and monotonous variation in regional-scale LE led to hypothesize that excess energy and water vapour brought through advection caused by pre-monsoon rainfall might have been recycled through rainfall to compensate for early part of monsoon rainfall at local-scale. However, long-term measurements and isotope analysis would be able to strengthen this hypothesis. This study would fill the key gaps in the global flux studies and improve understanding on local E–W exchange pathways, responses and feedbacks.

Appendix I

Turbulent heat fluxes (SH and LE) computation following Monin–Obukhov similarity theory (MOST).

Under neutral stability condition, the wind profile equation above a homogeneous vegetation canopy is thus, expressed as follows (Thom 1975).

$$u\left( z \right)=\frac{{{u_*}}}{k}~\ln \left( {\frac{z}{{{z_0}}}} \right)$$

(1)

$$\frac{{du}}{{dz}}=\frac{{{u_*}}}{{k~(z - d)}}~.$$

(2)

Here, z is the height of wind measurement above the zero plane displacement (m), z ₀ is roughness length (m), u _* is frictional velocity (m s⁻¹) and k is Von-Karman constant (=0.4). The parameters, d and z₀ depend upon the aerodynamic roughness of the surface over which measurements are made. The zero plane displacement (d) and roughness length (z₀) were computed from plant height (h) using simple empirical relations: d = 0.65 × h and z₀ = 0.1 × h given by Campbell and Norman (1998). The eddy viscosity was estimated from wind speed and roughness parameters using the logarithmic wind profile equation. Solution of logarithmic equation is simple under neutral stability condition. However, absolute neutral stability cannot be assumed always because very often, there is a thermal gradient in either direction in the vertical profile within the lower boundary layer. The atmosphere normally becomes stable during nighttime and unstable during daytime because of substantial heating during daylight hours especially in the sandy desert region. Under unstable or stable conditions, the stability parameter, sensible heat flux and eddy viscosity become auto-correlated. In the present study, stability corrected aerodynamic resistance was determined through iterative method based on MOST (Businger et al. 1971; Xu and Qiu 1997; De Bruin et al. 2000).

According to MOST, the vertical profiles of wind, temperature, and specific humidity for turbulent flows above a vegetation surface can be expressed as:

$${u_2} - {u_1}=\frac{{{u_*}}}{k}\left[ {\ln \left( {\frac{{{z_2}}}{{{z_1}}}} \right) - {\Psi _M}\left( {\frac{{{z_2}}}{L}} \right)+{\Psi _M}\left( {\frac{{{z_1}}}{L}} \right)} \right]$$

(3)

$${\theta _2} - {\theta _1}=\frac{{{\theta _*}}}{k}\left[ {\ln \left( {\frac{{{z_2}}}{{{Z_1}}}} \right) - {\Psi _H}\left( {\frac{{{z_2}}}{L}} \right)+{\Psi _H}\left( {\frac{{{z_1}}}{L}} \right)} \right]$$

(4)

$${q_2} - {q_1}=\frac{{{q_*}}}{k}\left[ {\ln \left( {\frac{{{z_2}}}{{{z_1}}}} \right) - {\Psi _Q}\left( {\frac{{{z_2}}}{L}} \right)+{\Psi _Q}\left( {\frac{{{z_1}}}{L}} \right)} \right].$$

(5)

In the above expressions, u = wind speed (m s⁻¹), θ = potential temperature (K), q = specific humidity (g kg⁻¹), z = height of measurement (m) above the zero plane displacement [(2.5 m−d) and (7.5 m−d)] in the present study. The subscripts, 1 and 2, represent lower and upper heights for the respective parameters in the profile measurement; θ _* = temperature scale, q _* = humidity scale and u _* = frictional velocity (wind scale); Ψ = universal stability function for momentum (subscript M), heat (subscript H) and vapour (subscript Q); L = Monin–Obukhov length which is expressed as follows:

$$L=\frac{{Tu_{*}^{3}}}{{gk{\theta _*}}}.$$

(6)

The eddy sensible and latent heat fluxes are expressed as:

$$SH= - \rho {C_p}{\theta _*}{u_*}$$

(7)

$$LE= - \rho \lambda {q_*}{u_*},$$

(8)

where, ρ = air density, T = mean absolute temperature, Cp = specific heat of air, and g = acceleration due to gravity, λ = latent heat of vaporization of water.

For computation of the stability functions (Ψ) for momentum (M), heat (H) and water vapour (Q) in this study, for unstable conditions, we adopted the version of the Businger–Dyer flux-relationships proposed by Dyer (1974), which read in integrated form as:

$${\Psi _M}=2ln\left( {\frac{{1+x}}{2}} \right)+ln\left( {\frac{{1+{x^2}}}{2}} \right) - 2ta{n^{ - 1}}\left( x \right)+\frac{\pi }{2}$$

(9)

$${\Psi _H}={\Psi _Q}=2ln\left( {\frac{{1+{x^2}}}{2}} \right),$$

(10)

where, $x={\left( {1 - 16\frac{Z}{L}} \right)^{\frac{1}{4}}}.$

For stable conditions, we used the equations proposed by Beljaars and Holtslag (1991):

$$- {\Psi _M}=\frac{{aZ}}{L}+b\left( {\frac{Z}{L} - \frac{c}{d}} \right)exp\left( { - \frac{{dZ}}{L}} \right)+\frac{{bc}}{d}$$

(11)

$$- {\Psi _H}={\left( {1+\frac{2}{3}\frac{{aZ}}{L}} \right)^{\frac{3}{2}}}+b\left( {\frac{Z}{L} - \frac{c}{d}} \right)exp\left( { - \frac{{dZ}}{L}} \right)+\left( {\frac{{bc}}{d} - 1} \right),$$

(12)

in which, a = 1, b = 0.667, c = 5 and d = 0.35 are constants taken as per Beljaars and Holtslag (1991).