Modelling Soil Carbon Content in South Patagonia and Evaluating Changes According to Climate, Vegetation, Desertification and Grazing

Abstract

In Southern Patagonia, a long-term monitoring network has been established to assess bio-indicators as an early warning of environmental changes due to climate change and human activities. Soil organic carbon (SOC) content in rangelands provides a range of important ecosystem services and supports the capacity of the land to sustain plant and animal productivity. The objectives in this study were to model SOC (30 cm) stocks at a regional scale using climatic, topographic and vegetation variables, and to establish a baseline that can be used as an indicator of rangeland condition. For modelling, we used a stepwise multiple regression to identify variables that explain SOC variation at the landscape scale. With the SOC model, we obtained a SOC map for the entire Santa Cruz province, where the variables derived from the multiple linear regression models were integrated into a geographic information system (GIS). SOC stock to 30 cm ranged from 1.38 to 32.63 kg C m⁻². The fitted model explained 76.4% of SOC variation using as independent variables isothermality, precipitation seasonality and vegetation cover expressed as a normalized difference vegetation index. The SOC map discriminated in three categories (low, medium, high) determined patterns among environmental and land use variables. For example, SOC decreased with desertification due to erosion processes. The understanding and mapping of SOC in Patagonia contributes as a bridge across main issues such as climate change, desertification and biodiversity conservation.

1. Introduction

Scientists and land managers of natural ecosystems acknowledge the importance of long-term monitoring systems for evaluating responses to disturbances (climate change or human activities) and providing baselines to evaluate potential changes [1,2]. In this context, since 2002, a long-term monitoring system (defined as a repeated field-based empirical measurements collected continuously and analyzed for at least 10 years) was established to monitor natural ecosystems and to produce scientific research focused on ecosystem function and ecosystem services, as well as on trends in biodiversity and the interactions between natural environments and land-use activities throughout southern Patagonia, Argentina [3].

Sustainable management of rangeland (rangeland can be defined as extensive areas of land that are occupied by native herbaceous or shrubby vegetation which are grazed by domestic or wild herbivores) for livestock production is an important economic activity in Southern Patagonia. Herbivores are known to be key drivers of soil processes in rangelands [4]. Changes in herbivore pressure (e.g., stocking rate) can have important consequences for ecosystem functioning [5]. Grazing can alter soil carbon stocks by changing the quality (e.g., dung, urine and litter inputs) and/or quantity (e.g., by causing compensatory regrowth in vegetation and/or by changing patterns of biomass allocation in the standing vegetation) of carbon that enters the soil. Grazing may also affect soil carbon stocks by altering rates of organic matter decomposition [4,6,7,8,9]. Soil organic carbon (SOC) content is important for ecosystem service provision, for example by supporting biodiversity [10], increasing soil aggregation, limiting soil erosion, and increasing water holding capacities [11]. Soil carbon also supports the capacity of the land to sustain plant and animal productivity and this potential depends on how rangelands are managed for livestock production [12]. Soil carbon is therefore a useful indicator for assessing the sustainability of livestock production on rangelands.

Rangelands are important economically and also culturally in Patagonia as rangelands provide the people of the region with a sense of place [13]. Given that maintenance of soil carbon is so important for the long term economic viability of rangelands it is surprising that little scientific research has focused on soil properties and on how these soil properties relate to grazing and land management more generally in these ecosystems. In Patagonia over the last 70 years, we have witnessed extensive degradation of once productive steppe ecosystems (desertification) [14]. As a result, stakeholders in Southern Patagonia developed a certification scheme for sustainable land management in the region [15]. However, all of the indicators are qualitative and quantitative indicators of rangeland condition are needed. SOC is one such possible indicator [16]. Data on carbon storage in forests, grasslands and shrublands at plant and stand level in Patagonia have been reported [6,7,17,18,19]. However, similar data for Patagonian rangelands are notable for their absence.

Long-term grazing intensity in arid and semiarid regions, such as southern Patagonia, may affect soil C stocks. The effects of grazing on soil C stocks likely also interact with other environmental variables that drive soil C [20]. Thus, variations in soil carbon pool may be influenced by climatic and topographic conditions. Understanding how these variables interact with grazing to alter soil carbon stocks at a regional scale is critical for understanding the impacts of land use decisions on both the sustainability of rangeland management and on atmospheric carbon concentrations and the climate.

The objectives in this study were (i) to model SOC stocks to 30 cm at a regional scale using climatic, topographic and vegetation variables; and (ii) to establish a regional baseline for SOC so that SOC can be used as an indicator of rangeland condition and therefore of the sustainability of land management. We hypothesize that (1) SOC would be lower where environmental conditions are harsh (low soil moisture conditions and high altitudes) and (2) that adverse environmental conditions would have a larger effect on soil carbon stocks than land use (stocking rate) at the regional scale in Patagonia.

2. Material and Methods

The study was conducted in of the PEBANPA (Parcelas de Ecología y Biodiversidad de Ambientes Naturales en Patagonia Austral-Biodiversity and Ecological long-term plots in Southern Patagonia) network of permanent plots [3]. There are 145 sites in the PEBANPA network of permanent plots (Figure 1), all of which were used in this analysis. Further detail of environmental conditions across Santa Cruz Province can be found in Peri et al. [3].

2.1. Soil Organic Carbon

For all 145 sites, we extracted data of SOC concentration (% C) and soil bulk density (BD) from the PEBANPA database (see Peri et al. [3] for details of the methodology). At each site, soil samples were collected from nine randomly selected points within a 20 m × 40 m quadrat using a hand auger (30 cm depth). Coarse root debris > 2 mm from soil samples had been removed by sieving. To reduce the number of chemical analyses we pooled individual soil samples into combined samples. From the nine samples collected within each quadrat, we created three composite samples so that each composite sample contained an equal proportion of soil from three auger holes (n = 3 for each site). The samples were finely ground to below 2 μm using a tungsten-carbide mill. Measurements of SOC concentration were derived from the dry combustion (induction furnace) method. Soil BD was estimated using the cylindrical core method (n = 3) by collecting a known volume of soil using a metal tube pressed into the soil (intact core), and determining the weight after drying. Knowing soil BD and depth of soil layers (0 to 30 cm) (Z), we applied the following equation to calculate the soil carbon stock:

$$\mathrm{SoilCstock}=\%\mathrm{C}\raisebox{1ex}{$\mathrm{g}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{g}$}\right.\times \mathrm{BD}\left(\raisebox{1ex}{$\mathrm{g}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{cm}}^3$}\right.\right)\mathrm{Xdepthtobedrock}\left(\mathrm{cm}\right).$$

(1)

2.2. GIS-Derived Independent Variables

Climatic, topographic, landscape and land-use variables were obtained for each sampling point (

Table 1. Explanatory variables used in soil carbon stock analysis.

Category	Description	Code	Unit	Data Source
Climate	mean annual temperature	AMT	°C	WorldClim (1)
	mean diurnal range	MDR	°C	WorldClim (1)
	isothermality	ISO	%	WorldClim (1)
	temperature seasonality	TS	°C	WorldClim (1)
	max temperature of warmest month	MAXWM	°C	WorldClim (1)
	min temperature of coldest month	MINCM	°C	WorldClim (1)
	temperature annual range	TAR	°C	WorldClim (1)
	mean temperature of wettest quarter	MTWEQ	°C	WorldClim (1)
	mean temperature of driest quarter	MTDQ	°C	WorldClim (1)
	mean temperature of warmest quarter	MTWAQ	°C	WorldClim (1)
	mean temperature of coldest quarter	MTCQ	°C	WorldClim (1)
	mean annual precipitation	AP	mm years−1	WorldClim (1)
	precipitation of wettest month	PWEM	mm years−1	WorldClim (1)
	precipitation of driest month	PDM	mm years−1	WorldClim (1)
	precipitation seasonality	PS	%	WorldClim (1)
	precipitation of wettest quarter	PWEQ	mm years−1	WorldClim (1)
	precipitation of driest quarter	PDQ	mm years−1	WorldClim (1)
	precipitation of warmest quarter	PWAQ	mm years−1	WorldClim (1)
	precipitation of coldest quarter	PCQ	mm years−1	WorldClim (1)
	global potential evapo-transpiration	EVTP	mm years−1	CSI (2)
Topography	elevation	ELE	m.a.s.l.	DEM (3)
	slope	SLO	%	DEM (3)
	aspect	ASPC	cosine	DEM (3)
	aspect	ASPS	sine	DEM (3)
	distance to water bodies	DWB	km	SIT Santa Cruz (4)
	distance to rivers	DR	km	SIT Santa Cruz (4)
Landscape and land-use	normalized difference vegetation index	NDVI	dimensionless	MODIS (5)
	net primary productivity	NPP	gr C m−2 year−1	MODIS (6)
	desertification	DES	degree	CENPAT (7)
	ecological area	EA	dimensionless	SIT Santa Cruz (4)
	stocking rate	SR	(ewe/ha/year)	SIT Santa Cruz (4)
	carrying capacity	RF	(ewe/ha/year)	SIT Santa Cruz (4)

⁽¹⁾ Hijmans et al. [23]; ⁽²⁾ Consortium for Spatial Information (CSI) [24]; ⁽³⁾ Farr et al. [25]; ⁽⁴⁾ SIT-Santa Cruz (http://www.sitsantacruz.gob.ar); ⁽⁵⁾ ORNL DAAC [26]; ⁽⁶⁾ Zhao et al. [27]; ⁽⁷⁾ Del Valle et al. [21].

). The methods used to generate the GIS-derived independent variables were also described in Peri et al. [3].

2.3. Modelling and Data Analyses

A pre-selection of variables was performed based on Pearson’s correlation indices obtained from paired analyses and considering the strength of the linear relationship (−1 to +1) and a p-value less than 0.05, with a confidence level of 95%.

For modelling, we used a stepwise multiple regression to identify which variables among these uncorrelated variables helped to explain SOC variation at landscape level. We employed a p value of <0.05 for the significance of each variable to be included into the model, analyzing the utility of the inclusion of the constant in the model, and used 100 steps for the final model selection. The model was evaluated through the standard error (SE) of estimation (the r²-adj), defined as the average of the difference between predicted versus observed values, and the mean absolute error (AE) defined as the average of the difference between predicted versus the observed absolute values (Statgraphics Centurion software, Statpoint Technologies, The Plains, VA, USA).

To test the model across different gradients we performed a calibration procedure, using the same database employed for the modelling (observed vs. modelled). The first test was carried out by analyzing the mean and absolute errors (differences between observed and modelled values of SOC expressed as kg m²). Secondarily, we tested the model performance by comparing SOC across different gradients of natural and human related variables: (i) vegetation types; (ii) stocking rates; (iii) soil covers (bare soil, shrubs, dwarf-shrubs, grasses, herbs, trees) (see further description of data and calculations in Peri et al. [27]).

With the SOC model, we obtained a SOC map for the entire Santa Cruz province (Argentina), where the variables derived from the multiple linear regression models were integrated into a geographical information system (GIS) using ArcMap 10.0 software [28]. For the SOC map, SOC values were assigned to three categories: low (0.01–4.47), medium (4.48–5.46), and high soil carbon stock (5.47–16.24 kg C m⁻²). The limits of each SOC class were defined to contain an equal quantity of pixels for the whole province. For each SOC class (low, medium and high) we calculated the mean values and standard deviation of 25 continuous variables (including climate, topographic and landscape variables, see Table 1) using data from the entire province. In addition, the mean values and standard deviation of SOC were also calculated for discrete variables of the animal stocking rate, forage receptivity, ecological area, and desertification.

3. Results

Across Santa Cruz province, SOC stock to 30 cm depth ranged from 1.38 to 32.63 kg C m⁻². Climate variables presented correlation indexes between 0.07 and 0.99. Some climate variables were greatly influenced by the landscape variables (e.g., precipitation of driest quarter was strongly negatively correlated with the NDVI index, 0.78). Landscape variables presented correlation indexes between 0.06 and 0.88. Potential evapo-transpiration was strongly and negatively correlated with the desertification index, 0.66. Finally, topography variables were the group with the lowest correlation indices (0.06 to 0.39). Most of the variables were highly correlated to SOC stock using the Pearson’s correlation index (

Table 2. Pearson’s correlation index used in soil carbon stock (SOC) analysis. (see Table 1 for variables definition).

Category	Variables	SOC
		Correlation	p-Value
Climate	AMT	−0.40	<0.001
	MDR	−0.32	<0.001
	ISO	0.35	<0.001
	TS	−0.51	<0.001
	MAXWM	−0.44	<0.001
	MINCM	−0.27	=0.001
	TAR	−0.44	<0.001
	MTWEQ	−0.10	=0.252
	MTDQ	−0.42	<0.001
	MTWAQ	−0.44	<0.001
	MTCQ	−0.30	<0.001
	AP	0.48	<0.001
	PWEM	0.41	<0.001
	PDM	0.54	<0.001
	PS	−0.42	<0.001
	PWEQ	0.42	<0.001
	PDQ	0.52	<0.001
	PWAQ	0.56	<0.001
	PCQ	0.39	<0.001
	EVTP	−0.50	<0.001
Topography	ELE	−0.06	=0.502
	SLO	0.24	=0.005
	ASPC	−0.16	=0.052
	ASPS	0.13	=0.118
	DWB	−0.15	=0.073
	DR	0.01	=0.868
Landscape and land-use	NDVI	0.60	<0.001
	NPP	0.53	<0.001
	DES	−0.46	<0.001

), where NDVI (0.600, p < 0.001) was the most correlated. The variables ASPC, DWB, ASPS, MTWEQ, ELE and DR were not significant correlated with SOC stock using the Pearson’s correlation index (Table 2).

The stepwise multiple regression selected three variables for the modelling: isothermality (ratio of average day variation in temperature divided by annual variability in temperature) (ISO, %), precipitation seasonality (PS, %) and normalized vegetation index (NVDI, dimensionless). These variables presented the best statistics, a high correlation with SOC stock, and low correlation among them (ISO × PS = −0.562, p < 0.001; ISO × NDVI = 0.384, p < 0.001; PS × NDVI = −0.471, p < 0.001). The inclusion of the constant in the model decreased the goodness of fit of the model; for this reason we decided to not include a constant in the modelling. The fitted model (r²-adj = 0.764; F = 156.1; SE = 4.08; AE = 2.71) explained 76.4% of variation in SOC values, and was expressed as:

SOC (kg m²) = 0.11 × ISO − 0.10 × PS + 12.03 × NDVI

(2)

The map of the adjusted SOC model showed a continuous decline from the northeast and central areas of Santa Cruz province where most forests and shrublands are growing to the south and southwest where rangelands dominate (Figure 2).

The characteristics of the climatic and topographic variables according to the SOC map developed for the entire study area and the different SOC map quantities (low, medium, high) determined patterns of change among environmental and land use variables (

Table 3. Mean (standard deviation) values of climatic, topographic and vegetation variables classified according to the soil carbon classes: low (0.01–4.47 kg C m⁻²), medium (4.48–5.46 kg C m⁻²) and high (5.47–16.24 kg C m⁻²). (see Table 1 for variables definition).

Variable	Total	Low	Middle	High
AMT	7.77 (2.40)	8.57 (2.66)	8.48 (2.06)	6.25 (1.59)
MDR	10.33 (0.65)	10.45 (0.67)	10.49 (0.54)	10.07 (0.67)
ISO	46.1 (1.5)	45.3 (1.1)	45.8 (1.1)	47.2 (1.5)
TS	4.47 (0.44)	4.67 (0.43)	4.59 (0.32)	4.16 (0.38)
MAXWM	19.56 (3.16)	20.77 (3.47)	20.50 (2.53)	17.42 (2.10)
MINCM	(-2.65) (2.20)	(−2.09) (2.43)	(−2.15) (2.12)	(−3.71) (1.54)
TAR	22.21 (1.76)	22.86 (1.73)	22.65 (1.38)	21.13 (1.60)
MTWEQ	5.67 (2.95)	6.03 (3.17)	4.58 (1.93)	6.42 (3.21)
MTDQ	9.81 (3.72)	12.46 (3.45)	10.06 (2.83)	6.90 (2.42)
MTWAQ	13.21 (2.84)	14.34 (3.10)	14.04 (2.28)	11.25 (1.84)
MTCQ	1.85 (2.11)	2.45 (2.34)	2.42 (1.93)	0.68 (1.43)
AP	245.92 (181.38)	251.45 (219.82)	222.28 (140.29)	262.68 (169.31)
PWEM	30.15 (18.90)	32.42 (22.31)	27.47 (14.82)	30.40 (18.15)
PDM	13.61 (12.52)	13.21 (15.48)	12.48 (9.88)	15.13 (11.30)
PS	24.41 (6.57)	29.26 (6.55)	23.58 (4.54)	20.35 (4.97)
PWEQ	79.81 (53.24)	84.29 (62.82)	72.67 (41.44)	82.36 (52.35)
PDQ	46.38 (41.07)	46.00 (51.11)	42.46 (32.03)	50.41 (36.41)
PWAQ	53.62 (42.92)	48.65 (54.03)	49.29 (33.29)	62.92 (36.76)
PCQ	67.33 (46.03)	71.97 (53.70)	62.10 (35.32)	67.56 (45.92)
EVTP	807.88 (101.56)	848.78 (107.69)	839.90 (78.57)	735.15 (71.02)
GAI	0.33 (0.36)	0.33 (0.45)	0.28 (0.30)	0.37 (0.28)
ELE	468.83 (383.85)	556.10 (411.60)	388.62 (366.35)	460.13 (348.21)
SLO	5.00 (5.76)	5.07 (5.99)	4.49 (5.00)	5.43 (6.14)
NDVI	0.21 (0.12)	0.13 (0.06)	0.17 (0.04)	0.32 (0.14)
NPP	1275.68 (684.62)	971.12 (306.50)	1124.78 (422.34)	1708.11 (899.45)

). Mean Annual Temperature (MAT) influenced SOC. SOC values are higher at lower temperatures compared to the average for the entire province (7.8 °C). Other related temperature variables followed the same pattern (TS, MAXWM, MINCM, MTDQ, MTWAQ, MTCQ). However, MTWEQ did not present a clear pattern of variation. Seasonal and daily variations of temperature (MDR and TAR) and isothermality (ISO) did not greatly influence SOC (Table 3). Rainfall (MAP) also influenced SOC. SOC values increased with precipitation. The correlation between SOC and other rainfall variables (PWEM, PDM, PWEQ, PDQ, PWAQ, PCQ) followed a similar pattern. The other studied climatic indices (EVTP and GAI) followed the combined patterns of the temperature and rainfall variables, where SOC values decreased with the evapotranspiration and aridity. SOC values were generally low in the mountain environments with SOC values generally increasing below 460 m above sea level (m.a.s.l.), and the topographic variable slope did not correlate with changes in SOC quantity. As normalized difference vegetation index (NDVI) and net primary productivity (NPP) increased so did SOC stock to 30 cm (Table 3).

Vegetation types also correlated with differences in SOC (3.8–5.5 kg C m⁻² shrubs, 5.9–6.8 kg C m⁻² grasslands, 12.1–12.3 kg C m⁻² forests), while animal stocking rate decreased SOC values along the studied gradient (10.9 kg C m⁻² in enclosures versus 4.6–6.7 kg C m⁻² with medium and high stocking rate) (

Table 4. Mean values (standard deviation, SD) and areas of soil carbon content (kg C m⁻²) at 30 cm depth sorted by discrete variables in Santa Cruz province, Patagonia, Argentina.

Variable	Category	Mean (SD)	Area (km2)
Desertification	No desertification	6.82 (3.38)	15,061
	Slight degraded	7.52 (2.47)	12,085
	Moderate desertification	6.52 (1.80)	34,135
	Moderate to severe desertification	5.11 (1.30)	84,011
	Severe desertification	4.72 (0.93)	63,502
	Very severe desertification	4.26 (0.61)	29,828
simple ecosystem classiffication	Humid magellanic grass steppe	9.64 (1.14)	6056
	Andean vegetation	7.51 (3.70)	15,815
	Dry magellanic grass steppe	7.29 (0.77)	11,796
	Mata negra thicket	6.35 (0.77)	28,374
	Sub-andean grassland	5.88 (1.74)	19,540
	Central plateau	4.62 (0.73)	131,911
	Shrub steppe San Jorge Gulf	4.59 (0.71)	11,990
	Mountains and plateaus	3.76 (0.87)	13,125
potential stocking rate (ewes/ha/yr)	<0.1	4.56 (0.46)	17,119
	0.1 a 0.20	4.59 (0.71)	84,363
	0.13 a 0.25	4.11 (0.68)	18,654
	0.16 a 0.30	3.76 (0.86)	13,166
	0.20 a 0.30	5.83 (1.22)	71,911
	0.3 a 0.4	7.40 (2.84)	27,493
	>0.4	9.71 (1.17)	5917
Actual Stocking rate (ewes/ha/yr)	0.07	4.62 (0.73)	131,758
	0.13	5.64 (1.55)	60,872
	0.18	4.60 (0.71)	11,886
	0.22	7.29 (0.76)	11,707
	0.3	7.52 (3.69)	15,738
	0.5	9.66 (1.14)	5886

). Finally, SOC decreased with desertification gradient (10.6 kg C m⁻² without presence, 8.3 kg C m⁻² low, 5.2 kg C m⁻² at moderate levels of desertification, and 4.4 kg C m⁻² at sites where desertification was pronounced) due to erosion processes.

The calibration of the model showed an average error of −0.01 kg m², and an absolute error of 2.72 kg m². When the performance of the model was tested across different natural and management related gradients, it can be observed that the error dispersion isn’t homogeneous (

Table 5. Model performance analysis using a calibration of the soil carbon content (kg C m⁻²) sorted by discrete variables: (i) vegetation types, (ii) stocking rate, (iii) soil covers (bare soil, shrubs, dwarf-shrubs, grasses, herb, trees).

Vegetation	N	Carbon	Modelled	Mean Error	Absolute Error
Shrub steppe	23	4.54	5.52	−0.98	2.27
Dwarf shrub steppe	28	3.98	4.77	−0.79	1.81
Shrub steppe	5	3.06	3.84	−0.78	1.82
Grass steppe	55	5.35	5.92	−0.57	2.15
Nothofagus antarctica forest	10	12.46	12.08	0.38	4.33
Nothofagus pumilio forest	12	13.73	12.29	1.44	4.12
Grassland	11	11.24	6.67	4.57	5.75
Nothofagus betuloides forest	1	20.25	12.34	7.91	7.91
Stocking rate (ewes/ha/yr)
<0.1	38	8.85	8.15	0.70	3.08
0.1–0.17	41	4.43	4.60	−0.17	1.71
0.18–0.25	30	5.79	6.76	−0.96	2.97
>0.25	36	7.41	7.18	0.23	3.27
% bare soil
<16.7	35	11.06	9.76	1.31	4.68
16.7–30.0	36	5.51	5.82	−0.31	2.10
30.1–48.0	39	5.11	5.57	−0.47	1.98
>48.0	35	4.97	5.45	−0.49	2.22
% shrubs
<1.0	35	6.91	7.11	−0.20	2.48
1–6.0	39	6.17	6.30	−0.13	2.43
6.1–21.3	35	4.32	5.25	−0.93	2.01
>21.3	36	9.03	7.81	1.22	3.96
% dwarf-shrubs
<0.3	35	7.85	8.01	−0.16	2.64
0.4–3.0	38	7.28	6.74	0.54	2.89
3.1–20.0	38	6.62	6.65	−0.03	3.24
>20.0	34	4.58	5.01	−0.43	2.02
% grasses
<10.9	37	8.18	8.84	−0.66	3.83
11.0–22.5	36	5.19	5.14	0.05	1.70
22.6–35.9	36	5.28	6.19	−0.91	2.16
>36.0	36	7.75	6.23	1.51	3.16
% Herbs
<1.8	35	4.83	4.77	0.06	2.01
1.9–4.3	38	4.86	5.53	−0.68	2.01
4.4–14.0	37	6.99	7.09	−0.11	2.60
>14.0	35	9.90	9.13	0.77	4.32
% trees
none	119	5.27	5.43	−0.16	2.33
0.1–50.0	16	11.50	11.84	−0.33	4.63
>50.0	10	14.72	12.38	2.34	4.29
Total general	145	6.61	6.62	−0.01	2.72

). In general, errors increased with SOC quantity. When vegetation type was considered, lower error values were observed in arid grasslands and shrublands, while greater errors were found in humid grasslands and forests.

When bare soil cover increased, the error in SOC predictions decreased. When herb and tree covers decreased, the error in SOC predictions also decreased. Prediction error didn’t vary systematically with stocking rate, or with shrub and/or grass cover (Table 5).

4. Discussion

Our model for SOC prediction was able to account for 76% of the variation of this soil property across the study area, with values ranging from 1.38 to 32.63 kg C m⁻². In the present study, SOC stock to 30 cm was mainly a function of climate and vegetation. As we have already shown [29,30], the prediction and mapping of soil carbon at the macro scale was possible using freely available geospatial data. The correlation between SOC and climate variables (isothermality and seasonality precipitation) may reflect the influence of climate variables on semi-arid ecosystem productivity, which are mainly related to water limitation. This highlights the importance of long term monitoring to know the processes that determine the magnitude of SOC variation, and to forecast how it may operate as climate and land use changes in the future. The effect of stocking rate on SOC was minimal in this analysis. However, grazing may indirectly affect SOC by modifying the type of vegetation cover [7]. In this study, vegetation cover, as represented by a Normalized Difference Vegetation Index (NDVI), also was a strong predictor of SOC in the fitted model. This was consistent with Kunkel et al. [31] who reported that NDVI-predicted soil carbon distribution in semi-arid montane ecosystem.

The characteristics of the climatic and topographic variables according to the SOC map developed for the entire study area and the different SOC map quantities (low, medium, high) determined patterns among the environmental and land use variables. Temperature (AMT) influenced SOC by increasing the quantity at lower temperatures compared to the average SOC stock for the entire province. Some models and observations suggest that high latitude forests and grasslands may behave as a C source in response to increased decomposition of soil organic matter resulting from temperature increases [32,33]. Thus, the temperature sensitivity of decomposing organic matter in soil partly determines how much carbon will be transferred to the atmosphere because of global warming. This is consistent with Peri et al. [29] who reported that soil respiration rates were correlated strongly to air and soil temperatures by evaluating seasonal dynamics in contrasting grasslands across gradients of climate (rainfall), long term grazing intensity (moderate and high stocking rates over the last 80 years) and land uses (silvopastoral system, primary forest and grassland) in Southern Patagonia. The interaction of climate change with C pools in high latitude ecosystems may be particularly important because climate change is expected to be greatest at high latitudes. For instance, in Southern Patagonia (Santa Cruz and Tierra del Fuego provinces), mean maximum annual temperature is predicted to increase by 2–3 °C by 2080 between 46° and 52°30′ SL [34], and this will have significant effects on Patagonian ecosystems.

Rainfall and other related rainfall variables (seasonality) also influenced SOC by increasing the quantity with precipitation. It has been demonstrated that increased variability in rainfall and soil water content significantly affected SOC in this grassland [35]. In the semiarid temperate steppe in northern China soil water availability was more important than temperature in regulating soil and microbial respiratory processes, microbial biomass and their responses to climate change [36]. The strong and direct relationship between rainfall and SOC may be related to the ANPP and mean soil water content. In the present work net primary productivity (NPP) influenced positively SOC quantity compared to the average for the province. There is evidence that ecosystem C inputs from ANPP can be directly affected by altered rainfall variability, independent of precipitation quantity [35]. Thus, precipitation constrains plant production and decomposition in arid ecosystems, with a greater response of plant production relative to decomposition [37]. However, it is probable that SOC is controlled by the complex interaction of environmental and biotic factors. At the regional scale, in this study, patterns of SOC were positively associated with mean annual precipitation and negatively correlated with mean annual temperature across a diverse range of soils and vegetation types. This has been also documented in grasslands in North America [38].

We found that SOC values were low in high altitude mountain environments (SOC values generally increasing below 460 m.a.s.l.). This may be due to changes in climatic variables, precipitation, temperature and vegetation types along altitudinal gradients that influence in consequence the quality, quantity and turnover of soil organic matter. For example, Mulugeta and Itanna [39] determined that soil carbon stocks and turnovers in various vegetation types was directly proportional to the mean annual precipitation and inversely proportional to the mean annual temperature prevailing along the elevation gradient.

Vegetation types also showed differences in SOC, being higher in forest than in shrublands. The C efflux from soils in these rangeland ecosystems was higher than in subtropical savanna grasslands of southern Texas, USA [40] and C effluxes measured in Patagonian forest and woodland ecosystems [29]. This highlights the importance of environmental conditions (mainly soil water availability and temperature), input of organic residues, soil microbial biomass, and soil properties on the magnitude on soil respiration among different ecosystems.

In this regional study, increased animal stocking rates decreased SOC values, consistent with Peri [7]. We believe that low SOC at high stocking rates may be due to both low vegetation cover (or high bare soil cover) and low ANPP. Also, Peri et al. [29] reported that litter cover, litter depth, and soil carbon concentration (C %) in the uppermost soil layers, in both the dry and humid Magellanic steppe areas, were lower under heavy long term stocking rates than sites under moderate grazing intensity. Bahn et al. [41] indicated that the degree to which soil CO₂ efflux is coupled to soil C content may be largely determined by the reductions of supply by removal of aboveground biomass through grazing. Grassland ecosystems with high soil organic matter may promote organic matter decomposition (microbial activity) by continuous addition of litter and root turnover, thereby increasing soil respiration rates.

Finally, the development of sustainable land management practices for Patagonian Rangeland could benefit from a paleo-ecological perspective. The results and underlying data presented here are heavily focused on the importance of domestic animal grazing as a key driver of SOC over the last 100 years, when sheep-farming by European settlers began. However, an applied palaeoecology perspective could address specific ecological and environmental questions highly relevant to nature conservation, allowing us to assess the naturalness, rarity and fragility of Patagonia’s current ecosystems. This in turn could provide a sound basis prioritizing conservation and ecosystem restoration investments [42]. For example, palaeoecological studies indicated that few vegetation types in western Europe are natural [43]. In northern Patagonia, Schäbitz [44] reported changes in vegetation in late Holocene due to a drying climate that favored semi-arid vegetation. In Santa Cruz province (Southern Patagonia), Horta et al. [45] showed also vegetation change from grass steppe to shrub steppe during the Middle Holocene, most likely related to climate changes and the subsequent consequences for human occupational dynamics. Thus while it is clear that current environmental conditions, including patterns of human land use provide a useful baseline against which land management can be assessed, it is also clear that additional research focused on the paleoecology of the region could enrich this baseline and facilitate the development of effective strategies for sustainable land management.

5. Conclusions

Soil carbon storage across sites was influenced by a large number of interacting variables the most important of which were climatic conditions and plant productivity. Best practice ecosystem management (grazing) can increase net carbon storage in grasslands, shrublands, wetlands, and forests, but economic incentives to maintain or increase soil C stocks are needed. Understanding the causes of variation and mapping of SOC in Patagonia is a first step which allows for assessment of the sustainability of land management at the local scale, which can help to not only increase resilience of rangeland regionally in Patagonia, but also help address issues at the global scale, such as climate change, desertification, and biodiversity conservation.

Better understanding the results of existing long-term studies such as soil carbon stocks, and also realizing the existence of permanent plots in Patagonia for future research (e.g., PEBANPA network) will hopefully contribute to solve regional ecological and socio-economic challenges in the sustainable use of our native ecosystems. However, research and administrative institutions, and farmers must cooperate and have a sustained commitment to finance and maintaining these large unique long-term plots and research platforms.