Mapping Changes in Carbon Storage and Productivity Services Provided by Riparian Ecosystems of Semi-Arid Environments in Northwestern Mexico

We analyze the importance of riparian ecosystems (RE) as critical areas for carbon storage and productivity in semi-arid regions of Northwest Mexico. We calculated the carbon storage by land cover and compared temporal trends of basal productivity (MODIS) and pre-monsoon productivity (Landsat) of RE, to other land cover types. We used land cover maps generated previously for the region (years 1993, 2002, and 2011), assigning values of carbon stored in aerial and root biomass, as well as organic carbon stored in the soil. To estimate productivity (proxy), time series were generated using the Normalized Difference Vegetation Index (NDVI) values of Landsat 4–5 TM and MODIS for each land cover type. We found that RE stores 93,147 tC/ha, about 1.5 times the estimated storage for oak forest (65,048 tC/ha). Productivity of RE was similar to highly productive land cover types, such as agriculture and oak forest, and higher than in the rest of the ecosystems of the region. We also found that changes from RE to agriculture and cultivated grasslands represented a decrease in productivity (p < 0.001). Finally, we report a gradual decrease in basal productivity (p = 0.0151) and pre-monsoon productivity (p = 0.031) in the RE. These results help us understand that changes in land use, intensive use of water, and climate can influence the ecosystem services of productivity and carbon storage offered by RE in semi-arid areas.


Introduction
Riparian ecosystems (RE) are key elements for the function and assemblage of habitats in arid and semi-arid environments [1][2][3].Moreover, numerous studies have found that many species (not only riparian species) in these areas depend on water, nutrients resources present in riparian habitats [4,5].RE resource richness is a consequence of the multiple ecological and hydrological processes such as high volumes of surface runoff, flooding regime, transport and deposition of organic matter, soil development, moisture retention, and biomass production [1][2][3]6].
Arid and semi-arid RE are often cataloged as large carbon reserves [7,8] due to their broad diversity and abundance of plant species, which can produce dense stands of individuals and large volumes of biomass, when compared to other arid ecosystems [9].
From an ecosystem services (ES) perspective, it is believed that riparian ecosystems of arid-semiarid zones (REAZ) provide a great variety and quantity of ES compared to the other ecosystems present in the same area [10].In addition to materials and goods obtained, the REAZ offer multiple support, regulation services, acting as filters of agricultural pollutants, water purifiers, erosion control, carbon capture and storage, flood control, among others, and cultural services, which contribute to the development of functional coupled human and natural systems [11][12][13][14][15].Despite the considerable provision of ES offered by the REAZ, there is a general lack of knowledge regarding key processes that affect the provision of these services.
In Mexico, only a limited quantity of ES has been subject to management or regulation in the REAZ.Due to their importance for economic-social development, ES related to water production and fertile land, are the most widely studied and regulated in REAZ in the country.However, services regarding primary productivity (support services) and carbon storage (regulation services) have been undervalued and scarcely studied.
Primary productivity is tightly related to ecosystem functioning and development, since it regulates material and energy exchange trough photosynthetic-respiration activity [16].It is also considered a support service, since it is the base for the "production" of other services.In this sense, when photosynthetically active radiation is converted into biomass and consumed by the livestock, or used as a fuel or as construction materials, productivity generates provision services [17].Energy that is transformed into biomass also provides the ES of carbon storage, promoting the regulation of climate [17].
The evidence of a gradual and significant increase in temperature over the last 100 years [18], derived from increased CO 2 levels in the atmosphere, has provoked the concern of the scientific and political community worldwide.In the face of this situation, decision-makers require accurate information regarding the current status and spatial distribution of sources, sinks and storage of carbon at global, regional, and even local scales [19].Interest has generally been focused on tropical and temperate forests, since those ecosystems represent the largest reservoirs of carbon in the world.However, it has been observed that carbon storage in the REAZ can match the levels of other ecosystems, such as the pine, oak [20,21], and tropical dry forests [22].
Due to the limited evidence provided by the scientific community [23], there is a general lack of information regarding productive capacity and carbon storage in the REAZ.Therefore, in the case of Mexico, there are no legal mechanisms to protect these ecosystems and ensure the preservation of their ES.
Previous studies have calculated that arid environments of Northwestern Mexico have undergone extensive conversions, primarily due to agricultural use [24].Therefore, it is imperative to improve our understanding regarding their capacity as and ES providers to be able to develop management guidelines according to their importance and prominence in arid landscapes.
The objective of the present study is to (1) provide evidence of the importance of the REAZ as critical areas for productivity and carbon storage services in a regional context, and (2) assess how these services can be affected by land and water use in semi-arid regions.For this, we developed the following objectives: (1) to estimate the carbon content in three different stocks (aerial biomass, root biomass, and in the soil) per land cover type present in an experimental REAZ; (2) determine the spatial distribution of carbon in three different years (1993, 2002, and 2011) using land cover maps of the study area [25]; (3) analyze NDVI values per land cover type in order to compare proxies of productivity between RE, and other ecosystems found in the same region; (4) quantify and analyze changes in carbon storage associated to changes in land use; and, finally (5) analyze the temporal trends in productivity per land cover type, as well as the behavior of these trends in RE.The present study will contribute to increase our knowledge regarding carbon dynamics, and vegetation function, in key ecosystems of semi-arid environments.

Study Area
We selected the sub-basin "El Cajoncito-Arroyo del Carrizo" as our study area, which includes the San Miguel (SMR) and Zanjón (ZR) river regions [26].It is located in North-Central Sonora (Mexico), approximately between latitudes 30 • 46 and 28 • 53 north and longitudes 110 • 33 and 110 • 45 west.The sub-basin approximate area is 9437 km 2 , and it covers 32% of the Sonora River basin.The El Cajoncito-Arroyo del Carrizo sub-basin provides most of the water for Hermosillo (the capital city of Sonora), since it feeds the Abelardo L. Rodríguez dam and multiple groundwater wells used to supply the city.Moreover, these river systems recharge aquifers that provide water for agricultural, public-urban, and fishing activities in at least eight municipalities in the State [27].
The mean elevation for the Zanjón River (ZR) region is 656 m, ranging from 1734 m in the Sierra de Cucurpe to 280 m in the lower reaches, where it converges with the San Miguel River (SMR).Both watersheds are subject to a bimodal precipitation pattern, with summer and winter monsoons.Precipitation in the ZR is the lowest in the entire basin of the Sonora River, with records of between 340 and 384 mm/year, occurring mostly during the summer monsoon.The annual mean temperature ranges between 21 • C and 22 • C [27].The dominant land cover types in this region are desert scrub, cultivated and induced grassland (mainly buffel-grass (Cenchrus ciliaris)), subtropical/succulent scrub, and agriculture (annual and perennial) [25].
The San Miguel River region has a mean elevation of 948 m, with a maximum of 2439 m in the Sierra Azul and a minimum of 225 m at the Abelardo L. Rodríguez dam.Annual accumulated precipitation ranges from 423 mm to 575 mm (most of it occurring during the summer monsoon), with an average annual temperature ranging between 19 • C and 22 • C [28].The dominant land cover types in this region are subtropical and desert scrub, oak/oak-pine forest, as well as natural/native grassland and riparian mesquite [25].
The presence of riparian vegetation is widespread in some sections of the SMR watercourse, mainly between Cucurpe and San Miguel de Horcasitas (Figure 1).Agricultural areas are found near main settlements like Cucurpe, Opodepe, Rayón, and San Miguel de Horcasitas.The most extensive agricultural fields are located at the riverbanks towards North Rayón.It is also worth noting that, despite the fact that topography makes it difficult to establish prairies of buffel-grass in the SMR sub-basin, extensive livestock production is also a significant activity.El Cajoncito-Arroyo del Carrizo sub-basin provides most of the water for Hermosillo (the capital city of Sonora), since it feeds the Abelardo L. Rodríguez dam and multiple groundwater wells used to supply the city.Moreover, these river systems recharge aquifers that provide water for agricultural, public-urban, and fishing activities in at least eight municipalities in the State [27].
The mean elevation for the Zanjón River (ZR) region is 656 m, ranging from 1734 m in the Sierra de Cucurpe to 280 m in the lower reaches, where it converges with the San Miguel River (SMR).Both watersheds are subject to a bimodal precipitation pattern, with summer and winter monsoons.Precipitation in the ZR is the lowest in the entire basin of the Sonora River, with records of between 340 and 384 mm/year, occurring mostly during the summer monsoon.The annual mean temperature ranges between 21 °C and 22 °C [27].The dominant land cover types in this region are desert scrub, cultivated and induced grassland (mainly buffel-grass (Cenchrus ciliaris)), subtropical/succulent scrub, and agriculture (annual and perennial) [25].
The San Miguel River region has a mean elevation of 948 m, with a maximum of 2439 m in the Sierra Azul and a minimum of 225 m at the Abelardo L. Rodríguez dam.Annual accumulated precipitation ranges from 423 mm to 575 mm (most of it occurring during the summer monsoon), with an average annual temperature ranging between 19 °C and 22 °C [28].The dominant land cover types in this region are subtropical and desert scrub, oak/oak-pine forest, as well as natural/native grassland and riparian mesquite [25].
The presence of riparian vegetation is widespread in some sections of the SMR watercourse, mainly between Cucurpe and San Miguel de Horcasitas (Figure 1).Agricultural areas are found near main settlements like Cucurpe, Opodepe, Rayón, and San Miguel de Horcasitas.The most extensive agricultural fields are located at the riverbanks towards North Rayón.It is also worth noting that, despite the fact that topography makes it difficult to establish prairies of buffel-grass in the SMR subbasin, extensive livestock production is also a significant activity.Social and economic development in both regions has been closely related to the use and exploitation of the RE.Provision of ground and surface water, as well as use of the fertile land on the floodplains, are the most important provision ES, and have allowed the development of agricultural and livestock production activities.The intensive use of the resources described above exerts a constant pressure on the RE, compromising its structural and functional integrity, as well as their capacity to provide ES.

Land Use and Cover Classification
This study used land cover maps generated by Méndez-Estrella et al. (2016) [25] (Figure 2), which were generated through supervised classification, using the CART (Classification and Regression Trees) model and the scheme of classes proposed by Anderson [29], considering a combination of physiognomic and life form criteria provided by diverse sources [30,31].and livestock production activities.The intensive use of the resources described above exerts a constant pressure on the RE, compromising its structural and functional integrity, as well as their capacity to provide ES.

Land Use and Cover Classification
This study used land cover maps generated by Méndez-Estrella et al. ( 2016) [25] (Figure 2), which were generated through supervised classification, using the CART (Classification and Regression Trees) model and the scheme of classes proposed by Anderson [29], considering a combination of physiognomic and life form criteria provided by diverse sources [30,31].Sampling (for classification training) for each of the land cover classes was made on the field, and also using orthorectified imagery with high spatial resolution.The training sample size ranged between 60 and 150 sites per land cover class defined in the scheme.
Once the training locations where collected, Landsat TM scenery (from Path 35-Row 39 and Path 35-Row 40) was selected for the years 1999, 2002, and 2011, assuring atmospheric and terrain correction [25].A total of two images per year (on both scene locations) were collected from the United States Geological Survey [25].
Using the TM scenes and a digital elevation model (DEM) obtained also from the United States Geological Survey, we derived an array of variables used to create the thematic classification maps [25].
Global precisions of the classifications were greater than 79%, the precision values of the user and producer are in the range of 52% to 100% and the kappa coefficient is greater than 75% [25].The classifications, therefore, comply with the acceptable limits of accuracy [32].

Carbon Storage by Land Cover Type
To identify critical areas for carbon storage we used land-use and land-cover maps of the area of interest, assigning each of the classes an average value of carbon, calculated from information collected in the field, or from previous studies [33].In order to obtain values for carbon storage per Sampling (for classification training) for each of the land cover classes was made on the field, and also using orthorectified imagery with high spatial resolution.The training sample size ranged between 60 and 150 sites per land cover class defined in the scheme.
Once the training locations where collected, Landsat TM scenery (from Path 35-Row 39 and Path 35-Row 40) was selected for the years 1999, 2002, and 2011, assuring atmospheric and terrain correction [25].A total of two images per year (on both scene locations) were collected from the United States Geological Survey [25].
Using the TM scenes and a digital elevation model (DEM) obtained also from the United States Geological Survey, we derived an array of variables used to create the thematic classification maps [25].
Global precisions of the classifications were greater than 79%, the precision values of the user and producer are in the range of 52% to 100% and the kappa coefficient is greater than 75% [25].The classifications, therefore, comply with the acceptable limits of accuracy [32].

Carbon Storage by Land Cover Type
To identify critical areas for carbon storage we used land-use and land-cover maps of the area of interest, assigning each of the classes an average value of carbon, calculated from information collected in the field, or from previous studies [33].In order to obtain values for carbon storage per land cover class, we used a combination of approaches ranging from the use of coefficients obtained from the literature to direct field sampling.
Our first approach to derive carbon storage values was through estimation obtained in previous studies, regarding biomass and carbon allocation in arid lands vegetation types [34][35][36][37][38][39][40][41][42][43].The second step was to complement the previous using datasets form the National Forestry Inventory (NFI), provided by the National Forestry Commission (CONAFOR, for its Spanish acronym).Specifically, we used dendrometric measurements of diameter at breast height (DBH), which is collected for every tree with DBH greater than 7.5 cm at each sample plot (conglomerate), conducted on land cover types present in our study site (Appendix A).Results obtained from the previous measurements were fed to specific allometric equations (Appendix B), to obtain estimates for carbon storage on this cover types.Finally, direct field sampling was carried out for riparian vegetation, where NFI estimation or data from the literature was not available.
Values of carbon stocks regarding (1) aerial biomass, (2) root biomass, and (3) organic carbon in the soil, were calculated and fed into the model through the following procedure: (a) Aerial biomass per land cover type was estimated using the NFI information for the study area (Appendix A).Each NFI sampling unit is called a "conglomerate" and, in arid zones, has an area of 1600 m 2 [44].The diameter at breast height (DBH) data was used to calculate biomass per species in each "conglomerate", applying the allometric equations reported in the literature [45,46].Each "conglomerate" sample size was between 7 and 39 trees (depending on land cover type), and the total biomass by conglomerate was estimated in ton/ha.(b) Aerial biomass of cultivated/induced grassland (not included in the NFI sampling sites), was obtained from regional studies regarding this variable [34], using similar sites to the land cover type of interest.(c) Aerial biomass of perennial and annual agriculture was obtained from the "Servicio de Información Agroalimentaria y Pesquera" database [35].(d) For the riparian ecosystems of the SMR (Riparian vegetation), which differs from Riparian mesquite (composed mainly by Prosopis, Cercidium, Acacia) due to the floristic composition (species of the genera Populus, Salix, Bacharis, etc.), no information was found in the NFI or in the literature.Therefore, eight field-sampling sites were established (with a sample size of 28 trees per site) using the NFI protocol, in which the DBH was measured (using only trees that had over 7.5 DBH) in order to enable the use of species-specific allometric equations (Appendix B).(e) For accumulated root biomass, we used previous works, which estimate shoot:root biomass ratios per vegetation type [36][37][38][39][41][42][43].(f) Finally, for the organic carbon stored in the soil, we obtained the reported carbon values for each soil type [36,40].
The information generated through allometric equations or reported in the literature is expressed in biomass per hectare, before being finally converted to ton/ha of carbon.This conversion was performed by multiplying the biomass by the standard factor of conversion to carbon, which is 0.47 according to the IPCC (2006) for most types of vegetal cover [47].In the case of agriculture, different conversion values were used according to the particular type of crop present, as described in Carvajal et al. [41].Appendix B presents the references for each calculation, or the bibliographic origins of the data that were used for the spatial modeling of critical areas for carbon storage.

Spatiotemporal Changes in Carbon Storage
Critical areas for carbon storage were obtained by applying the carbon storage and sequestration models, included in the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST 3.1.0)[20].
Once estimated, carbon stocks were assigned to each land cover type, assuming that spatiotemporal changes among the different vegetation types, (e.g., riparian vegetation transformed to agriculture) translates into carbon storage changes.Therefore, an ecosystems service model (InVEST 3.1.0)was applied to each of the classified thematic maps previously obtained (years 1993, 2002, and 2011), in order to estimate carbon gains and losses in our study site.For each year classified, a carbon storage map was obtained, with a spatial resolution of 30 m corresponding to that of Landsat TM 4-5 images.

Vegetation Indices as Proxies for Productivity
Vegetation indices are highly correlated to photosynthetic activity [16], and are often used to assess vegetation function.These indices can be derived from sensors mounted on satellite platforms, from where we can obtain images at specific temporal resolutions [16].We decided to use the MODIS and Landsat images to: (1) corroborate possible differences and trends in productivity with two sources of information; (2) complement the spatial resolution (with both MODIS and Landsat) of our analysis; and (3) to compare our results to those obtained in previous studies, with respect to the productivity trends found in the wetlands of the region [48,49].
The Normalized Difference Vegetation Index (NDVI) is frequently used to analyze the photosynthetic activity of vegetation, providing approximate information about the primary productivity, biomass, and the index of foliar area for extensive areas of the landscape [16,50].This index has also been used to observe changes in vegetation photosynthesis [51].For photosynthetically active vegetation, high reflectance in the near infrared and low reflectance in the red wavelength are characteristics utilized for the construction of the NDVI [52].This index is calculated by applying the following equation to each of the pixels of an image obtained through remote sensing: where NIR stands for near infrared reflectance and RED stands for red visible reflectance of the photosynthetically active radiation (PAR) spectral region.In this study, the parameters obtained based on NDVI were: (1) "basal productivity", using information from the MODIS sensor (MOD13Q1), which represents the proxy of vegetation productivity excluding the growing season [53], and (2) "pre-monsoon productivity", based on Landsat TM 4-5 satellite images (from 1988 to 2011), which represents the proxy of productivity of vegetation before the growing season.The previous analysis contributes to increase our understanding, regarding the importance riparian ecosystems, which could be key to maintain ecological processes at regional scales, even during drought periods.
By constructing a MODIS NDVI time series from 2001 to 2012, using the 16-day NDVI composite of the MODIS sensor (MOD13Q1), we obtained phenological parameters of the vegetation at the regional scale [53].This allowed us to understand key functions of the ecosystems analyzed, such as photosynthetic activity [16], and also helped us to analyze and predict the responses of the vegetation to changes in meteorological and management variables [54].In order to analyze long-term trends of the proxy represented by the basal productivity [53] of each of the coverage types in our study area, we excluded the productivity of the seasonally active vegetation (small integral) from the total productivity during the time of growth (large integral) [53].
An NDVI time series based on Landsat TM 4-5 satellite images (from 1988 to 2011) was also created.To accomplish this, we selected images taken prior to the summer monsoon for each year to analyze.The images were obtained from the Earth Explorer platform, managed by the United States Geological Survey [25] and chosen between March and April of each year.Two Landsat scenes were selected per year, due to the extent of our study area (Path 35-Row 39 and Path 35-Row 40).Once the images were obtained, information was extracted by coverage type that presented no changes over the period 1993-2011 [25].

Comparing Basal and Pre-Monsoon Productivity
Considering the hypothesis that the REAZ would be more productive than adjacent ecosystems, even in periods of the year where productivity is limited by the reduced availability of rainwater [1][2][3]6], it is important to determine annual differences in basal (MODIS 2001-2012) and pre-monsoon (Landsat 1993(Landsat -2011) ) productivity values between riparian vegetation and the surrounding land cover types.We used a Kruskal-Wallis test, to analyze general differences between land cover types and to compare land cover using Dunn's method [55].The data used was MODIS and Landsat NDVI values of the pixels corresponding to each land cover type.
To assess land use change impact on RE productivity, we compared NDVI values representing the growing season (large integral of MODIS) in RE which remained unchanged over the 2002-2011 period, versus those that underwent changes and were substituted by other land cover type.We focused primarily in changes between riparian vegetation and agriculture and cultivated/induced grassland and vice versa.For this analysis, we used the large integral of MODIS corresponding to the year 2011 only, in order to ensure that all of the transitions between land cover classes were complete.The data obtained was analyzed using a Mann-Whitney U rank sum test, to determine possible differences between productivity of preserved RE and sections of the river that have undergone changes.

Trends on Basal and Pre-Monsoon Productivity
Knowledge of trends in productivity over time provides us with an indicator of an ecosystem's condition facing climatic variations, and anthropogenic management [48,49,51].Previous studies have documented the progressive decline of productivity in the riparian ecosystems of the upper San Pedro River in Arizona [49], and in extensive areas of mangrove in the southeastern coast of the Gulf of California in Sonora, Mexico [48].Since our study is located between the two areas mentioned above, and to increase our understanding of ongoing processes on the region, we assess the temporal trends of the NDVI values per vegetation type, using a simple linear least squares regression analysis.Specifically, we analyze the proxies of a) productivity before the summer monsoon based on Landsat and b) basal productivity based on MODIS.

Temperature, Precipitation, Runoff and Water Depth Effects on Productivity Trends
To analyze the relationship between productivity trends for each land cover type and climatic variability, we obtained the temporal trends of (1) annual mean temperature, (2) mean temperature of spring-summer and autumn-winter, (3) annual accumulated precipitation, and (4) accumulated precipitation of spring-summer and autumn-winter.Climatic parameters for each region were obtained using daily records of three climatological stations from the National Meteorology Center [56], located within the ZR and five stations within the SMR [56].Then, the climatic parameters from each station were averaged and grouped according to the region (Table 1).A simple regression analysis was conducted for each parameter in of the regions (ZR and SMR).We also used data from the hydrological station "El Cajón", located on the course of the SMR [57] to observe annual accumulated runoff.To test correlation between climatological and hydrological variability and productivity trends (MODIS and Landsat) of riparian vegetation and riparian mesquite, a Pearson product moment correlation analysis was performed.
Finally, we used spatial information regarding the distribution of water depth, as derived by Méndez-Estrella et al. [25], and the Landsat TM derived NDVI pre-monsoon for years 2007, 2008, 2009, and 2010 to perform a non-linear regression analysis, to assess riparian vegetation productivity as a function of water depth.

Carbon Stocks and Allocation
Our findings show that the land cover types with greatest carbon storage in aerial biomass were riparian vegetation (33.5 TonC/ha), oak/oak-pine forest (11 TonC/ha), and riparian mesquite (4 TonC/ha).On the other hand, the land cover types storing the least carbon in their aerial tissue were: natural grassland, desert scrub, and cultivated grassland (prairies of buffel-grass) (with 2.85, 2.39 and 0.74 TonC/ha, respectively) (Table 2).Results from our analysis suggests higher carbon storage for riparian and riparian-mesquite root systems, than in the rest of the land cover types analyzed.In contrast, land cover types such as annual agriculture and cultivated grassland exhibit a lower quantity of accumulated carbon in the roots (Table 2).In subtropical scrub cover it has been reported that root systems have low storage carbon capacity [39].This is due to the fact that this vegetation type is distributed mainly on slopes and hills of medium elevation, where soils are shallow [44].Other vegetation types such as mesquite woodland and desert scrub are associated to deeper soils [44], allowing them to develop greater root biomass systems.Other studies have reported that root systems are the structural components with the highest proportion in the total biomass of mesquite woodland and desert scrub [9,58], which acts to increase the estimations of carbon stored below the soil, even compared to other coverage with greater aerial biomass (e.g., forest) (Table 2).
According to our analysis, oak/oak-pine forest, riparian vegetation, and subtropical scrub present high values for soil organic carbon, while grasslands, mesquite woodland, and desert scrub present lower values.It should be noted that the quantity of carbon stored in the soil is up to three times greater than the carbon stock present in the aerial biomass of the vegetation [59].
Comprising the three carbon stocks mentioned before, riparian vegetation, oak/oak-pine forest and the riparian mesquite present the highest stocks of carbon per unit area in the region.Nevertheless, considering the coverage of each vegetation type, subtropical scrub, forest (oak/oak-pine), mesquite woodland, and even desert scrub constitute a greater total storage of carbon compared to the riparian vegetation when considering the entire study area (Table 3).

Spatiotemporal Changes in Carbon Storage
Our results describe how spatial and temporal changes in the land cover reflect changes in carbon storage.Figure 4 shows the effect of the spatial change on the carbon stock, derived from exchanges among the landscape components in the period 1993-2011.
For the ZR, the reduction in carbon storage is related to agricultural area loss, substituted by classes such as desert scrub, mesquite woodland, and cultivated grassland, this change represents a decrease of between 4.491 and 20.401 TonC/ha.It should be noted that the agricultural biomass (mainly comprising annual crops) is constantly harvested, unlike the natural vegetation where carbon, although present in lower quantity, is stored for longer time periods.Likewise, the establishment of new prairies of buffel-grass in sites of natural vegetation, such as desert scrub and mesquite woodland represents a reduction in carbon storage of 6.175 and 15.91 TonC/ha, respectively.
On the other hand, changes related to riparian mesquite reestablishment as a consequence of agricultural field abandonment, lead to an increase of 14.09 TonC/ha in stored carbon (Figure 4).Establishment of perennial crops (such as grape, orange, and walnut) in sites of induced grasslands, mesquite woodland, or desert scrub, represents an increase in carbon storage of between 10.186 and 26.096 TonC/ha.The loss of cultivated grasslands and their substitution with mesquite and desert scrub represents a possible increase of 6.175 and 15.91 TonC/ha, respectively (Table 2).
Changes arising from the loss of land cover of anthropogenic origin and its substitution with "natural" vegetation have already been reported in another study [60], in which a relationship was found between the abandonment of agricultural fields on the peri-urban areas and the re-assignation of water from agricultural to exclusively urban use.

Spatiotemporal Changes in Carbon Storage
Our results describe how spatial and temporal changes in the land cover reflect changes in carbon storage.Figure 4 shows the effect of the spatial change on the carbon stock, derived from exchanges among the landscape components in the period 1993-2011.
For the ZR, the reduction in carbon storage is related to agricultural area loss, substituted by classes such as desert scrub, mesquite woodland, and cultivated grassland, this change represents a decrease of between 4.491 and 20.401 TonC/ha.It should be noted that the agricultural biomass (mainly comprising annual crops) is constantly harvested, unlike the natural vegetation where carbon, although present in lower quantity, is stored for longer time periods.Likewise, the establishment of new prairies of buffel-grass in sites of natural vegetation, such as desert scrub and mesquite woodland represents a reduction in carbon storage of 6.175 and 15.91 TonC/ha, respectively.
On the other hand, changes related to riparian mesquite reestablishment as a consequence of agricultural field abandonment, lead to an increase of 14.09 TonC/ha in stored carbon (Figure 4).Establishment of perennial crops (such as grape, orange, and walnut) in sites of induced grasslands, mesquite woodland, or desert scrub, represents an increase in carbon storage of between 10.186 and 26.096 TonC/ha.The loss of cultivated grasslands and their substitution with mesquite and desert scrub represents a possible increase of 6.175 and 15.91 TonC/ha, respectively (Table 2).
Changes arising from the loss of land cover of anthropogenic origin and its substitution with "natural" vegetation have already been reported in another study [60], in which a relationship was found between the abandonment of agricultural fields on the peri-urban areas and the re-assignation of water from agricultural to exclusively urban use.
Areas presenting the greatest reduction in carbon storage were those where riparian vegetation was lost, in a conversion to bare soil or to seasonal agriculture (decrease of between 44.197 and 49.992 TonC/ha).In contrast, areas with possible gains in carbon of up to 38.797 TonC/ha were also detected (Table 2 and Figures 2 and 3).Overall, change from riparian vegetation to bare soil and agriculture was 1987 and 220 ha, respectively.The previous represents a total carbon loss of 109,095 tons for our study area.Areas presenting the greatest reduction in carbon storage were those where riparian vegetation was lost, in a conversion to bare soil or to seasonal agriculture (decrease of between 44.197 and 49.992 TonC/ha).In contrast, areas with possible gains in carbon of up to 38.797 TonC/ha were also detected (Table 2 and Figures 2 and 3).Overall, change from riparian vegetation to bare soil and agriculture was 1987 and 220 ha, respectively.The previous represents a total carbon loss of 109,095 tons for our study area.

Differences in Productivity per Land Cover Type
We analyzed differences between productivity proxies for each land cover type to understand the relative importance of each ecosystem in function of its contribution to regional productivity.The results of multiple comparisons by treatment pairs procedure (using the Dunn's method), show that basal productivity (large integral minus small integral) of the riparian vegetation is equal to, or greater than, that of oak/oak-pine forest for years 2002, 2006, 2009, and 2011, and agriculture for 2001, 2002, 2004, 2009, 2010, and 2011 (Table 4, Figure 5a).

Differences in Productivity per Land Cover Type
We analyzed differences between productivity proxies for each land cover type to understand the relative importance of each ecosystem in function of its contribution to regional productivity.The results of multiple comparisons by treatment pairs procedure (using the Dunn's method), show that basal productivity (large integral minus small integral) of the riparian vegetation is equal to, or greater than, that of oak/oak-pine forest for years 2002, 2006, 2009, and 2011, and agriculture for 2001,2002,2004,2009,2010, and 2011 (Table 4, Figure 5a).For riparian mesquite, annual basal productivity during the period of analysis was generally lower than that of agriculture and riparian.When compared to oak/oak-pine forest, the riparian mesquite also exhibited lower basal productivity except during years 2003, 2006, 2009, 2011, and 2012 (no statistical differences during this years) (Table 4, Figure 5a).Annual basal productivity of riparian mesquite were consistently higher than those of natural grassland, mesquite woodland, desert scrub, and cultivated grassland (Table 4 with p < 0.05 and Q values > 4, Figure 5a) and the subtropical scrub, except in the latter for 2001, 2008, and 2010 (Table 4, Figure 5a).
Regarding the Landsat time series of pre-monsoon productivity, the results support those observed using MODIS-NDVI time series, since the riparian vegetation is among the three most productive types of land cover, together with agriculture and the oak/oak-pine forest (Figure 5b).Likewise, the riparian mesquite presents values (proxy of productivity) greater than those land cover types associated to more arid conditions, such as desert scrub, cultivated/induced grassland, and mesquite woodland (Figures 2 and 5b).It is important to observe that in years in which lower autumn-winter precipitation is recorded (1988, 1989, 1996, 1999, 2000, 2006, and 2009), there are increased differences in the proxies of productivity of the riparian vegetation and agriculture, compared to the other land cover types (Figure 5b).This is due to the fact that these two land cover types are less dependent on precipitation than the others present in the study area.

Trends in productivity
The NDVI time series, using MODIS (basal productivity 2001-2012) and Landsat (pre-monsoon productivity 1993-2011), shows a general negative trend in all land cover types in our study area (Figure 6).The MODIS data indicate that the observed trend is significant for six of the nine land cover types analyzed, with the exception of the natural grassland, cultivated grassland, and riparian mesquite (Figure 6a).However, the results show that the regression based on the NDVI of Landsat was only statistically significant for agriculture and riparian vegetation, suggesting a gradual decrease in NDVI values (Figure 6b).The greater temporal resolution of the observations made with MODIS provides us continuous information pertaining to the basal productivity of the different ecosystems, while Landsat represents, in this case, the productivity specific to a moment prior to the growth season which, in this region, generally coincides with the summer monsoon [61].

Trends in productivity
The NDVI time series, using MODIS (basal productivity 2001-2012) and Landsat (pre-monsoon productivity 1993-2011), shows a general negative trend in all land cover types in our study area (Figure 6).The MODIS data indicate that the observed trend is significant for six of the nine land cover types analyzed, with the exception of the natural grassland, cultivated grassland, and riparian mesquite (Figure 6a).However, the results show that the regression based on the NDVI of Landsat was only statistically significant for agriculture and riparian vegetation, suggesting a gradual decrease in NDVI values (Figure 6b).The greater temporal resolution of the observations made with MODIS provides us continuous information pertaining to the basal productivity of the different ecosystems, while Landsat represents, in this case, the productivity specific to a moment prior to the growth season which, in this region, generally coincides with the summer monsoon [61].

Decrease in Productivity and Its Possible Causes
In this study, we found a decreasing trend in annual accumulated precipitation across the entire study area, but this was only statistically significant in the ZR region (simple linear regression, p = 0.0178).Likewise, the increase in annual mean temperature was only statistically significant for the ZR region (simple linear regression, p = 0.0046) (Figure 7b,f).

Decrease in Productivity and Its Possible Causes
In this study, we found a decreasing trend in annual accumulated precipitation across the entire study area, but this was only statistically significant in the ZR region (simple linear regression, p = 0.0178).Likewise, the increase in annual mean temperature was only statistically significant for the ZR region (simple linear regression, p = 0.0046) (Figure 7b,f).
For the ZR region, which presents greater aridity than the SMR, the decrease in autumn-winter precipitation and the increase in temperature over the same period of the year were significant (simple linear regression, p < 0.05) (Figure 7d,h).This increase in air temperature, and precipitation reduction (Figure 7), have generated drought conditions repeatedly reported for the region (in the period between 1995 and 2011) [56,62] and could affect the productivity of all land cover types present in the region.On the other hand, while the reduction in precipitation was not significant for the SMR region (Figure 7a,c), the high exploitation of ground and surface water for sustaining agricultural and livestock production activities [28] could explain the reduction of surface runoff in the region [62], which could also be related to the decrease in productivity of the RE [49].
ISPRS Int.J. Geo-Inf.2017, 6, 298 15 of 26 For the ZR region, which presents greater aridity than the SMR, the decrease in autumn-winter precipitation and the increase in temperature over the same period of the year were significant (simple linear regression, p < 0.05) (Figure 7d,h).This increase in air temperature, and precipitation reduction (Figure 7), have generated drought conditions repeatedly reported for the region (in the period between 1995 and 2011) [56,62] and could affect the productivity of all land cover types present in the region.On the other hand, while the reduction in precipitation was not significant for the SMR region (Figure 7a,c), the high exploitation of ground and surface water for sustaining agricultural and livestock production activities [28] could explain the reduction of surface runoff in the region [62], which could also be related to the decrease in productivity of the RE [49].Despite the trends towards decreasing precipitation and increasing temperatures (mostly for the ZR), a significant correlation was only found between the decrease in the basal productivity of the riparian vegetation distributed in the SMR and an increase in the annual average temperature (Simple Despite the trends towards decreasing precipitation and increasing temperatures (mostly for the ZR), a significant correlation was only found between the decrease in the basal productivity of the riparian vegetation distributed in the SMR and an increase in the annual average temperature (Simple Linear Regression, p = 0.0468, Table 5).Likewise, for the riparian mesquite in the ZR region, there is a significant correlation between the decrease in basal productivity and the increase in average temperature of the autumn-winter period (simple linear regression, p = 0.0362, Table 6).
Table 5. Pearson correlation analysis between the proxy of the pre-monsoon and basal productivity of the riparian vegetation, accounting for climatic and surface hydrological parameters of the SMR.Agriculture, urban development, and livestock production are activities exerting great pressure on RE, since they use large quantities of ground and surface water.Ground water is obtained from wells and "norias" (shallow wells that do not reach the aquifer), while much of the surface water is obtained from water reservoirs known regionally as "represos", or dams.

Proxy of
Analyzing the spatial relationship between the static levels of water and the proxy of the pre-monsoon productivity (derived from Landsat images) of the riparian vegetation for different years of the study period, we observe that NDVI values exhibits an exponential decline with increasing water table depth (Figure 8).From this result, we can also observe that, after exceeding a threshold in the water table of approximately 10 to 12 m, photosynthetic activity declines considerably (Figure 8).In the SMR and ZR regions, close to 52% of the RE are located in areas with a water depth of between 1 and 12 m (according to previous studies) [25], the rest of the area is found at depths above this threshold (these sites might undergo functional and structural changes).

Productivity in Sites with Changes in Land Cover
The main changes undergone by riparian vegetation and riparian mesquite land cover types in the ZR and SMR are their conversion to agricultural areas or grasslands [25].In addition to a reduction in carbon storage (Table 2, Figure 4), substitution of riparian vegetation and riparian mesquite for other land cover types often results in a decrease in productivity.
Based on the results of a Mann-Whitney U rank sum test, in places where riparian vegetation has changed to agriculture, we were able to observe a significant decrease in productivity when the conversion occurs (p < 0.001, Mann-Whitney U statistic = 8480.5).Nevertheless, riparian areas that changed to agriculture have high NDVI values, depending on the crop type and the time of the year (perennial and annual crops during their productive cycle).Most of time, the opening of new agricultural fields implies the removal of natural vegetation and its substitution by annual or perennial crops, which for this semi-arid region of Mexico represents a decrease in productivity.
Abandonment of agricultural fields has been documented in this region [60]; the fields had been occupied mainly by riparian mesquite [25].Contrary to our expectations, this change also results in a decreasing productivity (p = 0.006, Mann-Whitney U statistic = 5438), which can be explained by the loss of vegetal coverage in the early stages of the replacement of agriculture by riparian vegetation (Figure 9b).Additionally, it is worth noting that new mesquite woodlands often present low values of foliar area.
Finally, changes from riparian vegetation to cultivated grasslands and vice versa also represent drastic changes in NDVI (Figure 9c,d).Substitution of RE for prairies of buffel-grass has been a common practice, above all in the region of ZR [63].In sites where these changes have been presented, a significant reduction in the values of NDVI is observed (p < 0.001, Mann-Whitney U statistic = 2379; Figure 9c).The change from cultivated grassland to riparian vegetation has also been reported, and has been explained from the perspective of the resilience of RE, since "it appears the conditions for the growth and maintenance of the prairies are not adequate to definitively substitute riparian vegetation, and for this reason, the land cover reverts to its previous state" [25].Recovery of the riparian vegetation is reflected in an apparent increase in productivity (p < 0.001, Mann-Whitney U statistic = 15,955; Figure 9d).

Productivity in Sites with Changes in Land Cover
The main changes undergone by riparian vegetation and riparian mesquite land cover types in the ZR and SMR are their conversion to agricultural areas or grasslands [25].In addition to a reduction in carbon storage (Table 2, Figure 4), substitution of riparian vegetation and riparian mesquite for other land cover types often results in a decrease in productivity.
Based on the results of a Mann-Whitney U rank sum test, in places where riparian vegetation has changed to agriculture, we were able to observe a significant decrease in productivity when the conversion occurs (p < 0.001, Mann-Whitney U statistic = 8480.5).Nevertheless, riparian areas that changed to agriculture have high NDVI values, depending on the crop type and the time of the year (perennial and annual crops during their productive cycle).Most of time, the opening of new agricultural fields implies the removal of natural vegetation and its substitution by annual or perennial crops, which for this semi-arid region of Mexico represents a decrease in productivity.
Abandonment of agricultural fields has been documented in this region [60]; the fields had been occupied mainly by riparian mesquite [25].Contrary to our expectations, this change also results in a decreasing productivity (p = 0.006, Mann-Whitney U statistic = 5438), which can be explained by the loss of vegetal coverage in the early stages of the replacement of agriculture by riparian vegetation (Figure 9b).Additionally, it is worth noting that new mesquite woodlands often present low values of foliar area.
Finally, changes from riparian vegetation to cultivated grasslands and vice versa also represent drastic changes in NDVI (Figure 9c,d).Substitution of RE for prairies of buffel-grass has been a common practice, above all in the region of ZR [63].In sites where these changes have been presented, a significant reduction in the values of NDVI is observed (p < 0.001, Mann-Whitney U statistic = 2379; Figure 9c).The change from cultivated grassland to riparian vegetation has also been reported, and has been explained from the perspective of the resilience of RE, since "it appears the conditions for the growth and maintenance of the prairies are not adequate to definitively substitute riparian vegetation, and for this reason, the land cover reverts to its previous state" [25].Recovery of the riparian vegetation is reflected in an apparent increase in productivity (p < 0.001, Mann-Whitney U statistic = 15,955; Figure 9d).

Carbon Storage in the REAZ
Our results from carbon stock estimation by area suggest that the REAZ has the capacity to store much more carbon than the rest of the land cover types analyzed.Nevertheless, the REAZ constitutes a small fraction of the landscape, when compared to the land cover types, such as oak/oak-pine forest or desert scrub, which might present a larger overall storage due to their extent.However, if we take into account the 59% of the riparian vegetation that has been substituted by agriculture in the region of the San Miguel River [24], the REAZ total values of carbon would be greater than those of cultivated grassland and riparian mesquite and close to those of the desert scrub (Table 3).Moreover, the importance of the REAZ resides not only on its overall capacity to store carbon but also in the fact that this ecosystem ranks among the highest carbon stocks in terms of its aerial and root biomass, as well as in its soils [7,8].To understand the importance of the REAZ in terms of carbon storage, it is necessary analyze the ratios per area (other land cover type/REAZ) with the rest of the land cover types in arid environments.From the previous, is worth mentioning that according to our results, conversion from riparian vegetation to any other land cover type in our study area represents a loss of carbon storage capacity.
Differences in carbon storage between ZR (16,976 TonC) and SMR (23,422 TonC) can be explained as a function of biophysical and management conditions of each region.The SMR presents greater annual precipitation and less extreme temperatures compared to the ZR [27,28].The SMR land cover types are related, to a lesser degree of aridity (e.g., oak/oak-pine forest and riparian vegetation), with greater carbon storage capacity.In contrast, the vegetation in the ZR, includes mostly desert scrub and mesquite woodland, which present lower values of carbon storage.Likewise, the ZR presents a larger area devoted to the development of economic activities related to intensive

Carbon Storage in the REAZ
Our results from carbon stock estimation by area suggest that the REAZ has the capacity to store much more carbon than the rest of the land cover types analyzed.Nevertheless, the REAZ constitutes a small fraction of the landscape, when compared to the land cover types, such as oak/oak-pine forest or desert scrub, which might present a larger overall storage due to their extent.However, if we take into account the 59% of the riparian vegetation that has been substituted by agriculture in the region of the San Miguel River [24], the REAZ total values of carbon would be greater than those of cultivated grassland and riparian mesquite and close to those of the desert scrub (Table 3).Moreover, the importance of the REAZ resides not only on its overall capacity to store carbon but also in the fact that this ecosystem ranks among the highest carbon stocks in terms of its aerial and root biomass, as well as in its soils [7,8].To understand the importance of the REAZ in terms of carbon storage, it is necessary analyze the ratios per area (other land cover type/REAZ) with the rest of the land cover types in arid environments.From the previous, is worth mentioning that according to our results, conversion from riparian vegetation to any other land cover type in our study area represents a loss of carbon storage capacity.
Differences in carbon storage between ZR (16,976 TonC) and SMR (23,422 TonC) can be explained as a function of biophysical and management conditions of each region.The SMR presents greater annual precipitation and less extreme temperatures compared to the ZR [27,28].The SMR land cover types are related, to a lesser degree of aridity (e.g., oak/oak-pine forest and riparian vegetation), with greater carbon storage capacity.In contrast, the vegetation in the ZR, includes mostly desert scrub and mesquite woodland, which present lower values of carbon storage.Likewise, the ZR presents a larger area devoted to the development of economic activities related to intensive use of soil and water, such as large-scale grape cultivars and extensive livestock production that promotes the use of exotic grasses [25].Spatial information regarding the location of loss/gain of carbon storage capacity has been assessed in previous studies [20,33,64].However, the value of assessing potential carbon storage capacity in arid environments, and especially in the REAZ, has been underestimated [64].In order to mitigate the effects of greenhouse gas accumulation, it is of great importance to provide accurate information to stakeholders, managers and policy-makers regarding the current status, dynamics, and spatial distribution of the sources and sinks of carbon [65].Tools such as those utilized in this study can facilitate the ES mapping process, linking the spatial information of land cover and use to biophysical values obtained from official databases, local studies, or field samples [20].

Productivity Differences and Trends
Productivity proxies were derived using the NDVI, since the index is highly correlated to the Leaf Areas Index [66], canopy coverage, biomass, chlorophyll content, photosynthesis, respiration, etc. [16].High NDVI values represent productive vegetation coverage, while a decrease through time of NDVI, in a specific area, might suggest a steady degradation of the ecosystem biophysical conditions [67].
Given the time period analyzed (2001-2012), we found that the REAZ presents an equal, or greater, productivity than (1) agriculture, which is supplied with water and nutrients, and (2) oak/oak-pine forest, which is distributed in areas with more precipitation and cooler summer temperatures [56].However, compared to the rest the land cover types analyzed, basal productivity of riparian vegetation was always greater for all the years analyzed (Table 4, Figure 5a).The greater basal and pre-monsoon productivity of the RE (riparian vegetation and riparian mesquite) could be explained as a consequence of the constant flow of materials, water, and energy, in the systems.Our results show that agriculture, riparian vegetation, oak/oak-pine forest, and riparian mesquite land cover types present the higher values of NDVI in relation to the rest of the classes analyzed.These results provide evidence of the importance of the RE as areas of high productivity in the arid zones of Northeastern Mexico, particularly in regions such as the ZR, in which RE adjoin other less productive vegetation types, such as desert scrub, mesquite woodland, and cultivated grassland (Figures 2 and 5a,b).
Our analysis of trends shows a significant decrease in the productivity for the REAZ over time (both Landsat and MODIS sensors suggested this).In the regional context, other studies that used NDVI data from Landsat have found a progressive decrease in pre-monsoon productivity in mangroves of the Gulf of California in the period 1990-2010 [48], and in agriculture and riparian vegetation in the San Pedro River (USA) between 1984 and 2012 [49].The tendency towards a decrease in the productive capacity of the vegetation appears to be a regional phenomenon and this was also found in our analysis based on MODIS and Landsat data (Figure 6a).The gradual decrease on NDVI values in other REAZ of the region, such as the San Pedro River (USA), has been explained as a response to the reduced flow of surface water, increased air temperature, and increased depth of water table levels [49].The negative trend in productivity of the riparian vegetation and riparian mesquite could be the consequence of a synergy between the set of climatic conditions and the state of water management in the region that can give rise to the decrease in surface runoff and a possible decrease in phreatic levels, as shown in previous studies [25].
Water availability and management are considered the most prominent drivers for productivity in the REAZ.In our present analysis we found an inverse relationship between groundwater depth and REAZ productivity, demonstrating the phreatic level modifications could explain, in part, the trends mentioned before.There are around 1544 subterranean water concessions in the region, of which 898 correspond to the SMR and 646 to the ZR.The volume of water under concession for the region, according to the Public Register of Water Rights of the National Water Commission (REPDA, by its Spanish acronym), is approximately 144.9 Mm 3 /year, while the annual mean recharge is 147.3Mm 3 /year [27].Water recharge and extraction in the aquifers is, thus, very similar and most of the wells are found close to the course of the rivers, which could increase the water table depth, mostly in areas with greater agricultural activity [25].On the other hand, a total of 387 "represos" have been reported in SMR alone [68], covering a total area of 185 ha and retaining the runoff that would otherwise reach the watercourses.Most of the dams are constructed without authorization, since CONAGUA only has around 80 surface water concessions registered in both of the sub-basins [69].Human activities present in the region (deviation of watercourses and detention of surface water flow) have reduced the availability of water for other purposes, such as urban and ecological uses.This might alter the structure and functioning of the REAZ.With a reduction on surface and ground water availability, riparian plant species of arid and semi-arid zones make structural adjustments, reducing their LAI and canopy coverage [70], a phenomenon that can be identified remotely through a reduced NDVI [49].These adjustments seek to maintain the function of the RE when water is scarce, but if the lack of water is maintained, the system can undergo a transition from the strictly riparian vegetation formed by genera such as Populus and Salix, towards plant communities composed by facultative species such as Prosopis, Tamarix, Celtis, and other xerophytic trees [1][2][3].

Conclusions
There is a greater quantity of carbon stored in the central and northern sections of the SMR, a region that has been less dynamic in terms of the development of agricultural and livestock production activities.This region also features vegetation types with greater values of the productivity proxies (NDVI: greater integral-lesser integral of MODIS and NDVI: pre-monsoon of Landsat).In contrast, the ZR shows lower quantity of carbon stored per hectare.This is partly due to the introduction of large fields of induced grasslands in areas where riparian vegetation, mesquite woodland, and desert scrub were present.
Abandoned agricultural areas can be converted over time into carbon sinks, which represents an exchange of costs and benefits between a social contingency (land abandonment) and an ecological benefit (recovery of vegetation).
Despite the fact that other vegetation types of greater coverage in the region (such as the subtropical scrub oak/oak-pine forest or mesquite woodland) represent a greater carbon storage, the RE feature the land cover types with the greatest capacity for carbon storage per unit area (even greater than that of the oak/oak-pine forest).Likewise, basal and pre-summer monsoon productivity values in the REAZ are equal, or higher (some years), than land cover types such as oak/oak-pine forest or agriculture, which are considered to be highly productive.
A possible synergy between the increase in temperature, mainly of autumn-winter, the decrease in precipitation, and of surface runoff, as well as changes in land use and the intensity of water use, could have a negative influence on productivity in the REAZ.
This study provides evidence to enable the following conclusions about productivity: 1.
There is a decrease in the time of the basal and pre-monsoon productivity of the REAZ.This is partly due to the drier conditions observed during the previous decades in the region, partly provoked by ENSO dominance in recent years [71].

2.
Productivity of riparian vegetation decreases with the increased depth of the static levels of subterranean water.

3.
Productivity decreases in sites where the REAZ has transformed to other land cover types.

4.
In sites where there is an intensive use of water and soil, a transition is seen between strictly riparian (genera Populus and Salix) and facultative (genera Prosopis, Acacia and Celtis) species.
The present work improves our understanding of how land use change dynamics, intensive water use, and changes in meteorological conditions can have an impact on the ES of productivity and carbon storage offered by the REAZ.

Figure 1 .
Figure 1.Study area location.The region of the San Miguel-Zanjón Rivers is located in the central part of the state of Sonora, Mexico.Social and economic development in both regions has been closely related to the use and exploitation of the RE.Provision of ground and surface water, as well as use of the fertile land on the floodplains, are the most important provision ES, and have allowed the development of agricultural

Figure 1 .
Figure 1.Study area location.The region of the San Miguel-Zanjón Rivers is located in the central part of the state of Sonora, Mexico.

Figure 2 .
Figure 2. Classification of land cover used in this study [25].

Figure 2 .
Figure 2. Classification of land cover used in this study [25].

Figure 3 .
Figure 3. Spatial distribution of the carbon storage in the region of San Miguel River (SMR)-Zanjón River (ZR).

Figure 3 .
Figure 3. Spatial distribution of the carbon storage in the region of San Miguel River (SMR)-Zanjón River (ZR).

Figure 4 .
Figure 4. Changes in carbon storage in the SMR-ZR region.

Figure 4 .
Figure 4. Changes in carbon storage in the SMR-ZR region.

Figure 5 .
Figure 5. Accumulated Normalized Difference Vegetation Index (NDVI) values of MODIS, representing the proxy of the basal productivity (a) and NDVI values of Landsat, representing the proxy of productivity of the vegetation before the summer monsoon (b).

Figure 5 .
Figure 5. Accumulated Normalized Difference Vegetation Index (NDVI) values of MODIS, representing the proxy of the basal productivity (a) and NDVI values of Landsat, representing the proxy of productivity of the vegetation before the summer monsoon (b).

Figure 6 .
Figure 6.Trends in productivity (proxy) from the NDVI time series of MODIS (a) and Landsat (b), representing the proxies of the basal (large integral-small integral) and pre-monsoon productivity, respectively, for each land cover type (p values in bold, show significance).

Figure 6 .
Figure 6.Trends in productivity (proxy) from the NDVI time series of MODIS (a) and Landsat (b), representing the proxies of the basal (large integral-small integral) and pre-monsoon productivity, respectively, for each land cover type (p values in bold, show significance).

Figure 7 .
Figure 7. Trends in total annual precipitation and annual average temperature for the SMR (a,e) and ZR (b,f).Accumulated precipitation and the autumn-winter average temperature for the SMR (c,g) and for the ZR (d,h).

Figure 7 .
Figure 7. Trends in total annual precipitation and annual average temperature for the SMR (a,e) and ZR (b,f).Accumulated precipitation and the autumn-winter average temperature for the SMR (c,g) and for the ZR (d,h).

Figure 8 .
Figure 8. Spatial relationship between water depth and pre-monsoon NDVI values of the riparian vegetation for different years of the study period.

Figure
Figure Spatial relationship between water depth and pre-monsoon NDVI values of the riparian vegetation for different years of the study period.

Figure 9 .
Figure 9.Comparison of the proxy of total productivity during the growing season (greater integral of MODIS) among different land cover types that remained unchanged during the period 2002-2011, relative to the same classes that underwent changes from riparian vegetation to agriculture or cultivated/induced grassland, and vice versa.

Figure 9 .
Figure 9.Comparison of the proxy of total productivity during the growing season (greater integral of MODIS) among different land cover types that remained unchanged during the period 2002-2011, relative to the same classes that underwent changes from riparian vegetation to agriculture or cultivated/induced grassland, and vice versa.

Table 1 .
Meteorological stations used to analyze the trends in temperature and precipitation between 1988 and 2010.

Table 2 .
Carbon stocks per land cover type.

Table 3 .
Values of total carbon considering the total area of each land cover type for 2011.

Table 4 .
Results of all pairwise multiple comparison procedures (Dunn's method) between basal productivity (NDVI-MODIS) of riparian vegetation and the rest of land cover classes (for the period from 2001 to 2012).Cells in gray indicate no significant differences between treatments (p > 0.05).

Table 4 .
Results of all pairwise multiple comparison procedures (Dunn's method) between basal productivity (NDVI-MODIS) of riparian vegetation and the rest of land cover classes (for the period from 2001 to 2012).Cells in gray indicate no significant differences between treatments (p > 0.05).

Pre-Monsoon Productivity of the Riparian Vegetation (Landsat)
* Statistically significant at p < 0.05.

Table 6 .
Pearson correlation analysis between the proxy of the pre-monsoon and basal productivity of the riparian mesquite, accounting for climatic parameters of the ZR.

of Pre-Monsoon Productivity of the Riparian Mesquite (Landsat)
* Statistically significant at p < 0.05.