The relationship between ecosystem services and human modification displays decoupling across global delta systems

The ties between a society and its local ecosystem can decouple as societies develop and replace ecosystem services such as food or water regulation via trade and technology. River deltas have developed into important, yet threatened, urban, agricultural and industrial centres. Here, we use global spatial datasets to explore how 49 ecosystem services respond to four human modification indicators, e.g. population density, across 235 large deltas. We formed bundles of statistically correlated ecosystem services and examined if their relationship with modification changed. Decoupling of all robust ecosystem service bundles from at least one modification indicator was indicated in 34% of deltas, while 53% displayed decoupling for at least one bundle. Food-related ecosystem services increased with modification, while the other bundles declined. Our findings suggest two developmental pathways for deltas: as coupled agricultural systems risking irreversible local biodiversity loss; and as decoupled urban centres externalising the impact of their growing demands. The direct relationship between societal development and local ecosystem services breaks down at relatively minor levels of human modification of large river delta landscapes, according to a statistical analysis of 235 deltas.

N ature provides humanity many benefits, yet as societies develop they increasingly turn to global markets and engineered solutions to meet their needs 1 . Societal development can therefore weaken the ties between people and their local ecosystems, and disconnect resource extraction from resource use, a process called decoupling 1,2 . Ecosystem services (ES) represent the contribution of ecosystems to human welfare and Earth's life support system 3,4 . The decoupling of societies from their local ES has numerous consequences. For example, international trade means that food 5 and fish 6 production and prices are often based on global demands rather than local stocks. Decoupling economic growth from natural resource consumption in many richer countries may appear sustainable 7 , yet can externalise environmental costs to the poorest people and nations 8 . Decoupling can make local resource governance more problematic 6 , and can reduce resilience, making it more difficult for systems to deal with shocks 9 such as extreme climate events. In extremis, decoupling can lead to system-wide degradation, as hypothesised for past civilizations such as the Roman Empire where expanding trade networks allowed urban centres to grow beyond their local limits, threatening the whole system with collapse under broader food production problems 10 .
We propose an empirical approach to identify decoupling using the characteristic response curves 11 between human societal development and local ES. To do this, we build on the ecosystem service cascade framework: ecosystem properties and functions supply potential ES, which become services when used, and then provide benefits to society (Fig. 1a) 12 . These benefits can enhance societal development, which in turn modifies natural systems and the processes that enable ecosystem service supply, improving some ES, while degrading others (Fig. 1b i. and ii.). Human modification can improve some ES through, for example irrigation and fertiliser inputs to increase food supply 4 . Modification can equally replace, pollute, or otherwise degrade local ES 4 . However, the direct relationship between societal development and local ES can be broken in wealthier societies which can develop beyond the supply capacity of the local ecosystem by replacing or importing ES (Fig. 1b iii.). This response was proposed between agricultural ES and land-use intensity, ES initially increasing before declining beyond a certain level of modification 4 . Such change points have also been identified in the relationships between ES, population and economic indicators in urban areas 13 . We posit that we can identify decoupling using a change point in the response curve between societal development and local ES, which indicates the point at which systems persist independent of their local ES. The shape of this response curve indicates the vulnerability of ES to modification and allows us to infer thresholds of modification beyond which ES decouple. Fig. 1 Graphical representations of decoupling and its mechanisms. a Decoupling of the ecosystem service cascade. The traditional ecosystem service cascade flows from ecosystem property or function to benefit 12 ; these ecosystem services (ES) can maintain and promote human societal development which can then enhance or degrade the ecosystem properties on which ES are based. Decoupling disconnects societal development from this cascade, as local ES are replaced with external inputs. b Characteristic response curves of ES to the modification of natural systems by societal development. i. Positive response of ES to modification in a coupled system. ii. Negative response of ES to modification in a coupled system. iii. Weakening of responses (or feedbacks) between modification and ES beyond a change point in a decoupled system. c Mechanisms through which the modification of river deltas can impact local ES. The mechanisms and impacts on ES of an increase in population, the disruption of water flow through dams and reservoirs, human appropriation of net primary productivity (HANPP), and built infrastructure. Impacts can affect both supply, the potential ES available from nature; and service, the realised ES used by people. Positive and negative signs indicate, respectively, an increase or decrease in the component. Different types of ES may respond, and decouple, differently to modification. The different categorisations of ES could therefore potentially inform their responses. The foundational conceptual framework of ES, based on the Millennium Ecosystem Assessment (MEA), divides ES into provisioning goods, regulating flows, services supporting the others, and cultural services 14 . The Intergovernmental Platform on Biodiversity and Ecosystem Services framework updates these contributions into the material, nonmaterial, regulating, and the nature that supports them; as well as intrinsic, instrumental and relational valuations 15 . We might expect provisioning or material ES, such as food, to respond positively to modification, while regulating services respond negatively. ES can also be represented, as above, as a cascade from ecosystem property, to supply, to service 12 ; different positions along this cascade may decouple at different levels of modification. However, using a priori categorisations can introduce problems-such as which categorisation to follow 16 or how to classify ES fitting into multiple categories, and within category trade-offs may weaken or cancel the response of combined ES 17 . Using multiple approaches to group ES helps overcome the limitations of any individual categorisation, allowing a more robust assessment of the response of ES to modification and in turn, decoupling.
The development or degradation of one ecosystem service can affect others 17,18 , meaning that synergies and trade-offs, or positive and negative correlations between ES, are critical in understanding the diverse impacts of modification 19,20 . For example, the loss of a wetland habitat would simultaneously degrade flood control, water filtration, and fish stocking ES 21 . In this way, the effects of modification can interact, escalate and spread beyond the individual ecosystem service indicator assessed. Data-driven unsupervised methods are widely used to group ES 22,23 . These approaches potentially provide more generalisable results by incorporating the correlations between individual ES to establish ES bundles 17 , sets of consistently associated ES cooccurring in space and time 22 . Using bundles and examining their synergies and trade-offs can therefore provide further insights on the effect of modification on ES.
River deltas present a globally important system, well-suited to assess decoupling given their range of ES and gradient of modification. Deltas are rich in ES with fertile land for agriculture 24 , as well as rich fisheries 25,26 . Internationally rare and valuable ecosystems, such as mangroves, freshwater swamps, and saltwater marshes offer important habitats and ES 25,27 . The unique supply of terrestrial, coastal and marine ES have allowed deltas to become hotspots of human development, from the cradles of civilization to modern mega-cities, currently housing 340 million people 24,28 . This means they reflect a gradient of modification from relatively pristine, natural systems to some of the most modified places on Earth 29 . Increasing human modification of deltas endangers both the ES they provide, and their long-term sustainability 29,30 . Degraded flood protection and water purification ES may require engineered solutions, which can be unsustainable and, in many deltas, economically impossible 31 . Given the significance of deltas as urban, industrial and agricultural centres, decoupling in these systems has local and global implications 30 .
The term decoupling can refer to the disconnection between deltas and their contributing basins, here we use it to refer to the disconnection between human modification and local ES within delta systems. Human modification impacts delta ES in several ways (Fig. 1c), we analyse four of its most important indicators: population density, upstream influences, ecosystem flows and built infrastructure. Population density represents human occupancy and ES demand, identified as the principal predictor for global modification 32 , and rising in increasingly urbanised deltas 33 . Upstream influences, including water and sediment flow are critical to delta functioning-interception by dams and reservoirs affects water quality 34 , fish recruitment 35 , floodplain habitats and fluvial state 36 . We measure upstream influences using flow disruption, the additional time water takes to reach the river mouth due to dams and reservoirs 37 . Humans also exploit and modify delta ecosystem flows. This modification can be measured by human appropriation of net primary productivity (HANPP), the proportion of natural productivity used by harvest, burning, and other land use 38 . Finally, built human infrastructure in deltas is increasingly replacing the ecosystems supplying many ES 39 . We measure this using human footprint 40 , an integrative indicator, which we use to combine measurements of the urban area, roads, railways and electrical infrastructure. In sum, modification in delta systems may challenge ecosystem service persistence, and broader delta sustainability.
Here, we examine the extent to which societies have become decoupled from their local ES across delta systems due to the modification of natural processes. We investigate: (i) the degree to which ES are decoupled from indicators of modification across delta systems, (ii) the characteristic response curves of ES to modification and how they vary between modification indicators; and (iii) the synergies and trade-offs between ecosystem service bundles and how sensitive they are to modification. To do this we created a dataset of all 235 large (>10 km 2 ) delta distributary networks worldwide. We then bundled 49 robust observed and modelled global delta ecosystem service indicators (Supplementary Table 1) using an unsupervised approach. This approach, hereafter called unsupervised bundling, groups ecosystem service indicators based on their correlation structure and the consistency in their clustering over an ensemble of six unsupervised analyses. We found three robust bundles of ES which we identified as primarily 'intactness', 'productivity' and 'food' related. Strong synergies and trade-offs were present within and between bundles, which weakened as modification increased excepting the persistent trade-off between 'food' and 'intactness'. We contrasted our unsupervised approach with the response curves of individual ES, and bundles created using the MEA framework and ecosystem service cascade and found largely supportive results. We then examined the characteristic response curves of each ecosystem service bundle with our four indicators of human modification: population density, flow disruption, HANPP and human footprint. Decoupling was indicated in eight of twelve response curves, occurring with built infrastructure for every bundle, and evident at relatively low levels of modification.

Results
Delta ecosystem services indicate decoupling as modification increases. For at least one modification indicator 53% (125 of 235) of deltas indicated decoupling (i.e. were beyond a change point in the response curve between modification and ES) for one or more of the robust unsupervised bundles, while 34% (79 of 235) indicated decoupling for all three (Fig. 2, Table 1). The deltas over which we explored this decoupling represented broad gradients of modification and ecosystem service supply (see Supplementary Notes 1-4 for geographical variations in ecosystem service supply and modification, Supplementary Note 5 for statistical summary and Supplementary Note 6 for comparison of ecosystem service diversity with modification). Decoupling was identified by classifying smoothed regression curves. We found five characteristic response curves between ecosystem service bundles and modification indicators, varying from linear to saturating to parabolic, the latter two containing a change point indicating decoupling (Fig. 2; for individual ES see Supplementary Note 2 and for other bundling methods see Supplementary Note 7). Almost 25% of deltas showed decoupling at relatively low levels of human footprint across the unsupervised, MEA and ecosystem service cascade methods (Fig. 2, Table 1 and Supplementary Note 7). Within the bundles, individual ES often showed weaker response curves, but the better fitting (r > 0.1) curves were consistently of the same characteristic shape as the overall bundle (73 of 80 responses, Table 2, Supplementary Note 2).
Unsupervised bundling identified four bundles across all deltas and ES globally, which we then assessed for decoupling. Three of these bundles were generally robust, meaning the ES included in each bundle were consistent across the six clustering methods used (Table 3), and broadly consistent as modification increased (Supplementary Note 8). We labelled these bundles 'intactness', 'productivity', and 'food'. 'Intactness' was the most consistent bundle (ES indicators within were grouped with this bundle 85% of the time), containing biodiversity intactness indicators, soil quality indicators, and the absence of flood risk and water quality deficit. 'Productivity' was the second most consistent bundle (82% consistency), made up of net primary productivity (NPP), carbon storage, and terrestrial and marine species richness indicators.  Table 3 Table 3 for the composition of each bundle. The rightmost column indicates the percentage of deltas showing decoupling for that ecosystem service bundle for one or more modification indicators. The bottom row indicates the percentage of deltas showing decoupling for that modification indicator for one or more ecosystem service bundles. The bottom right number (53.2%) indicates the total percentage of deltas that show decoupling at least once. All thresholds marked indicate a response curve with a pseudo-r > 0.1 (see Fig. 2).
'Food' (76% consistency) included indicators for crop, fish and aquaculture production. The final bundle 'other' was the least robust (50% consistency, less than the 57% average found when random data was clustered), including water-related indicators such as sediment flux, discharge and wetland area alongside oil area and pasture. Given this inconsistency in composition and response, we restrict our analysis to the three robust bundles (see Supplementary Note 7 for more discussion of bundle consistency and appropriateness). Response curves were typically weaker when ES were grouped by the supervised methods (Supplementary Note 7). However, applying unsupervised clustering within these frameworks gave more consistent responses with our unsupervised bundles ( Supplementary Figures 60 and 61). Across the methods, 'food' indicators grouped and responded similarly. For the MEA categories, regulating 'soil and water' and supporting 'intactness' share similar responses to unsupervised 'intactness'; similarly 'richness and productivity' and 'carbon and water' share similar responses to 'productivity'. For the ecosystem service cascade, the property 'intactness' and supply 'soil' responded similarly to unsupervised 'intactness'; 'richness' and 'productivity and water' with unsupervised 'productivity'. Finally, the distinct service 'regulating' bundle showed very weak relationships to modification, mirroring the weak relationships of its constituent ES to modification. Where decoupling occurred, it appears later along the gradient of modification for service-based indicators than supply or property. Based on the similarity of these results, our further analyses are based on unsupervised bundling.
Responses to modification vary with type of modification and ecosystem service. The type and gradient of the response curve varied between the different ES and modification indicators. We found that modification had the greatest negative impact on the 'intactness' bundle, which displayed negative relationships across each of the different modification indicators (Fig. 2). 'Intactness' showed a consistent linear decline as population density and HANPP increased, but halved before decoupling with human footprint in a negative saturating relationship. The strength of these relationships was variable (pseudo r-values from 0.11 < r < 0.56), but weaker for human footprint (r = 0.18) where the relationship decouples; and flow disruption (r = 0.11) which showed little obvious pattern throughout (similar or weaker results were found using another indicator for dam capacity per watershed, Supplementary Note 9). 'Productivity' showed positive then decoupled relationships across all modification types, decoupling at a population density from 10 to 100 people per km 2 , and between 30 and 40% HANPP. The fits of the regression curves for the 'other' bundle, and of the individual ES within it, were very weak (0.00 < r < 0.16), limiting any interpretation of these responses. 'Food' consistently increased with population density and HANPP, an inverse of 'intactness', and displayed the strongest overall fit to the regression curve (0.12 < r < 0.67). 'Food' decoupled with human footprint, but at almost double the level as the other bundles.
Synergies and trade-offs between individual and bundles of ecosystem services. The bundles produced were consistent in the correlations between their constituent ES (Fig. 3a, b). We used pairwise correlations to infer synergies and trade-offs. We found strong positive correlations, i.e. synergies, between the 'food', 'productivity' and 'other' bundles (Spearman's ρ > 0.4, illustrated by blue in Fig. 3), which all had strong negative correlations, or trade-offs, with 'intactness' (ρ < −0.42, illustrated by red). The pairwise ecosystem service correlations were very consistent within the bundles-excepting weak or non-existent trade-offs of forest cover and wetlands with 'intactness' indicators, and generally weaker correlations within the 'other' bundle (Fig. 3a).
Correlations were typically strongest at the lowest stages of modification, weakening as modification increased though broadly remaining in the same direction. The negative 'food' and 'intactness' correlation was the only one that remained strong even at the highest levels of modification (ρ < −0.39, Fig. 3c, Supplementary Note 10).

Discussion
We found evidence for the decoupling of ES across four gradients of human modification in global deltas. Every ecosystem service  We used an ensemble of six clustering methods to identify four bundles of ecosystem services. The number beside each ecosystem service indicates how many times it was clustered in that bundle, out of six methods. The percentage robustness for each bundle is the average times each ecosystem service a bundle contains was selected within it (see Method and Fig. 5 for more details). We also created bundles using ecosystem service cascade and MEA frameworks (see Supplementary Note 7).  bundle displayed decoupling from at least one modification type in 34% of large deltas, while one or more bundles indicated decoupling in 53%. Decoupling often occurred at relatively low levels of modification. The bundle 'intactness' displayed negative responses to modification, while 'food' had a positive response, and all bundles decoupled with built infrastructure. Our results were reasonably consistent between unsupervised and supervised bundling, and between individual indicators and their unsupervised bundles; we will therefore mainly discuss the results of the unsupervised groups. Intactness had a strong trade-off with the other bundles, which shared synergies with one another. As modification increased, bundles appeared less stable around the decoupling threshold but were broadly consistent, and while synergies and trade-offs weakened, the trade-off between 'intactness' and 'food' was persistent across the modification gradients.
We empirically show trends of decoupling across different modification indicators and ecosystem service types. These nonlinear response curves fit the proposed feedbacks between society, nature and ES 1,2 . Such decoupling could be forced by poor environmental conditions limiting ecosystem service supply, or by societies choosing to emphasise different land uses and instead import or engineer services. Implications will depend on the type of ecosystem service and its replaceability, and will vary locally and globally. However, we should remember that we are indicating decoupling statistically, across delta systems, at a snapshot in time. Multiple response curves for individual systems are likely, and other mechanisms may explain the response curves found. However, a change in response with modification, for whatever mechanism, is the most likely explanation for the patterns found. As an example, if monotonic relationships continued, it would imply, paradoxically, that the most populous deltas have the lowest productivity per capita. Decoupling may also emerge from different development pathways at different velocities. Such legacies can lock-in deltas to for example agricultural production 30 requiring increasing artificial inputs to maintain output 9 . Rapid modification may not allow for a full replacement of ES, or a full awareness of the unpredictable hazards they can mitigate. Finally, the directionality of decoupling is uncertain-as illustrated by Fig. 1, human modification and ES to some extent drive one another. The presence of ES allows the migration and settlement which in turn leads to modification. While our methods cannot establish causality, decoupling does indicate that beyond a threshold local ES no longer drive modification, as modification is able to increase while ES decrease or stay constant. Longitudinal studies could better elucidate both the directionality of decoupling, and the mechanisms that drive it.
The implications of decoupling, and the synergies and tradeoffs found, vary between the bundles. The 'food' bundle was a combination of agricultural and aquaculture services, including water withdrawal, which is predominately used by the agricultural sector 41 . 'Food' did not decouple with population, continuing to increase at higher population densities. This coupling indicates that the capacity for agricultural development remains in the majority of deltas, even with the increased competition for land accompanying high population densities. However, such development may risk delta sustainability given the strong trade-off between the 'food' and 'intactness' bundles, which persists through the modification gradient. These results mirror the predominantly negative associations generally found between provisioning services and other ES 17 . On one hand, continued coupling presents more predictable responses to modification, which could aid management. On the other, without decoupling, 'food' production can continue increase to the extent where 'intactness' services collapse. Agricultural production has already been linked to flooding, subsidence, pollution 42 and salinisation 43 and in turn a reduction in human wellbeing in several deltas 44 . In contrast to the population, 'food' became decoupled from the human footprint in 9% of deltas. This decoupling matches the inverted U-relationship proposed between land-use intensity and agricultural services 4 . As land-use intensity increases agriculture is supplanted by urban land uses, as is currently occurring in Southeast Asian deltas 45 . This change is permitted by the replaceable nature of 'food' services, which can be easily imported, in comparison to 'intactness' services. This developmental pathway from natural ecosystems to agriculture to urban land uses may then circumvent further local negative impacts from intensive agriculture. Conversely, it weakens the warning feedback signal local ecosystems can provide 9 , externalises the impacts of these populous, growing urban centres, and may lead to wider food security issues given the importance of deltas to global food production 46 .
'Intactness', in contrast to the other bundles, showed consistent negative responses to modification-two-thirds of this bundle are lost at the highest levels of modification. While some 'intactness' ES, such as flood control and water purification, can be substituted for more expensive engineered solutions 31 , others, such as biodiversity and soil quality services regenerate very slowly 47,48 and are more difficult to replace. The consistent negative associations of 'intactness' with population density and HANPP show their continued impact even in highly modified deltas. However, decoupling already occurs at low levels of the human footprint, suggesting that more fragile biodiversity ES may be lost with relatively little built infrastructure. These disproportionate impacts of even limited land-use change on ES concur with other studies 49,50 . In contrast, contrary to typical findings 17 , 'intactness' initially shows strong trade-offs with the other bundles, even 'productivity', which includes many indicators that measure or are strongly related to biodiversity 51,52 . This trade-off could be explained by the majority of high 'intactness' deltas being less accessible and climatically suitable for human development given their very low productivity. Minimum amounts of ES are required for a delta to develop; even if this very development eventually degrades these services. The decoupling of 'intactness' could be interpreted as a positive sign for deltas-suggesting that a minimum set of 'intactness' related ES can be maintained as modification increases. Yet individual deltas beyond the decoupling threshold display widely different outcomes; the more modified Scandinavian deltas retained near-maximum 'intactness' ES while Southeast Asian deltas showed amongst the lowest. The sustainability and replaceability of this decoupled level of 'intactness' is a critical question in delta futures.
The 'productivity' bundle decoupled with relatively little modification. This bundle included primary productivity and species richness indicators-the nature and habitats that support ES supply. The initially positive responses we observed between 'productivity' and modification, again reflect that some supply of these ES is necessary for deltas to develop. However, beyond the decoupling threshold, this positive relationship breaks down. For HANPP, this suggests that artificial means, e.g. fertilisers or biotechnology, are required to maintain agricultural output for the most modified deltas 53 , which may explain why 'food' services decoupled later than others. However, the many species richness indicators included in 'productivity' will, as with 'intactness', be more difficult to replace and therefore potentially lead to decline in other ES as modification continues. The loss of natural productivity has major implications beyond individual delta systems, distorting biogeochemical carbon and water cycling and diminishing associated climate and water regulation.
The remaining ES indicators, clustered in the 'other' bundle, showed individually weak response curves from which decoupling could not be identified, and clustered inconsistently due to their weak correlations with one another. Surprisingly, many of these indicators were water-related ES, such as sediment flux, wetlands, and discharge, which play a vital role in deltas. We may have expected decoupling of such 'water' ES with modification through their replacement through virtual water trade 54 or engineering, or the different demands of changing land uses 30 . Similarly, wetland area showed little to no response to modification despite their widely reported importance in providing ES and vulnerability to human influence 55 . This perhaps reflects the absence of wetland quality from the indicator-many wetlands are highly modified, and artificial wetlands are increasing in area 56 . Equally surprisingly, upstream flow disruption had weak relationships with water indicators and indeed all bundles, despite the mechanistic links with water, sediment flows and services such as agriculture and fisheries 35,57 . These weak relationships may reflect issues with the indicators, or the complex interrelations of upstream influences and human modification on delta systems. This complexity is illustrated by recent findings that despite severe sediment flow reduction in many deltas due to upstream modification 58 , the global net delta area is increasing 59 . Therefore, while the modification of delta hydrology is extremely important to delta sustainability 29 , it appears to have no simple relationship with ES.
We used a breadth of indicators to as fully as possible represent the modification types and ES and related properties and processes important to deltas. The nature of the available indicators means their interpretation requires care. Some indicators, such as forest and wetland area represent multiple ES, perhaps accounting for their absence of clear synergies and trade-offs; repeating the analyses with these in different categories did not substantially affect our results (Supplementary Note 11). Other indicators, such as flood regulation and the prevention of soil loss, used where direct ecosystem service metrics were unavailable, are impossible to separate from technological solutions or external influences. Indicators for cultural or relational ES were also unavailable, meaning these understudied ES, which may have shown different responses, could not be included explicitly. The selection of indicators will also affect the bundling process. When combining ecosystem service indicators with varying responses, e.g., within the ecosystem property and supply categories, internal trade-offs can hide individual responses. We used unsupervised bundling to overcome this problem by grouping correlated ES. While clustering will always produce a least robust cluster, here creating a bundle ('other') with some less logical components, our ensemble approach has the benefit of measuring this robustness. Using an unsupervised approach with pre-existing ecosystem service frameworks can combine the advantages of unsupervised bundling with greater logical consistency within bundles. We found identifiable bundles and consistent response curves across each unsupervised method.
Further, we should be cautious against drawing conclusions about individual deltas and their decoupling thresholds. Global datasets, often based on modelling, may be less reliable in remote areas with fewer data points, sometimes coarse resolutions may over or underestimate values in small deltas, and temporal mismatches could affect the responses seen. The distributary networks we based our analysis on will have different local contexts and influences, and their delineation will also affect results. Decoupling thresholds will evidently rely on the modification indicators. Here, human footprint had several higher-thanexpected scores, potentially weakening the response curves-for example, overglow from night lights and the equal weighting of different road sizes can overestimate rural infrastructure 40 in Arctic deltas such as the Kuparak. Finally, responses may be affected by overlaps between ecosystem service and modification indicators, given that ES depend on human interactions with ecosystems. These overlaps are clearest between HANPP and agricultural ES. However, such overlaps would be expected to lead to consistent responses; the presence of response curves indicating decoupling, despite this, if anything reinforces our findings. There are also overlaps within the modification indicators: population represents a basis for other modification 32 , while HANPP and human footprint both include land-use change, and dams are linked to infrastructure and agriculture 57 . Despite these factors, the distinct responses and different synergies and tradeoffs that develop show that each type of modification has a distinct effect on ES.
We present a global analysis of ecosystem service modification in deltas and demonstrate multiple, mutually supporting lines of evidence for the decoupling of ecosystem service bundles with built infrastructure therein. We observed two main delta development pathways-towards agriculture or built infrastructure. Given intensive agriculture in deltas is typically for global trade more than local consumption 45 , decoupled urban deltas are perhaps benefitting from a 'coerced' coupling of agricultural deltas only possible with ever-increasing artificial inputs 9 . In both cases, local needs are decoupled from ecosystem service supply, complicating resource management and threatening local and wider sustainability. Our analyses present several exciting directions for future research. We can use longitudinal data to examine and identify mechanisms responsible for the decoupling of ES from modification in individual areas. These mechanisms, and ecosystem service management could be further informed by adopting a 'bright spots' approach 60 and assessing the similarities in those deltas which maintain ES at higher modification levels. Finally, we can match these mostly supply-based services to their users, identifying the winners and losers from modification and rehabilitation. Observing ecosystem service response to modification may help us predict the impacts of societal development in deltas and elsewhere. Such response curves have the potential to inform us where and to what extent development can be combined with sustainable ecosystem service supply, and where ecosystem service degradation can occur 61 .

Method
Study system. We mapped all large deltas globally, i.e. those over 10 km 2 , as these were more likely to feature a variety of ES relevant to delta systems. We visually identified delta distributary networks, using high spatial resolution satellite imagery, and the then most comprehensive global delta dataset 62 , additionally including Arctic, island and inland deltas. We then manually delineated the area of the 235 large deltas found using the watershed either side of the furthest distributary 63,64 ; see Supplementary Methods 1 for protocol. Although the distributary network may not correspond with the sedimentary processes that define a delta, it presents an internally consistent dataset by providing a definable, observable, and mappable area containing the vast majority of the largest global deltas.
Delta modification and ecosystem service data. We measured delta modification using four pre-existing datasets: population density, flow disruption, HANPP, and human footprint (see Supplementary Note 3). Population density (people per km 2 ) was based on 2010 census-data, at a 30 arc-sec (~1 km) resolution 65 . HANPP, the percentage of net primary productivity appropriated by humans was calculated for 2000 using vegetation modelling and land use statistics, at a 5 arc-minute (~10 km) resolution 38 . Flow disruption, the increased residence time of water within a watershed due to reservoirs, was modelled on global discharge, watershed and large dam data, at a resolution of 0.5 degrees (~55 km) 37 . Finally, from human footprint, an index of cumulative human pressure, we selected those components not directly overlapping our other indicators: built-up environments, electric power infrastructure, roads and railways, measured for 2009 at a 30 arc-sec (~1 km) resolution 66,67 .
To provide a broad picture of ecosystem service relationships in deltas, and given the gaps in global data availability, we use a mixture of ecosystem property, supply and service indicators. We selected 49 ecosystem service indicators (Table 3,  Supplementary Table 1) based on their global coverage and reliability (i.e. inclusion in peer-reviewed papers or creation by recognised agency). Where multiple indicators were available for the same ecosystem service, we selected among them based on recency and resolution, or included those showing appreciably different aspects of an ecosystem service, e.g. crop area, output, and economic value.
All indicators were based on publicly available datasets, using the most recent data available. We calculated the average value for each indicator for each delta. For raster data we took the average value of the centroids of every cell within the area of each delta. For vector data we intersected the data with each delta area and calculated the proportion covered. Where necessary (e.g. indicators involving intact species richness) we calculated proportions. We rescaled indicators where Fig. 4 Dichotomous key to classify the relationship between each modification indicator and ecosystem service bundle. The relationship between modification indicator and ecosystem service bundle was visualised using a LOESS regression (Fig. 2), then the regression curve described was classified with this key. Decoupled relationships are indicated by the red boxes with asterisks. Fig. 5 Workflow to identify ecosystem service bundles. The numbered boxes distinguish the individual ecosystem service indicators. We first established the optimal number of clusters for our methods to resolve to using consensus clustering, in this case four. Next, we used six clustering methods to create bundles of ecosystem services (ES). Co-occurring ES have the same colour. We then created an algorithm to relabel co-occurring ES, so bundles containing similar ES were given the same label. Lastly, we built the final bundles from the ES which clustered within a particular bundle most frequently. We calculated the robustness of these bundles by taking the average number of times each constituent ecosystem service was found in its bundle-if every ecosystem service within a bundle was assigned to it by six of the six methods, its robustness would be 100%. necessary so that larger values indicated a greater 'service', and then standardised them for comparison. This centres the indicators around zero, so we subtracted the minimum value of the indicator from each value to make zero the lowest value for visualisations. We describe the indicators, their sources, and additional processing information in Supplementary Table 1.
Identifying ecosystem service and modification relationships and decoupling. We assumed decoupling when there was a shift in the relationship between ES and modification-either a reversal, or to no relationship; which can be described as parabolic or saturating response curves. To identify the shape of the response curve we used a locally estimated scatterplot smoothing (LOESS) regression, using a span of 1-meaning the whole dataset affects each predicted point, resulting in the smoothest curve. We then used a dichotomous key to classify the regression line for each relationship into one of seven types (Fig. 4). For the population density modification type, we applied a log transformation given the right skew. We created a pseudo r-value for each regression (based on the sum of squares of the data divided by the sum of squares of the residuals) to show how well it fit the data. Decoupled deltas were identified as those beyond the point where the relationships shifted from positive to negative or vice versa. We also used breakpoint analysis to identify shifts between one stable regression relationship to another; using the 'strucchange' package 68 . Breakpoint analysis partitions each modification gradient into several regression models by identifying at which points the gradient should be divided to minimise the residual sum of squares and Bayesian information criterion 68 .
Constructing ecosystem service bundles. We clustered the large number of ecosystem service indicators into bundles (see Fig. 5 for workflow) to reduce the overlap between many of the indicators and the unevenness in data availability between different ecosystem service types. We used supervised and unsupervised approaches to group the ecosystem service indicators. In the supervised approach we divided the indicators into the MEA provisioning, regulating, and supporting groups; and the ecosystem service cascade property, supply and service categories. The unsupervised approach used an ensemble of six clustering methods: k-means, k-medoids (partitioning around medoids), and two variations each of random forest and hierarchical clustering (see Supplementary Methods 2 for explanation). Ensemble methods increase classification accuracy and give a metric of clustering robustness, while reducing the subjectivity of algorithm choice. These methods require a defined number of clusters to resolve to. We used consensus clustering 69 which performs multiple runs (1000) over subsets of the data (0.8), then used the 'elbow' method to indicate where diminishing returns in cluster consensus began (Supplementary Methods 2.1).
We applied the six clustering methods on the 49 ecosystem service indicators across all the deltas. This produced six sets of four ecosystem service groups made up of correlated ES. The labels each method assigned to clusters were often arbitrary (i.e. cluster 1 from one method may share the ES of cluster 4 from another). Therefore, we developed an algorithm to relabel the clusters, which iterated through combinations of cluster labels to maximise the number of matching clusters for each ecosystem service. We then assigned each ecosystem service to the bundle it appeared in the most times, assigning them to the most logical bundle when there was an even split between bundles. This became the final bundle set. The number of times an ecosystem service was assigned to a bundle (with a maximum of six out of six methods) indicates how strongly it associates with that bundle, while the average of these values for each bundle indicates its overall robustness. We finally assembled bundles within the MEA and ecosystem service cascade categories using k-medoids to cluster the ES.
Assessing ecosystem service synergies and trade-offs. We used a pairwise Spearman's rank correlation coefficient (ρ) to show positive and negative correlations indicating spatial synergies and trade-offs between ES. We stratified the deltas using the breakpoint analysis method to order each modification indicator against an evenly spaced series (1:235), setting an arbitrary number of four breaks to subset deltas into four stages of modification. We then repeated the tests for Spearman's ρ to assess synergies and trade-offs for deltas of each stage of modification. We finally created ecosystem service bundles using k-medoids clustering to assess how bundle composition shifted at each stage of modification. We used k-medoids clustering to provide consistent results across the stages and indicate the relative strength of clusters.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The delta shape files created for this project are available from https://doi.org/10.5281/ zenodo.6346472 70

Code availability
Spatial analyses were conducted in QGIS 3.10 71 . Clustering and statistical analyses were performed in R 3.6.3 72 . The clustering algorithm was implemented in MATLAB R2016a.