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:
radiation balance features
vegetation cover characteristics
meteorological 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 (3):
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:
where
The value of \(\delta_s\) can be calculated using equation:
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.
Depending on the longitude, the solar time (\(S_t\)) was determined according to the relation:
Since the module uses the global radiation measured per horizontal surface as input, the value of \(I_s\) is calculated based on the relationship:
where
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):
The \(w_b\) constants for each Landsat satellite are given in the table:
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):
The reflected shortwave radiation flux is calculated as follows:
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:
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:
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:
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:
where fraction of vegetation cover is calculated as follows:
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):
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:
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:
Normalized Difference Moisture Index (NDMI) can be used for surface moisture estimation. NDMI can be calculated as follows:
Soil Adjusted Vegetation Index (SAVI) is next widely used vegetation index. It is similar to MSAVI. SAVI can be calculated as follows:
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):
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\):
The air pressure \(P\) for height \(z\) is calculated using the relation:
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:
where \(E_z\) is the pressure of saturated water vapour at height \(z\) calculated using the Magnus-Tetens equation (Buck 1981):
Similarly, the saturated water vapour pressure is calculated for the surface:
Water vapour pressure deficit is calculated as follows:
The air density for the height \(z\) is calculated by the relation:
Latent heat is calculated by the relation:
The psychrometric constant \(\gamma\) is calculated by the relation:
The gradient of saturated water vapour \(\Delta\) is calculated by the relation:
where