Calculation

Using the SEBCS for QGIS plug-in, it is possible to calculate several characteristics related to energy exchange at the surface and energy transformation into individual heat fluxes. It also calculates meteorological features and vegetation indices that have some connection with the land microclimate. The following features can be calculated:

  • meteorological features

  • vegetation cover characteristics

  • radiation balance features

  • surface aerodynamic features and boundary layer (atmospheric) stability

  • heat balance features

Radiation balance features

The calculation of the radiation balance features includes the balance of short-wave and long-wave radiation on the active surface and the resulting total net radiation. The calculation procedure includes the calculation of the albedo, surface reflectance and atmospheric emissivity. The SEBCS for QGIS includes surface geometry issues in the radiation balance features calculation. The relationship of the individual radiation balance components is summarized by the equation (1):

(3)\[Rn=Rs_{\downarrow}-Rs_{\uparrow}+Rl_{\downarrow}-Rl_{\uparrow}\]

The individual radiation balance components calculation procedure is described in the following text.

Incoming short-wave radiation flux

The incoming (incident) short-wave radiation flux (irradiance) is the main source of energy for evaporation and energy transformations on the Earth’s surface. The flux is usually measured at meteorological stations as the energy of radiation incident on a horizontal surface (global radiation), but the amount of energy reaching the surface is significantly related to its shape and the geometry of the incident radiation. The SEBCS for QGIS module calculates the flux of incident shortwave radiation as a function of surface shape and radiation geometry. The computational approach summarized by Kumar (1997) was used as follows. The value of the incident short-wave radiation flux on the surface was calculated using equation:

(4)\[Rs_{\downarrow} = I_s \cdot \cos i\]

where

(5)\[\begin{split}\cos i = \sin \delta_s (\sin Lat \cos \beta_s - \cos Lat \sin \beta_s \cos a_{w}) \\ + \cos \delta_s \cos H_s (\cos Lat \cos \beta_s + \sin Lat \sin \beta_s \cos a_{w}) \\ + \cos \delta_s \sin \beta_s \sin a_{w} \sin H_s\end{split}\]

The value of \(\delta_s\) can be calculated using equation:

(6)\[\delta_s = 23.45 \sin \left(\frac{360(284+N)}{365}\right)\]

The hour angle corresponds to the condition where \(H_s = 0\) at noon (when the Sun is on the meridian). Each hour plus represents a change in angle of +15°; each hour minus represents a change in angle of -15°. By convention morning values are positive and afternoon values are negative.

(7)\[H_s = (12 - S_t)\cdot 15\]

Depending on the longitude, the solar time (\(S_t\)) was determined according to the relation:

(8)\[S_t = GMT + \frac{24 \cdot Long}{360}\]

Since the module uses the global radiation measured per horizontal surface as input, the value of \(I_s\) is calculated based on the relationship:

(9)\[I_s = \frac{R_{s\downarrow konst}}{\sin \alpha_s}\]

where

(10)\[\sin \alpha_s = \sin \delta_s \sin Lat + \cos \delta_s \cos Lat \cos H_s\]

Surface reflectance (albedo)

The calculation of the broadband reflectance or albedo is based on empirical approach. An empirical relation was used for the calculation, which calculates the broadband albedo based on the spectral reflectance of the surface for individual spectral bands according to the relation (Tasumi et al. 2008):

(11)\[\alpha = \displaystyle\sum_{b=1}^{7} (\rho_{s\_b} \cdot w_b)\]

The \(w_b\) constants for each Landsat satellite are given in the table:

Table 2 Values of \(w_b\) for particular spectral bands according to Liang et al. (2001 and 2003) and Tasumi et al. (2008) for Landsat 4, 5 and 7 and according to Olmeo et al (2017) for Landsat 8 and 9.

Spectral band

Blue

Green

Red

NIR

SWIR1

SWIR2

Landsat 8, 9

0.246

0.146

0.191

0.304

0.105

0.008

Landsat 4, 5, 7

0.254

0.149

0.147

0.311

0.103

0.036

Alternatively, for other data sources, an empirical approach based on the use of vegetation indices NDVI and MSAVI can be used (Duffková et al. 2012):

(12)\[\begin{split}\alpha = 0.08611 + 0.89472 \cdot MSAVI + 5.55866 \cdot MSAVI^2 -0.1183 \cdot NDVI\\ - 1.9818 \cdot MSAVI^3 - 4.5034 \cdot MSAVI \cdot NDVI - 11.463 \cdot MSAVI^2 \cdot NDVI\\ + 7.46145 \cdot MSAVI \cdot NDVI^2 + 5.2994 \cdot MSAVI^2 \cdot NDVI^2\\ + 4.76657 \cdot MSAVI^3 \cdot NDVI - 2.3127 \cdot MSAVI^3 \cdot NDVI^2\\ - 3.4274 \cdot MSAVI \cdot NDVI^3\end{split}\]

The reflected shortwave radiation flux is calculated as follows:

(13)\[Rs_\uparrow = \alpha \cdot Rs_\downarrow\]

Incoming long-wave radiation flux

The incoming longwave radiation emitted by the atmosphere is calculated from the air temperature measured at \(z\). The calculation is based on the Stefan-Boltzmann law:

(14)\[Rl_{\downarrow} = \varepsilon_{ac} \sigma {T_{a\_K}}^4\]

The emissivity of the atmosphere is calculated based on the air temperature at the \(z\) level and the amount of water vapour in the air. According to Brutsaert (1982) it is calculated:

(15)\[\varepsilon_{ac} = 1.24 \left( \frac{e_a \cdot 10}{T_a + 273.16} \right)^{\frac{1}{7}}\]

This approach is largely approximate and can be used for clear-sky weather conditions using the calculated \(\varepsilon_{ac}\) value. A more appropriate approach is to measure the value of \(Rl_{\uparrow}\) directly.

Emitted long-wave radiation flux

The longwave radiation emitted by a surface is determined by the temperature of the surface. The calculation is based on the Stefan-Boltzmann law:

(16)\[Rl_{\uparrow} = \varepsilon \sigma {T_{s\_K}}^4\]

The surface emissivity is calculated based on an empirical approach of emissivity determination using the NDVI Treshold Method - \(NDVI^{THM}\) (Sobrino et al. 2004). The method uses the NDVI index (Normalized Difference Vegetation Index, Tucker 1979). For the emissivity determination, the range of index values is divided into three categories:

  • \(NDVI < 0.2\) - In this case the surface is considered as bare ground and the emissivity values are derived from the reflectance values in the red region of the spectrum

  • \(NDVI > 0.5\) - In this case the surface is fully covered by vegetation and a typical emissivity value of \(\varepsilon = 0.99\) is determined.

  • \(0.2 ≤ NDVI ≤ 0.5\) - In this case the surface can be considered as a mixture of bare soil and vegetation cover.

The relationship between emissivity and surface cover is shown by the following relationship:

(17)\[\varepsilon = \varepsilon_v P_v + \varepsilon_s (1 - P_v) +d\varepsilon\]

where fraction of vegetation cover is calculated as follows:

(18)\[P_v = \left(\frac{NDVI - NDVI_{min}}{NDVI_{max} - NDVI_{min}}\right)^2\]

where \(NDVI_{min} = 0.2\) and \(NDVI_{max} = 0.5\). As a result, the emissivity can be calculated based on the empirical relation as follows (Sobrino et al. 2004):

(19)\[\varepsilon = 0.004 P_v + 0.986\]

Vegetation cover characteristics

The calculation of vegetation cover characteristics includes the estimation of vegetation spectral indices and leaf area index.

Normalized Difference Vegetation Index (NDVI) is one of the most widely used spectral vegetation indices. NDVI provides information on vegetation cover, its quality and possibly also its quantity (biomass). NDVI can be calculated as follows:

(20)\[NDVI=\frac{R_{NIR}-R_{RED}}{R_{NIR}+R_{RED}}\]

Modified Soil Adjusted Vegetation Index (MSAVI) has similar uses to the NDVI. Unlike the NDVI, there is no oversaturation of values with higher vegetation cover. The MSAVI can also be used to estimate vegetation height. MSAVI can be calculated as follows:

(21)\[MSAVI=0.5\cdot(2R_{NIR}+1-\sqrt{(2R_{NIR}+1)^{2}-8\cdot(R_{NIR}-R_{RED})}\]

Normalized Difference Moisture Index (NDMI) can be used for surface moisture estimation. NDMI can be calculated as follows:

(22)\[NDVI=\frac{R_{NIR}-R_{SWIR1}}{R_{NIR}+R_{SWIR1}}\]

Soil Adjusted Vegetation Index (SAVI) is next widely used vegetation index. It is similar to MSAVI. SAVI can be calculated as follows:

(23)\[SAVI = \frac{NIR - R}{NIR + R + L} (1 + L)\]

where \(L\) is soil brightness correction factor defined as 0.5 to accommodate most land cover types. SEBCS for QGIS calculates leaf area index (\(LAI\)) according to Jafaar & Ahmad (2016):

(24)\[\begin{split}LAI_1= \begin{cases} 11 \cdot SAVI^3 & SAVI > 0;\ SAVI \leq 0.817\\ 6 & SAVI > 0.817 \end{cases}\end{split}\]
(25)\[\begin{split}LAI_2 = \begin{cases} -\frac{\ln \frac{0.61-SAVI}{0.51}}{0.91} & SAVI >0;\ SAVI\leq 0.61\\ 6 & SAVI > 0.61 \end{cases}\end{split}\]
(26)\[LAI = \frac{LAI_1 + LAI_2}{2}\]

The SEBCS_lib library contains some additional methods of \(LAI\) calculation. See documentation.

Meteorological features

SEBCS for QGIS uses several variables to calculate aerodynamic quantities and heat fluxes, either constants (see Abbreviation list) or input variables such as water vapour pressure etc.

The calculation of heat fluxes by SEBCS for QGIS uses the vertical gradient of the meteorological characteristics such as temperature gradient, humidity gradient or change in air velocity. Because the models usually work with a large spatial extent, the calculation typically uses the gradient within the atmospheric boundary layer, i.e. between the Earth’s surface and the mixing layer, which is set by default to \(z = 200 m\).

The calculation of air temperature for height \(z\) above the surface assumes an adiabatic change in air temperature \(\Gamma\):

(27)\[T_{z} = T_{a} + \Gamma(Z - Z_{st})\]

The air pressure \(P\) for height \(z\) is calculated using the relation:

(28)\[P = 101.3 \left( \frac{293-\Gamma \cdot (DMT + z)}{293} \right)^{5.26}\]

The characteristics of the humidity and its gradient are based on the assumption of a constant value of relative humidity in the vertical gradient. The water vapour pressure is then expressed by the relation:

(29)\[e_z = \frac{E_z \cdot Rh}{100}\]

where \(E_z\) is the pressure of saturated water vapour at height \(z\) calculated using the Magnus-Tetens equation (Buck 1981):

(30)\[E_z = 0.61121 \cdot \exp \left( \frac{17.502 \cdot T_z}{240.97 + T_z} \right)\]

Similarly, the saturated water vapour pressure is calculated for the surface:

(31)\[E_s = 0.61121 \cdot \exp{\left(\frac{17.502 \cdot T_s}{240.97 + T_s}\right)}\]

Water vapour pressure deficit is calculated as follows:

(32)\[VPD = E_z - e_z\]

The air density for the height \(z\) is calculated by the relation:

(33)\[\rho = \frac{353.4}{T_z + 273}\]

Latent heat is calculated by the relation:

(34)\[\lambda = 2501 - 2.3723 \cdot T_z\]

The psychrometric constant \(\gamma\) is calculated by the relation:

(35)\[\gamma = \frac{c_p \cdot P}{\lambda \cdot 0.622}\]

The gradient of saturated water vapour \(\Delta\) is calculated by the relation:

(36)\[\Delta = 45.03 + 3.014 T + 0.05345 T^2 + 0.00224 T^3\]

where

(37)\[T = \frac{T_z + T_s}{2}\]