Spatio-temporal characterization of landscape fire in relation to anthropogenic activity and climatic variability over the Western Himalaya, India

ABSTRACT

The frequent occurrence of a forest fire can be detrimental to the functioning of ecologically fragile regions such as the Western Himalaya. This study assesses the spatio-temporal variability and trends of fire occurrence in forests, shrubland, and grassland in relation to anthropogenic activities and climate oscillations over Uttarakhand and Himachal Pradesh; the most fire-affected states in the Western Himalaya. The Pre-monsoon (March to June) temporal distribution of fire occurrence is bimodal with the first and second fire peaks occurring in the last week of April and late May/early June, respectively. The length of fire season was found to have reduced by ~ 10 days in 2010s compared to 2000s with the increasing burned area concentration over a shorter period (~110 days). The Mann-Kendall Tau (τ) and Sen’s-slope suggested an increasing trend (τ = 0.3, Sen’s slope = 46.78 km² year⁻¹, p = 0.02) of burn area over Uttarakhand where most fire activity occurs. The Getis-Ords (Gi*) statistic was derived by using Fire Radiative Power (FRP) data to infer the spatial distribution of hot and cold clusters which were found to be greatest (Gi* = 1–10) in evergreen needleleaf forest which also accounted for a marginally higher mean per-pixel FRP than other land cover types. Analyses of the role of climatic variables on fire activity revealed moderate positive and negative correspondence between the burned area and average maximum-temperature (R = 0.57, p = 0.01) and total precipitation (R = – 0.53, p = 0.02), respectively. Assessment of the variation of fire activity with global climatic oscillation indices found the Niño-3.4 index to have the strongest, albeit modest, relationship (R = 0.51) with burned area anomalies. This study highlights the increase in fire intensity in the fragile forests of the Western Himalayan region and the requirement for further research into the role of climate and anthropogenic activities on fire occurrence in order to best preserve regional biodiversity and develop an effective forest management strategy.

KEYWORDS:

• Fire radiative power
• burned area
• regional climate
• global climatic Oscillations
• western Himalaya

1. Introduction

Landscape fire occurrence is influenced by both anthropogenic activity and climate variability and is primarily concentrated around tropical and subtropical biomes (Moritz et al. 2012). As a global phenomenon, wildfire occurrence results from the interaction of climate, weather, physiography, fuel characteristics, and anthropogenic activity (Flannigan et al. 2009). Fire is an integral part of many ecosystems (Fernandes, Rego, and Rigolot 2011; Stott, Goldammer, and Werner 1990) but currently, climate change and increasing human activities have bolstered fire frequency which can have an adverse impact on the ecosystem and environmental functioning (Flannigan et al. 2009; Whitlock et al. 2010). Along with alteration of the land surface characteristics, fires also release large quantities of aerosols and trace gases into the atmosphere which alters atmospheric composition and reduce air quality (Ramanathan et al. 2005; Crutzen et al. 1979; Keywood et al. 2013; Crippa et al. 2016). The majority of global fire activity occurs in the tropics, with savannas being the most affected land cover type (Dwyer, Gregoire, and Malingreau 1998; Maurin et al. 2014) although large-scale fires frequently occur throughout the tropical forests of South America and Southeast Asia (Roberts 2001).

The spatio-temporal variability of fire severity can have long-term consequences on flora and fauna, including forest structure, composition, diversity, and productivity (Parashar and Biswas, 2018), and the physical and chemical composition of the soil (Knicker 2007). The surface characteristics (e.g. elevation, topography, slope), environmental (e.g. forest type, fuel type, loading, and moisture, etc.), and climatic or weather (precipitation, temperature, wind, humidity, solar radiation, etc.) conditions, directly and indirectly influence fire behavior (Rothermel 1983) and fire regime (Negi and Kumar 2016). The decrease in precipitation during the dry season reduces the surface and fuel moisture which increases fire risk which, during periods of prolonged drought, can result in larger and more destructive landscape fires (Mandel et al. 2014; Randerson et al. 2006). Furthermore, regional climatic fluctuations due to climatic oscillations, such as El Niño and La Niña, impact local, regional to global fire regimes (Mason et al. 2017; Alencar, Nepstad, and Diaz 2006) and fire severity (Simard and Haines 1985; Skinner et al. 2006; Fuller and Murphy 2006). Along with environmental connection, fire severity and fire regime are tightly coupled with anthropogenic activities such as increasing population density (Marlon et al. 2008), the rapid growth of wildland-urban interface (Radeloff et al., 2018), land clearance (Bowman et al., 2020) and human modifications (Haas, Calkin, and Thompson 2013). Whilst nearly 90% of the fire ignitions in the United States are due to human activity, climate change has increased the length of the fire season and the size of the burned area (Rogers et al. 2020; National Research Council 2010). The situation is predicted to get worse under a warmer climate and increased population pressure, with increased fire occurrence, greater burned area, and longer fire seasons across the globe (Flannigan et al. 2006).

Fire activity has been increasing in many countries in South and Southeast Asia over the last two decades, particularly in India, Pakistan, Indonesia, and Myanmar (Vadrevu et al. 2019). A total of 520,861 active-fire detections were reported between 2003 and 2017 across India (Sannigrahi et al. 2020) which are now considered a major cause of forest degradation in central India (Chandra and Bhardwaj 2015), western Himalaya (Dobriyal and Bijalwan 2017), and the Northeast Himalaya (Puri et al. 2011). The Western Himalaya is a fragile biodiversity hotspot where the landscape is composed of evergreen broad and needle leaf forest, deciduous broadleaf forest, mixed forest, shrubland, and grassland (Reddy et al. 2015; Roy et al. 2015). The prevalence of coniferous forest cover (~16% of the forested area) (Dobriyal and Bijalwan 2017) plays a vital role in maintaining primary productivity in the region through the deposition of needles (Singh and Parida 2019; Singh et al. 2019) which also serves as the main source of flammable fuel (Dobriyal and Bijalwan 2017). During the pre-monsoon season (March to June), higher temperatures and reduced precipitation (Singh and Mal 2014) increase the severity of fires over recent decades which are responsible for widespread forest degradation and biodiversity loss (Dobriyal and Bijalwan 2017; Prabhakar, Somanathan, and Mehta 2006; Roy 2005).

Remote sensing data and geospatial techniques are key tools in the identification, monitoring, and assessment of landscape fires (Chuvieco 2009; Reddy et al. 2019; Bar et al., 2020). Detection of the burned area using remotely sensed data is typically carried out using surface reflectance measurements and, in some cases thermal observations. Burned area algorithms often discriminate between burned surfaces, where the over-story vegetation has been removed to varying degrees and replaced with char and ash, and unburned surfaces using spectral reflectance or vegetation indices (VIs) (Roy et al. 1999; Stroppiana et al. 2002; Giglio et al. 2010). A number of operation datasets such as MOD14GN, MOD14A1, MOD14A2, MOD14CMQ, MYD14CMQ, MCD14ML, MCD64CMQ, and MCD64A1 have been developed using Moderate Resolution Imaging Spectroradiometer (MODIS) data which provide burned area measurements at spatial resolutions of 250–1000 m and daily temporal resolution (Giglio, Schroeder, and Justice 2016; Chuvieco et al. 2016; Urbanski et al. 2018). On the other hand, active-fire detection utilizes thermal observations in the middle and longwave infrared to detect pixels that contain one or more fires that are actively burning at the time of satellite overpass. Unlike burned area products, active-fire algorithms detect thermal anomalies during each overpass, providing information on the diurnal variation in fire activity, and can do so in near real-time. In addition to the direct detection of fires and fire-affected areas, remote sensing data are also a key resource in the assessment of the fire risk (Matin et al. 2017; Bar et al., 2020), fuel-type mapping (García et al. 2017), burn severity assessment (Meng et al. 2017) and in the analysis of vegetation recovery (SParks et al. 2018; Szpakowski and Jensen 2019).

Several studies have investigated the dynamics of fire behavior and fire risk in the Western Himalaya region (Joseph, Anitha, and Murthy 2009; Babu, Roy, and Prasad 2016; Gupta et al. 2018; Bar et al., 2020) but these studies have been limited to a few years (~ 3–10) of analysis and do not explicitly assess the role between regional and global climatic variation and fire intensity (e.g. burn area). This study for the first time investigates the spatial and temporal dynamics of landscape fire in the Western Himalaya region over the past 19 years (2001–2019) and assesses the role of regional and global climatic oscillations on fire activities. The overarching objectives of this study are to: (i) characterize the spatial distribution and temporal dynamics of landscape fire, (ii) assess the extent to which regional and global climatic oscillations regulate the fire regime, and (iii) investigate the influence of anthropogenic activity, such as land management practices, on fire occurrence using human modification (HM) and population density data.

2. Materials and methods

2.1 Study area

The Western Himalayas is a region rich in biodiversity but fire is a common anthropogenic and naturally driven phenomenon (Dobriyal and Bijalwan 2017)⁠. This study was restricted to two states in Western Himalaya: Uttarakhand and Himachal Pradesh (Figure 1), the most fire-prone states of the region, where the dominant forest types are evergreen broadleaf, needle leaf forest, and deciduous broadleaf forest (Reddy et al. 2015). The extension of the region is 75° 34ʹ 12” E to 81° 1ʹ 48” E longitude and 28° 42ʹ 0” N to 33° 15ʹ 36” N latitude. Surface fires are a frequent occurrence in the Chir pine (Pinus roxburghii Sarg) forest in northern regions of the Western Himalaya due to their foliage, which is highly flammable and has a large surface-to-volume ratio, and the mountainous terrain (Brown, Bhattacharyya, and Shah 2011). The topography of the region is characterized by steep mountain ridges and narrow valleys, with an elevation between 173 and 7764 m (Figure 1).The climate of the region ranges from tropical to temperate (Singh and Mal 2014; Parida, Pandey, and Patel 2020)⁠ with a warmer summer from March to June (during pre-monsoon) when most fire activity occurs. The average maximum temperature varies from 15°C to 40°C (peaks in May and June) whilst the annual precipitation ranges from 600 to 2000 mm with most falling between July and September (Mohd et al. 2017; Kumar, Shanu, and Jahangeer 2017).

Figure 1. (a) Land use and land cover (LULC) map of 2005 of Uttarakhand and Himachal Pradesh extracted from (Roy et al. 2015) and (b) ALOS World 3D (30 m) digital surface model (DSM)

2.2 Datasets

Fire activity was characterized using the MODIS burned area (MCD64 at 500 m spatial resolution; Giglio et al. 2018) and active fire (MYD14A1 at 1 km spatial resolution; Giglio, Schroeder, and Justice 2016) and the 250 m MODIS Fire (Climate Change Initiative) CCI product (Chuvieco et al. 2018). Land use and land cover (LULC) data, regional climatic data (maximum temperature and precipitation), and climatic oscillation indices (i.e. Niño-3.4, Niño-4, and Dipole mode index) were used which are summarized in Table 1. All the spatial data were processed in World Geodetic System 84 (WGS84) coordinate system.

Table 1. Data used to accomplish the vegetation fire studies over different land use and land cover type (LULC) in India

Data used	Product ID	Sensors and characteristics	Unit	Data Sources
Fire products	MCD64A1 (MODIS)	500 m, monthly	Km^2	LPDAAC https://lpdaac.usgs.gov/
	Fire CCI 5.1	250 m, monthly		ESA https://www.esa-fire-cci.org/
LULC	Decadal LULC, 2005	100 m, decadal		ORNL DAAC https://daac.ornl.gov/
DSM	ALOS World 3D (AW3D30)	30 m	meter	JAXA https://www.eorc.jaxa.jp/ALOS/en/aw3d30/
FRP	MYD14A1 (MODIS)	1 km and 8 days composite	Megawatt	ORNL DAAC https://daac.ornl.gov/
Maximum temperature (Tmax) and Precipitation (Preci)		1 °× 1° (Tmax) and 0.25 °× 0.25° (Preci), daily	°C and millimeter	Indian Metrological Department https://mausam.imd.gov.in/
SST indices	Niño-3.4	5°N to 5°S and 170°W to 120°W, monthly (1° × 1°)		NCDC, NOAA https://www.ncdc.noaa.gov/
	Niño-4	5°N to 5°S and 160E to 150 W, monthly (1° × 1°)
	DMI	Western (50°E to 70°E and 10°S to 10°N) and Eastern (90°E to 110°E and 10°S to 10°N), monthly (1° × 1°)		ESRL, NOAA https://www.esrl.noaa.gov/psd/
Gridded Population Density	GPWv411	2005,2010 and 2015 30 arc-second (~1 km)	Persons/km^2	https://sedac.ciesin.columbia.edu/data/collection/gpw-v4
Global Human Modification (GHM)	gHM-2016	30 arc-second (~1 km)		10.1111/gcb.14549

2.2.1 Land use land cover (LULC) and Digital Surface Model (DSM) data

A decadal LULC dataset (Roy et al. 2015) for the year 2005 was utilized for categorizing fires into forest, shrubland, or grassland fires. This LULC dataset contains 19 land cover classes at a spatial resolution of 100 m and, was derived from a combination of Landsat-5 Enhanced Thematic Mapper (ETM+) and India Remote Sensing satellites (IRS) Linear Imaging self-Scanning System (LISS-III) data set. The product, which is explained in detail in Roy et al. (2015), has been validated and found to have an overall accuracy of 94% and a kappa coefficient of 0.94. The regional land cover consists of six cover types (Figure 1A) comprising forest (i.e. Evergreen Needleleaf (ENF), Evergreen Broadleaf (EBF), Deciduous Broadleaf (DBF), and Mixed forest (MF)), Shrubland, and Grassland. The 30 m World 3D Topographic dataset (version 2.2), derived using Advanced Land Observing Satellite (ALOS) data (Takaku et al. 2007), was used to investigate the spatial distribution of fire intensity with elevation (Figure 1B).

2.2.2 Burned area products

Burned area spatial and temporal dynamics in Uttarakhand and Himachal Pradesh were assessed using two burned area products: MCD64A1 and Fire CCIv5.1. The MODIS 500 m burned area product (MCD64A1, collection 6) detects burned area using a hybrid approach that utilizes temporal composites of a vegetation index and active-fire detections to identify fire-affected pixels. Validation using Landsat burned area datasets indicates that the MCD64A1 product has low errors of omission (0.37 km²) and commission (0.23 km²) (Giglio et al. 2018). The MODIS-derived global burned area dataset, Fire CCI v5.1 was also used in this study which developed as part of the Climate Change Initiative (Chuvieco et al. 2018) makes use of the 250 m spatial resolution red and near-infrared wavebands of MODIS onboard Terra in addition to active-fire information to detect burned area. This product has higher errors of omission and commission (0.4 km² and 0.4 km², respectively) than the MCD64A1 product but the improved spatial resolution facilitates the detection of small (< 100 ha) burned area (Chuvieco et al. 2018). Monthly burned area data covering the main fire season (March–June) was downloaded for both products for each year (19 years) to assess the total pre-monsoonal burn area over the region.

2.2.3 Fire radiative power data

Fire Radiative Power (FRP, MW) retrievals, which describes the rate of radiative energy release during combustion, from the Aqua MODIS Thermal Anomalies product (MYD14A1) were employed to investigate the influence of land cover type on fire characteristics. MYD14A1 is a level-3 composite product that provides daily active-fire detections and FRP retrievals at 1 km² spatial resolution over 8 days. The FRP data between March and June from 2003 to 2019 were extracted for further analysis. FRP is retrieved in the collection of 6 product using the middle-infrared (MIR) radiance method present in equation-1 (Wooster et al. 2005).

(1) $FRP=\frac{P_a.\sigma .\epsilon }{.{\epsilon }_{MIR}}\left(L_f-L_b\right)$(1)

where $L_f$and $L_b$ are the MIR spectral radiance of the active fire and ambient background pixels (W m⁻2 sr⁻¹ µm⁻¹) respectively, σ is the Stefan-Boltzmann constant (5.67x10⁻⁸ J s⁻¹ m⁻² K⁻⁴), α denotes a sensor-specific empirical constant (α = 3.0 × 10 − 9 W m⁻² sr⁻¹ μm⁻¹ K⁻⁴), ε and ε_MIR are the broadband and middle-infrared spectral emissivity’s, respectively, and $P_a$ is the pixel area (m²). FRP retrievals have been shown to be related to the fuel consumption rate (Smith and Wooster 2005) and smoke and trace gas emissions (Freeborn et al. 2008).

2.2.4 Gridded human population and global human modification (HM) data

Previous research suggests that anthropogenic activity contributes to ~ 90% of the fire ignitions in the study region (Dobriyal and Bijalwan 2017). Gridded population density (GPWv411) data for 2005, 2010, and 2015, estimated from extrapolation of census data (Morton 2019), and a global Human Modification (gHM) dataset from 2016 (Kennedy et al. 2019) were used to investigate the relationship between fire occurrence and the scale of anthropogenic activity. The gHM index, which varies between 0 and 1, provides a cumulative measure of human land modification using metrics including population density, built-up area, transport and utilities (road, railways, and power lines) networks, nighttime lights, and mining locations. Both datasets are available in 30 arc-second (~1 km) spatial resolution.

2.2.5 Climate data

The majority of fires in the region are caused by anthropogenic activity (Dobriyal and Bijalwan 2017) whilst environmental and climatic variables contribute to creating the conditions conducive for fires to spread. The role of regional climatic conditions on fire occurrence was assessed using the Indian Meteorological Department (IMD) precipitation data (0.25°×0.25°) and maximum temperature (1°×1°) dataset, downloaded from National Data Center (NDC), Pune. The gridded daily precipitation data, discussed in detail by Pai et al.(2014), was produced by the National Climate Center using 7000 rain-gauge stations data. Inverse Distance Weighted (IDW) interpolation was used to generate a gridded dataset with the weights based on the distance between the station and the grid point (IMD 2005). Daily maximum temperature (°C) and precipitation (mm) data during the pre-monsoon (March to June) were extracted to assess their influence on regional fire occurrence.

2.2.6 Global climatic oscillation indices

The regional weather and climatic cycles are governed by macro-scale global climatic oscillations (Mason et al. 2017). It is well established that in the last few decades, climatic extremes have enhanced fire activity around the world (Doerr and Santín 2016; Mason et al. 2017). In this study, the relationship between burned area and the Niño-3.4, Niño-4, and Dipole Mode Index (DMI) climate anomaly indices were investigated. The Niño-3.4 and Niño-4 are indices that characterize the average Sea Surface Temperature (SST) anomalies. The Niño-3.4 represents the average state of equatorial Pacific SSTs from the dateline to the South American coast (5°N – 5°S and 170°W – 120°W) whereas Niño-4 covers central equatorial Pacific (5°N to 5°S and 160°E – 150°W) (Trenberth and Stepaniak 2001). Typically, these indices calculation uses average SST over the region then a consecutive 5 monthly running mean method was performed to compute the monthly anomaly. The Niño-3.4 has been used to define El Niño (> 0.4) and La Niño-3.4 (< −0.4) events whilst the Niño-4 region typically has less variation than other Niño regions (Trenberth 2020). The Indian Ocean Dipole (IOD) is a tropical ocean-atmosphere phenomenon that leads to climatic variability, particularly over South-East Asia and the Indian continent (Ashok, Guan, and Yamagata 2001), East Africa, and Australia. The DMI (e.g. quantitative measures of IOD phenomenon) is defined as the SST differences between the western (from 50°E to 70°E and 10°S to 10°N) and eastern (from 90°E to 110°E and 10°S to 10°N) regions of the Indian Ocean (Saji et al. 1999)⁠. The monthly anomaly of these indices was downloaded and the median pre-monsoonal (March–June) anomaly was calculated to assess its relationship to fire anomalies over the same period.

2.2.7 Field observations of fire activity

A ground survey, consisting of 25 burn sites, was carried out between 28^th May to 2^nd June (2019) over parts of Uttarakhand Himalaya to document the site and effects of fire on the landscape. According to the Forest Survey of India (FSI), Nainital, Almora districts are the most fire-affected areas of the south-eastern parts of Uttarakhand (ISFR. 2017). The field sites were selected based on the fire risk maps from the FSI over Nainital and Almora districts. The survey plots were defined using three categories, (i) burnt (burnt a few days prior to the campaign), (ii) burning (actively smoldering or flaming), and (iii) unburned. The field survey was conducted during the peak fire period and documented the type, condition, and thickness of the surface fuel layer, effects of the fire, fire spreading mechanism (surface or crown fire), and fire management. The thickness of the fuel was measured by the graduated scale and information on the fire spread was estimated using photographs.

2.3 Methods

All the spatial datasets were initially processed in Python-2.7 environment using the “ArcPy” package available in ArcGIS-10.3. The basic statistical analysis and inferential statistical methods were conducted using the R-programming platform (R Core Team 2020). All the maps and figures were made in R using different packages (e.g. raster, rgdal, rasterVis, ggplot2, ggpubr, plot3D, zpy, spdep, spgwr, and tmap, among others). The conceptual flowchart of the methodology is presented in Figure 2.

Figure 2. Methodological flow of spatio-temporal variability of forests, shrubland, and grassland fire and its relationship with regional climatic and anthropogenic factors and global climatic oscillations

2.3.1 Temporal analysis of fire

The focus of the study is to analyze the spatial and temporal dynamics of fire occurrence in the Western Himalaya and its relationship with regional climatic conditions and anthropogenic activity. The analysis of the temporal trends in the burned area over the past 19 years was carried out using the Mann-Kendall (MK) nonparametric statistical test which is more robust to outliers and sudden breaks in nonhomogeneous time-series (Önöz and Bayazit 2003). Sen’s (Sen 1968) slope estimator was applied to estimate the magnitude of burn area trend. The “YuePilon” method (Yue et al. 2002) was adopted which reduces the influence of autocorrelation over time-series. The length of the fire season and distribution pattern are important indicators of fire regime and complex fire dynamic (Rogers et al. 2020; National Research Council 2010). The length and distribution pattern of fire during pre-monsoon was analyzed using a normal probability function over two-time periods (2001–2010 and 2011–2019). An equality test between fire distribution in 2001–2010 and 2011–2019 was determined using the nonparametric Kolmogorov–Smirnov test (K–S test, Marsaglia, Tsang, and Wang 2003).

2.3.2 Spatial fire frequency and burn area fraction

The pre-monsoonal (March–June) burn area of Fire-CCI was used to calculate the annual fire occurrence in a 1 × 1 km grid over a 19-years’ (2001–2019) time span. The burn area fraction (%) within a 1 km² grid was calculated using the 250 m Fire-CCI burn area data and the average (2001–2019) burn area fraction determined using the annual burn fraction. These data were used to characterize the gradient of fire intensity across the landscape.

2.3.3 FRP Hot/Cold spots mapping

FRP magnitude has been found to vary with the land cover type (Vadrevu et al. 2015; Chuvieco et al. 2018) and this was investigated in this study using the mean per-pixel FRP (2003–2019). The spatial patterns of the FRP were investigated in relation to land cover type by using the Getis-Ords Gi* spatial statistic to derive hot and coldspot clusters (Getis and Ord 2010). When applied to FRP retrievals, the Getis-Ords function describes the spatial autocorrelation of higher and lower “intensity” active fires. A detailed description of Getis-Ords Gi* is available in Reddy et al. (2019).

2.3.4 Correlation and regression analysis

The pre-monsoonal burn area was correlated and regressed with selected climatic variables (i.e. maximum temperature and precipitation). The Pearson correlation coefficient (R) (Nahler 2009) statistical function was applied between climatic variables and burn area to quantify the linear association (R) ranging from – 1 to +1. Linear least square regression was employed to compute the contribution (regression coefficient or β) of anthropogenic (population density and human modification), regional (maximum temperature and precipitation), and global climatic indices in fire intensity (burn area).

2.3.5 Standard anomaly estimation

The pre-monsoonal (median of March–June) global climatic oscillations indices (Niño-3.4, 4 and DMI) are analyzed in relation to pre-monsoonal burn area standard anomaly (SA) over 19-years (2001–2019). The SA was estimated using equation 2(2) $\mathrm{S}\mathrm{A}=\frac{{\mathrm{x}}_{\mathrm{i}-\bar{\mathrm{x}}}}{\delta }$(2) .

(2) $\mathrm{S}\mathrm{A}=\frac{{\mathrm{x}}_{\mathrm{i}-\bar{\mathrm{x}}}}{\delta }$(2)

where $x_i$ is defined as the individual observation of burn area, $\bar{x}$ is the mean burn area (2001–2019) of the whole observation and $\delta$ represents the standard deviation.

3. Results

3.1 Spatial and temporal variation in fire

Over the past 19 years (2001–2019), the average annual pre-monsoonal burned area in Uttarakhand and Himachal Pradesh is 867 km² (sd = 780 km²) and 120 km² (sd = 124 km²), respectively, which equates to 1.6% and 0.2% of the land surface area of each state, respectively. There was a significant inter-annual variation in the burned area (Figure 3A-B) which ranges between 3 and 451 km² and 50–3083 km² in Himachal Pradesh and Uttarakhand, respectively. The overall trend (Table 2 and Supplement 1) indicates that burned area is increasing by 2.6 and 46.7 km² year⁻¹ in Himachal Pradesh and Uttarakhand, respectively. DBF and EBF are the most fire-prone land cover types in Himachal Pradesh and account for 42 and 32% of the total burned area, respectively. In Uttarakhand, ENF and DBF account for 43 and 41% of the total burned area. In Uttarakhand and Himachal Pradesh, the average fire size is small at 0.30 and 0.24 km² whilst the largest burned area detected was 100.4 and 38.3 km², respectively.

Table 2. Trends of burn area (including forest, shrubland, and grassland) over the period 2001–2019 as per the Fire_CCI fire data in Uttarakhand and Himachal Pradesh. * indicates over 90% confidence level (CL)

Name of the state	MK Tau (τ)	Two-sided p value	Sen’s slope (km^2 year^−1)	Linear trend (km^2 year^−1)	p value
Himachal Pradesh	0.14	0.3	2.64	4.67	0.25
Uttarakhand	0.43*	0.09	46.78*	56.53*	0.02

Figure 3. Interannual variability of (a) pre-monsoonal total and (b) land cover wise (DBF: Deciduous Broadleaf Forest, EBF: Evergreen Broadleaf Forest, ENF: Evergreen Needleleaf Forest, Grass: Grassland, MF: Mixed Forest and Shrub: Shrubland) burn area (km²) over Himachal Pradesh (left panel), and Uttarakhand (right panel) from 2001 to 2019 using Fire-CCI product

Fire occurrence among different land cover types over western Himalaya is shown in Supplement 2A and 2B and their temporal variation in Figures 3C and Figures 3D. It is clear from Figure 3D, that there has been a large increase in burned area in Uttarakhand between 2011 and 2019. Notably, ENF cover has an 87% increase in the burned area relative to that which occurred between 2001 and 2010. A similar inter-annual burn area pattern has also been noticed in the MCD64A1 burn area dataset (Supplement 3). Burned area is typically low in shrublands and grasslands with DBF being the only land cover type with a moderate degree of fire activity in both regions. Most fires occur in the dry season between March and June where daily total rainfall is < 50 mm and the average maximum temperature is 31.5°C. The nonparametric K-S test revealed that the distribution of burned area in 2001–2010 and 2011–2019 was not equal or significantly dissimilar. The absolute maximum difference of cumulative probability density is 0.19 (p < 0.01) (Figure 4a). The decadal daily fire distribution (Figure 4b) shows a bimodal peak in the burned area during pre-monsoon which results from fires occurring in the last week of March and the last week of May to the 1^st week of June. The length of fire season, defined as decadal daily start-up days to recession of fire time (shown as the normal distribution curve of Figure 4b), is ~ 10 days shorter in 2011–2019 compared to 2001–2010 although the burned area is concentrated over a shorter time period (Figure 4c).

Figure 4. (a) Kolmogorov–Smirnov or K-S test between 2001–2010 and 2011–2019 daily fire distribution where d and p value denote as max absolute difference of cumulative probability density and significant level, (b) probability density curve of daily fire events distribution, and (c) boxplot between fire occurrences of two time slots. The blue dot in boxplot indicated the mean of the distribution

Analysis of the fire frequency between 2001 and 2019 (Figure 5A) indicates that fire occurrences in the region are largely sporadic with the highest fire frequency of 8 over a 19-year period. Only 0.1% of burned area detections were detected eight times over the 19-year period and these fires predominantly occur in ENF and DBF. Analysis of the average burned area fraction within 1 km² area (Figure 5B) indicates that 50% of the fire patches having a burned area fraction > 30% (sd = 25%) and these are predominantly found in southern Uttarakhand. There appears to be a little variation in fire frequency with altitude and slope with 70% of the burned area occurring at altitudes between 550 and 1750 m and slopes of 10–40° (Figure 5C).

Figure 5. (a) Spatial fire frequency and (b) average burn area fraction (%) over Himachal Pradesh and Uttarakhand from 2001 to 2019. (c) Altitude, ground slope and land cover (DBF: Deciduous Broadleaf Forest, EBF: Evergreen Broadleaf Forest, ENF: Evergreen Needleleaf Forest, Grass: Grassland, MF: Mixed Forest and Shrub: Shrubland) wise fire frequency distribution

3.2 FRP variation with land cover type and intensity distributions

The spatial distribution of the hot and cold clusters derived using the z-score from the Getis-Ords function (Figure 6a) applied to the mean pixel FRP indicates that the spatial patterns of the FRP in Himachal Pradesh are more irregular and isolated than that found in Uttarakhand. In the latter, the spatial distribution of FRP displays more spatial clustering and a greater range in FRP magnitude with a number of cold (–2 to – 4) and hot (2–10) clusters. Areas where the z-score centered on zero were related to fewer fires with large uncertainties (p value > 0.6). Histograms of the Gi* or Z score (Supplement 4) for different land cover types indicate that ENF and shrublands have the widest range in z-score and the highest frequency of hot clusters. The altitude and ground slope wise distribution of FRP indicated that fire events were concentrated in ~ 2000–3800 m altitude and 45–65° slope zone had higher energy (Figure 6b). There is a little variation in the average FRP with slope angle (ranges between 10 and 40°) suggesting that most of the forest fires are surface rather than crown fires which are supported by observations from the field assessments (Supplement 7). Figure 6c presents a box plot showing the range in FRP magnitude with land cover type for each state which indicates that fires in Himachal Pradesh typically have a greater range and higher FRP magnitude for most land cover types than Uttarakhand (Supplement 5A and 5B). The FRP magnitude for different land cover types is broadly similar but with mixed forest and ENF have higher FRP.

Figure 6. (a) Hotspot and clodspot spatial analysis over Uttarakhand and Himachal Pradesh by using the mean FRP data over the 16 years (2003–2019), wherein Getis-Ord Gi* statistical function was applied. The Z-score represents the magnitudes of spatial clustering (hot and clodspot) of FRP. The Z-score is statistically significant at p value < 0.1 when Z-score exceeds ± 1. (b) Slope and Altitude wise FRP distribution. (c) Forest type-wise (DBF: Deciduous Broadleaf Forest, EBF: Evergreen Broadleaf Forest, ENF: Evergreen Needleleaf Forest, Grass: Grassland, MF: Mixed Forest and Shrub: Shrubland) FRP variability over Himachal Pradesh and Uttarakhand from 2003 to 2019. The blue dot represents the mean FRP (MW)

3.3 Gridded population density, global human modification (GHM), and spatial coverage of fire

In this section, the spatial association between the human footprint and fire intensity (burn area) was investigated over Himachal Pradesh and Uttarakhand. The demi-decadal mean population density in the region increased from 148 to 173 people per km² between 2005 and 2015 with the highest population densities found at lower elevations (Figure 7a). According to human modification intensity, increased human activity is found along the southern and central parts of both states (Figure 7d) which is also the most fire-prone areas (Figures 5 - Figure 6a). The spatial correspondence between the burn area fraction (%) and population density (Figure 7b-c) and human modification intensity (Figure 7e-f) at varying altitudes (m) and slopes (°) is shown as a contour plots. About 40–60% of burn area pixels were associated with moderate population density (200–1000 persons/km²) at lower elevations (< 1000 m) and moderate slopes (10°-50°). The human modification ranged between 0.15 and 0.8 in western Himalaya, but areas of higher burned area fraction corresponded to 0.2–0.5 human modification intensity. The correspondence between mean burn area fraction (%) and human modification and population density is with low coefficients of determination (R²) of 0.12 and 0.13, respectively. Comparatively, the human modification (β coefficient = 39.58%) exhibited a significant (p < 0.01) higher contribution in fire occurrences than population density (β coefficient = 0.089%) (Supplement 6).

Figure 7. (a) Gridded population density (persons per km²) of 2015, (b–c) association of burn area fraction (%) and population density (2015) in varying altitude (m) and slope (°). (d) Human modification (HM) intensity (2016), value of HM ranges from 0 (unmodified) to 1 (maximum modification). (e–f) Association of burn area fraction (%) and intensity of human modification in gradient of altitude and ground slope

3.4 Field-based fire characterization

The mountainous physiography of the Western Himalaya makes it more difficult to control and regulate fire events. The photographs of the field observation were presented in Supplement 7 and the corresponding coordinates and altitude of these sites were tabulated in Supplement 8. The ground observation suggested a good deposition of inflammable dry Chir pine leaves (~ 10–12 inches thickness measured by graduated scale), which helps to ignite (Supplement 7b). Creeping is a mostly dominated type of fire (Supplement 7 c–d) where an ample supply of ground fuel support for spreading of the fire line. To immediately control or put out the fire, the local authority has used the rudimentary method (Supplement 7d). Due to ground or creeping fire, the lower parts of trunks are partially or fully burnt (Supplement 7e–f).

3.5 Climate variability and fire occurrence

A daily average burn area from March to June (i.e. 60–180 days of the year) computed over the study area and related with the average maximum temperature (°C) and total precipitation (mm) are presented in Figure 8A. The maximum temperature increases from March to June which coincides with periods of time when precipitation is low, during which time most fire activity occurs. A statistical analysis between the pre-monsoonal total burn area (March–June), cumulative precipitation, and average maximum temperature for both regions is shown in Figure 8B. The correlation between the burned area and average maximum temperature (R = 0.57, p = 0.01) and precipitation (R = −0.53, p = 0.2) is moderate. The average maximum temperature is positively related to the burn area whilst the total precipitation is negatively related to the burn area. However, it is evident that there is considerable variability in rainfall and the average maximum temperature with the total burned area (Supplement 9, Supplement 10). The relationship between inter-annual pre-monsoonal burn area standard anomaly (SA), and the global oscillation indices (Niño-3.4, Niño-4, and DMI) was shown in Figure 9. The association between the burned area and the climate anomalies exhibits a moderate correlation with Niño-3.4 (R = 0.51 and p = 0.03) (Figure 9A) whilst the Niño-4 (R = 0.37) and DMI have a weak agreement (R = 0.1) (Figure 9B and C).

Figure 8. (a) Association among daily mean burn area (km²), maximum temperature (°C) and precipitation (mm). (b) statistical relation of burn area and maximim temperature and precipittaion over Himachal Pradesh and Uttarakhnad

Figure 9. Association and linear relationship between (a) Niño-3.4, (b) Niño-4, and (c) DMI and burned area standard anomaly from 2001 to 2019 in western Himalayan (Uttarakhand and Himachal Pradesh). The boxplot present the axis-wise distribution of fire, Niño-3.4, Niño-4, and DMI anomlay

4. Discussion

The spatial characterization of landscape fire in the Western Himalaya revealed the average pre-monsoonal burned area has increased by 72.94 km² year-¹ with the greatest rate of the increase occurring between 2011 and 2019. The burned area also displayed significant inter-annual variation with peaks in burned area evident in 2005, 2008, 2009, 2012, 2016, and 2018 along with a cohesive agreement with maximum temperature and precipitation (Supplement 9). Other studies have also found an increasing trend in the burned area although in different regions of western Himalaya (Dobriyal and Bijalwan 2017; Srivastava and Garg 2013). This study found that Uttarakhand accounts for 88% of the fire incidence in the Western Himalaya with most fires occurring in ENF and DBF landcover. The increasing trend in the burned area is also found in the active-fire detections which has been increased by 218 fire counts over the time period. Whilst this is a low number relative to many fire-affected regions, they could pose a threat to the biodiversity of the Western Himalaya (Chitale and Behera 2019). Analysis of the temporal distribution of fire occurrence revealed the length of the fire season reduced slightly (almost 10 days) in 2011–2019 compared to 2001–2010 with fire events being more concentrated. The decadal daily fire occurrences depicted a bimodal distribution. The first peak was in late March due to the huge accumulation of fresh litter, and the second peak between late May and 1^st week of June due to higher temperature, low humidity, and extreme dryness (Chauhan, Gautam, and Negi 2018; Singh et al. 2014). The highest concentrations of the burned area were found in the central and southern parts of both states over 550 to 1750 m altitude.

Analysis of MODIS FRP retrievals found similar mean FRP magnitudes and distributions with land cover type. The dominant tree species in ENF is Chir pine (Pinus roxburghii Sarg) whose needles are highly flammable due to their resin which accelerates the ignition and fire intensity. However, the low tree density (Supplementary 7) and lack of ladder fuels to instigate crown fires (Brown, Bhattacharyya, and Shah 2011) typically result in low-intensity surface fires. The median FRP (28 MW) in ENF is marginally higher than that found in other land cover types which suggests fire occurrences were surface rather than a crown as the latter typically result in higher FRP magnitudes (Wooster and Zhang 2004).

The southern and central (lower elevated) districts of Uttarakhand and Himachal Pradesh have higher population density (200–1800 persons/km²) and increased anthropogenic activities in comparison to the northernmost higher elevated regions (Sati 2014), and in these regions, the fire occurrences were greatest. These findings are well corroborated with the earlier findings by Babu, Roy, and Prasad (2016) who concluded that high fire risk was co-located in the elevation ranges from 173 to 2000 m. Most forest fires (around 90%) in this region result from human-made (intentional or accidental) (Chauhan, Gautam, and Negi 2018; Dobriyal and Bijalwan 2017; Liu, Goodrick, and Heilman 2014) but environmental and climatic variables also play a very important role (Srivastava and Garg 2013). Here we found moderate to week linkages between regional climatic variables (temperature and precipitation) and the burned area which is in line with the study by Vadrevu et al. (2019) but over the south and south-east Asia. Other studies have also found limited agreement between fire occurrences and regional climate (Srivastava and Garg 2013) which may be masked by the other environmental and climatic factors (such as topography, forest type, fuel characteristics, wind, relative humidity, and solar radiation, among others) and strong anthropogenic influence on fire intensity.

Analysis of the relationship between a number of climate indices and burned area revealed some correspondence with warmer (El Niño) and cooler (La Niña) periods with anomalies in fire occurrences in western Himalaya. Peaks in the burned area were found to occur in some of El Niño years (e.g. 2005, 2006, 2016, and 2018) although there is some inconsistency which was also observed by Vadrevu et al. (2019) over the south and southeast Asia. The highest correspondence was found between the Niño-3.4 index whereas the Niño-4 and DMI had a weak correspondence with fire intensity which has been attributed to regional affecting factors (Srivastava and Garg 2013). These findings are in line with the literature that the El Niño correlates with forest fires (Fuller and Murphy 2006; Skinner et al. 2006; Simard and Haines 1985) but in a different region. This tendency seems to be a response to the increase in fires favored by the high temperature, deficit rainfall, and drought, which are also the crucial characteristics of El Niño years (Alencar, Nepstad, and Diaz 2006).

This study used burned area and active-fire data to investigate the spatial and temporal dynamics of fire activities and the role of the local meteorology and climate in these dynamics over western Himalaya. The moderate spatial resolution of the fire observation (250 m and 500 m) may omit smaller and less intense fires that are more difficult to detect, particularly when burning under forest canopies (Pereira et al. 2004). Future investigations should assess the utility of higher spatial resolution burned area (e.g. from Sentinel-2) and active-fire data (e.g. from Visible Infrared Imaging Radiometer Suite) to monitor landscape fire in the region to determine the level of omission when using moderate resolution data. Higher spatial resolution meteorological data, currently to 1° (temperature) and 0.25° (precipitation), may improve understanding the impact of local meteorology on fire intensity and how this links with anthropogenic activity. The decadal LULC data have limited availability and, to perform a forest type and other landcover classification a good number of ground samples were required that is why the study was only used LULC of 2005 which could reduce the precision of categorical burn area segregation in extreme years.

5. Conclusions

The space-time dynamics of western Himalaya’s forest fire landscape has peculiar characteristics. The present study assesses the spatio-temporal variability of forest, shrubland, and grassland fire in the two most fire-affected states of western Himalaya (Uttarakhand and Himachal Pradesh) over a 19 year period and its relation to regional climatic variability (maximum temperature and precipitation) and global climatic oscillations (Niño-3.4, Niño-4 and DMI). The overall trend in the burned area estimated from the Fire-CCI and MCD64A1 datasets is one of the increasing fire occurrences which presents a risk to the fragile forest environment of the region. Whilst the trend in the burned area is increasing there is considerable variability both inter-annually and within forest type. Analysis of MODIS FRP retrievals indicates little variation in mean FRP with land cover type which may be indicative of the fires being surface rather than crown fires due to the lower tree cover density and lack of ladder fuels. However, the spatial clustering using Getis-Ords found most hotspots are associated with ENF which may result from the prevalence of surface needles in these forests. Most burned area is found at lower elevations on the southern slopes of western Himalaya which falls in an intermediate region with respect to the level of human modification and population density. The regional climatic conditions (maximum temperature and precipitation) were positively and negatively related to the total burned area although further analysis is needed to disentangle the role of climate and anthropogenic factors on fire occurrence and to better understand the impact of fire on this fragile mountainous environment.

Highlights

1. Forest fire burn area has increased by 72.94 km² yr⁻¹ over Western Himalaya.

2. Length of the fire season in 2010–2019 was ~10 days shorten in comparison to 2001–2010.

3. Temperature has greater control in the fire landscape than precipitation.

4. Fire occurrences were moderate to weakly in correspondence with regional climatic and human activities.

5. Creeping was a major fire type to induce tree mortality and forest distraction.

Acknowledgments

The authors are thankful to the European Space Agency (ESA), Earth Data of NASA for Fire CCI, MCD64A1 fire datasets, and IMD climatic data. First author thanks to UGC for assisting NET-JRF fellowship (3289/(SC NET-JAN 2017) for doctoral research. We are also thankful to Nanital District Forest Officer (DFO) for their eminent help during fieldwork.

CRediT authorship contribution statement: S.B and B.R.P.: Conceptualization, Investigation, Methodology, Software, Analysis, Writing and editing–original draft. A.C.P. and P.A.: Conceptualization, Supervision, Methodology, Software, Writing – review and editing. G.R and J.D.: Conceptualization, Supervision, Data support, Writing – review and editing.

Disclosure statement

The authors declare no potential conflict of interest.

Supplementary material

Supplemental data for this article can be accessedhere.

Additional information

Funding

This work was supported by the University Grants Commission [3289/(SC NET-JAN 2017].