Modelling Maritime Pine ( Pinus pinaster Aiton) Spatial Distribution and Productivity in Portugal: Tools for Forest Management

: Research Highlights: Modelling species’ distribution and productivity is key to support integrated landscape planning, species’ afforestation, and sustainable forest management. Back-ground and Objectives: Maritime pine ( Pinus pinaster Aiton) forests in Portugal were lately affected by wildﬁres and measures to overcome this situation are needed. The aims of this study were: (1) to model species’ spatial distribution and productivity using a machine learning (ML) regression approach to produce current species’ distribution and productivity maps; (2) to model the species’ spatial productivity using a stochastic sequential simulation approach to produce the species’ current productivity map; (3) to produce the species’ potential distribution map, by using a ML classiﬁcation approach to deﬁne species’ ecological envelope thresholds; and (4) to identify present and future key factors for the species’ afforestation and management. Materials and Methods: Spatial land cover/land use data, inventory, and environmental data (climate, topography, and soil) were used in a coupled ML regression and stochastic sequential simulation approaches to model species’ current and potential distributions and productivity. Results: Maritime pine spatial distribution modelling by the ML approach provided 69% ﬁtting efﬁciency, while species productivity modelling achieved only 43%. The species’ potential area covered 60% of the country’s area, where 78% of the species’ forest inventory plots (1995) were found. The change in the Maritime pine stands’ age structure observed in the last decades is causing the species’ recovery by natural regeneration to be at risk. Conclusions: The maps produced allow for best site identiﬁcation for species afforestation, wood production regulation support, landscape planning considering species’ diversity, and ﬁre hazard mitigation. These maps were obtained by modelling using environmental covariates, such as climate attributes, so their projection in future climate change scenarios can be performed.


Introduction
Nowadays, access to big spatial data on climate and other environmental variables has fostered the use of powerful techniques from artificial intelligence and spatial statistics, such as machine learning (ML) and geostatistical modeling, which, coupled with geographic information systems (GIS) allow for the construction of simulated maps for the species' habitat suitability and productivity under the impacts of climate change [1][2][3][4][5][6][7][8][9][10][11].
Indeed, statistical modelling techniques such as classical regression (CR), generalized linear models (GLM), algorithmic modelling based on machine learning (ML), e.g., Bayesian networks (BNs), maximum entropy (MaxENT), and classification and regression trees (CART) have become increasingly popular [12]. For instance, some of these statistical modelling techniques were used to predict current and future suitable habitats and productivity for Maritime pine populations (Pinus pinaster Aiton) in Spain [1], the impact of climate change on the potential distribution of Mediterranean pines [2], and the current and future conflicts between Eucalypt plantations and high biodiversity areas in the Iberian Peninsula [3]. In Portugal, a bioclimatic modelling approach (MaxENT) was used to study the influence of environmental variables explaining the presence of strawberry tree (Arbutus unedo L.), a typically adapted species to the Mediterranean region, under climate change scenarios [4].
Regarding forest species' productive potential mapping, several methodological approaches have been essayed in various studies [6][7][8][9][10][11]. In Portugal, forest species' productive potential maps for present and future climate change scenarios are available, for each one of the seven forest management regions the country, as divided in [13]. Nevertheless, each of these regional maps were produced with different methodological approaches (e.g., using ecological-cultural characteristics and/or edaphic-climatic characteristics, or bioclimatic indices or productivity estimation), thus no consistency exists when considering the whole country.
Yet, in a recent study [11], current species' potential maps were produced for Portugal using the species ecological envelopes methodology, proposed by DGRF (Direção Geral dos Recursos Florestais) [14], and afterwards, projected into the future under climate change scenarios. Moreover, in that study a Bayesian ML approach was used to explore the most influential environmental variables on current species' distribution and productivity for the species essayed (Maritime pine and Eucalypt). The findings in that study indicated that the species' ecological envelope methodology was robust and encouraged further investigation on fitting statistical ML models for the species' distribution and productivity for Maritime pine in Portugal.
Maritime pine forests provide a significant amount of the wood harvested in Portugal [15] as raw material for the wood-based industry. According to the last National Forest Inventory (NFI 2015) this species' forests area represents 23% of Portuguese forests (713 × 10 3 ha). However, the impact of climate change in Maritime pine forest area must not be neglected. Indeed, several studies have consistently reported a future decrease of potential areas of the Maritime pine, particularly in south of the Tagus river and the inland of Portugal, associated with a species' migration trend from south to north and from the interior to coastal areas [11].
Maritime pine forests are mostly established by natural regeneration and managed in medium to long revolutions (e.g., more than 35 years) [15][16][17]. Maritime pine is a light pioneer species and viable seed production occurs later in stands of 15 to 20 years that ensures desirable natural regeneration [18].
In Portugal wildfires are currently a major issue and a serious challenge for forest management. In fact, during the 1950s and 1960s, the country observed a total burnt area of about 5000 ha year −1 . Between 1975 and 1980, the total burnt area increased remarkably, reaching about 43 [19,20].
Undeniably, wildfires' severity and recurrency, e.g., within a 10 to 15 year interval, have had a huge impact in Maritime pine forests, making the species' post-fire recovery uncertain [21]. For instance, during the past decades Maritime pine's post-fire area has been mainly converted into scrub, pastures, and Eucalypt plantations [19]. Moreover, forest owner's confidence levels in sustaining Maritime pine management significantly dropped afterwards. In fact, a study in the central inland of Portugal noticed that these stands became less stable (inventory data 1991-1996 and 2007-2010) and under-stocked, with large amounts of small-diameter poles, and enlarged tree size variability, showing a decline in their management [22].
More recently, the industrial development of bioenergy production from biomass residues created another industrial destination for forests' raw material which traditionally had no commercial value [23][24][25]. The opportunity to value forest residues resulting from silvicultural operations (e.g., brushing, cleaning, pruning, thinning, and harvesting) are an inducement to forest management and fire severity mitigation.
It is well known that fire hazard mitigation can be attained by a combination of different associated effects, such as forest type-appropriate management, together with stands' age planification, density, and spatial landscape pattern arrangements (i.e., patch size, forest fragmentation) [26][27][28][29]. Indeed, what contributes most to fire severity are the species' flammability and combustibility, spatial arrangement, vertical structure (high/low and close/open), and land topographic characteristics (slope and aspect) [27,28].
To that end, the definition of a forest type implies considering, among others, straightforward afforestation planning at the landscape level. The stand's vertical structure can be managed through corresponding rearrangement by different ages, heights, or eligible density. The area control method by Davis and Johnson [30] allows for forest production regulation, introducing the age diversity factor to the landscape. Therefore, this method works as a tool on fire hazards mitigation by creating a stand age class progression (e.g., different succession stages) with different fuel loads, as well as vertical fuel continuities between harvesting compartments. Over time, the stands' age class progression (e.g., different succession stages), due to their establishment and growth, results in a vertical structural change, from short and open (in the establishment phase), to short and close, and, eventually, tall and close (in the harvest phase). Finally, when the compartment is harvested (e.g., clear-cutting), a break in the fuel load is observed [26]. It is important to stress that application of the area control method requires previous knowledge of the area under productivity regulation [9,26]. Thus, the availability of cartographic support on species' current and potential distributions and productivity is key as a supporting tool for integrated landscape planning, species' afforestation, and sustainable forest management [26].
The present study focused on producing Maritime pine distribution and productivity maps by using a ML modelling approach. In a previous study [11], the authors used a predefined set of environmental variables to produce an ecological envelope map for the species' potential distribution. Afterwards, a Bayesian ML approach highlighted that the species' current distribution was mainly determined by precipitation-related variables, and that elevation-and temperature-related variables were very important to differentiate the species' current productivity. However, these last variables were not the ones used to obtain the species' ecological envelope map.
Therefore, in the current study, the research question was to verify if the set of environmental variables used in species' ecological envelope map were also good regressors in the species' distribution model (SDM) and species' productivity model (SPM). Additionally, species' spatial productivity modelling was also performed using an alternative methodological approach (sequential stochastic simulation). Hence, a full discussion about high and low Maritime pine's productivity spots was made. Furthermore, the fitted species' distribution model (SDM) and species' productivity model (SPM) may be used to make projections of the future under climate change scenarios.
The specific aims of this study were: (1) firstly, to model the species' spatial distribution and correspondent productivity using a ML Regression approach; (2) in a second step, to model the species' spatial productivity through an alternative approach based on a sequential stochastic simulation process; (3) to produce a species' potential distribution cartography and the definition of the ecological envelope's thresholds by using a ML classification technique; and, finally (4) to identify present and future key factors of the species' afforestation and optimized management. To that end, land cover and forest inventory, used as hard data, together with soft data attributes such as climate variables, topography, and soil information, were used to model species' current and potential spatial distribution and productivity.

Study Area
The study area is Mainland Portugal (36.9636 • N 9.4944 • W and 42.1543 • N 6.1892 • W), located in the Iberian Peninsula (southwestern Europe), between Spain, to the east, and the Atlantic Ocean to the west. The annual temperature ranges from about 7 • C in the northern and central inland highlands to about 18 • C in the south coast. The average annual precipitation has the highest values in the north and northwest regions and the lowest values in the south and inland. The study area includes two climate classifications: (i) hot-summer Mediterranean in the south and central inland and lower elevations (below 800 m); and (ii) temperate Mediterranean in the north and northwest region of the country and higher elevations (above 800 m) [31].

Hard Data-Land Cover and Forest Inventory
The Portuguese Land Cover and Land Use (LCLU) 1995 thematic map (COS1995) [32] was used to extract the Maritime pine's forest cover layer. The COS1995 map has 1:25,000 scale, corresponding to the minimum cartographic unit of one ha. Furthermore, a five-level hierarchical classification with 89 classes segmentation at the most detailed level allowed a detailed analysis regarding this species' forest cover ( Figure 1). The COS1995 was available through the Web Feature Services (WFS) [33]. potential spatial distribution and productivity.

Study Area
The study area is Mainland Portugal (36.9636° N 9.4944° W and 42.1543° N 6.18 located in the Iberian Peninsula (southwestern Europe), between Spain, to the e the Atlantic Ocean to the west. The annual temperature ranges from about 7 °C northern and central inland highlands to about 18 °C in the south coast. The annual precipitation has the highest values in the north and northwest regions lowest values in the south and inland. The study area includes two climate classifi (i) hot-summer Mediterranean in the south and central inland and lower ele (below 800 m); and (ii) temperate Mediterranean in the north and northwest regio country and higher elevations (above 800 m) [31].

Hard Data-Land Cover and Forest Inventory
The Portuguese Land Cover and Land Use (LCLU) 1995 thematic map (CO [32] was used to extract the Maritime pine's forest cover layer. The COS1995 m 1:25,000 scale, corresponding to the minimum cartographic unit of one ha. Furth a five-level hierarchical classification with 89 classes segmentation at the most d level allowed a detailed analysis regarding this species' forest cover ( Figure  COS1995 was available through the Web Feature Services (WFS) [33] The National Forest Inventory 1995 was used to assess pure Maritime pine productivity, defined as forests where the dominant species represents more th and the ground cover is higher than 10% [17,33,34] (Figure 1).  The National Forest Inventory 1995 was used to assess pure Maritime pine stands' productivity, defined as forests where the dominant species represents more than 75% and the ground cover is higher than 10% [17,33,34] (Figure 1).

Soft Data-Climate, Topography, and Soil Attributtes
The climate data sets (Table 1) for the period between 1960 and 1990 was downloaded from the WorldClim site version 1.4, with a 1 km × 1 km grid [35]. The elevation data (Table 1) was derived from the Shuttle Radar Topography Mission 1 Arc-Second Global (SRTM1N22W016V3) with a spatial resolution of 30 m [36]. The soil data (Table 1) was derived from the ESDBv2 Raster Library in a 1 km × 1 km grid [37][38][39]. Aspect-The orientation of slope measured clockwise in degrees from 0 to 360, where 0 is north-facing, 90 is east-facing, 180 is south-facing, and 270 is west-facing.
Soil WRBFU Soil codes from the international soil classification system for naming soils and creating legends for soil maps.
Climate variables (Table 1) followed the original WorldClim symbology, wherein the first variables from Tmax to BIO11 are temperature-related variables and variables from BIO12 to BIO19 are precipitation-related variables. The topographic-related variables considered were elevation, slope, and aspect. The soil-related variables used the soil codes according to the international classification system (Appendix A, Table A1-soil codes legend for the soil data in this study).

Species' Distribution and Productivity Modeling-A Machine Learning Regression Approach
The species' presence was extracted from COS1995 classes, in 1 km × 1 km grid, as a coverage percentage in each pixel (Pb95; species' presence occurs when pixel percentage coverage is higher than 0). The same procedure was used for the correspondent bioclimatic attributes , topography and soil variables (Table 1).
To model the pure Maritime pine stands' productivity, the site index was used, since it is the most widely accepted site productivity indicator [40]. The classical site index SI 35 (i.e., the stand dominant height at the age of 35 years) fitted by Páscoa et al. [41] was evaluated for each inventory plot. It should be stressed that this site index model was fitted for the species and country using forest inventory data from 1995, as well. The extracted climate attributes, elevation, and soil data sets were projected in each plot location (Table 1).
WEKA (Waikato Environment for Knowledge Analysis) open-source software [42,43] was subsequently used for modelling purposes. WEKA allows for big data pre-processing, and the supported algorithms are organized under Classify, Cluster, Associate, and Select Attributes (SA). The SA was used as the regressor for covariates selection in the subsequent ML analysis.
The search method "Rankers" under the "Correlation Attribute Evaluation" provides a selection ranking technique for the predictors' inclusion in the final regression model, by measuring the correlation (Pearson) between the predictor and the class. Thevalues were expressed in relative terms for comparability purposes, as follows: Afterwards, the seven best fitted variables were retained and used to model the species' distribution (Pb95) and correspondent productivity (SI 35 ). Additionally, a set of environmental variables were used for species' ecological envelope delimitation, accordingly: BIO5, BIO6, BIO7, BIO12, E, WRBFU (Table 1).
ML models were fitted, under Classify, for the species' distribution (Pb95; species' presence occurs when pixel percentage coverage is higher than 0) and productivity (SI 35 ). The variables Pb95 and SI 35 are continuous variables, thus three ML regression trees algorithms were tested, namely: M5P, Random Forests (RF), and Random Tree (RT) [44,45]. WEKA software provides k-fold cross-validation to test the algorithm's accuracy, by dividing the data set into k subsets. One of these subsets is randomly selected as the test set, and the other k-1 subset constitutes the training sample, in each run. The models' fitting performance is assessed using statistics such as the pseudo-coefficient of determination (R 2 ) [46], the mean absolute error (MAE) and the root mean squared error (RMSE). Under Classify, the species distribution (P; species presence and absence, that is P = 1 when Pb95 > 0 and P = 0 otherwise) was fitted. As the variable P is a discrete (nominal) variable, then three ML classification tree algorithms were tested, namely: J48, random forests (RF), and random tree (RT) [44,45]. Statistics such as Cohen's kappa (kappa), MAE and RMSE statistics were computed for the models' fitting performance assessments.

Species' Productivity Modelling-A Sequential Gaussian Simulation Approach
The sequential Gaussian simulation (SGS), a stochastic sequential simulation algorithm, was used. This geostatistical approach starts with the definition of the univariate distribution normal scored values. Indeed, the first step is the normal score transformation of the measured values [47]. Furthermore, normal scores at 2 km × 2 km grid node locations were sequentially simulated using simple kriging (SK) with zero mean [48]. The SI 35 spatial structure was performed through variography. Finally, the obtained outputs went through a back-transform function, allowing the results to be obtained in the original grade values. The Space-Stat software v.4.0.18 from Biomedware was used for computation purposes [48,49]. It is important to highlight that when multiple sequences of simulations are designed (100 different scenarios were performed), it is possible to obtain more accurate probabilistic maps. The mean image (MI) together with a representation of the spatial uncertainty (the 100 simulations standard deviation was used as the uncertainty measurement), allows the spatial productivity clusters robustness to be ascertained, and therefore constitutes a powerful tool for productivity evaluation and future management policies [11,50].

Species' Ecological Envelope-Potential Distribution
The species' potential distribution map was based on the methodology introduced by the DGRF (Direção Geral dos Recursos Florestais) [14] concerning the definition of the species' ecological envelopes. Maritime pine's ecological envelope is based on literature information [51] and was defined using the following five variables: (i) temperature range (T max Aug-T min Jan <26 • C); (ii) temperature limits (T max Aug <29.9 • C); (iii) precipitation (BIO12 > 850 mm); (iv) elevation (E < 800 m); and (v) lithology (different of limestone).
In this study, the variables' threshold aimed to ensure that at least 75% of plots (1995), as the dominant species [33], are inside of the defined envelope. The soft covariates used for the species' ecological envelope definition were selected through the ML regression tree algorithm J58. Additionally, to safeguard consistency in data usage in all produced maps, the soil data was used to set the lithology constraint as follows: WRFBU-soils differing by luvisols calcic, cambisols calcarean, and fluvisols calcic.
Afterward, the covariates' thresholds were obtained fitting a ML classification tree model. The maps of the selected variables were constructed using the new thresholds and a subsequent Boolean reclassification (0-not suitable, 1-suitable). A map algebra approach was used for the species' potential distribution map computation (i.e., the ecological envelope) as follows: . Finally, the final map was obtained, reclassifying into four suitability classes: 3-excellent (without temperature range, temperature limits, and precipitation restrictions), 2-good (with one restriction), 1-regular (with two restrictions), and 0-unsuitable (null).
The obtained map for the species' potential distribution was then compared with the empirical maps for forest regions in Portugal by Alves et al. [52] and the Maritime pine occupation map by Oliveira [17] (Figure 2). The coverage of the Maritime pine plots in 1995 ( Figure 1) was also assessed using its suitability class. of simulations are designed (100 different scenarios were performed), it is possible to obtain more accurate probabilistic maps. The mean image (MI) together with a representation of the spatial uncertainty (the 100 simulations standard deviation was used as the uncertainty measurement), allows the spatial productivity clusters robustness to be ascertained, and therefore constitutes a powerful tool for productivity evaluation and future management policies [11,50].

Species' Ecological Envelope-Potential Distribution
The species' potential distribution map was based on the methodology introduced by the DGRF (Direção Geral dos Recursos Florestais) [14] concerning the definition of the species' ecological envelopes. Maritime pine's ecological envelope is based on literature information [51] and was defined using the following five variables: (i) temperature range (T max Aug-T min Jan <26 °C); (ii) temperature limits (T max Aug <29.9 °C); (iii) precipitation (BIO12 > 850 mm); (iv) elevation (E < 800 m); and (v) lithology (different of limestone).
In this study, the variables' threshold aimed to ensure that at least 75% of plots (1995), as the dominant species [33], are inside of the defined envelope. The soft covariates used for the species' ecological envelope definition were selected through the ML regression tree algorithm J58. Additionally, to safeguard consistency in data usage in all produced maps, the soil data was used to set the lithology constraint as follows: WRFBU-soils differing by luvisols calcic, cambisols calcarean, and fluvisols calcic.
Afterward, the covariates' thresholds were obtained fitting a ML classification tree model. The maps of the selected variables were constructed using the new thresholds and a subsequent Boolean reclassification (0-not suitable, 1-suitable). A map algebra approach was used for the species' potential distribution map computation (i.e., the ecological envelope) as follows: Finally, the final map was obtained, reclassifying into four suitability classes: 3-excellent (without temperature range, temperature limits, and precipitation restrictions), 2-good (with one restriction), 1-regular (with two restrictions), and 0-unsuitable (null).
The obtained map for the species' potential distribution was then compared with the empirical maps for forest regions in Portugal by Alves et al. [52] and the Maritime pine occupation map by Oliveira [17] (Figure 2). The coverage of the Maritime pine plots in 1995 (Figure 1) was also assessed using its suitability class.
The maps produced in the study (current species distribution, species productivity, and species potential distribution) were evaluated together, aiming for an evaluation of the best sites for species settlement.

Species' Distribution and Productivity Modelling-A Machine Learning Regression Approach
Two types of models (ME and MR) were used to predict species grid percentage coverage (Pb95), species presence (P), and site index (SI 35 ) ( Table 3). The ME models used as regressor variables were the same variables as the ecological envelope. The MR models used as regressor variables were the seven best variables selected by the search method "Rankers" according to their relative importance (RI; Table 2). It should be noted that all variables are continuous, with exception made to the variables P and WRBFU which are discrete (nominal variables). The relative importance (RI) of regressor variables for the Maritime pine distribution (Pb95) in the MR model was BIO13 > BIO19 > BIO16 > BIO12 > BIO18 > BIO17 > BIO2, and for Maritime pine presence (P) in the MR model, The best performance was achieved by the "random forest" algorithm. The selected model for the species' distribution (Pb95-species grid percentage coverage and P-presence/absence) as well as for the species' productivity (SI 35 ) was the ME model fitted by the algorithm RF (random forest), a model using conceptual environmental variables as regressors. The ME models' fitting efficiency for Maritime pine distribution (Pb95 and P), were 69% and 65%, respectively. The ME model for Maritime pine productivity (SI 35 ) explained only 43% of the observed variability ( Table 3). The ecological envelopes' variables proved to be robust in defining the species' presence and corresponding productivity with potential areas. The observed versus residuals scatterplot shows a slightly higher performance than the ME models ( Figure 3).

Species' Productivity Modelling-A Sequential Gaussian Simulation Approach
The previous structural evaluation-variography-revealed no evidence of spatial anisotropy and, therefore, the omnidirectional variogram was used for the subsequent stochastic simulation technique ( Figure 5). The species' spatial productivity map produced by the SGS approach highlights the highest productivities of timber (for fully-mature 35 year old trees) in the coastal center regions followed by the northern and central montane and submontane regions ( Figure 6).

Species' Ecological Envelope-Potential Distribution
The summary data for the species' presence (Pb95 > 0) provides a few highlights about the ecological envelope variables range (Table 4) and the new computed thresholds through the ML classification tree algorithm, J48 (Table 5). A comparative assessment of these values allows one to conclude that the new defined thresholds are inside the range of the species' presence and, therefore, consistent with the literature values (Table 5; Figure 7a,b). As expected, a broader range is observed when considering the species' presence (Table 4) than when considering the species' ecological envelope range, which sets the species' potential distribution area (Figure 7c). As a result, the species' potential distribution area was assessed in 60% of the country's area and 78% of the species' forest inventory plots (1995), as a dominant species, were found inside this ecological envelope (Figure 7c,d). The species' spatial productivity map produced by the SGS approach highlights the highest productivities of timber (for fully-mature 35 year old trees) in the coastal center regions followed by the northern and central montane and submontane regions ( Figure 6). The species' spatial productivity map produced by the SGS approach highlig highest productivities of timber (for fully-mature 35 year old trees) in the coastal regions followed by the northern and central montane and submontane regions (Fig   (a) (b)

Species' Ecological Envelope-Potential Distribution
The summary data for the species' presence (Pb95 > 0) provides a few highlights the ecological envelope variables range (Table 4) and the new computed thresholds th the ML classification tree algorithm, J48 (Table 5). A comparative assessment of these allows one to conclude that the new defined thresholds are inside the range of the s presence and, therefore, consistent with the literature values (

Species' Ecological Envelope-Potential Distribution
The summary data for the species' presence (Pb95 > 0) provides a few highlights about the ecological envelope variables range (Table 4) and the new computed thresholds through the ML classification tree algorithm, J48 (Table 5). A comparative assessment of these values allows one to conclude that the new defined thresholds are inside the range of the species' presence and, therefore, consistent with the literature values (Table 5; Figure 7a,b). As expected, a broader range is observed when considering the species' presence (Table 4) than when considering the species' ecological envelope range, which sets the species' potential distribution area (Figure 7c). As a result, the species' potential distribution area was assessed in 60% of the country's area and 78% of the species' forest inventory plots (1995), as a dominant species, were found inside this ecological envelope (Figure 7c,d).     Soils different of LVcc, CMca, and FLca Legend: Temperature-related variables BIO7-Temperature annual range (°C); BIO6-Minimum temperature of the coldest month (°C); BIO5-Maximum temperature of the warmest month (°C); Precipitation-related variables BIO12-Annual precipitation (mm); Topographic-related variables E-Elevation (m); and Soil related variables WRFBU-soils different of luvisols calcic (LVcc), cambisols calcarean (CMca), and fluvisols calcic (FLca).

Species' Afforestation and Management
The National Forest Inventory statistics  show the noticeable decrease in Maritime pine areas over the last 50 years (Figure 8a). Besides, as the species' area decreased, a long-term damaging effect is observed in the stand age structure distribution of the past 20 years (1995-2005-2015) (Figure 8b-d). In 1995 (Figure 8b), Maritime pine forest had a regular age structure distribution. But in 2005 (Figure 8c), an important decrease in mature and old stands (>40 years) is already observed. The large wildfire which occurred in 2003 had a significant impact on the new stands' increase by natural regeneration (<10 years). This trend was worse in 2015 (Figure 8d), with a massive increase in stands established by natural regeneration (<10 years) in opposition to the mature and old stands' disappearance (>40 years). Again, 14 years after the 2003 wildfire, another large and extremely severe wildfire occurred in 2017. As a result, the young stands' burnt area (<10 years; without viable seed production) will hardly recover by natural regeneration, so conversion to scrubland or other species is expected. From 1995 to 2015, the increase of uneven-aged stands was also observed, due to the absence of management (human desertification, forest owners ageing and lack of forest investment due to the high fire hazard) (Figure 8b,d). large and extremely severe wildfire occurred in 2017. As a result, the young stands' burnt area (<10 years; without viable seed production) will hardly recover by natural regeneration, so conversion to scrubland or other species is expected. From 1995 to 2015, the increase of uneven-aged stands was also observed, due to the absence of management (human desertification, forest owners ageing and lack of forest investment due to the high fire hazard) (Figure 8b,d).

Discussion
The results obtained in this study revealed that: (1) Maritime pine spatial distribution (Pb95 and P) modelling by ML approach provided a 69% and 65% fitting efficiency, respectively, while only a 43% fitting efficiency was obtained with the species' productivity (SI 35 ) modelling; (2) the species' spatial productivity map corresponding to the SGS mean image (MI) shows the spatial distribution for high and low Maritime pine productivity spots in concordance with the results obtained when using the ML approach. Together, the approaches of the evaluation stresses two important points: (i) the species' productivity is optimized in the northern coast (Atlantic influence) followed by the northern and central montane and submontane regions, and (ii) the obtained results overlap, given the two used techniques, validating the distribution of the associated spatial uncertainty and signaling a less robust assessment of the Maritime pine's productivity in the south of the country; (3) the new obtained thresholds, through ML, were consistent with the literature [17,52], and refined the species' potential distribution map. The species' potential area was assessed in 60% of country's area with 78% of the species' forest inventory plots (1995) found inside this ecological envelope; and (4) the change in the Maritime pine stands' age structure observed in the last decades put the species' recovery by natural regeneration at risk.
The species spatial productivity map produced by the SGS approach uses the site index SI 35 as a forest productivity indicator. This site index is highly linearly correlated to the site index SI 50 (fitting efficiency of 90%) used in the PBRAVO growth and yield model for the species in Portugal. So, correspondence between the site index SI 35 classes and timber productivities (e.g., mean annual increment at the age of 45 years) is possible. Low productivities (SI 35 < 15 m) correspond to around 5 m 3 ha year −1 or lower and high productivities (SI 35 > 18 m) correspond to around 8 m 3 ha year −1 or higher, which is in accordance with other studies and authors [56][57][58].
The species' potential distribution map (ecological envelope) identifies the coastal strip in the north and center penetrating inland up to elevations of 700-900 m (e.g., smSA) as the Maritime pine distribution in Portugal. However, elevations from 0-400 m correspond to the most favorable sites (e.g., bSA, bMA, bAM, bsM). Current knowledge calls attention to the fact that the species does not well tolerate intense and prolonged colds and snow [17,52]. The species has severe limitations to growth in elevations above 800 m due to wind and snow, showing increased susceptibility to pests and diseases [17,52]. It was found that the species occurs in a variety of soil types, but mostly in cambisols molic, podzols haplic and regosols dystric [39]. Indeed, the species grows well in poor soil with light textures, preferably siliceous (e.g., podzols) [17,39].
Comparing the current Maritime pine distribution with the species' potential distribution map, only 39% of species' area stands in unsuitable sites compared to 51% in excellent and good sites. Regarding the productivity map, most of the species' current distribution is in intermediate productive sites (56%), with high productive sites representing only 17%. When comparing the potential and productivity maps, a narrow concordance of only 28% between potential and productivity classes was obtained. Since forest productivity depends both on the environmental factors inherent to the site (climate, topography, and soil) and on management-related factors, it is possible that non-concordant classes (potentiality vs. productivity) can be explained to some extent by site micro-environmental conditions and/or species' management (active/absent). Furthermore, the species' range expansion in Portugal during the twentieth century (mainly due to human activity), and the extensive gene flow among populations related to this expansion, indicates that there is no genetic structure for the Portuguese populations of this species, meaning that the species is genetically homogeneous in the country [59].
Species management can be supported by applying the intensive management stand prescription, round wood production-oriented, as defended by Louro et al. [60], that uses artificial regeneration by plantation. In intense natural regeneration areas, a fully stocked stand prescription, pulp wood production-oriented, argued by Oliveira [17], is well suited. Moreover, the later prescription saves site preparation and plantation costs, being the most appropriate for the existing naturally regenerated Maritime pine stands of Portuguese private forest areas. It should be stressed that the two above recommendations were considered the most economically efficient in a study conducted for the Portuguese central inland species [58].
Forest planning and management at the local level is argued to be mandatory [26]. Firstly, the Forest Management Unit (FMU) must be organized into administrative units (units). Secondly, the Maritime pine forest must be identified by using the latest official land cover map and updated with the latest annual burnt area maps. Thirdly, the method for area control can be applied, following the species' productivity map in defining the management compartments' number and size, necessary to wood production regulation (annual constant production). The number of compartments is defined by the harvesting age, and the size of each compartment considering the site class productivity aiming to ensure a constant wood production in each year. Fourth, the silvicultural schedule is applied to each compartment, supported by the species' potential distribution map, and will allow after harvesting the decision of whether to regenerate or to convert to another species or land use [9,26]. If the decision is to regenerate the Maritime pine stand, then measures to ensure natural regeneration sustainability should be considered.
In sum, the maps for the species' current and potential distribution and productivity produced in this study allows us to identify the best sites for Maritime pine establishment, support this species' wood production regulation and landscape planning, to decide about species' diversity introduction (other products and services), and for fire hazard mitigation (fragmentation/compartmentalization). Furthermore, these maps were produced for the Maritime pine considering the entire country and with a higher resolution than others previously produced [7,13,61], as more detailed input data was used. Besides, the use of a ML approach to model species' distribution model (SDM) and species' productivity model (SPM) allowed us to validate the regressors' variables used as covariates and to assess the resulting maps accuracy and precision, nonexistent information in previous studies. Moreover, it is worth noting that the seven-best selected regressor variables, either for the SDM and SPM, obtained in this study were consistent with the findings of a previous study [11]. Indeed, emphasizing that precipitation-related variables are more important in explaining the species' distribution while a combination of temperature-and precipitationrelated variables are important to explain species' presence. On the other hand, elevationand temperature-related variables are more important to explain the species' productivity.
Additionally, in a similar study for Maritime pine in Spain, the best models selected for the species current distribution and productivity were the ML classification tree algorithm random forest and the ML regression tree algorithm by random forest, respectively [1]. These authors obtained an overall accuracy of 73% for the SDM and 60% for SPM, which is slightly higher than in the present study. This may be explained as in that study, more detailed data regarding the lithostratigraphy, geology and soil physical properties were used (e.g., bulk density of fine earth fraction (<2 mm); percentage of silt in soil; probability of occurrence of R horizon within 200 cm; depth to bedrock (R horizon) up to 200 cm; cation exchange capacity) [1].
Finally, further investigation is already being undertaken to compare maps' agreement, for the present and for future projections under climate change scenarios, using the species' distribution model (SDM) and species' productivity model (SPM) fitted in this study and using the ML maximum entropy modelling by the MaxENT software.

Conclusions
Modelling species' distribution and productivity are key to supporting integrated landscape planning, species afforestation, and sustainable forest management. Maritime pine forest area has decreased remarkably in the past years and requires urgent measures optimize recovery. The potential distribution map (ecological envelope) and the productivity map, obtained through a coupled stochastic-deterministic approach, are key for decision making on identifying the best sites for species afforestation, to apply the best management practices, and to improve species productivity. Post-fire areas in the best sites for the species should also be evaluated regarding the species' natural regeneration capacity.
Since the species' current and potential distributions and associated productivity were modeled using environmental covariates, such as climate attributes, their projection in future climate change scenarios is an achievable goal. This prospective view supports the design of future integrated planning, regarding fire hazard mitigation, and promoting forests multifunctionality and ecosystems services, under climate change scenarios.  Data Availability Statement: All the data used in this study will be available upon request to the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest. Table 2. Relative importance analysis to obtain a ranking regarding the significance of each regressor variable to the species distribution (Pb95 and P) and productivity (SI 35 ) modelling in ML tree classification and regression.