Ecosystem Services Mapping Uncertainty Assessment : A Case Study in the Fitzroy Basin Mining Region

Ecosystem services mapping is becoming increasingly popular through the use of various readily available mapping tools, however, uncertainties in assessment outputs are commonly ignored. Uncertainties from different sources have the potential to lower the accuracy of mapping outputs and reduce their reliability for decision-making. Using a case study in an Australian mining region, this paper assessed the impact of uncertainties on the modelling of the hydrological ecosystem service, water provision. Three types of uncertainty were modelled using multiple uncertainty scenarios: (1) spatial data sources; (2) modelling scales (temporal and spatial) and (3) parameterization and model selection. We found that the mapping scales can induce significant changes to the spatial pattern of outputs and annual totals of water provision. In addition, differences in parameterization using differing sources from the literature also led to obvious differences in base flow. However, the impact of each uncertainty associated with differences in spatial data sources were not so great. The results of this study demonstrate the importance of uncertainty assessment and highlight that any conclusions drawn from ecosystem services mapping, such as the impacts of mining, are likely to also be a property of the uncertainty in ecosystem services mapping methods.


Introduction
Ecosystem services (ES) are the benefits that ecosystems provide to human beings, bridging the connection between ecological and social systems [1][2][3][4].This concept is increasingly used in both theoretical and applied research [5,6] such as spatial planning [7,8], environmental economics [9] and risk management [10].ES play a crucial role for supporting ecosystem-based management and decisions [11,12] and the mismanagement of ES can potentially have large economic impacts and may affect thousands of people [13,14].In order to support the decision-making, accurate and reliable methods and outputs for ecosystem services assessments are increasingly demanded [15,16].
Among the various ecosystem services assessment methods, ecosystem services mapping has received increasing attention [17,18].Mapping is a significant tool for identifying and quantifying the spatial characteristics of ecosystem services [19].Mapping can provide researchers and decision makers with spatially explicit information to help make feasible policies and strategies [20,21].There has been a growing number of publications on ES mapping including significant reviews [17,18,22,23] and case study applications from different perspectives and scales [24,25].A range of mapping methods and integrated mapping tools have been developed, such as the widely used InVEST (Integrated Valuation of Ecosystem Services and Trade-offs) [26] and ARIES (Artificial Intelligence for Ecosystem Services) [27].Furthermore, a growing number of ES mapping research groups and projects have been established by different research organizations and governments, such as the Ecosystem Services Partnership and Natural Capital Project [19].
Although the number of ES mapping studies has increased, there are still multiple challenges due to models' limitations, data availability and other relevant issues [17,28].For instance, different mapping methods may produce completely different outputs for the same ecosystem service [29,30], while data availability may limit the choice of mapping methods and affect map accuracy [31,32].Consequently, there are many uncertainties in ecosystem services mapping, which are complicated and expressed in different ways.These uncertainties make accurately assessing and explicitly conveying ecosystem services values and spatial distributions for researchers and decision-makers challenging.A crucial step forward is conducting an assessment of uncertainty to quantify the robustness and reliability of a study, test key hypotheses and identify hedging opportunities [33].
Although it is recognized that uncertainties are significant for ecosystem services mapping and modelling, only a small fraction of the literature has assessed uncertainty and only a few studies have focused specifically on uncertainty [16,[33][34][35][36].This contrasts with the literature for other spatially explicit modelling in related disciplines such as landscape ecology and hydrology [37][38][39][40][41].However, the assessment of uncertainty in ES mapping is growing and there are various definitions and perspectives used in analysing uncertainties [16,33,35,36].
Existing studies on uncertainty and uncertainties derivations in ES assessment generally consider three types of uncertainty [33,35].Firstly, there are uncertainties associated with characterizing ecological functions and processes.Ecosystem functions are complicated and for many mechanisms there is still not a clear and sufficient understanding [20,42].Secondly, uncertainties arise from subjective preferences of researchers and other participants.For instance, as the purpose of ES modelling may vary from theoretical predictions, making practical decisions or communication with stakeholders, the preferences of researchers for a specific model and parameterization will result in different forms of uncertainty [43,44].Also differences in expert experience or the subjective judgments about relative ecosystem services can lead to uncertainties [45,46].Finally, practical uncertainty may result from the modelling process employed such as the modelling methods and data.Differences between diverse modelling methods used for the same ecosystem service may lead to differences in outputs [35].
The objective of this study was to assess the impact of uncertainty on ecosystem services mapping.We used the Fitzroy Basin mining region as a case study to map hydrological ecosystem service water provision with InVEST.Uncertainty was modelled by altering various inputs and parameterizations of the ES model.A baseline was constructed using the best available data and parameterization using the seasonal water yield model of InVEST.We then assessed three types of uncertainty: (1) spatial data sources, (2) modelling scales (temporal and spatial) and (3) parameterization and model selection.We compared outputs from the baseline model to the uncertainty scenarios using a range of indices.We conclude by discussing the implications of uncertainty for ecosystem services mapping and potential for future research on ecosystem services mapping uncertainty.

Study Area
The study area was the Fitzroy Basin (Figure 1), one of the largest (142,665 km 2 ) basins in the Great Barrier Reef catchments and also one of the largest region in Australia [47].There are six sub-catchments drained by six rivers including the Isaac-Connors, Mackenzie, Comet, Nogoa, Fitzroy and Dawson rivers.The region has a dry and wet season with distinctive climatic characteristics and the mean rainfall for the basin is approximately 630 mm .Most of the study area is relatively flat except the north-eastern area with low hills.The main land uses are grazing, cropping and forest, as well as mining.The urban and towns mainly service the coal mining industry and agriculture [48].The Fitzroy Basin catchment is of great strategic significance; both environmental as well as mining.The urban and towns mainly service the coal mining industry and agriculture [48].The Fitzroy Basin catchment is of great strategic significance; both environmental and economic [49,50].On one hand, minerals production provides enormous economic stimuli for the industry and for local community development.On the other hand, the environmental and ecological protection have attracted more and more attention from researchers and public, especially concerns associated with hydrology such as water supply and water quality and the impacts on Australia's Great Barrier Reef World Heritage Marine Park.

Baseline Scenario and Ecosystem Service Models
This study focused on the hydrological ecosystem service-water provision-characterized by water yield at both the pixel and catchment scales.We set a Baseline scenario based on the InVEST seasonal water yield model to characterize the water provision service of the study area.This scenario was used as reference which was based on the most feasible ES model with the best available data, precise ES model and appropriate parameterization for local conditions.For example, the seasonal water yield model was used rather than an annual model as it can better reflect the real climate conditions caused by the distinctive wet and dry season of the study area as opposed to an annual water yield model which is also commonly used.The uncertainty scenarios represent realistic alternative ways of modelling ES for the study area.

Seasonal Water Yield Model
The seasonal water yield model in InVEST is a spatially explicit raster based method for characterizing water yield based on three outputs: quick flow, local recharge and base flow.Quick flow characterizes short term stream flow generated from storm rainfall [51], which is calculated based on the runoff curve number method, developed by USDA (United States Department of Agriculture) Natural Resources Conservation Service (the runoff curve number predicts runoff or infiltration from excess rainfall for classified hydrological soil groups).Base flow characterizes stream flow emanating from groundwater.Finally, local recharge is water which moves from the surface to recharge groundwater, which in turn can contribute to potential base flow.Details of the model processes can be found in Table 1 [26,52,53].

Baseline Scenario and Ecosystem Service Models
This study focused on the hydrological ecosystem service-water provision-characterized by water yield at both the pixel and catchment scales.We set a Baseline scenario based on the InVEST seasonal water yield model to characterize the water provision service of the study area.This scenario was used as reference which was based on the most feasible ES model with the best available data, precise ES model and appropriate parameterization for local conditions.For example, the seasonal water yield model was used rather than an annual model as it can better reflect the real climate conditions caused by the distinctive wet and dry season of the study area as opposed to an annual water yield model which is also commonly used.The uncertainty scenarios represent realistic alternative ways of modelling ES for the study area.

Seasonal Water Yield Model
The seasonal water yield model in InVEST is a spatially explicit raster based method for characterizing water yield based on three outputs: quick flow, local recharge and base flow.Quick flow characterizes short term stream flow generated from storm rainfall [51], which is calculated based on the runoff curve number method, developed by USDA (United States Department of Agriculture) Natural Resources Conservation Service (the runoff curve number predicts runoff or infiltration from excess rainfall for classified hydrological soil groups).Base flow characterizes stream flow emanating from groundwater.Finally, local recharge is water which moves from the surface to recharge groundwater, which in turn can contribute to potential base flow.Details of the model processes can be found in Table 1 [26,52,53].
a i,m = mean rain depth on a rainy day at pixel i on month m n i,m = number of events at pixel i in month m P i,m = monthly precipitation for pixel i in month m CN i = curve number for pixel i E 1 = exponential integral function QF stream,m = monthly quick flow of stream pixel in month m QF i,m = monthly quick flow of non-stream pixel i in month m QF i = annual quick flow of pixel i S i = potential maximum soil moisture retention after runoff begins Local Recharge L avail,i = max(γL i , 0) L i = annual local recharge of pixel i P i = annual precipitation of pixel i P i.m = monthly precipitation of pixel i AET = actual evapotranspiration PET = potential evapotranspiration ET 0,i,m = reference evapotranspiration for month m L avail,i = the available recharge to pixel i p ij ∈ [0, 1] is the proportion of flow from cell i to j K c,i,m = monthly crop factor α m = the fraction of upslope annual available recharge that is available in month m β i = the fraction of the upgradient subsidy that is available for downgradient evapotranspiration γ = the fraction of pixel recharge that is available to downgradient pixels Base Flow B sum,j L sum,j −L j i f j is a nonstream pixel B sum,i = L sum,i ∑ j p ij i f j is a stream pixel B i = base flow of pixel i L sum,i = cumulative upstream recharge of pixel i j ∈ {cells to which cell i pours} 1 for pixels adjacent to stream 2 for pixels not adjacent to the stream Water 2018, 10, 88 5 of 25

Data Sources and Processing
The inputs for the Baseline scenario using the seasonal water yield model include spatial GIS layers and biophysical parameters.Among the spatial layers, monthly precipitation and reference evapotranspiration were provided by the climate database of Australian Bureau of Meteorology [26].The Digital Elevation Model (DEM) layer was derived from the SRTM-derived 3 Second Digital Elevation Model Version 1.0 [54], which provides better quality and consistency compared to the original SRTM data [55].Four hydrologic soil groups (HSG A, B, C and D) were classified based on the CSIRO dataset of saturated hydraulic conductivity [56] for assessing infiltration capability of different soil, which represent low to high runoff potential from HSG A to D. The values for the Land use land cover (LULC) were reclassified from the QLUMP (Queensland Land Use Mapping Program) land use classification maps [57].The QLUMP land use mapping methodology integrated satellite images together with field survey and local expert knowledge and provides the most reliable and accurate LULC for the region.We reclassified the mined land area based on a previous research [58,59] into four types including bare land, industrial, grazing and water, where the biophysical parameters were the same as the QLUMP classifications (Table 2).All the spatial layers were processed in ArcGIS [60] and resampled to the same resolution of 90 m based on a majority for the categorical data and average values for the quantitative data.
The parameters, runoff curve number, monthly rain events, threshold flow accumulation and crop factor K c values were derived from a range of sources including previous studies and weather station data.The curve number is related to land use and soil attributes such as soil type and infiltration capability [61].The runoff curve numbers (Table 2) were determined by an Australian developed hydrological software CatchmentSIM [62] and related literature [63,64].The monthly rain events values were the numbers of rain events higher than 0.1 mm [26], which were derived from the long term weather station observation data in the study area [65].The threshold flow accumulation is the number of upstream cells that must flow into a cell before it is considered part of a stream, used to classify streams in the model.According to the InVEST handbook [26], previous local research [66] and the stream network dataset from local government [67], we assumed a flow contribution area of 1 km 2 .The monthly K c (Table A1, Appendix A) were calculated from remote sensing data of fPAR (fraction of Photosynthetically Active Radiation).The calculation method uses the following equations [26,[68][69][70][71], where FC is the Fraction Cover and LAI describes the Leaf Area Index (Equation ( 1)).The fPAR data was extracted from the Australia fPAR dataset version 5 [72].

Uncertainty Assessment Overview
Besides the Baseline scenario, six uncertainty scenarios were developed, where a single input for the ES model was altered for each scenario (Table 3).The alternative uncertainty scenarios represent reasonable input parameterizations or modelling methods, however, less precise and/or with greater uncertainty.For example, they might represent the type of data that might be used in data poor regions.Through the six uncertainty scenarios, we assessed three types of uncertainty including: (1) spatial data sources, (2) modelling scales (temporal and spatial) and (3) parameterization and model selection.

Digital Elevation Model (DEM)
This scenario utilized a DEM layer derived from 1:100,000 contours map sheets surveyed by the AUSLIG (Australian Surveying and Land Information Group) [73].Topography is the key factor for determining stream networks and other hydrological characteristics of a catchment [74,75].The scenario aimed to characterize uncertainties caused by different data sources and data precision.This source is different to the SRTM derived DEM used in the baseline as it is less precise and has a coarser map scale than the baseline.In addition, it only delineated the original natural topography for the whole basin and does not include anthropogenic changes such as those associated with mining disturbance.

Land Cover Map
This scenario aimed to assess the uncertainties associated with land use and land cover classification system mapping.For this scenario the National Dynamic Land Cover (NDLC) [76] was used instead of QLUMP.The NDLC Land cover was reclassified for the modelling of ecosystem services to fifteen land cover types (Table 2).The NDLC dataset is a national-scale land cover dataset for Australia, produced by Geoscience Australia and the Australian Bureau of Agriculture and Resource Economics and Sciences.Compared with the QLUMP map used in the baseline, there are a greater number of detailed vegetation cover classifications, which are useful to identify changes in vegetation cover and extent, however, as it is a national scale map it is less accurate thematically and spatially.The K c values for each land cover (Table A2, Appendix A) and the other parameters were calculated using the same methods as the Baseline scenario.

Climate Data Temporal Scale
This scenario substitutes monthly data for precipitation and evapotranspiration with annual data to assess the impact of the temporal scale of the climate data.Under this scenario, we still utilized the seasonal water yield model but assumed that only annual climate data were available.It is common for water yield modelling to be based on annual data.But for areas where the climate conditions have distinctive seasons, using annual data and models may provide inaccurate assessment of water yield.For example, high intensity rainfall over short periods can saturate soils resulting in greater overland flow and less infiltration than the same amount of rainfall over longer periods.
The monthly inputs into the seasonal water yield model were calculated by dividing annual data into twelve equal values.This scenario is similar to the Annual model scenario in that it considers only annual climate data, however, the annual model uses a different method for modelling water yields and includes different input data i.e. plant available water fraction and root depth.Thus, we could assess differences associated with ignoring seasonal variability without the assessment being confounded by using a different model with different inputs.This scenario tested the impact of spatial scale by resampling all the baseline data by 10.This scenario aimed to assess the impact of spatial scale on ecosystem services mapping and modelling.Assessing spatial scale-dependencies of spatially explicit phenomena is one of the central questions in landscape ecology and geographic research [77,78].

Annual Model
The annual water yield model was applied in this scenario instead of the seasonal water yield model.As described above in Section 2.4.3, it is common for water yield modelling to use only annual data.In fact, earlier versions of InVEST (i.e., version 3.2.0)only included the annual water yield model.The annual water yield model is based on the Budyko curve and annual average precipitation.The Budyko curve simulates the empirical relation between evapotranspiration, potential evapotranspiration and precipitation [79,80].The basic process of the water yield model is to calculate the water run off for every pixel as the precipitation minus evapotranspiration.It assumes the water yield can reach specified sites through relative flow pathways.Then the sub-watershed level water yield is summed up from the pixel-scale results.The equations for annual water yield and the parameters are listed in Table 4 [26,52,81,82].While the outputs of annual model include pixel-scale water yield map and catchment-scale annual total water yield volumes, these pixel-scale results are only for calibration purpose but not for scientific interpretation or decision-making [26] and thus only the catchment-scale values were analysed under this scenario.For the spatial layer inputs, the root restricting depth layer was used as a proxy for soil depth due to data availability.The soil depth and plant available water content were all obtained from the soil attributes of ASRIS (Australian Soil Resource Information System) datasets.The annually precipitation and evapotranspiration layers were derived from the database of Australian Bureau of Meteorology.For parameters (Table A3, Appendix A), the annual K c value for each land use type was calculated based on the monthly K c value by InVEST K c calculator [26], where the monthly K c values were also calculated from fPAR values.The root depth values were based on advice from experts and the literatures [52,83] due to a lack of field data.

K c
This scenario used K c values derived from literature (Table A4, Appendix A), which represents a common situation when local field data are scarce.The literature-based monthly K c value were determined on the basis of the plants' K c values at different growth stage from FAO data [52] with the Australian growing season being integrated.Of the K c values derived for different land covers, the K c value of cropping land were defined based on the value of several major crop types [57], such as cereal, cotton, hay and silage, while the horticulture was characterized by the value of citrus, vine and nuts fruits.

Comparison of Scenario Outputs
This study compared the spatial distribution and catchment scale values of quick flow, base flow and local recharge to analyse the effect of uncertainty on water provision services.Firstly, the spatial patterns of water provision were analysed based on quick flow and base flow maps and value percentile distributions were calculated in the 'Raster' package of R 3.2.0[84].Secondly, we looked at differences in particular regions associated with the main land uses (i.e., mining).The land type boundaries of four dominant land use including grazing, forestry, cropping and mining were extracted based on the QLUMP land cover classification.A simple Regional Value Change (RVC) index (Equation (2)) was used to quantify the differences to the baseline.Finally, this study calculated the annual total water yield volume to compare the annual catchment-scale differences associated with each of the scenarios.This was the only time the annual model was included in the assessment.For all scenarios except for the annual model, we used the ArcGIS raster calculator and zonal analyse tool to calculate the annual total volume of quick flow and local recharge and then summed them up to calculate the annual total water yield.
where M uj and M bj are the mean values of an ecosystem service in j land cover under the uncertainty scenario u and Baseline scenario b, respectively.This index is normalized with a value range from 0 to 1, where higher value indicates greater differences.

Spatial Distributions of Water Provision Service
In this section, we first analysed the quick flow outputs by comparing the Baseline scenario to the results for the uncertainty scenarios.Under the Baseline scenario, the quick flow maps (Figure 2) showed that the highest pixel value was 1662 mm, while the median value was only 59 mm (Figure 3).The values between the first quartile and the third quartile ranged from 28 mm to 92 mm.The highest value was located on the stream pixels in the northeast edge of the region.Besides the stream network, there were large apparent differences in spatial distribution across the region.Large areas with high quick flow values were found in the north-eastern and eastern edge and scattered across the cropping land and the mining footprint.Besides these high value areas, other regions had comparatively lower values.
A qualitative comparison between the Baseline and uncertainty scenarios (Figures 2 and 3) showed that the spatial distribution of quick flow for all the uncertainty scenarios apart from spatial scale were almost the same as the Baseline.High values were mainly distributed along the stream network and in the north-eastern part, the values of other areas were comparatively low.The stream quick flow of Spatial scale scenario was uniquely different to the Baseline and other uncertainty scenarios in terms of spatial distribution (Figure 2) and overall values (Figure 3).
The relative spatial differences of pixel values between uncertainty scenarios and Baseline scenarios for quick flow are shown in Figure 4.The DEM scenario showed the largest differences; especially for the stream pixels, with the maximum positive value of 1472 mm and the minimum negative value of −1473 mm, however these differences were confined to a few pixels.The Land cover map scenario showed much more heterogeneous differences i.e. higher in cropping areas while lower in the mining footprint area.The difference associated with the Climate data temporal scale scenario were less than the above two scenarios with values ranging from −127 mm to 123 mm, which mainly occurred in the eastern part of the region.For the Spatial scale scenario maps, the pixel values in most of the region were higher, among which the highest value was 1467 mm.In contrast to the other scenarios, there was no difference between the Baseline and K c scenario for quick flow.
Water 2018, 10, 88 11 of 25 of the region were higher, among which the highest value was 1467 mm.In contrast to the other scenarios, there was no difference between the Baseline and Kc scenario for quick flow.In terms of distribution of the pixel values for the quick flow maps, the maximum values were similar to each other and the Baseline scenario with the highest value of 1666 mm for the Climate data temporal scale scenario and the lowest value of 1649 mm for the Spatial scale scenario (Figure 2).However, the value for the 25th percentile and 75th percentile was much lower compared with the maximum value and varied between scenarios (Figure 3).The range between the 25th percentile and 75th percentile for the Spatial scale scenario was the highest at 52 mm to 656 mm respectively.While the Land cover map scenario was the lowest with values of 17 mm to 70 mm for the 25th and 75th percentile.The DEM, Climate data temporal scale and Kc scenarios had similar ranges to the baseline.
The base flow maps (Figure 2) also showed large differences in the spatial patterns between scenarios.The spatial distribution of pixel values for the DEM scenario and Spatial scale scenario overall were approximately the same as the Baseline scenario, where high values were mainly located at the eastern edge, mining footprints and central-southern area.Though at fine scales there was much greater heterogeneity in the values of the Spatial scale scenario.The Climate data temporal scale and Kc scenarios had only very few areas with comparatively high values.Besides the differences in spatial pattern, pixel value distributions also showed differences (Figure 3).A maximum value of 1103 mm for the Land cover map scenario was 162 mm higher than that of the Baseline scenario.While the maximum values of other scenarios were all much less than the 940 mm of the Baseline scenario.The value ranges between the 25th and 75th percentile were much lower than the maximum values, among which the first quartile values of all scenarios were 0; except for the Spatial scale scenario with a value of 9.77 mm.The values of Kc scenario and Climate data temporal scale scenario were so small that even the third quartile values were 0. The differences in pixel values for base flow are shown in Figure 4.The DEM scenario and Land cover map scenario had unique patterns of differences compared to the other scenarios.The Climate data temporal scale scenario and Kc scenario showed smaller values at the eastern boundary (coast), mining footprint and central southern area.The Spatial scale scenario tended to have higher pixel values across the region with the highest value of 554 mm.In terms of distribution of the pixel values for the quick flow maps, the maximum values were similar to each other and the Baseline scenario with the highest value of 1666 mm for the Climate data temporal scale scenario and the lowest value of 1649 mm for the Spatial scale scenario (Figure 2).However, the value for the 25th percentile and 75th percentile was much lower compared with the maximum value and varied between scenarios (Figure 3).The range between the 25th percentile and 75th percentile for the Spatial scale scenario was the highest at 52 mm to 656 mm respectively.While the Land cover map scenario was the lowest with values of 17 mm to 70 mm for the 25th and 75th percentile.The DEM, Climate data temporal scale and K c scenarios had similar ranges to the baseline.
The base flow maps (Figure 2) also showed large differences in the spatial patterns between scenarios.The spatial distribution of pixel values for the DEM scenario and Spatial scale scenario overall were approximately the same as the Baseline scenario, where high values were mainly located at the eastern edge, mining footprints and central-southern area.Though at fine scales there was much greater heterogeneity in the values of the Spatial scale scenario.The Climate data temporal scale and K c scenarios had only very few areas with comparatively high values.Besides the differences in spatial pattern, pixel value distributions also showed differences (Figure 3).A maximum value of 1103 mm for the Land cover map scenario was 162 mm higher than that of the Baseline scenario.While the maximum values of other scenarios were all much less than the 940 mm of the Baseline scenario.The value ranges between the 25th and 75th percentile were much lower than the maximum values, among which the first quartile values of all scenarios were 0; except for the Spatial scale scenario with a value of 9.77 mm.The values of K c scenario and Climate data temporal scale scenario were so small that even the third quartile values were 0. The differences in pixel values for base flow are shown in Figure 4.The DEM scenario and Land cover map scenario had unique patterns of differences compared to the other scenarios.The Climate data temporal scale scenario and K c scenario showed smaller values at the eastern boundary (coast), mining footprint and central southern area.The Spatial scale scenario tended to have higher pixel values across the region with the highest value of 554 mm.

Value Comparisons for Particular Land Uses
To quantify the impact of uncertainty scenarios for different land uses, the RVC index was utilized to calculate the base flow and quick flow differences between the baseline and uncertainty scenarios.The RVC index value was normalized with the values ranging between 0 and 1 and a higher value indicating a greater difference.Figure 5 showed the RVC indices of different scenarios for the four main land use types which occupied more than 99% of the total area and the mean RVC value of all land use types (including all the main and small covers).For the quick flow, the RVC values of Spatial scale scenario were the highest in all the different areas, among which the maximum value was 0.80 in forestry.The highest value of the other scenarios in forestry was only 0.03 for the Climate data temporal scale scenario, which was much lower than the value for the Spatial scale scenario.For other areas not including forestry, the RVC values of Land cover map scenario were the second highest after the Spatial scale scenario and much higher than the other scenarios.The lowest value was 0 in K c scenario for all of the four land uses.

Value Comparisons for Particular Land Uses
To quantify the impact of uncertainty scenarios for different land uses, the RVC index was utilized to calculate the base flow and quick flow differences between the baseline and uncertainty scenarios.The RVC index value was normalized with the values ranging between 0 and 1 and a higher value indicating a greater difference.Figure 5 showed the RVC indices of different scenarios for the four main land use types which occupied more than 99% of the total area and the mean RVC value of all land use types (including all the main and small covers).For the quick flow, the RVC values of Spatial scale scenario were the highest in all the different areas, among which the maximum value was 0.80 in forestry.The highest value of the other scenarios in forestry was only 0.03 for the Climate data temporal scale scenario, which was much lower than the value for the Spatial scale scenario.For other areas not including forestry, the RVC values of Land cover map scenario were the second highest after the Spatial scale scenario and much higher than the other scenarios.The lowest value was 0 in Kc scenario for all of the four land uses.For the base flow for all of the four land use areas, the RVC values of Kc and Climate data temporal scale scenarios were much higher than the values of the other three scenarios, among which the highest value was 0.98 for mining footprint area under Climate data temporal scale scenario.The value distribution of the other three scenarios for the mining area ranged from 0.02 to 0.17, which were much lower than those for cropping, forestry and grazing area.The Land cover map scenario produced the lowest RVC value of 0.02 in mining areas while the lowest value for the other land uses were found in the DEM scenario.The RVC values of all land uses under Kc and Climate data temporal scale scenarios were also much higher than for the other scenarios.

Catchment-Scale Annual Total Water Yield
Annual total water yield represents water provisioning at the whole catchment scale.The water yield amount for the Annual model scenario was a direct output from the annual water yield model, while the value for other scenarios using the Seasonal water yield model were the sum of quick flow and local recharge.Figure 6 represents the annual total water yield volume of all seven scenarios.The total volume of water yield for the Baseline scenario was 1.69 × 10 10 m 3 and was less than the values For the base flow for all of the four land use areas, the RVC values of K c and Climate data temporal scale scenarios were much higher than the values of the other three scenarios, among which the highest value was 0.98 for mining footprint area under Climate data temporal scale scenario.The value distribution of the other three scenarios for the mining area ranged from 0.02 to 0.17, which were much lower than those for cropping, forestry and grazing area.The Land cover map scenario produced the lowest RVC value of 0.02 in mining areas while the lowest value for the other land uses were found in the DEM scenario.The RVC values of all land uses under K c and Climate data temporal scale scenarios were also much higher than for the other scenarios.

Catchment-Scale Annual Total Water Yield
Annual total water yield represents water provisioning at the whole catchment scale.The water yield amount for the Annual model scenario was a direct output from the annual water yield model, while the value for other scenarios using the Seasonal water yield model were the sum of quick flow and local recharge.Figure 6 represents the annual total water yield volume of all seven scenarios.The total volume of water yield for the Baseline scenario was 1.69 × 10 10 m 3 and was less than the values 5.12 × 10 10 m 3 of Spatial scale scenario and 2.8951 × 10 10 m 3 of Annual model scenario.While the total volume from Climate data temporal scale (1.49 × 10 10 m 3 ), DEM (1.62 × 10 10 m 3 ), Land cover map (1.56 × 10 10 m 3 ), K c (1.47 × 10 10 m 3 ) scenarios were less than the Baseline scenario with the K c scenario producing the lowest water volume.

Impact of Uncertainty on Water Provision
The Baseline scenario outputs represented the current water provision service condition of the study area.As previously observed in a previous study area [85], the average water flow volume of the Fitzroy River and Isaac river sub-catchments were much higher than other sub-catchments, which was in accordance with the general spatial patterns in our study.Similarly, erosion research in this area [86][87][88] showed water quantity spatial distributions comparable to our study with values of both the quick flow and base flow high in the eastern sections.The reason for this spatial distribution is due to higher elevation and topographical heterogeneity and high annual precipitation in these regions [74,[89][90][91] which generate higher stream flow and local recharge.In the mining area which was mainly occupied by bare pits, coal stockpile and industrial plants, the poor permeability (high runoff curve number) resulted in high quick flow [63], while the low evapotranspiration affected the recharge to base flow.In addition, the comparatively high quick flow of cropping land was also a result of hydrologic conditions associated with high runoff curve numbers, which are always accompanied by large amount of soil loss [59,92].
The scenarios analysis demonstrated how uncertainties derived from the mapping and modelling methods affects ecosystem service modelling outputs and that the sensitivity of different outputs (i.e., base flow versus quick flow) varied with the uncertainty scenario.The uncertainties in quick flow associated with plant evapotranspiration (Kc scenario) were minor compared to the baseline.While changes to topography under the spatial scale and DEM scenario lead to different stream directions and water flow values.Of the uncertainty factors assessed, a coarsening of spatial scale (Spatial scale scenario) had the greatest effect on both overall values, spatial distribution of values and even the characterization of stream direction which corresponds with findings of other

Impact of Uncertainty on Water Provision
The Baseline scenario outputs represented the current water provision service condition of the study area.As previously observed in a previous study area [85], the average water flow volume of the Fitzroy River and Isaac river sub-catchments were much higher than other sub-catchments, which was in accordance with the general spatial patterns in our study.Similarly, erosion research in this area [86][87][88] showed water quantity spatial distributions comparable to our study with values of both the quick flow and base flow high in the eastern sections.The reason for this spatial distribution is due to higher elevation and topographical heterogeneity and high annual precipitation in these regions [74,[89][90][91] which generate higher stream flow and local recharge.In the mining area which was mainly occupied by bare pits, coal stockpile and industrial plants, the poor permeability (high runoff curve number) resulted in high quick flow [63], while the low evapotranspiration affected the recharge to base flow.In addition, the comparatively high quick flow of cropping land was also a result of hydrologic conditions associated with high runoff curve numbers, which are always accompanied by large amount of soil loss [59,92].
The scenarios analysis demonstrated how uncertainties derived from the mapping and modelling methods affects ecosystem service modelling outputs and that the sensitivity of different outputs (i.e., base flow versus quick flow) varied with the uncertainty scenario.The uncertainties in quick flow associated with plant evapotranspiration (K c scenario) were minor compared to the baseline.While changes to topography under the spatial scale and DEM scenario lead to different stream directions and water flow values.Of the uncertainty factors assessed, a coarsening of spatial scale (Spatial scale scenario) had the greatest effect on both overall values, spatial distribution of values and even the characterization of stream direction which corresponds with findings of other studies such as [15,75].For base flow, spatial patterns and values changed with all uncertainty scenarios except for the DEM scenario.Differences due to the Climate data temporal scale, K c and Land cover map scenarios show that uncertainty associated with precipitation, evapotranspiration, plant factors and land classification schemes affect modelling outcomes.These factors have also been recognized as key variables affecting relevant ecosystem services in other studies [24,93,94].A critical and interesting result of this study was from the Climate data temporal scale scenario which showed that ignoring distinctive wet and dry seasons by using an annual model are likely to produce model outputs that are less likely to reflect true values.
Another interesting observation from the analysis were that literature-based parameterization resulted in very different modelling outputs than direct measurements.For example, for the K c scenario, the literature-based K c of forest was between 0.5 and 1.2 [52] while the remote sensing measured K c values ranged from 0.56 to 0.80.Reason for this difference is that the vegetation characteristics of the forests in the study area are characterized by sparse vegetation in contrast to the less open forests from the locations in the literature [95,96].In addition, the literature-based K c value were from more developed locations with less vegetation and thus had lower K c values of 0.2 for industrial land and 0.3 for urban residential area [26,52].
Of all the scenarios, the Spatial scale scenario had the greatest overall impact.Spatial scale affects climate, land cover, plant parameters and the DEM.The large apparent impact of spatial scale scenario for both quick flow and base flow demonstrates, as shown in many studies, that scale is a critical driver of uncertainty in spatially explicit modelling [15,39,97].As well as impacts on spatial patterns of quick and base flow, this study also developed and utilized the RVC index to quantify water provision changes under the uncertainty scenarios for specific land uses.The RVC analysis of quick flow showed that the Land cover map and Spatial scale scenarios had higher average RVC values than other scenarios, which are likely to be caused by the spatial variability in the runoff curve numbers and DEM [98,99].In contrast for the base flow, the average RVC values for K c and Climate data temporal scale scenario had much higher values than for other scenarios.Similar trends were found with the RVC assessment in the qualitative analysis of spatial patterns with distinct RVC values in different areas.
The final assessment based on the annual total water yield of a catchment represents the water provision ability and hydrological conditions of the catchment [100][101][102].Previous studies in the region on the annual water yield volume were mainly from gauged stream flow during wet season [85,87].These studies provide water yield volumes based on distinctive weather conditions for specific years as opposed to average climate condition used in our modelling.For example, Jones et al. gauged total stream flow volume were 2.16 × 10 10 m 3 in 2012/13 and 4.96 × 10 10 m 3 in 2013/14 [85], which are similar to the water yield value for the Baseline scenario of 1.69 × 10 10 m 3 .However, under different uncertainty scenarios such as spatial scale the values were much higher.The differences between the uncertainty scenarios and gauged total stream flow volume were in the same order of magnitude suggesting that even the outputs from the uncertainty scenarios are likely to provide some useful information.

ES Mapping Uncertainty in the Context of the Fitzroy Basin and Mining
This study used the water provision mapping in the Fitzroy Basin mining region as a case study to quantify uncertainties in ES mapping.The differences in outputs associated with each uncertainty scenario reflect the result of using models, parameterization and data that could reasonably be used for mapping ES.Our study found relatively large differences in model outputs due to uncertainty, demonstrating the need to identify and categorize uncertainty sources when conducting analyses [33,35].The consideration, identification and categorization of uncertainties are likely to be different and specific to an application due to differences in data sources, catchment characteristics etc. [103][104][105].The choice of uncertainty methods assessment selection are an important consideration as uncertainty assessment can be complex and time-consuming [33].This study applied scenario analysis in an uncertainty assessment, which is common and considered to be an effective approach [35,38,41].However, as spatial data and modelling for ES mapping is complex with almost an infinite source of uncertainty, it is difficult to characterize all sources of uncertainties in the spatial attributes of ES.Spatial data are more complex than parametric values and different observation scales may lead to different outcomes and the accuracy of each pixel may be distinctive [33].
The application of uncertainty analysis such as for assessing mining impacts is difficult from the perspective of decision makers.Some researchers have argued that most uncertainty analysis do not assist decision-makers due to the non-practical outputs [33,106,107].Our study offers spatially explicit maps and statistical diagrams, which provide decision-makers and stakeholders with explicit information on how these uncertainties are distributed.
The information provided by ES studies such as ours need to be interpreted with caution due to the inherent uncertainty associated with modelling methods.Our study provides information on water provision conditions such as hotspot areas and annual volumes that could potentially improve hydrological ecosystem management for management of Great Barrier Reef catchments and for assessing impacts associated with economically valuable minerals production [59,108,109] as water provision could affect either the reef environment or the large water demands from the mining industry in the region [110,111].There needs to be and there is a large focus on planning for water resources and regional ecological sustainable development in the Fitzroy Basin [50,86,112].However, management also need to consider uncertainty and the limitations associated with modelling outputs used to inform decision making.

Conclusions
This study mapped water provision ES service, quantifying uncertainties within the ES mapping process in one of the largest mining resources region in Australia.This study demonstrated the diversity of uncertainties associated with ES mapping and the many ways in which uncertainty can affect mapping outputs from spatial patterns to total water provision.The results showed that varying mapping scales had the highest impact on uncertainty.However, most forms of uncertainty led to changes in ES mapping outputs i.e. data source, parameters and mapping methods.This study provides a useful demonstration of the impact of uncertainty that is likely to be present in many ES mapping applications, prompting greater need for further research on uncertainty in ES mapping.

Figure 1 .
Figure 1.Location and land use map in the study area.

Figure 1 .
Figure 1.Location and land use map in the study area.
water yield for each pixel AET xj = Annual actual evapotranspiration on the pixel x with LULC j P x Annual precipitation on the pixel x R xj = The ratio of potential evapotranspiration to precipitation w x = Modified dimensionless ratio of plant accessible water storage to expected precipitation AWC x = The volumetric (mm) plant available water content ETo x = The reference evapotranspiration from pixel x K xj = Plant (vegetation) evapotranspiration coefficient ETo x = Reference evapotranspiration K c ( x )= The evaporation factor for each LULC

Figure 2 .
Figure 2. Quick flow and Base flow maps for Baseline and uncertainty scenarios.

Figure 2 .
Figure 2. Quick flow and Base flow maps for Baseline and uncertainty scenarios.

Figure 4 .
Figure 4. Relative difference between Baseline and uncertainty scenarios.Figure 4. Relative difference between Baseline and uncertainty scenarios.

Figure 4 .
Figure 4. Relative difference between Baseline and uncertainty scenarios.Figure 4. Relative difference between Baseline and uncertainty scenarios.

Figure 5 .
Figure 5. RVC values of four dominant land use area.

Figure 5 .
Figure 5. RVC values of four dominant land use area.

Table 1 .
Inputs and description of Seasonal Water Yield Model.

Table 2 .
Land use and land cover classification for hydrological soil groups (HSG), land types and relevant runoff curve numbers.

Table 3 .
Baseline and uncertainty scenarios with alternative spatial data source or modelling scales or parameterization and model tested for a particular scenario in bold.

Table 4 .
Inputs and description of Annual Water Yield model.

Table A2 .
Monthly K c values for each land classification under Land cover map scenario.

Table A3 .
Biophysical parameters for Annual model scenario.

Table A4 .
Monthly K c values for each land classification under K c scenario.