The Extent of Infrastructure Causing Fragmentation in the Hydrocarbon Basin in the Arid and Semi-Arid Zones of Patagonia (Argentina)

: Fragmentation is a disruption in the connectivity of landscapes. The aims of this paper are (i) to quantitatively assess the fragmentation rates in three landscape units located in a hydrocarbon basin, and (ii) to model their behavior between 2001 and 2013 using landscape metrics at di ﬀ erent scales of resolution. The following metrics were selected using principal component analysis (PCA): The Clumpiness Index (CLUMPY), patch density (PD), perimeter-area fractal dimension (PAFRAC) and e ﬀ ective mesh size (MESH). Results from our investigations pointed out that hydrocarbon activity increased the fragmentation at the sites. In particular, the CLUMPY index increased in all three landscape units, the average of PD decreased from 60 to 14 patches per 100 hectares, whereas the mean of MESH was quite constant, however, due to oil production, it decreased mainly in the coastal valleys. Finally, the PAFRAC also decreased at sites with oil production, being more evident in the plateau and coastal canyons. As a whole, outputs from our analyses clearly pointed out that the monitoring of landscape fragmentation trends in arid and semi-arid zones can be successfully achieved using metrics derived from satellite spectral information.


Introduction
Oil resources are geographically unevenly distributed throughout the world. The concentration of reserves does not coincide with the main consumption areas, being that, in general, the concentration of reserves is in developing countries, whereas, most of the industrialized countries are importers of oil. The main oil consumers are as follows: The United States of America, Japan and Germany, while the main producers are Saudi Arabia, the United States of America and the countries of the North Sea (such as Norway and the United Kingdom). In terms of the ranking of producing countries, Argentina is positioned in 25th place [1]. In Argentina, the oil exploitation started in 1907 in the San Jorge Gulf basin (Chubut and Santa Cruz provinces). This basin has a total area estimated at around 200,000 km 2 and currently produces around 37,000 m 3 of oil per day, obtained from more than 6000 oil wells (9000 are inactive) [2].
The oil industry has significant economic importance in both producing and consumer countries, being that it is a crucial element of national and international economic policies. According to the cycles that characterize this economic activity, the oil industry had a production crisis that declined during the years of 2009 to 2014. Future exploration projects will focus on obtaining oil through non-conventional techniques in mature fields where the application of assisted recovery techniques is Sustainability 2019, 11, 5956 3 of 20 that reduces variable redundancy and extracts uncorrelated information represented by a new set of orthogonal (uncorrelated) variables [29][30][31][32].
In the hydrocarbon basin under investigation in this paper, the main anthropogenic disturbances that alter landscape patterns and generate new patches and corridors are mainly sheep farming and the petroleum industry. In particular, the latter generates roads, extraction facilities, machinery service areas, and pipelines built during the different stages of exploration, exploitation and transportation, thus producing specific disturbances such as the removal of vegetation cover and soil compaction [33]. In extra-Andean Patagonia, sheep were introduced at the beginning of the 20th century through a system of farms with sizes ranging from 10,000 to 50,000 ha [34,35]. Sheep grazing is considered one of the largest consumers of natural vegetation at a global level [36,37] since it generates changes in vegetation communities [38] and affects soil structure and resources.
In 2003, Latin America entered the Commodities Consensus [39], that underlines the consolidation of a new economic, political and ideological order, sustained by the boom in international prices of raw materials and consumer goods demanded by wealthier nations and emerging powers, which generates undoubted comparative advantages in economic growth, but at the same time, produces new asymmetries and deep inequalities in Latin American societies. Oil exploitation has been characterized by cyclical behaviors with ups and downs, as for example, in the boom of the last decade and the crisis currently ongoing today. In 2015, the price of an oil barrel fell by 41% in relation to the previous year, and according to the International Energy Agency, the price remained low throughout the whole 2016 year.
The aim of this study is to quantitatively evaluate landscape fragmentation in the hydrocarbon basin located in the arid and semi-arid zones of Patagonia (Argentina) over a period of 12 years. To this aim, the smallest package of metrics that best explain the fragmentation process in the study area has been herein selected, analyzed and discussed. In particular, we sought to quantify the fragmentation rates in three landscape units and to explain their behavior between 2001 and 2013, using landscape metrics at different scales of resolution. This enabled us to quantify the effects of fragmentation and support a mitigation strategy to reduce the environmental impacts produced by the oil industry.

The Study Was Performed in the Hydrocarbon Basin Located in the Chubut and Santa Cruz Provinces
The hydrocarbon basin of San Jorge Gulf covers an area of approximately 170,000 km 2 , and is located in Patagonia, Argentina ( Figure 1). The most important city in the region is Comodoro Rivadavia with a total of 186,583 habitants [40]. This region is divided into three landscape units that from east to west are as follows: The western valleys, the plateaus and the coastal valleys. The western valleys have an area of 532 km 2 and are located in the phytogeographic region of the Central District, Subdistrito Chubutense [41], which is the most widespread unit in Central Patagonia. The plateaus have an area of 628 km 2 , among them, the main ones are as follows: Pampa del Castillo and Pampa Salamanca, with an average height of approximately 750 m, constituting one of the most important topographic features of the zone. It has a southwest-northeast orientation, from where many canyons originate, sloping towards the Atlantic. The coastal valleys have a surface area of approximately 300 km 2 , the boundary is the supra-tidal line on San Jorge Gulf (Mar Argentino), the northern and western boundaries are at the maximum level of the plateaus of Castillo and Salamanca; the southern limit is at latitude 45 • 48 S. In this region the climate is semi-arid and cold. The mean annual temperature is 13 • C and precipitation is concentrated during the coldest months of the year (June and July). The average precipitation is 247.50 mm (1981-2010) [42].

Methods
The analysis herein performed was based on the following steps: (i) Selection of sampling sites with and without hydrocarbon activity in the landscape units ( Figure 2), (ii) digitization of landscape elements, (iii) calculation of metrics, (iv) running PCA and selection of the most representative landscape metrics; (v) re-sampling of the satellite images; and (vi) analysis of the obtained results. A summary of these methods is shown in Figure 3.

Methods
The analysis herein performed was based on the following steps: (i) Selection of sampling sites with and without hydrocarbon activity in the landscape units ( Figure 2), (ii) digitization of landscape elements, (iii) calculation of metrics, (iv) running PCA and selection of the most representative landscape metrics; (v) re-sampling of the satellite images; and (vi) analysis of the obtained results. A summary of these methods is shown in Figure 3.

Methods
The analysis herein performed was based on the following steps: (i) Selection of sampling sites with and without hydrocarbon activity in the landscape units ( Figure 2), (ii) digitization of landscape elements, (iii) calculation of metrics, (iv) running PCA and selection of the most representative landscape metrics; (v) re-sampling of the satellite images; and (vi) analysis of the obtained results. A summary of these methods is shown in Figure 3.    A grid was drawn in vector format to randomly select five sampling sites for each landscape unit with hydrocarbon activity and (three) control sites without hydrocarbon activity.
High resolution multispectral SPOT 5 satellite images (19 February 2013) were used along with Landsat 7 Enhanced Thematic Mapper Plus (ETM +) of medium spatial resolution (19 December 2001). The data pre-processing consisted in rectification, co-registration and resampling of the different spatial resolutions (10, 30, 60 m) made using the closest neighbor (with an error less than 0.5) [43]. Both of the two images were rectified to the plane coordinate system POSGAR 1994.
A mesh of 77 polygons were placed on the satellite images in a vector format with a resolution of 2500 ha. Five random sites were selected in each of the landscape units where both the linear (roads, tracks, pipelines and seismic lines) and nonlinear elements (oilfields) were digitized. Then the transformation of the vector files to raster images was carried out. The rasterization of the initial maps was generated using a binary code: Anthropogenic activity (linear + nonlinear elements) and the matrix (natural vegetation) ( Figure 4). The process was carried out using QGIS v. 2.14 [44]. A grid was drawn in vector format to randomly select five sampling sites for each landscape unit with hydrocarbon activity and (three) control sites without hydrocarbon activity.
High resolution multispectral SPOT 5 satellite images (19 February 2013) were used along with Landsat 7 Enhanced Thematic Mapper Plus (ETM +) of medium spatial resolution (19 December 2001). The data pre-processing consisted in rectification, co-registration and resampling of the different spatial resolutions (10, 30, 60 m) made using the closest neighbor (with an error less than 0.5) [43]. Both of the two images were rectified to the plane coordinate system POSGAR 1994.
A mesh of 77 polygons were placed on the satellite images in a vector format with a resolution of 2500 ha. Five random sites were selected in each of the landscape units where both the linear (roads, tracks, pipelines and seismic lines) and nonlinear elements (oilfields) were digitized. Then the transformation of the vector files to raster images was carried out. The rasterization of the initial maps was generated using a binary code: Anthropogenic activity (linear + nonlinear elements) and the matrix (natural vegetation) ( Figure 4). The process was carried out using QGIS v. 2.14 [44]. Some studies suggest that the most suitable and representative metrics may differ at class and landscape levels [28], therefore, in this study, a large number of metrics were obtained and calculated at patch, class and landscape levels. This set gave measurements of the landscape spatial structure in terms of form, complexity, connectivity and fragmentation. The landscape metrics were calculated using the FRAGSTATS 4.2 program. It is a computer software program designed to compute a wide variety of landscape metrics for categorical map patterns [32], which is the program most used to calculate landscape metrics [45]. The cartography used as input for the calculation of the metrics was a rasterized image reclassified into a binary code. The landscape metrics were computed using the four-cell rule with nearest neighbor and those used for PCA. The list of landscape metrics used in the study is shown in Appendix A.
A correlation analysis was performed to assess whether there are any differences between the metrics calculated at the class and landscape levels. Nine metrics were selected at the landscape level and eleven at the class level. The selection was based on some previous studies [27,29]. In order to obtain the minimum set of metrics that most adequately describe the landscape pattern, PCA was performed. Following this, a correlation analysis of each metric was undertaken at both class and landscape levels to evaluate if the spatial scale affects the metric.
To evaluate the fragmentation, 20 landscape metrics were analyzed and after a redundancy analysis, only four were selected: Patch density (PD), effective mesh size (MESH), the Clumpiness Index (CLUMPY) and the fractal dimension of the perimeter-area (PAFRAC).
The minimum package of metrics obtained in this study includes those that express form, like PAFRAC, aggregation such as CLUMPY and PD, and contagion and interspersion such as MESH. Form metrics measure the geometric complexity of a landscape, as well as the influence of the interaction between the shape and size of the patch on ecological processes. In more detail, the PD measures the number of patches present in 100 hectares [32], and it is characterized, as all of the selected indicators, by the fact that it is not affected by the scale of the analysis [28]. The MESH is based on the probability that two points randomly selected in a region are connected within a unit without encountering a physical barrier [46,47]. The MESH is a measure of fragmentation, it may be influenced by the value of the metric at the patch level [48] but it has a high capacity for determining temporal changes [49] and is one of the most frequently used [38,49,50]. The PAFRAC is a metric that Some studies suggest that the most suitable and representative metrics may differ at class and landscape levels [28], therefore, in this study, a large number of metrics were obtained and calculated at patch, class and landscape levels. This set gave measurements of the landscape spatial structure in terms of form, complexity, connectivity and fragmentation. The landscape metrics were calculated using the FRAGSTATS 4.2 program. It is a computer software program designed to compute a wide variety of landscape metrics for categorical map patterns [32], which is the program most used to calculate landscape metrics [45]. The cartography used as input for the calculation of the metrics was a rasterized image reclassified into a binary code. The landscape metrics were computed using the four-cell rule with nearest neighbor and those used for PCA. The list of landscape metrics used in the study is shown in Appendix A.
A correlation analysis was performed to assess whether there are any differences between the metrics calculated at the class and landscape levels. Nine metrics were selected at the landscape level and eleven at the class level. The selection was based on some previous studies [27,29]. In order to obtain the minimum set of metrics that most adequately describe the landscape pattern, PCA was performed. Following this, a correlation analysis of each metric was undertaken at both class and landscape levels to evaluate if the spatial scale affects the metric.
To evaluate the fragmentation, 20 landscape metrics were analyzed and after a redundancy analysis, only four were selected: Patch density (PD), effective mesh size (MESH), the Clumpiness Index (CLUMPY) and the fractal dimension of the perimeter-area (PAFRAC).
The minimum package of metrics obtained in this study includes those that express form, like PAFRAC, aggregation such as CLUMPY and PD, and contagion and interspersion such as MESH. Form metrics measure the geometric complexity of a landscape, as well as the influence of the interaction between the shape and size of the patch on ecological processes. In more detail, the PD measures the number of patches present in 100 hectares [32], and it is characterized, as all of the selected indicators, by the fact that it is not affected by the scale of the analysis [28]. The MESH is based on the probability that two points randomly selected in a region are connected within a unit without encountering a physical barrier [46,47]. The MESH is a measure of fragmentation, it may be influenced by the value of the metric at the patch level [48] but it has a high capacity for determining temporal changes [49] and is one of the most frequently used [38,49,50]. The PAFRAC is a metric that informs us about the complexity of the form (it is binary and provides one for shapes with simple perimeters and two for complex perimeters). It is a sound metric that is not influenced by the scale of the analysis [51]. Finally, CLUMPY isolates the configuration and area components, giving an effective index of fragmentation calculated using an adjacency matrix, which shows the frequency at which different types of patch pairs appear side by side on the map [52].
To evaluate the effect of the pixel size of the Landsat and SPOT satellite images, spatial resampling was performed using the nearest neighbor method. This data processing was carried out using QGIS v. 2.14 and GRASS software (resampled tool). The spatial resolution was 10, 30 and 60 m for both of the images used. Following this, these classified images were incorporated into the FRAGSTATS program to calculate landscape metrics.

Results
The results are presented in two sections: (i) Analysis of landscape metrics, and (ii) analysis of the temporal change in fragmentation. informs us about the complexity of the form (it is binary and provides one for shapes with simple perimeters and two for complex perimeters). It is a sound metric that is not influenced by the scale of the analysis [51]. Finally, CLUMPY isolates the configuration and area components, giving an effective index of fragmentation calculated using an adjacency matrix, which shows the frequency at which different types of patch pairs appear side by side on the map [52].

Analysis of Landscape Metrics at the Class Level
To evaluate the effect of the pixel size of the Landsat and SPOT satellite images, spatial resampling was performed using the nearest neighbor method. This data processing was carried out using QGIS v. 2.14 and GRASS software (resampled tool). The spatial resolution was 10, 30 and 60 m for both of the images used. Following this, these classified images were incorporated into the FRAGSTATS program to calculate landscape metrics.

Results
The results are presented in two sections: (i) Analysis of landscape metrics , and (ii) analysis of the temporal change in fragmentation .  Based on the obtained results, the CLUMPY metric was selected at the class level which explains 93.44% of the total variance ( Figure 6, Table 1). Based on the obtained results, the CLUMPY metric was selected at the class level which explains 93.44% of the total variance ( Figure 6, Table 1).  At the landscape level, PCA showed that factors 1 and 2 explain more than 95% of the total variance of the metrics analyzed ( Figure 7, Table 2). From this analysis, the metrics selected were as follows: PD, MESH and PAFRAC, which are the ones that best explain the total variability in this case.  The correlation between the levels of analysis (class vs. landscape) was high and significant in all the analyzed metrics ( Figure 8).  At the landscape level, PCA showed that factors 1 and 2 explain more than 95% of the total variance of the metrics analyzed ( Figure 7, Table 2). From this analysis, the metrics selected were as follows: PD, MESH and PAFRAC, which are the ones that best explain the total variability in this case.  At the landscape level, PCA showed that factors 1 and 2 explain more than 95% of the total variance of the metrics analyzed ( Figure 7, Table 2). From this analysis, the metrics selected were as follows: PD, MESH and PAFRAC, which are the ones that best explain the total variability in this case.  The correlation between the levels of analysis (class vs. landscape) was high and significant in all the analyzed metrics ( Figure 8).  The correlation between the levels of analysis (class vs. landscape) was high and significant in all the analyzed metrics ( Figure 8). When resampling the satellite images, the ratio of the PD calculated with Landsat and SPOT images was significant and higher than 69% for the spatial resolutions of 10, 30 and 60 m. The PAFRAC showed a significant relationship higher than 82% for all the scales analyzed, and eventually, the MESH showed a relationship of 71% with a spatial resolution of 10 m, while when the resolution was 30 and 60 m the ratio was 95% and 74%, respectively ( Figure 9). When resampling the satellite images, the ratio of the PD calculated with Landsat and SPOT images was significant and higher than 69% for the spatial resolutions of 10, 30 and 60 m. The PAFRAC showed a significant relationship higher than 82% for all the scales analyzed, and eventually, the MESH showed a relationship of 71% with a spatial resolution of 10 m, while when the resolution was 30 and 60 m the ratio was 95% and 74%, respectively ( Figure 9).

Analysis of the Temporal Change in Fragmentation
Between 2001 and 2013, fragmentation was greater at the sites with hydrocarbon activity than without ( Figure 10). The average density of the patches in the landscape units went down from 60 to 14 patches per 100 hectares at the sites without any hydrocarbon activity between the years 2001 and 2013, whereas the mean density of patches increased from 309 to 444 patches at the sites with activity. The mean mesh size (MESH) remained constant between years in all three landscape units, whereas at the sites with oil production they decreased, mainly in the coastal valleys. The PAFRAC decreased in all three landscape units with hydrocarbon activity in the period analyzed, this decrease being more evident in the plateaus and the coastal valleys. On the other hand, this metric did not show a common pattern in the landscape units at the sites without any activity. CLUMPY showed an increase in all three landscape units between 2001 and 2013 at the sites with hydrocarbon activity. This increase was more evident in the coastal valleys. This metric did not show a common pattern of behavior in the landscape units at the control sites over time.

Analysis of the Temporal Change in Fragmentation
Between 2001 and 2013, fragmentation was greater at the sites with hydrocarbon activity than without ( Figure 10). The average density of the patches in the landscape units went down from 60 to 14 patches per 100 hectares at the sites without any hydrocarbon activity between the years 2001 and decreased in all three landscape units with hydrocarbon activity in the period analyzed, this decrease being more evident in the plateaus and the coastal valleys. On the other hand, this metric did not show a common pattern in the landscape units at the sites without any activity. CLUMPY showed an increase in all three landscape units between 2001 and 2013 at the sites with hydrocarbon activity. This increase was more evident in the coastal valleys. This metric did not show a common pattern of behavior in the landscape units at the control sites over time.

Discussion
Understanding and assessing trends in landscape fragmentation in arid and semi-arid zones can be achieved using metrics derived from spectral information. However, selecting and interpreting the minimum set of metrics that represent change most effectively is a challenge, since all metrics have limitations that restrict their use and application [52]. Multivariate analysis enables a reduction in the number of metrics by showing the redundancy between them. In this study, the metrics that best expressed the pattern of change in a Patagonian hydrocarbon basin were reduced to four from fourteen; several papers have focused on explaining the reduction of metrics for avoiding redundancy between them (e.g., [28][29][30][31]). The minimum package of metrics obtained in this study includes those that express form, like PAFRAC, aggregation such as CLUMPY and PD, and contagion and interspersion such as MESH. Form metrics measure the geometric complexity of a landscape, as well as the influence of the interaction between the shape and size of the patch on ecological processes. Aggregation metrics measure the tendency of patch types to be spatially aggregated which refers to the texture of the landscape. The contagion and interspersion metrics are based on adjacent patch types, considering information about border segments [46]. Numerous studies in the last two decades have reported on the behavior of landscape metrics; however, few have focused on changes in resolution [28]. Li and Wu [53] have shown that changes in the level of analysis could affect the behavior of the metrics. In this study, we show a high relationship between the metrics analyzed at both class and landscape levels (r 2 > 0.97) and the scales (r 2 > 0.69). This could be due to the fact that a binary classification was analyzed between an anthropogenic activity and a matrix without any disturbance, the latter being dominant in the landscapes irrespective of the scale.
Our results showed that the landscape units in the Patagonian hydrocarbon basin have undergone an important fragmentation process between 2001 and 2013. This loss of connectivity was more evident in the coastal valleys than in the rest of the landscape units. The first oil deposits in the region were located near to the coast at the beginning of the 20th century, in the costal valleys that descend to the sea from the upper terraced levels. The loss of connectivity was a consequence of hydrocarbon activity, since the applied metrics showed a tendency to fragmentation when compared with the control sites without any activity. This loss of connectivity, intensified in the analyzed period, facilitated the serious desertification process, caused by the mismanaging of the sheep grazing that occurred in the last century. The advance of desertification in the arid and semi-arid zones of Patagonia is one of the main socioecological and environmental problems present. The fragmentation of the plant cover matrix in the hydrocarbon basin under study can affect ecosystem functional and structural attributes such as increased runoff, diminished above-ground net primary productivity (and, consequently, secondary production), and increased physiognomic changes (e.g., replacement of grasslands by shrub-lands). This last aspect was already observed and measured in arid and semi-arid Patagonia were plant communities showed changes in land cover; in the western valleys the shrub-lands increased their coverage demonstrating the process of expansion of thickets of shrub plants in environments dominated by grasses or other herbaceous plants.
The probability that two animals located at two different points within the study area can be found without having to cross a barrier, such as a road or an urban area [49], decreased as a result of the oil activity; the MESH had a value close to 2500 ha in the control sites and of 1957 ha in the disturbed sites. The PD increased from 37 to 376 patches per 100 hectares at sites with hydrocarbon activity and it was observed that over time the forms had been simplified, as shown by PAFRAC, which decreased from a value of 1.38 to 1.31, possibly due to the increase in infrastructure and roads required by the petroleum industry. In the control sites, the temporal difference of the shape was not so marked, showing that there have been no modifications in the landscape configuration by the petroleum industry. CLUMPY showed that the patches are more aggregated in sites with hydrocarbon activity, which increased during the period analyzed. In the control sites, the aggregation between patches presented a random pattern. It should be noted that CLUMPY isolates the configuration component from the area component, thus giving an effective index of fragmentation that is not confused with changes in the area [21]. Between 2003 and 2004 there was an increase in the international price of a barrel of crude oil, and a consequent increase in oil exploitation, which is reflected in the opening of new locations, roads and tracks. This is due to the oil boom that has occurred since 2003 in the hydrocarbon basin in the Chubut province, as part of neo-extractivism and an era of well-being and prosperity that has had repercussions on the configuration of the regional landscape [54].

Conclusions
As a whole, this study showed that landscape is a complex system, with dynamics, spatial configuration, structure and functionality that are the result of the interaction between natural, economic and socio-cultural factors. In particular, in this paper, landscape fragmentation rates were quantitatively assessed in the hydrocarbon basin, located in the Chubut and Santa Cruz provinces of the arid and semi-arid zones of Patagonia (Argentina). It is important to measure and document fragmentation in the landscape, to support mitigation strategies, sustainable planning and policy development for hydrocarbon activity.
The fragmentation rates were herein evaluated over a period of twelve years, in three landscape units, using PCA to select the smallest package of metrics that best explain the fragmentation process. The minimum package of metrics obtained in this study includes those that express form, like PAFRAC, aggregation such as CLUMPY and PD, and contagion and interspersion such as MESH. The advantages of the methodology herein adopted is the use of a supervised classification to evaluate and categorize the impacts produced by the oil industry. The disadvantage is mainly linked with the binary approach which does not allow the inclusion of different land use and land cover types. As a future project, it would be interesting to include information on the combination of different land uses and land covers and, moreover, further research should be extended to consider the fragmentation of the landscape from an interdisciplinary perspective including the relationship between cultural and natural landscapes. Appendix A Table A1. Definition and description of FRAGSTATS metrics [29].

Landscape Metric Description Range Units Formula
Total Area (TA) TA equals the sum of the areas (m 2 ) of all patches of the corresponding patch type, divided by 10,000 (to convert to hectares); that is, total class area.
TA ≥ 0 Ha. TA = n j=l aij Clumpiness Index (CLUMPY) Equals the proportional deviation of the proportion of like adjacencies involving the corresponding class from that expected under a spatially random distribution.
Number of like adjacencies (joins) between pixels of patch type (class) i based on the double-count method. g ik = Number of adjacencies (joins) between pixels of patch types (classes) i and k based on the double-count method. P i = Proportion of the landscape occupied by patch type (class) i.

Aggregation Index (AI)
Aggregation index is calculated from an adjacency matrix, which shows the frequency with which different pairs of patch types (including like adjacencies between the same patch type) appear side-by-side on the map. 0 AI 100 Percent. AI = g ii maxg ii (100) g ii = Number of like adjacencies (joins) between pixels of patch type (class) i based on the single-count method. max-g ii = Maximum number of like adjacencies (joins) between pixels of patch type (class) i based on the single-count method.