Effects of Land Use Changes on Ecosystem Services Value Provided By Coastal Wetlands: Recent and Future Landscape Scenarios

Wetlands are very sensitive to land use changes which alter supply and quality of their ecosystem services. This study analyzes the effect of land cover changes on the spatial and temporal patterns of the value of ecosystem services provided by coastal wetlands in northwest Mexico. These ecosystems have a high degree of naturalness, but are at risk because of land use changes promoted to reactivate regional development. Remote sensing and Markov chains modeling were used to estimate change trends, together with a value transfer approach for the ecosystem service valuation. The results indicated that the total flow of ecosystem service value tended to increase (18 million dollars (2007 USD)), presumably biased by the highest worldwide value assigned to saltmarsh/unconsolidated bottom, which increased in area by 8% during the study period. The most notable transition probability was observed between natural wetlands, highlighting littoral and saltmarsh as the classes with the highest probability of change over time. The southern part of the study area is the most susceptible to change, where unconsolidated bottom and forested mangrove (saltmarsh) are prevalent. Therefore, we can argue that the conservation of these coastal environments should be a priority in future land use management.


Introduction
The interest in the ecosystem services and the importance they have increased since the publication of the Millennium Ecosystem Assessment [1], with the report of global decline of 15 of the 25 listed ecosystem services, which anticipates a large and negative impact on future human welfare. Since then, one of the demands from the MA was to increase research on measuring, modeling and mapping ecosystem services, resulting in a positive response to increase the efforts to identify, quantify and value ecosystem services [2,3]. In fact, identifying the economic value of these services is essential in revealing their societal value, providing common metrics to facilitate comparisons among attributes and divergent scenarios in policy assessment [4]. These values differ significantly depending on the socioeconomic and geographic contexts, assigning higher values to wetland ecosystem services in countries with higher incomes, or valuing lower the hectare of a large wetland compared with the same but smaller wetland type [5].
The provisioning and quality of ecosystem services depend on the integrity and functioning of ecosystems, but the environment is systematically transformed, sometimes to increase production in agriculture, to increase space for human settlements or for resource extraction, with immediate or short time positive impacts on the local economy. However, it causes a decline in the provision of ecosystem services, which sooner or later will create a cost that must be paid by the population to compensate for the loss.
Although the impact of human activities is common to all natural ecosystems worldwide, coastal wetlands have been particularly affected because they receive the cumulated impact from all activities in the watershed. Although there are many drivers of change, Land Use/Land Cover (LULC) changes have been responsible, directly or indirectly, for the loss of 40% of all coastal wetlands and are the main source of pollutants discharged into coastal environments [6]. Moreover, in many parts of the world, local natural beauty, the availability of resources or the high productivity of natural wetlands are rapidly being converted into tourism facilities, salt evaporation or aquaculture ponds, housing developments, roads, ports, hotels and golf courses, thereby resulting in reduced or modified biodiversity, altered functional processes and the diminishing provision of ecosystem services to society, causing substantial social and environmental costs [7].
Mexico is not an exception to this global trend because the most important tourism infrastructure and technical agriculture have been developed along the coastal zone, modifying land use patterns and resulting in degradation and loss of natural ecosystems and their functions. Consequently, a decline in the provision of ecosystem services to local communities occurs [8,9].
In northwest Mexico, the coastal zone of Sinaloa maintains ecologically important wetland zones with a high degree of naturalness, but particularly in southern Sinaloa, they are at risk mainly as consequence of governmental development plans to reactivate the regional economy and promote tourism and agriculture by creating facilities and infrastructure such as resorts, golf courses, hydrolectrical dams and irrigation and drainage channels. Therefore, to guarantee the provision of ecosystem services from wetlands, policy interventions are needed to steer these processes and mitigate their negative impacts on ecosystems and society [10]. Such policy interventions require accurate assessments of the relevant ecological, socioeconomic, and political developments as well as of their potential future implications. LULC change analyses have been indicated as one of the high priority concerns for research and for the development of strategies for sustainable management [11]. In this context, it is important to produce information related to landscape changes and their effects on the value of ecosystem services provided by wetlands in this zone. This study focuses on landscape changes, but it also contributes an update to the current status and trends of wetlands, making it possible to compare and evaluate the output offered by alternative investments and highlighting the importance of wetlands as natural capital. Therefore, the main objectives of this study were to assess the dynamics of LULC based on recent Earth Observation data and forecast future changes in Ecosystem Service Values (ESV) provided by wetlands by modeling future LULC scenarios.

Study area
The Southern Coastal Zone of Sinaloa (SCZS) is integrated by four municipalities (Mazatlan, Concordia, Rosario and Escuinapa, from north to south) on the west coast of Mexico, with an approximately 120 km coastline and a total area of approximately 2500 km 2 ( Figure 1).
Despite that Rosario does not have coastline, the municipalities of the SCZS share similar environmental conditions, including diverse coastal wetlands such as mangroves, saltmarshes, estuarine and fluvial systems. The climate is warm and humid, with rainfall in the summer (928 to 1,457 mm) and an annual average temperature varying between 22 and 26°C [12].
The most important wetland ecosystems are represented by the Estero de Urias and the Presidio River at Mazatlan, the Huizache-Caimanero lagoon system and the Baluarte River in Rosario, and diverse lagoons and saltmarsh ecosystems integrated into the Marismas Nacionales (Sinaloa) system, which also includes the Teacapan Estuary in Escuinapa ( Figure 1).
There are more than 35 urban and rural communities in the region, and all of these communities have less than 5,000 inhabitants, with the exception of Mazatlan, Villa Union, El Rosario, and Escuinapa (all connected by federal highways). Among them, Mazatlan is the most important city, with a population above 430,000 people, 87% of the total municipal population, and infrastructure to support tourism, industrial and fishing harbor, and other economic activities [12].
The population growth in Mazatlan (5% annual rate) and a population density of 173 people per km 2 produces high pressure on natural land cover, and, because of the lack of appropriate land use planning and measures for sustainable development, environmental stress is also growing. The remaining municipalities have population densities below the state average (48 people per km 2 ), but productive developments (aquaculture, agriculture, tourism) are growing, which also adds pressure to the local natural resources.
Despite the importance of coastal ecosystems and after several attempts, there are currently no natural protected areas in the region, except for the federally protected oceanic islands and mangrove systems.

Methods
To analyze LULC changes and forecast trends for the next two decades in both land cover and ecosystem service values, the methodological approach used here included three steps: (i) the production of LULC thematic maps from 2000 and 2010; (ii) an estimation of change trends using Markov chains; and (iii) an analysis of the annual ESV flow based on wetland type and change tendencies.

LULC classification
In this study, the land cover data sets were produced from multi-spectral Landsat Thematic Mapper (TM) imagery (path/row: 31/44), acquired from the USGS Global Visualization Viewer (http://glovis.usgs.gov/), corresponding to the dry season in May 2000 and March 2010. Both images were clear and nearly free of clouds, previously projected to UTM coordinates (zone 13Q, WGS84) before limit the study area to the boundaries of municipalities and the physiographic province of the Pacific Coastal Plain (PCP) by a masking process. The image processing, GIS development and output of the final coverage maps were conducted using IDRISI Taiga software and ESRI ArcGis 9.3.
Each image was updated with vector features such as rural and urban areas, channels, rivers and shrimp farms, which were digitized on-screen from QuickBird images from 2002-2010 available from Google Earth. These features were added as a mask in every spectral band, eliminating the corresponding pixels for further classification processes. The modified images were later independently classified into six informational classes of natural wetlands as defined by Berlanga et al. [13] using a supervised method with the maximum likelihood algorithm. Once the classification process was concluded, the previously defined polygon vectors were added and reclassified to produce a total of eight LULC categories ( Table 1).
The validation of the output maps was conducted based on the proposal by Pontius et al. [14]. Each classified map was compared with a reference map to generate an observed sample matrix (cross-tabulation) that represents the absolute frequencies n ij of the i-th category in the classification map (rows) and the j-th category in the reference map (columns). The values in the diagonal n jj correspond to concordances in both the rows and columns classes. From this, a population matrix representing the proportions of the different categories in the classified map and the reference map (p ij ) was created: where N i is the observed frequency of the i-th category in the classified map (total of the i-th row). The population matrix allows the estimation of unbiased statistics, as the quantity disagreement (Q) and allocation disagreement (A) parameters, calculated for each class and overall. The quantity disagreement parameter is derived from the differences between the proportions of categories in the classified map and the reference map, taking value of zero when the categories have the same proportion in both maps. The allocation disagreement parameter occurs when differences in the spatial allocation of the categories in both maps are present, and its value depends on the number of pixels in the classified map that need to be reallocated to increase the agreement with the reference map [14].
In addition to the Q and A disagreement values, the standard Kappa index (K') was also estimated. This index accounts for the expected agreement due to random spatial relocation of the categories in the classified map given the proportions of the categories in both maps, regardless of the size of the quantity disagreement parameter. The index takes values in the range from 0 to 1, when the observed agreement is equivalent to the statistically expected random agreement and when the agreement is perfect, respectively. All formulas of these indicators of agreement and disagreement are fully explained in [14].

Modelling LULC changes
This process was achieved using Markov chains to represent a spatial transition model, which is widely used to model LULC changes [15]. The Markov chain analysis follows a stochastic process in which distinct LULC categories are considered as states of the chain. A chain is defined as a stochastic process having the property that the value of the process at time t, X t , depends only on its value at time t-1, X t-1 , and not on the sequence of values X t-2 , X t-3 ,…,X 0 , which the process passed to arrive to the stage X t-1 [16,17]. The Markov chain can be expressed as: where the conditional probabilities P ij represents single-step transition probabilities or simply transition probabilities of the Markov chain; and P ij represents the probability estimated during the transition process from state i to state j in a given period of time [16]. P ij can be estimated from the values observed in a change detection matrix by tabulating the observed area of data i and j (n ij ) and dividing it by the sum of the number of times (N i ) that state i has occurred [17]: The 2000 and 2010 land-cover maps were contrasted with a crosstabulation change detection matrix, and the statistical independence of LULC was tested with the chi-squared (X 2 ) test [17]. Later, the transition probabilities P ij were calculated to output a single step transition probability matrix, a discrete-time Markov chain (DTMC), and this was generalized to an n-step transition probability matrix with the Chapman-Kolmogorov for Markov chains [16]: In matrix notation: The matrix of n-step transition probabilities is obtained by multiplying the matrix of one-step transition probabilities by itself (n -1) times n n P P = ) ( [16]. We projected the matrix P for 2020 (n = 1) and 2030 (n = 2). The analysis was completed using the Idrisi Taiga Markov module [18].

Estimation of ESV
A wide array of market and non-market methods has been used to economically estimate the ecosystem services values of wetlands [19][20][21][22]. Based on the data set provided by Ghermandi et al. [23], which uses value references from empirical studies, Valdez et al. [24] applied a value transfer method to derive economic values for ecosystem services provided by coastal wetlands in Sinaloa, validating the method through a meta-regression model. Here, those values were used to obtain estimates of ESV for 2000, 2010 and for the 2020-2030 projection.
The ESVs for each land cover type and the total flow for the two study years (2000 and 2010) were calculated by multiplying the value coefficients obtained from Valdez et al. [24] by the respective cover extent. The change in the ESV was estimated by summarizing the value of each LULC type from 2000 and by calculating how each type changed in the following decade (2010). Based on the projected areas of the Markov chains model (ha), we calculated the total flow in the horizon years (2020 and 2030) using the same value coefficients as in previous years (2000 and 2010) with regard to the natural wetland categories.

Overall trend in the LULC
Once the imagery classification processes were completed, the accuracy of the output maps was validated to define the extent of each land cover category. From the confusion matrix, the proportion of correctly classified pixels was estimated at 0.83 (83%). Consequently the total disagreement was 0.17, and most of this, around 80%, corresponded to Q, while A rise up to 20%. This implies that the classification errors are mainly associated with the overestimation/underestimation of the areas of each category, rather than the distribution of the categories on the map. The classes with major allocation disagreement were Land cover (Lco) and Agriculture (Agc) reaching 97% and 100%, respectively Land cover Land cover: tropical forests, secondary succession. Cover symbols: Littoral (Lit), Coastal lagoon (Cla), Saltmarsh/unconsolidated bottom (Sun), Saltmarsh/forested mangrove (Smn), Riverine (Riv), Aquaculture ponds Aps), Urban (Urb), Agriculture (Agc) and Land cover (Lco).  ( Figure 2). Conversely, the Littoral category (Lit) was classified without error and the disagreement in the Coastal lagoon category (Cla) was mostly due to A (>60%). Considering all the land use and cover categories, the standard Kappa index was K'= 0.77, which corresponds to a substantial agreement on the [25] classification scale.
Based on those results, the output maps were assumed to be a reliable interpretation of the landscape of the study area in both years analyzed. The LULC maps for the dry season in years 2000 and 2010 are displayed in Figure 3. The dominant covers for both dates were land cover (Lco) and agriculture (Agc), but regarding wetlands, as defined by Berlanga et al. [13] in their classification system, the most important were coastal lagoon (Cla), saltmarsh unconsolidated bottom (Sun) and saltmarsh forested mangrove (Sme).
Most of the cover classes displayed surface gains and losses, reflecting an spatial dynamics (Figure 4a), but at the final balance, saltmarsh/unconsolidated bottom and coastal lagoon, were the wetland covers that most positively contributed to the landscape transformation, with net gains of 1990 ha and 1718 ha, respectively, while saltmarsh/forested mangrove and littoral classes had a net loss of 1567 and 529 ha (Figure 4b). The agriculture cover decreased considerably, losing a net area of 8,634 ha during the study period. Table 2 shows the LULC change matrix from 2000 to 2010, and it is clear that significant change (13% of the total area) happened during the 10-year period, with the main changes attributed to interactions between agriculture soils and other land uses or coverages (Lco). Also the natural wetland classes (Lit, Cla, Sun and Smn) contributed to the Lco increasing with around 3500 ha or about 30% of the total gain for this cover from 2000 to 2010.
Although the highest proportion of change involved non-wetland inland covers, littoral and the only forested wetland class (Smn) suffered a considerable decrease of about 31% and 14% respectively. By contrast, coastal lagoon and saltmarsh/unconsolidated bottom classes increased about a tenth of their former area, while the artificial wetland Aps also gained around 300 ha.

LULC-Markov analysis
The simple-step transition probability matrix P, resulting from the contrast of the 2000 and 2010 study area LULC patterns, is shown in Table 3, displaying statistical dependence, X 2 = 1.3 × 10 6 (P < 0.05). The resistance to change from one class to another can be observed on the diagonal of the matrix, representing the fraction of pixels that maintain the same LULC in the initial and final images. In this time period, littoral and saltmarsh (both unconsolidated bottom and forested mangrove) were the wetland classes most susceptible to change, together with agriculture. By contrast, coastal lagoon, riverine and aquaculture pond wetlands, together with urban category, were the most resistant to change. It highlights that lower resistance was observed among natural wetlands.
According to Table 3, the Lco class, involving different land use and coverages, was the most transformed, gaining surface from agriculture and littoral classes (P>0.30). Just regarding wetlands, saltmarsh/ forested mangrove class was the most affected among wetland classes, being transformed into saltmarsh/unconsolidated bottom class, with a transformation probability of 0.23, the highest if littoral to land cover change is disregarded; at the same time, the saltmarsh/unconsolidated bottom class decreased and transformed into the coastal lagoon class and the saltmarsh/forested mangrove class with probabilities of 0.22 and 0.01, respectively. The littoral class is also transformed into land cover with a probability of 0.33.

Extent of LULC categories
Based in the calculated transitional probabilities, the most extensive wetland categories during the study years (2000, 2010 and the predicted 2020-2030) are represented by the coastal lagoon and saltmarsh (both unconsolidated bottom and forested mangrove) classes, amounting together around 21% of the area, increasing 1.1% in the whole study period. The land cover class, which represents the natural vegetation and other land cover classes, was always above 50% of the total study area, excepting in 2000 when occupied 48.1%. In addition, it is expected that for 2030, the agriculture class will loss an equivalent of 10% of the total area, while urban cover and aquaculture ponds will increase less than 1% of the same extent, probably growing on agriculture land.
All the categories will maintain their trend observed for the period 2000-2010, with the exception of the saltmarsh/unconsolidated bottom class that maintained an increasing trend, but at the end it is expected a marginal reduction (193ha) during the 2030 projected year ( Table 4).

Variation of ESV
The estimation of the total Environmental Services Value (ESV) flow was obtained based only in the temporal patterns associated with natural wetland cover types (Table 5). It was calculated an increase from about $223 million (2007 USD dollars) up to $241 million when the ESV was projected to 2020. However, in the next projected decade (2030), the total flow decreases by approximately 2 million (2007 USD). The ESV of the saltmarsh classes (unconsolidated bottom and forested mangrove) accounts for 90% of the total estimated ESV, followed by the coastal lagoon (7%), and riverine (1%) classes. The total ecosystem service value increases over time, amounting to approximately 18 million (2007 USD).
The increase in the total ESV during the study period is presumably a result of the highest value assigned worldwide to the saltmarsh/ unconsolidated bottom cover. However, the analysis projects a significant decline (approximately 12%) in the ESV provided by the saltmarsh/forested mangrove class, representing potential high costs for the local inhabitants.

Discussion and Conclusions
The ecosystem service valuation issue is far from solved, mostly because there is not agreement on the ES concepts and classification, although the Millennium Ecosystem Assessment (MA) scheme is currently the most widely used despite its inconsistencies [26]. While the concepts evolve and criteria are unified, this classification system has been used for landscape management and valuation purposes, but more studies on mapping, modelling and spatial dynamics of the SE are       required, to anticipate changes, understand processes and planning [1].
Depending on the extent or characteristics of the study area, in many cases, applications from remote sensing and GIS provide the only economically feasible way to gather regular land cover information to achieve the above mentioned goals [27,28]. Here, we use the potential information provided by multi-temporal Landsat TM data analyses to accurately map and analyze trends that can be used as inputs in the ES valuation, for further management processes regarding wetlands in southern Sinaloa.
Together with other analytical tools that allow some level of forecast, it is possible to produce indicators on the direction (from one class to another) and extent (ha) of LULC change in the future with Markov chains that output transition probabilities, a measure of the possible change. Here, this methodology forecast that littoral and saltmarsh wetland classes (both unconsolidated bottom and forested mangrove) are those with the highest probability of change in the near future. These results are consistent with findings from Robles et al. [24] for the northern coast of Nayarit, who found that saltmarsh, was the type of wetland most susceptible to change, particularly by effect of shrimp farm growing.
The total flow of ESV by wetland category was estimated based only on the land cover that resulted from the projection derived from Markov chains. The results indicate that the rates of change in the total ecosystem service value increase during 2000 to 2010 period, with the total flow of ESV increasing by 18 million (2007 USD), mostly as consequence of the increase in area of the unconsolidated saltmarsh cover. This positive trend in the total ESV flow contrasts with the trends found in other studies. For example, Kreuter et al. [29] and Zhao et al. [30] conducted assessments on ecosystem service values, and both studies found that the total annual ecosystem service value declined over the time and this reduction is attributed mainly to the effect of urban sprawl on land cover.
In this sense, our study suggests that even when there has been an urban growth during the study period, there is still no influence of this trend on the ESV provided by natural wetlands, and changes occur among wetland covers that supply different ES and their related values. In addition, it is important to note that the increase in unconsolidated saltmarsh is most likely a consequence of mangrove loss, but unlike other wetlands, saltmarsh is highly valued in other latitudes but not in Mexico, where this cover is considered unproductive, salty and bare soil. As a consequence, increases in this cover could be overvaluing the estimation. Future research to assess the ecosystem services value for this cover at a regional level is necessary to adjust the total flow. In fact, the use of estimates derived from Camacho-Valdez et al. [24] was a good starting point to assess the values of ecosystem services in the study area, but also is required to conduct a primary study that focus on the evaluation of saltmarshes.
Based on the landscape trends observed from 2000 to 2010, a change in total flow is expected in the future, and in the projection for the year 2030, the ESV provided by the saltmarsh/unconsolidated bottom cover could decrease by approximately 200 ha, mostly as consequence of the increase in tourism infrastructure, which is projected to encompass approximately 40,000 hotel rooms and the services they require. It is also expected that the population will grow with urban expansion, leading to a decline in natural wetland covers and therefore, in the ecosystem services they provide.
Focusing only on the ESV estimations for each wetland, the reduction of the saltmarsh/forested mangrove category over the time was the most important change in the study period, although changes affecting the value of mangrove coverage are due to interactions between natural wetlands instead of changes in land use. The provision of ecosystem services and the ESV of this wetland type could be altered if the negative trend continues. However, despite the relevance of this wetland cover, its contribution to the total flow is not significant (5% on average). Thus, its variation has a relatively small impact in monetary terms, as measured in the present approach.
The importance of ESV has been widely recognized as a useful tool to enhance land use planning [7]. Unfortunately, the underlying uncertainties and constraints in an ecosystem valuation using value transfer make it difficult to encourage policy makers to incorporate such valuation in environmental management studies, even when ESV estimation continues to be an efficient element to guide decisions for policy makers. Some of the limitations related to the use of value transfer in decision making could be addressed through greater interactions between researchers and policy analysts [31].
Although the described limitations exist, research using more accurate ESV methods that incorporate the spatial dimension, together with the best economic indicators at a regional level, are urgently needed to encourage policy makers to use ESV as an additional tool in their efforts to balance sustainable land use and regional development conditions [32].