Assessment of Post-Fire Vegetation Recovery Using Fire Severity and Geographical Data in the Mediterranean Region (Spain)

Abstract

Wildfires cause disturbances in ecosystems and generate environmental, economic, and social costs. Studies focused on vegetation regeneration in burned areas acquire interest because of the need to understand the species dynamics and to apply an adequate restoration policy. In this work we intend to study the variables that condition short-term regeneration (5 years) of three species of the genus Pinus in the Mediterranean region of the Iberian Peninsula. Regeneration modelling has been performed through multiple regressions, using Ordinary Least Squares (OLS) and Geographic Weight Regression (GWR). The variables used were fire severity, measured through the Composite Burn Index (CBI), and a set of environmental variables (topography, post-fire climate, vegetation type, and state after fire). The regeneration dynamics were measured through the Normalized Difference Vegetation Index (NDVI) obtained from Landsat images. The relationship between fire severity and regeneration dynamics showed consistent results. Short-term regeneration was slowed down when severity was higher. The models generated by GWR showed better results in comparison with OLS (adjusted R² = 0.77 for Pinus nigra and Pinus pinaster; adjusted R² = 0.80 for Pinus halepensis). Further studies should focus on obtaining more precise variables and considering new factors which help to better explain post-fire vegetation recovery.

1. Introduction

Wildfires are one of the most important environmental problems at present. In the European context, Spain is one of the countries which registers the highest fire incidence in number and surface burned [1,2,3]. These events cause a great amount of damages (i.e., degradation of soils, water and biodiversity), which may further lead to economic and even human costs. From an environmental point of view, wildfires cause alterations of landcover, loss of carbon reserves, and changes in the soil composition and its hydrogeomorphologic behavior [4]. Nevertheless, Mediterranean vegetation is quite adapted to fire disturbance. Mediterranean species have post-fire ecological strategies, like the resprout ability, the seed bank persistency, or the growth or dispersal ability [4]. In this sense, the analysis of Large Forest Fires (LFF) acquires a special interest because the effects they provoke are devastating. For this reason, the application of an adequate restoration policy requires the exhaustive study of the physical environment and dynamic evolution of the affected area [5].

Over the last few years, several studies have underlined the importance of remote sensing in analyzing ecological dynamics following fire and studying post-fire vegetation regeneration [5,6,7]. Compared to field surveying, satellite images offer a less expensive alternative and provide broader information of burned areas by obtaining biophysical variables of wildfires [7].

Medium-resolution optical sensors such as MSS, TM, ETM+, and OLI of Landsat series enables the monitoring of vegetation in burned areas for more than 40 years [8,9,10]. There are also researches studying vegetation regeneration using other sensors such as AVIRIS [11], AVHRR [12], MODIS [13], SPOT-VEGETATION [14], or RADAR [15]. According to Hirschmugl et al. [16], time series analysis is the most used approach in the monitoring of forest disturbances and degradation, as well as regeneration processes with the objective to analyze spectral variations in the forest cover [8,11,17].

The factors which define the vegetation regeneration rate after a wildfire are multiple, and their complete identification or modeling is difficult. These factors are related to wildfire characteristics, environmental conditions, and the life history of plant species [18]. Some studies applied to northern regions [10,12,19,20] and Mediterranean regions [8,11,21] have pointed out some of the main drivers of vegetation regeneration following forest fires. Fire severity levels, topography (elevation, slope, and orientation), post-fire climate, or vegetation cover class are the most used in regeneration estimates. Other studies have emphasized the influence of solar radiation on the water availability for vegetation growth [8] or the effects of applying different restoration models [9]. Several studies have revealed that the influence of environmental factors on regeneration can vary across vegetation types [8,20,22].

In order to evaluate fire severity, defined as the magnitude of the ecological change produced by fire [6], different methods have been proposed. Some studies have successfully applied vegetation indices such as the Normalized Difference Vegetation Index (NDVI) [23], due to the relationship between the amount of vegetation consumed and fire severity [7]. On the other hand, specific indices have been developed which record with greater spectral contrast the fire effects as the Normalized Burn Ratio (NBR) [24] subsequently modified [25,26]. Another method based on parametrized variables to estimate fire effects in vertical strata is the Composite Burn Index (CBI) [25], widely and effectively applied with good results (correlation between CBI and dNBR R² = 0.83) [27] and modified by De Santis and Chuvieco [28].

To assess vegetation regeneration using time series, several methodologies have also been proposed. The most common is to monitor the vegetation state from spectral indices. Among the most used indices are the Normalized Difference Vegetation Index [17,27,29], the Regeneration Index [11,14,19], the Normalized Difference Infrared Index [30], and the Soil-Adjusted Vegetation Index [31]. Spectral Mixture Analysis has also effectively been applied in that context [8].

Although several studies have evaluated fire severity and vegetation regeneration using remote sensing data, few have integrated both aspects in a single study. Several studies have identified the significant influence of fire severity [21,32] and environmental factors, such as topography and climate [8,17], in vegetation regeneration. Therefore, the need to integrate both analyses to improve the predictive models of these dynamics has arisen. In recent years, more studies have been undertaken to investigate the factors that determine post-fire regeneration patterns with satisfactory results [19,20].

In this study, we hypothesize that fire severity and environmental variables such as vegetation type, meteorology, and topography determine the post-fire vegetation regeneration. Therefore, regeneration patterns on burned surfaces will vary between areas presenting different severity levels. Furthermore, meteorological post-fire conditions and topography will have different impacts on the regeneration of different forest covers.

The general objective will be to model the short-term vegetation regeneration (five years after fire) in three large forest fires in the Mediterranean region of the Iberian Peninsula, by knowing the fire severity and the interacting environmental variables. In this case, we will model the forest regeneration of the genus Pinus, which is widely extended in Mediterranean forests. In addition, Pinus communities are one of the most affected by fires due to the high content of resins that promote both fire start and spread [33].

Specific objectives of the study area include: (1) To model the evolution of post-fire vegetation to obtain regrowth patterns through spectral indices; (2) Generate environmental variables involved in regeneration and (3) Identify the variables which are most relevant in explaining the short-term regeneration using multiple regression models and to estimate the short-term regeneration of Pinus species.

2. Materials and Methods

2.1. Study Area

The study is based on three large fires that occurred in the summer of 1994 (Figure 1) in the municipalities of Castrocontrigo (Province of León), Uncastillo (Province of Zaragoza), and Moratalla (Province of Murcia). These areas have been selected according to several criteria:

(i)

The burned area exceeds the 500 ha lower limit to be considered a large forest fire in Spain [34].

(ii)

Burned areas where post-fire restoration activities were not implemented, and therefore did not alter the natural regeneration of the burned area.

(iii)

Burned areas which have not been burned in subsequent fires since 1994.

These areas are located in the Mediterranean biogeographic region [35]. Topographic and climatic characteristics are quite different; the elevation varies between 1100 m in the highest levels in the Castrocontrigo and Uncastillo fires to 1400 m in the Moratalla fire. On the other hand, mean annual rainfall differs between Moratalla fire (ranging from 300 to 400 mm) and Castrocontrigo or Uncastillo fires (recording from 600 to 700 mm), while mean annual temperatures range from 10 °C in Castrocontrigo, 15 °C in Uncastillo, to 17.5 °C in Moratalla [36]. In the summer season, the scarcity of precipitation and high temperatures, characteristic of the Mediterranean climate, increase the risk of fire.

All study areas are dominated by anthropogenic coniferous forests, mainly with species of the genus Pinus along with certain deciduous species of the genus Quercus. In the shrub strata, sclerophyll species such as Rosmarinus, Thymus, or Juniperus species have been identified. Nevertheless, there are latitudinal and altitudinal differences between the three burned areas, which generate particularities that are described in the Map of Vegetation Series [35]. The area of Castrocontrigo fire is placed in the Supramesomediterranean series dominated by Pinus pinaster and Pinus sylvestris among Quercus pyrenaica. The Uncastillo area is located on the limit between the Supramesomediterranean and Mesomediterranean series. Therefore, most of this area is covered by different species of pines (Pinus halepensis, Pinus nigra, and Pinus sylvestris). The area of Moratalla is in the Mesomediterranean series, an area with a more arid climate than the previous ones, where Pinus halepensis is dominant.

As can be seen, the selected study areas represent a broad gradient of environmental conditions representative of the Mediterranean ecoregion. Therefore, the models derived from the analysis of these variables can be applied to a wide variety of Mediterranean environments.

Table 1. Fire data of the studied forest fires [37].

Municipal Term of Origin	Start Date	Total Forest Area (ha)	Wooded Area (ha)	Non-Wooded Area (ha)	Non-Forest Area (ha)	Cause
Castrocontrigo	22/08/1994	1150	500	650	0	Unknown
Uncastillo	16/07/1994	6589	4849	1740	2049	Lightning
Moratalla	04/07/1994	24,817	16,262	9555	2828	Powerline

shows the forest fires data related to total burned area, type of vegetation affected, cause, and date.

2.2. Datasets

In this research, several data sources have been used to generate the spatial information used to model the vegetation regeneration. Satellite images have been obtained for the respective fires and spatial data have been collected to generate the environmental variables.

Landsat TM and Landsat ETM+ images were obtained for the pre- and post-fire period between 1994 and 1999 from the EOLI-SA server of the European Space Agency (ESA) (

Table 2. List of images used.

Castrocontrigo		Uncastillo		Moratalla		Sensor
Date	Days since Fire	Date	Days since Fire	Date	Days since Fire
08/07/1994	−45	06/07/1994	−18	29/06/1994	−6	TM
28/08/1994	6	23/08/1994	30	15/08/1994	38	TM
16/08/1995	359	10/08/1995	382	18/07/1995	375	TM
01/08/1996	709	11/07/1996	718	21/08/1996	774	TM
05/08/1997	1078	14/07/1997	1086	07/07/1997	1094	TM
08/09/1998	1477	21/10/1998	1550	07/05/1998	1398	TM
01/07/1999	1773	29/08/1999	1910	21/07/1999	1845	ETM+

).

Images radiometric resolution is 8 bits with visible spectral information (RGB), near infrared (NIR), and short wave infrared (SWIR1 and SWIR2). Images spatial resolution is 30 m, so all the information used in the study will be adapted to that resolution. The selected scenes belong to the months of July and August, with the exception of images which have been chosen from the end of June, beginning of September, or October. The selection of similar dates was prioritized in order to avoid the influence of seasonal differences in spectral radiation [19]. In addition, the choice of images with zero cloud in the study area was prioritized so as not to interfere in the area of interest [27].

A Digital Elevation Model (DEM) was acquired with a spatial resolution of 5 m from the National Geographic Institute to generate the topographic variables. Meteorological data were obtained from raster models generated from a set of data of average monthly precipitations and average maximum and minimum monthly temperatures for a period of 60 years (1950–2010). Raster models have been generated by interpolating the weather station dataset recorded in the network of AEMET (Agencia Estatal de Meteorología) onto a regular grid (0.1 × 0.1 degree). The grid was developed following an improved version of the interpolation technique described by Brunetti et al. [38]. This improvement combines a radial weight with a Gaussian shape with an angular weight [39,40]. The vegetation cover map was derived from a combination of two vectorial datasets: the Second National Forest Inventory (INF2) at a scale of 1:50,000 [41] and the Corine Land Cover (CLC) 1990 project at a scale of 1:100,000 [42]. This vegetation cover map only includes forest areas because non-forest areas have been removed due to their lack of interest to our study.

Fire severity was measured using the GeoCBI proposed by De Santis and Chuvieco [28]. It is a method based on the Composite Burn Index [25,43], designed to estimate the damage caused by fire, taking into account five forest strata: soil, herbaceous (<1 m), shrubs (1–5 m), trees (5–20 m), and trees (>20 m). The GeoCBI incorporates two new variables for each plant stratum (percentage of change in the leaf area and vegetation covered area), and considers litter and fuel consumption, changes in soil colour, foliage alteration, canopy mortality, and char height. The severity variable was stratified into five levels to facilitate understanding: 0 = Not burned, >0 to <1 = Low, 1 to <2.5 = Medium, >2.5 to <2.85 = High, and >2.85 to 3 = Very high.

2.3. Methods

The methodology for modelling post-fire regeneration using time series requires a series of processes summarized in Figure 2.

2.3.1. Image Processing

Satellite image acquisition is subject to errors as a result of alterations in satellite motion, sensing data mechanism, and atmosphere interference. A multitemporal study like this requires radiometric homogenization and geometric correction due to the differences in the conditions in which images were acquired [11].

(A) Geometric correction

This process involves transforming the pixels coordinates in ESA image series so that they match precisely with those of the reference image. In this work, we have taken images from U.S. Geological Survey (USGS) as reference, because they have the greatest quality processing. The WGS84 UTM 30N Coordinate System was established for the whole series. The correction was carried out in ArcGIS 10.2.2 by taking 30 control points distributed evenly throughout the image in order to reduce the geometric error [44]. A second-degree polynomial fit and a cubic resampling were applied, which considers the values of the closest 16 pixels and allows for greater accuracy.

(B) Atmospheric correction

In this step, the digital numbers (DN) are transformed to physical magnitudes, which allow the use of images obtained by different sensors [44]. The DN of each band was transformed to reflectance values ($\rho$), following the method based on the dark object proposed by Chávez [45].

$${\rho }_k=\frac{d^2 \times \mathsf{\pi }\times a_{1,k} \left({\mathrm{DN}}_k - {\mathrm{DN}}_{\min ,k}\right)}{E_{0,k} \times sin {\theta }_e \times {\tau }_{k,i}}$$

(1)

where $\rho$ is the reflectance value for each pixel, d is the corrector factor of the distance Earth-Sun, which is calculated from the Julian day, τ_k,i is the average transmissivity (the values used are those proposed by Gilbert et al. [46] for each band: 0.7–0.78–0.85–0.91–0.95–0.97), E_0,k is the solar irradiance on the top of the atmosphere whose most suitable values are currently available in landsat.usgs.gov/esun; the other parameters, a_1,k (multiplicative coefficient of conversion to radiance) and θ_e (sun elevation angle) are obtained from the image metadata.

(C) Topographical shading correction

Topographic correction solves the errors in the reflectance values caused by solar illumination differences due to relief variations. From the resampled DEM to the same images spatial resolution (30 m), the incidence angle (γ) is calculated. The cosine of this angle allows us to measure illumination variations [44]. The topographic shading correction was carried out under the criterion of c-Teillet [47] for each band.

$${\rho }_{h,i}= {\rho }_{h,i} \left(\frac{ cos {\theta }_i }{cos {\gamma }_i }\right)$$

(2)

where ${\rho }_{h,i}$ indicates the reflectance of a pixel i in a horizontal terrain, ${\rho }_i$ the corresponding slope, θ_i the sun azimut angle of the scene, and γ_i the incidence angle calculated corresponding to that pixel.

2.3.2. Use of Spectral Vegetation Index (NDVI)

Several studies have supported the utility of spectral indices derived from Landsat to study the post-fire regeneration dynamics of vegetation [9,27,48]. The spectral indices help identify the state of the vegetation along time series due to the enhancement in the sensitivity of the different spectral bands. One of the most used indices is the Normalized Difference Vegetation Index proposed by Rouse et al. [23], which is calculated by combining the near infrared (NIR) and red (R) bands. The combination of its normalized difference formulation and the use of the higher absorption and reflectance regions of chlorophyll make this index robust in different conditions [9]. In contrast to other vegetation index such as Enhanced Vegetation Index (EVI), the NDVI is easier to calculate and has been widely applied to study vegetation regeneration [14,21,49,50,51].

For this study, it is interesting to model the short-term regeneration of vegetation. For this reason, the NDVI for the complete series (1994–1999) has been calculated to study the regeneration patterns and to generate the variables to be included in the model.

2.3.3. Environmental Variables

The environmental variables which were used to correlate the Pinus species recovery, measured by the NDVI in year 5 after fire, have been defined based on accurate data availability.

A vegetation cover map of burned areas was generated by rasterizing the vegetation cover map in vector format. This map was overlaid with severity maps in order to represent the spatial distribution of Pinus species on burned areas with the same spatial resolution as satellite images. The Pinus species affected by fires were Pinus pinaster and Pinus nigra in the Castrocontrigo and Uncastillo fires, and Pinus halepensis in the Uncastillo and Moratalla fires. Based on the physiological characteristics of each species, as well as the survival and growth capacity of the seeds after fire [33,52], two models have been carried out: one for Pinus halepensis, with a greater regenerative ability [53], and other for Pinus nigra and Pinus pinaster.

Topographic variables (elevation, slope, and aspect) were generated from the DEM mosaic for each fire, resampled to a resolution of 30 m, and projected in the Reference System WGS84 UTM 30N. In order to generate solar orientation variables, the Meng et al. [20] methodology was used. The cosine and sine aspect variables were calculated to provide information on surfaces with northerly and easterly orientations, respectively. Sine and cosine transformations are done so as to convert a nonlinear variable to a linear variable. These variables incorporate information on the time of solar exposure and how it impacts regeneration

With regard to climatic variables, anomalies in precipitation and temperature during the year of the fire (1994) and the successive one (1995), as well as the mean of the anomalies for the 1994–1999 period, were calculated. It is assumed that excess or deficiency in precipitation, as well as abnormally high or low temperatures, may act as limiting factors for vegetation regeneration [20]. Climatic anomalies were calculated for each month at the pixel level. Values were standardized considering the mean and the standard deviation of the 1950–2010 series. Finally, anomalies were obtained by defining thresholds using the first and third quartile of the new Z values.

$$Z = \frac{X - \mu }{\sigma }$$

(3)

where X is the pixel value (precipitation values in mm and temperatures in °C), μ is the mean, and σ is the standard deviation. Before calculating the anomalies, extreme values in the datasets were removed (the lower and upper values of 1 and 99 percentiles, respectively).

2.3.4. Post-Fire Regeneration Assessment

Regression analysis is used to model, predict, or explain complex phenomena. In this case, the influence of fire severity and certain environmental variables on the short-term regeneration rate of selected Pinus species is explained by multiple regression models (

Table 3. List of variables for regression analysis.

Variable		Units	Description
Dependent	NDVI YEAR + 5	Values between −1 and 1	It shows vegetation greenness 5 years after fire.
Explanatory
Post-fire effects	NDVI YEAR + 1	Values between −1 and 1	Vegetation greenness immediately after fire.
Fire severity	CBI	Values between 0 and 3	A CBI value = 0 indicates unburned surface and CBI = 3, an area of very high severity.
Topography	Elevation	Metres
	Slope	Grades
	Northerly	Values between −1 and 1	Obtained by the orientation cosine. Values near 1 indicate North orientation.
	Easterly	Values between −1 and 1	Obtained through the orientation sine. Values near 1 indicate East orientation.
Climatic anomalies	Precipitations	Z value	Total anomalies in the wet season (months that concentrate 80% of annual precipitation: October–May) in 1994 and 1995 and the mean of the anomalies for 1994–1999 period.
	Minimum temperature on January	Z value	Total anomalies in January in 1994 and 1995 and the mean of the anomalies for 1994–1999 period.
	Maximum temperature on July	Z value	Total anomalies in July in 1994 and 1995 and the mean of the anomalies for 1994–1999 period.

NDVI: Normalized Difference Vegetation Index.

).

First of all, a regression analysis was carried out with Ordinary Least Squares (OLS) to diagnose if the selected variables for the models are suitable. However, spatial data show properties that fail to carry out the assumptions and requirements of statistical methods such as OLS regression. This global regression technique assumes the existence of spatial stationarity. It means the model has stable predictive capacity throughout the study area. In this case, the explanatory variables show nonstationary relationships. In addition, spatial autocorrelation in the NDVI values for the 1994–1999 period has been detected by the Moran Index [54]. As an alternative, some authors propose to use regression models that incorporate the variation in geographic space, such as the Geographically Weighted Regression (GWR) [55]. For this reason, this method will be used to estimate post-fire regeneration values.

Pixels used in regression analysis were obtained from a random sample select on the set of burned pixels, using the post-classification functions. Border pixels were previously removed because they were not representative for characterizing the burned surfaces, using a 3 × 3 sum space filter.

$$n={\left(\frac{Z_{\alpha /2}\times \sigma }{E}\right)}^2$$

(4)

where α is the confidence level, in this case 5%, Z_α/2 is the corresponding critical value (1.96), σ is the standard deviation of the dependent variable values, and E is the margin error in estimation. Final samples of 15,814 pixels for Pinus halepensis model and a second sample of 7129 pixels for Pinus nigra and Pinus pinaster model were selected.

3. Results

3.1. Regeneration Evolution According to Fire Severity

First analysis has focused on the evolution of the short-term regeneration of the selected Pinus species according to the severity level. According to fire severity level, Pinus halepensis evolution on the one hand, and other Pinus (nigra and pinaster) on the other, from the mean NDVI values, has been represented (Figure 3). This separation is based on the unequal germination capacity of the species, being higher in Pinus halepensis [53].

Results show a significant decrease in the post-fire values and a later progressive recovery, being faster in the early years in the case of Pinus nigra and pinaster, while in Pinus halepensis, recovery is slower. Fire severity level affects Pinus regeneration significantly, being higher NDVI values associated with lower fire severity. Pre-fire values are different between fire severity levels due to the influence of other environmental variables according to the zone. In addition, areas which were drier before fire may have burned more intensely. Some fluctuations in the short-term regeneration may be caused by variations in the images acquisition dates after summer. This could explain the fall in NDVI values in Pinus nigra and Pinus pinaster between 1997 and 1998.

3.2. Multiple Linear Regression with Ordinary Least Squares

An exploratory regression analysis with OLS allowed for the elimination of some variables initially defined because of multicollinearity problems (spatial autocorrelation) or the low significance for the predictive models.

The average of the climatic anomalies between 1994 and 1999 were dismissed because of their high multicollinearity with post-fire climatic anomalies (1994 and 1995). Moreover, the low spatial resolution of the meteorological data implies a reduced variability of precipitation and temperature data, and there is a high multicollinearity between maximum and minimum temperatures. Therefore, two models were created for each group of species, with one model having maximum temperature, while the other model used minimum temperature. The resulting models were then evaluated and the model with best fit selected. Additionally, the explanatory variables included are the most consistent, that is, those that were significant a greater number of times (>80) when executing the regression with OLS 100 times (

Table 4. Frequency of explanatory variables at the 95% level of significance.

Variable	Pinus halepensis	Other Pinus: P. nigra and P. pinaster
NDVI YEAR + 1	100	100
Fire severity	100	100
Elevation	100	93
Slope	100	82
Northness	2	0
Eastness	0	3
Prec. anomalies +0 and +1	100	100
Tmax anomalies +1	93	85
Tmin anomalies +1	82	83

). In this case, the orientation variables have not been significant and have been removed.

Table 5. Frequency Ordinary Least Squares (OLS) results for Pinus halepensis in the 95% confidence interval.

Variable	Coefficient	Standard Error	t-Statistic	Robust Prob.	VIF
Intercept	0.175275	0.003849	45.536378	0.000000 *	-
Slope	−0.000772	0.000040	−19.111591	0.000000 *	1.326365
Elevation	0.039808	0.002215	17.970745	0.000000 *	2.611230
NDVI+1	0.828052	0.007227	114.571406	0.000000 *	1.489753
Severity	−0.016329	0.001308	−12.486857	0.000000 *	1.173556
Prec. anomalies +0 and +1	0.025462	0.000942	27.018506	0.000000 *	6.863740
Tmax anomalies +0 and +1	−1.104539	0.105490	−10.470555	0.000000 *	6.022914

Multiple R²: 0.702618; Adjusted R²: 0.702505; AICc: −58,550.823000; Joint F-Statistic [e]: 6224.47229; Prob. (>F), (6, 15807) degrees of freedom: 0.000000 *; Koenker (BP) Statistic [f]: 2117.39545; Prob. (>chi-squared), (6) degrees of freedom: 0.000000 *; Jarque-Bera Statistic [g]: 23,935.24457; Prob. (>chi-squared), (2) degrees of freedom: 0.000000 *.

and

Table 6. OLS results for other Pinus (P. nigra and P. pinaster) at 95% confidence interval.

Variable	Coefficient	Standard Error	t-Statistic	Robust Prob.	VIF
Intercept	0.162278	0.007830	20.725448	0.000000 *	-
Slope	−0.000311	−0.000179	−17.37394	0.008237 *	1.373015
Elevation	0.035865	0.005188	6.912623	0.000000 *	1.293531
NDVI+1	0.528326	0.011330	46.630112	0.000000 *	2.385877
Severity	−0.039403	0.003988	−9.881250	0.000000 *	1.234551
Prec. anomalies +0 and +1	0.118088	0.004990	23.665083	0.000000 *	2.248585
Tmax anomalies +0 and +1	0.077906	0.005814	13.399919	0.000000 *	3.996377

Multiple R²: 0.457279; Adjusted R²: 0.456745; AICc: −14,906.559792; Joint F-Statistic [e]: 1000.12498; Prob. (>F), (6, 7122)) degrees of freedom: 0.000000 *; Koenker (BP) Statistic [f]: 285.86905; Prob. (>chi-squared), (6) degrees of freedom: 0.000000 *; Jarque-Bera Statistic [g]: 892.04419; Prob. (>chi-squared), (2) degrees of freedom: 0.000000 *.

contain the first regression models results. The models that have shown the best results have been those in which the minimum temperatures have not been included because they were redundant (high value of the Variance Inflation Factor—VIF).

Multiple R² and Adjusted R² values show a good fit for Pinus halepensis model (Table 5). In order to assess in detail each explanatory variable, the coefficient, t-statistic, robust probability, and VIF were used.

Each explanatory variable coefficient shows the relationship between each explanatory and dependent variable. In this case, slope, severity levels, and Tmax anomalies have a negative influence. High slope slows down rooting vegetation while abnormally high temperatures can cause thermal stress and limit regeneration. As anticipated, short-term regeneration rate was lower when severity was higher. Nevertheless, the severity seems to influence less strongly because sample pixels used coincided with high severity values. Consequently, the influence of a low or medium severity on regeneration is not collected, leaving this variable weak against others. On the other hand, the variable with a greater positive impact is the state of the vegetation after fire (NDVI+1). That is, short-term regeneration rate was higher when post-fire vegetation greenness was higher. There is also a positive influence on elevation and precipitation. In the Mediterranean region, above-average rainfall can help the growth of vegetation by increasing soil moisture. The elevation could be related to lower temperatures avoiding thermal stress in summer.

Standard errors allow assessment of the coefficients obtained. Low standard errors across all variables indicate coefficients are consistent. On the other hand, the robust probability or p-value showed that all variables are statistically significant and, therefore, important for the regression model. VIF provided information on the possible redundancy in explanatory variables. In this case, the relatively high values for climatic anomalies suggest removing one of the two variables could increase the model fit.

As regards to the regression model for Pinus nigra and pinaster (Table 6), the model shows lower Multiple R² and Adjusted R² values due to differences in the vegetative cycle between modelled species. Again, variables with inverse relationships are fire severity and slope, whereas the post-fire vegetation state (greeness) is the variable with more positive influence. Certain above-average rainfall and temperatures higher than usual show a positive coefficient that could be due to compensations in evapotranspiration. Standard errors are reduced and VIF values confirm the absence of variable redundancy in this model.

However, in both models, the significance of Koenker statistic (p-value less than 0.05 for a 95% confidence level) indicated biased standard errors due to heteroscedasticity. In addition, the significant p-value in the Jarque-Bera statistic showed that the residual values deviated slightly from a normal distribution. Finally, an analysis of the residuals was carried out to study spatial autocorrelation using the Moran Index. In both, autocorrelation has been positive, with values around 0.3. Therefore, it was suitable to run the Geographically Weighted Regression (GWR) analysis.

3.3. Multiple Linear Regression with Geographically Weighted Regression (GWR)

Regression with GWR allows for the improvement of the adjustments and for the neutralizing of the spatial dependence in residual values. In comparison with a global regression, the coefficients in GWR are functions of spatial location [56]. Regression analyses were carried out using a defined kernel with a two-square function in which the bandwidth was determined by an optimal number of neighbors. The optimal number of neighbors has been 500, where the value of the Akaike Information Criterion (AIC), which is a relative quality indicator of a model, is the lowest. In final models, anomalies in maximum temperatures were removed because of local multicollinearity problems.

Multiple R² and adjusted R² values obtained in the models using GWR show a significant improvement over the OLS model. In this case, the variables used provide a model explaining 80% of Pinus halepensis regeneration and 78% for Pinus nigra and Pinus pinaster (

Table 7. GWR results at 95% confidence interval.

Variable	Pinus halepensis		Other Pinus (P. nigra and P. pinaster)
	Coefficient	Standard Error	Coefficient	Standard Error
Intercept	0.177311	0.066181	0.256385	0.122476
Slope	−0.000587	0.000408	0.019479	0.189559
Elevation	0.000038	0.000068	−0.258428	0.052855
NDVI+1	0.738801	0.081561	0.503855	0.050379
Severity	−0.007752	0.011347	−0.003933	0.014675
Prec. anomalies +0 and +1	0.031667	0.035612	0.188268	0.022894

Pinus halepensis: Multiple R²: 0.802229; Adjusted R²: 0.801694; AICc: −108.228365; Other Pinus (P. nigra and P. pinaster): Multiple R²: 0.782399; Adjusted R²: 0.774602; AICc: −557.066965.

). It is noticed that explanatory variable forces and influences have changed over OLS models for Pinus nigra and Pinus pinaster. Higher slope shows a slightly positive coefficient, as in the study by Meng et al. [20] for spruce forests. In contrast, the elevation has a negative coefficient that could be related to temperatures too low, which can limit regeneration. However, forces and influences are practically maintained in the Pinus halepensis regression model. Higher slope and higher fire severity have negative influence, while increasing the elevation and the anomalies in precipitations have a positive influence.

However, the regression models do not have an equal predictive capacity for all areas, so there are regional variations (Figure 4). Models show a better fit for Pinus halepensis regeneration in Moratalla and Pinus pinaster in Castrocontrigo. These spatial variations could be related to the sampling point distribution by a smaller number of neighbors, because in Uncastillo sampling points are distributed in a more dispersed way, and by the local conditions effects (microclimate, edaphic composition).

Finally, the relative importance of the GWR explanatory variables were analyzed from the t-statistic (Figure 5). As expected, the post-fire vegetation state (NDVI+1) is the most important variable to the model. In contrast, fire severity was found to have a weak influence on the other variables in both models. Elevation is the variable that displayed the most difference among species, and this was attributed to species adaptation to temperatures decreasing with elevations. This decrease can be positive to avoid abnormally high temperatures in summer in Moratalla (excess evapotranspiration), but can act as a restriction in Castrocontrigo in case of too low temperatures (reduced plant activity).

Moreover, analysis of the residual values also shows better results in GWR than in OLS. The Moran Index values are much closer to the expected values with values below 0.05, also showing a lower variance and greater probability of random distribution (p value and Z score).

3.4. Validation

In order to validate the models, we used the root-mean-square error (RMSE) which measures the average of the squared errors (errors among NDVI values on 1999 and the values predicted by OLS and GWR) (

Table 8. Root-mean-square error (RMSE) from Regression models.

Error	Pinus halepensis		Other Pinus (P. nigra and P. pinaster)
	OLS	GWR	OLS	GWR
RMSE	0.00144253	0.000833933	0.007218565	0.003186855

). Results show that GWR models have smaller error, which implies a better adjustment than OLS models. Taking into account that NDVI values fluctuate among −1 and 1, errors are low. Therefore, predictions obtained using GWR could be considered valid.

4. Discussion

Several studies have researched post-fire vegetation regeneration based on growth patterns according to species, or on the relationship among fire severity and regeneration dynamics. Nevertheless, other previous studies have supported the importance of delving into the analysis of factors which determine regeneration [10,19,20,21]. Therefore, in spite of the difficulties to obtain information to generate certain environmental variables, in this study it has been considered important to model the regeneration considering post-fire climate and topography.

The results show a clear relationship between fire severity levels and regeneration rate values, which are slower when fire severity level is higher [19,27]. Meng et al. [20] and Ireland & Petropoulos [19] agree on the importance that solar radiation levels have on vegetation growth rate. In Mediterranean contexts, north aspect and high slopes present higher regeneration rates due to less evapotranspiration and higher humidity content. However, this phenomenon has not been analyzed because orientation variables have not been significant for the model in the cases studied. With regards to elevation, it has been shown that it could act positively to avoid high temperatures and periods of post-fire drought [53] in more arid environments, such as Moratalla. This is in agreement with Meng et al. [20] for coniferous forests in Sierra Nevada (USA).

With regards to the influence of post-fire weather, a positive relationship between short-term regeneration and above-average precipitation in the months after fire has been found. This relationship was also found for Pinus halepensis and Pinus nigra in Catalonia [21,22] and Pinus halepensis in Ayora in Valencia [8]. Nevertheless, climatic data have a lower spatial resolution compared to Landsat images and topographic data. Consequently, multicollinearity problems have not allowed for a robust analysis of the temperatures or precipitation impact on regeneration. In addition, it would be interesting to include drought indices, which would require greater spatial and temporal resolution data.

Modeling results could also be improved by incorporating new variables that capture local conditions such as landscape configuration, local microclimate, and hydrological processes which can determine vegetation recovery after fire [4,19]. In the study carried out by Chu et al. [10], water soil content showed a positive relationship with the regeneration of the larch forest, being the second conditioning factor in the recovery of this species after fire. Moreover, Röder et al. [8] concluded that water availability, in relation to slope and edaphic composition, is the most important factor that limits vegetation recovery.

In addition, it should be noted that short-term regeneration measured, through the NDVI, represents the relative vegetation cover. A direct measure of Pinus forests regeneration could be obtained also considering the vegetation structure using LiDAR techniques [57].

5. Conclusions

Our study provides advances in the analysis of the impact of fire severity and environmental variables on the short-term vegetation regeneration in Mediterranean regions.

Regeneration measurement for different species of Pinus from the NDVI has been related to the fire severity levels with good results. Severity degree measured by the CBI indicates that short-term regeneration was slowed down when severity was higher. In addition, the immediate NDVI values after the fire allowed us to verify that when the less damage produced in the vegetation, the short-term regeneration is greater.

From the multiple linear regression models generated, the explanatory capacity of the environmental variables of topography (elevation and slope) and post-fire climate (anomalies in precipitation) in post-fire vegetation recovery after fire has been verified/tested. The impact that each has on regeneration is closely related to the environmental characteristics to which each species is adapted. Thus, the elevation can be a driving factor for the Pinus nigra regeneration linked to temperatures too low or, in contrast, can favor Pinus halepensis growth, by avoiding high temperatures in summer in the southeast of the Iberian Peninsula.

In contrast to OLS, it has been confirmed that GWR is an important and valid local regression technique to explore spatial heterogeneity in the relations of explanatory variables. With this methodology, it has been possible to model the regeneration obtaining high adjustment values, with adjusted R² values of 0.80 for Pinus halepensis and 0.77 for Pinus nigra and Pinus pinaster. The models’ improvement should focus on the generation of more precise environmental variables and considering new factors that can increase it explanatory power: lithologic characteristics, alteration of edaphic composition, solar radiation, time elapsed since the last fire, vegetation physical characteristics, etc.

Moreover, the proposed method is an approximation to the modelling of the short-term regeneration for the Iberian Peninsula but exportable to other territories with input variables transformation. The results obtained are useful in improving knowledge about the factors which determine the post-fire regeneration patterns of a forest ecosystem under different environmental and climatic conditions. Therefore, these advances could help decision-makers in determining which areas vegetation will not regenerate naturally after large fires and thus requires the implementation of specific restoration programs.