High potential for loss of permafrost landforms in a changing climate

The presence of ground ice in the surface layers of permafrost is central to the geomorphic and hydro-ecological functioning of high-latitude and altitude permafrost landscapes (Rowland et al 2010, Kokelj et al 2014, Liljedahl et al 2016). Involved geomorphic processes give rise to permafrost landforms, such as ice-cored pingo mounds (Mackay 1988), ice wedges that form distinctive polygonal patterned ground (Washburn 1980), and slowly creeping rock glaciers in mountains (Barsch 1996). Ice-rich permafrost landforms are responsive to climatic perturbations as indicated by the intensified permafrost degradation during the past decades (Liljedahl et al 2016, Farquharson et al 2019). Relict permafrost landforms found in temperate climates (e.g. pingo remnants and ice-wedge casts in western and central Europe (Vandenberghe and Pissart 1993) imply that distributions of permafrost landforms have changed before and are likely to change in future climates (Washburn 1980, Grosse and Jones 2011, Liljedahl et al 2016).

Projected increases in air temperatures and precipitation with Arctic amplification (Hoegh-Guldberg et al 2018) are anticipated to intensify permafrost warming (Biskaborn et al 2019) and associated thermokarst processes, i.e. thawing and settling of ice-rich permafrost (Jorgenson et al 2006, Grosse et al 2016, Douglas et al 2020). Landscapes with current thermokarst features or susceptibility to future thermokarst are estimated to cover up to 20% of circumpolar permafrost regions (Olefeldt et al 2016). Rising air and ground temperatures pose a potential constraint for pingo persistence across northern Siberia (Grosse and Jones 2011), Alaska (Jones et al 2012), and Canada (Mackay 1988), for example. Degradation of ice-wedge polygons has been observed across the Arctic, including continuous permafrost regions (Jorgenson et al 2006, 2015, Grosse et al 2016, Liljedahl et al 2016, Steedman et al 2017, Fraser et al 2018, Farquharson et al 2019). Moreover, downslope creeping of rock glaciers due to ground ice deformation has shown increased surface velocities related to increased air temperatures in some regions (Kääb et al 2007, Jones et al 2019). The degradation of ground ice in permafrost landforms can induce thermokarst or mass-movement hazards (Kääb et al 2007) and thereby threaten infrastructure development (Rowland et al 2010, Kanevskiy et al 2017, Hjort et al 2018, Streletskiy et al 2019).

Permafrost landforms bear local-scale geomorphic (Humlum 2000, Knight et al 2019, Ward Jones et al 2019), hydrological (Liljedahl et al 2016, Brighenti et al 2019, Nitzbon et al 2019) and biological (Lara et al 2018) significance, but also affect global systems through alterations in atmospheric (Schuur et al 2015) and aquatic (Serikova et al 2018) greenhouse gas fluxes from permafrost regions. These landforms also contribute to local biodiversity (Jorgenson et al 2006) by serving as highly specialized ecosystems for flora (Jorgenson et al 2015) and fauna (Fountain et al 2012). In some regions they also provide a freshwater supply for human use (Jones et al 2018).

At present, the information on current and potential future distributions of permafrost landforms at the circumpolar scale is sparse. This is despite an increasing amount of circumpolar to global permafrost extent simulations conducted both with process-based (e.g. Guo and Wang 2016) and empirical/statistical approaches (Chadburn et al 2017, Aalto et al 2018, Obu et al 2019). The lack of spatially detailed knowledge on geomorphic features poses a limitation to the assessments of climate change impacts on permafrost environments. We aim to reduce these knowledge gaps by applying statistical distribution modelling. Statistical modelling methods have been used in geomorphological contexts (Luoto et al 2010, Aalto et al 2014, Rudy et al 2017, Fewster et al 2020). In our previous work (Aalto et al 2018), we have demonstrated that a statistical approach produces similar permafrost extents with process-based models (e.g. Guo and Wang 2016) and other empirical models (Chadburn et al 2017, Obu et al 2019), and have relatively small prediction errors. The advantage of the approach over physically based models which require extensive parameterization of physical soil properties, for example, is that statistical modelling allows for predicting landform occurrence at high spatial resolution across the circumpolar region using readily available geospatial data on climate, soil and topography. Thus, our approach does not consider various processes affecting heat fluxes or ground ice accumulation in permafrost but is suitable for examining how different environmental factors affect the distribution of permafrost landforms across the circumpolar region.

Our objective is to model, for the first time, the potential environmental spaces for the occurrence of pingos, ice-wedge polygons and rock glaciers across the Northern Hemisphere. These spaces encompass conditions where climate, topography and soil properties are suitable for landform presence. In the modelling we assumed the landforms to be in equilibrium with climate conditions for the considered period, that is, we did not make assumptions concerning whether the landforms are degrading or still growing. Based on the statistical modelling results, we present the circumpolar coverages of potential environmental spaces for the studied landforms during a recent period (1950–2000) and predict their regional changes for two future periods (2041–2060 and 2061–2080) using representative concentration pathway RCP4.5 (stabilization scenario) and RCP8.5 (rising radiative forcing pathway) climate-change scenarios. Moreover, the results offer insights into the fine-scale distribution of various ground ice types and help inventory permafrost landforms in unmapped regions.

2.1. Permafrost landform observations

Observations of pingos, ice-wedge polygons and rock glaciers across the Northern Hemisphere (i.e. presence data, figure 1(a)) made since late 20th century were compiled from available inventories, published studies and geomorphological maps (see section S1 (available online at stacks.iop.org/ERL/15/104065/mmedia) in supplementary material). Ice-cored pingo mounds (figure 1(b)) form in permafrost environments by injection and freezing of pressurized water in the soil, and results in massive-ice formation and upheaving of overlying surface material (Mackay 1988). Two primary pingo types are often distinguished: hydrostatic, predominantly occurring in low-lying continuous permafrost regions, and hydraulic, which are more typical to discontinuous permafrost and depend on water moving under a hydraulic gradient (Mackay 1988, Grosse and Jones 2011, Jones et al 2012). Distinct patterned ground features, ice-wedge polygons (figure 1(c)), occur nearly ubiquitously in unconsolidated deposits across the Northern Hemisphere permafrost domain and wedge ice is the most common form of ground ice in continuous permafrost (Bernard-Grand'Maison and Pollard 2018). The polygon shape results from the cyclical growth of ice wedges in thermal contraction cracks formed in susceptible soils during cold spells in wintertime (Washburn 1980, Kokelj et al 2014).

Zoom In Zoom Out Reset image size

Rock glaciers are bodies of poorly sorted debris and ice that move due to deformation of internal ice (Barsch 1996, Berthling 2011, figure 1(d)) in permafrost environments. They occupy most mountain regions in both continental and maritime climates including high-Arctic conditions (Evans 1993) and high elevations in warmer climates in, e.g. the Middle East (Gorbunov 2013). We included both active (flowing or creeping downward) and inactive (stagnant) rock glaciers that contain ice (Barsch 1996) regardless of whether they were primarily influenced by talus slope dynamics or by glacier dynamics (moraine-derived rock glaciers) (see section S1, supplementary material) but excluded relict landforms. In addition to presence data, the applied modelling required information on environments where the landforms are not present. Thus, we compiled datasets with 1349 absence grid cells for pingos, 494 for ice-wedge polygons and 1549 for rock glaciers (section S1, supplementary material).

2.2. Environmental data

2.3. Statistical modelling

Based on the predicted environmental spaces for current and future time periods, we classify three pathways of suitability for landform occurrence; 1) persisting conditions with suitable topographic and soil conditions for landform occurrence and climate remaining favourable, 2) environmental space loss, in which climatic conditions become unsuitable, and 3) new environmental spaces with potential for occurrences outside current distribution (i.e. suitable topography, soil and newly formed climatic conditions). Finally, we examine the factors affecting landform occurrence and evaluate the reliability of the predictions.

3.1. Regional changes in environmental spaces

Based on ensemble modelling, the predicted current coverage of environmentally suitable spaces is 2.04 × 10⁶ km² for pingos, 3.12 × 10⁶ km² for ice-wedge polygons and 1.40 × 10⁶ km² for rock glaciers (figure 2, union of orange and red zones). Under the moderate 2041–2060 climate-forcing scenario (RCP4.5), the total areas of potential environmental spaces for the studied permafrost landforms are predicted to reduce to 1.44 × 10⁶ km² (by 29%), 2.47 × 10⁶ km² (21%) and 1.08 × 10⁶ km² (23%), respectively (figure 2, supplementary tables 2 and 3). By 2061–2080, assuming the RCP8.5 scenario, on average over 50% of recent environmental spaces may be lost.

Zoom In Zoom Out Reset image size

The largest losses of environmental spaces for pingos and ice-wedge polygons are concentrated on northwestern Russia and the middle reaches of major Arctic rivers (figures 2(a) and (b)). In the zone of discontinuous permafrost, suitable conditions for rock glacier occurrence show a distinctive shift towards higher elevations. As a result, areal losses are projected in most areas of current occurrences (figure 2(c)). In continuous permafrost regions, the projected losses are generally small.

New environmental spaces are projected to emerge where air temperatures and precipitation attain suitable levels for occurrence and other factors are favourable. For pingos and ice-wedge polygons, new suitable areas (circumpolar totals of 0.67 and 0.33 × 10⁶ km²) predominantly situate across Canada's Arctic Archipelago and in the Taymyr Peninsula (figures 2(a) and (b)). Pingo environments show a poleward shift, while suitable conditions for ice-wedge polygons are predicted in diverse landscapes amidst the recent occurrences. However, suitability for recent pingos, also, is predicted in many locations on Canadas's Arctic islands even though not visible in the small-scale map (figure 2(a)). The largest areal increases in potentially suitable conditions for rock glacier occurrence (0.15 × 10⁶ km²) appear in northern Alaska, high-Arctic Canada, eastern Greenland and Far-Eastern Russia (figure 2(c)). The general trajectory of the changes is towards colder regions of continuous permafrost. Based on the climate-change scenarios, both air temperature and rainfall show distinct increases at grid cells with observed landform occurrences (figure 3). Snowfall is projected to decrease at rock glacier sites, in contrast to the increase in rainfall (figures 3(c) and (d)).

Zoom In Zoom Out Reset image size

3.2. Climatic controls and sensitivity of permafrost landforms

According to the analyses of variable importance and response shapes, the optimal conditions for pingos and ice-wedge polygons at a circumpolar scale are similar (figures 4(a)–(c), supplementary figure 3) and characterized by flat topography with moderately high soil moisture, low rain- and snowfall (<∼300 mm), and TDD (°C-days) of less than ∼2000. For ice-wedge polygons, FDD also have a notable contribution indicating that their occurrence requires at least ∼2500 °C-days. Notably, figure 3 suggests that FDD will drop below 2000 at many pingo and ice-wedge polygon sites regardless of the climate-change scenario.

Zoom In Zoom Out Reset image size

Rock glacier occurrence increased on coarse-sediment slopes with low wetness index, TDD less than ∼1200 °C-days, and FDD between ∼1000 and 4500 °C-days (figure 4(d)). Increased probability of occurrence towards minimum FDD values in the dataset suggests that rock glaciers exist across a wide range of cold climates although most landforms in the dataset have relatively high FDD values (figure 3(b)). Each landform shows a clear decline in probability of occurrence when TDD exceeds ∼1000 °C-days (figures 4(b)–(d)). The response of rock glaciers to increasing TDD is especially sharp. An even more abrupt threshold is found between rainfall and both pingos and ice-wedge polygons.

3.3. Statistical and spatial evaluations

To assess model fit and predictive performance we used a split-sample approach; the models were calibrated using 70% of the data and evaluated against the remaining 30% of observations. Based on 100-fold cross-validation, our approach predicted well the permafrost landform occurrences (figure 5). Ensemble method yielded the highest classification accuracy (AUC and TSS) in all cases and thus superior performance. Rock glacier models had the highest R², AUC and TSS values, yet respective ensemble averages in pingo and ice-wedge polygon models also indicated excellent predictive performance. Moderate between-model variability was observed. GBM and RF best explained the variation in each landform occurrence having R² values at or above 0.90 (figure 5).

Zoom In Zoom Out Reset image size

In key regions of landform distribution, the models largely agreed (figure 6). High inter-model variability was encountered in poorly sampled regions, such as the Tibetan Plateau where GLM and GAM predicted notably larger ice-wedge polygon distributions than GBM and RF, which both yielded more constrained environmental spaces (figure 6(d)). The lower evaluation statistics of GLM and GAM in pingo modelling (figure 5(a)) may explain the model disagreement in many regions (figure 6(a)). In cases like these, the ensemble method facilitates controlling for prediction uncertainty by allowing for a conservative model agreement-based assessment. Despite superior model performance, ensemble models failed to identify a few known regions for observed occurrence. For example, only GBM and RF predicted pingo occurrence potential in Svalbard and eastern Greenland (figure 6(a)). Evaluations of spatial uncertainty such as presented here help identifying regions and environmental conditions where additional investigations should be made and thereby improve future model performance.

Zoom In Zoom Out Reset image size

The predicted losses of environmental spaces are largely associated with regions of projected near-surface permafrost thaw in the near future due to warming climate (e.g. Slater and Lawrence 2013, Chadburn et al 2017), but not exclusively, as areas of continuous permafrost in, e.g. northern Siberia may also become unsuitable for pingos and ice-wedge polygons. Based on our analyses, the environmental suitability for these landforms will be foremostly constrained by increasing TDD and rainfall. These findings suggest that warmer and wetter summers, and the projected thickening of the active layer (Aalto et al 2018, Peng et al 2018) can melt pingo ice (Grosse and Jones 2011) and ice wedges close to the ground surface in these regions, as observed elsewhere in the Arctic (e.g. Fortier et al 2007, Liljedahl et al 2016, Bernard-Grand'Maison and Pollard 2018, Farquharson et al 2019). First studies considering thermokarst-inducing processes in a numerical model indeed project substantial surface permafrost degradation and widespread landscape collapse even for the currently still cold North-East Siberian permafrost region under a strongly warming climate scenario (RCP8.5) (Nitzbon et al 2020), which is largely in line with our findings. Alongside increasing air temperatures, increased rainfall has been shown to stimulate permafrost thaw (Douglas et al 2020).

Importantly, predicted losses do not imply complete thaw of permafrost at depth but rather that the permafrost landforms can be lost due to conditions becoming unfavourable for their occurrence. However, ground ice in the landforms can be still protected from thaw by ecosystem properties, such as insulation by thick layers of peat (Shur and Jorgenson 2007) or cooling effect of cold air circulation in blocky rock glacier surfaces (Jones et al 2019). Local hydro-topographical and ecological properties can also determine whether prolonged degradation leads to drainage or wetting of landscape, and thus the potential release of greenhouse gases (Jorgenson et al 2015, Liljedahl et al 2016, Nitzbon et al 2020).

Rainfall was the most important single variable delineating suitable conditions for pingos and ice-wedge polygons at the circumpolar scale, while TWI and FDD foremostly affected rock glacier occurrence potential. Even though these were central factors, the results implied that local topography and soil properties hold notable roles in constraining finer scale environmental suitability for the landforms. Based on the findings, it is suggested that increases in the amount of rainfall to above ∼200 mm (figures 4(b) and (c)) can severely limit pingo and ice-wedge polygon occurrence potential. In the case of ice-wedge polygons, TDD had nearly as an important of a contribution suggesting that air temperature is the central factor for their occurrence. Pingos also showed a secondary relation to TDD yet rainfall was more dominant in their case.

The effect of rainfall is due arguably to the increased conduction of heat into the ground (Jorgenson et al 2006, Iijima et al 2010) but also the introduction of advective heat (Kane et al 2001). Rainfall also has been shown to induce abrupt thaw processes in ice-rich permafrost terrain (Kokelj et al 2015). Owing to the nonlinear response shape, the found effect also applies the other way; suitable conditions can emerge if currently very dry areas begin to receive higher rainfall (>∼100 mm, figures 4(b) and (c)) in the future, provided that air temperature and other conditions are suitable. Some of rainfall's contribution could be attributed to its moderately strong correlation with TDD in the modelling data (supplementary figure 4). Notwithstanding, the analysis of variable importance suggested that rainfall had a greater relative effect than TDD, especially for pingos (figure 4(a)).

The responses of pingos and ice-wedge polygons with rainfall, and rock glaciers with air temperature, were strong and supported by each method. The projections of future climate, however, include uncertainty, especially in mountainous areas (Hijmans et al 2005, Hoegh-Guldberg et al 2018). The interpolated climate surfaces cannot fully account for orographic precipitation patterns or wind-induced snow transport, for example. Ice-wedge polygons' response to snowfall indicated a lowered occurrence probability with increasing snowfall (supplementary figure 3) in line with the observations by Kokelj et al (2014). However, climatic predictors averaged over a series of decades cannot represent future dynamics of snow depth-controlled frost cracking under wintertime cold spells or other extreme air temperature, precipitation or disturbance (i.e. fire) events, which play a role in the formation and degradation of especially ice-wedge polygons (Jorgenson et al 2006, Jones et al 2015, Liljedahl et al 2016, Kanevskiy et al 2017). Based on experimental evidence (Allard and Kasper 1998, Christiansen 2005), thermal contraction cracking occurs during prolonged cooling periods with ground temperatures below −20 to −15 °C. In this study, the moderate importance of FDD and snowfall for ice-wedge polygons suggests that higher temporal resolution of climate predictors (monthly to sub-monthly) and a more explicit representation of snow cover variation are needed to estimate thermal cracking potential.

Recent studies show that ice-wedge polygons formed in centennial to millennial time scales can be lost by melting of their uppermost portion in a few decades (Jorgenson et al 2006, Liljedahl et al 2016, Fraser et al 2018, Farquharson et al 2019). Large pingos with thick insulating overburdens are to a degree resilient to short-term changes in climate (Mackay 1998), but evidence of collapsed and paleo-pingos in concurrent non-permafrost environments confirms their climate sensitivity over longer time scales (Vandenberghe and Pissart 1993, Jones et al 2012). Degradation of ice-wedge polygons at local scales demonstrates large spatial variation owing to complex interactions between, e.g. ground ice, wetness conditions and microtopography (Kanevskiy et al 2017, Steedman et al 2017, Fraser et al 2018, Nitzbon et al 2019, 2020). Notwithstanding, some have argued that future thermokarst rates at broader scales mainly depend on the magnitude of regional climate change (Olefeldt et al 2016). Degradation of any of the studied permafrost landforms can also occur due to their natural lifecycles without external forcing (Washburn 1980, Mackay 1988).

Few quantitative or spatial estimates exist to compare with the predicted coverages. Mackay (1972) estimated that ice-wedge polygons covered 2.6 × 10⁶km² of the Northern Hemisphere, which falls relatively close to our ensemble prediction for 1950–2000 (3.1 × 10⁶ km²). On a regional level, our ice-wedge polygon predictions for the main occurrence areas are in line with a recent Canada-wide modelling of wedge-ice content (O'Neill et al 2019), yet regional discrepancies are visible (section S3, supplementary material).

The discrepancies between the models (figures 4 and 5) are argued to have stemmed from different initial model assumptions and parameterizations (Thuiller et al 2009). For example, only GBM and RF captured a decline in rock glacier occurrence probability when FDD approached zero, i.e. milder freezing-season conditions (figure 4(d)). GLM and GAM overlooked this threshold and extended to too warm and flat areas, such as periglacial environments more characteristically affected by extensive solifluction processes, e.g. in northern Scandinavia (Aalto et al 2017). The wide distribution may also have been partly attributed to the inclusion of moraine-derived rock glaciers (section 1, supplementary material) which are easily mixed with debris-covered glaciers, and thus may extend below lower permafrost boundary limits. Moreover, some of the compiled observations might have been relict rock glaciers, and thus possibly lacked permafrost, as suggested by relatively high TDD and FDD for many objects in the dataset (figures 3(a) and (b)). The discrimination of intact rock glaciers from these landforms is often challenging, especially from optical remote sensing data from which a large part of observations was mapped.

Considering the projected new regions with suitable environmental spaces for landforms, we argue that the dependency of landform growth on local geomorphic and hydrologic conditions and glacial history means that they are likely to occur in more constrained environments than estimated in this study. Prediction of pingo formation potential, for example, is difficult owing to its dependency on special groundwater pressure conditions (governed by underlying substrate properties) and dynamic processes such as lake drainage (Mackay 1998). Here, coarse sediments showed moderate importance and realistic response shapes with the landforms. Due to relatively few observations of soil properties within the permafrost domain available to produce the soil data layers (Hengl et al 2017), the data used here likely have a limited ability in depicting landform-scale variation in grain size proportions. Moreover, although we argue that the ∼1 km resolution is sufficient to delineate topo-climatic conditions for most rock glaciers, there remains a need to account for site-specific debris availability, lithology and melt-water dynamics, or groundwater conditions and sedimentary history in the cases of pingos and ice-wedge polygons. Not considering all constraining factors for landform development may cause modelling to overestimate the extent of new suitable conditions, and subsequently the presented proportions of lost environmental spaces are arguably conservative.

Ultimately, we model suitable conditions for landform occurrences in recent and future climates, but cannot presently address the physical processes, such as ground-ice degradation or aggradation and associated heat fluxes, behind landform development using our approach. The vulnerability of permafrost to atmospheric warming is highly complex and depends on multiple dynamic processes acting in the surface water, snow, vegetation and soil layers and related positive and negative feedbacks (Shur and Jorgenson 2007, Nitzbon et al 2019). Moreover, the thermal inertia of ice-rich ground encumbers permafrost warming (Grosse et al 2011). Ground ice and associated phase changes of water (i.e. the latent heat effect), as well as the content of unfrozen water in frozen soils (Hu et al 2020), are key constituents of the near-surface permafrost response to climate forcing that may not be uniform in different initial ground temperature conditions (Riseborough 1990). Not considering such transient effects lowers our ability to assess landform occurrence potential in the future, given that the used approach assumes equilibrium state between landforms and climate. Process-based models could be appropriate for assessing the physical responses and related feedbacks of climate change because of their explicit parameterization of heat fluxes and transient effects in the surface layers (e.g. Harris et al 2009). However, they are often conducted at coarse spatial resolutions, which limits their applicability in assessing permafrost landforms at a relevant scale in a hemispheric extent.

Statistical modelling allows for producing reasonably accurate information on the suitable environments for permafrost landforms at a high resolution. This information can be used to link locally resolved detailed land surface processes to globally relevant developments due to climate change. Predictions can also facilitate initial landform mapping in understudied regions although their accuracy at pixel resolution (∼1 km²) is to be taken cautiously due to the documented prediction errors (figure 5). In addition to likely profound geomorphic changes, ground ice degradation in the areas of projected lost environmental conditions may greatly alter hydrological and biogeochemical processes. In addition, the eventual disappearance of permafrost landforms represents a loss in local geodiversity, and thereby affects biodiversity (Antonelli et al 2018). From a human activity point of view, ground ice degradation and loss of permafrost landforms can affect freshwater availability (Jones et al 2018, 2019) and infrastructure integrity (Raynolds et al 2014, Kanevskiy et al 2017). Using permafrost landforms as proxies for massive ice, our modelling yielded novel information on fine-scale spatial variation of high ground ice content across permafrost landscapes. It should be noted that pingos occupy a very small proportion of landscapes relative to ice-wedge polygonal terrain or even large rock glaciers, and thus the potential amount of ground ice expected in a grid cell is different for each landform.

From this study the following conclusions could be drawn:

• The suitable environmental conditions for permafrost landforms at the circumpolar scale are closely related to climate and prone to shift along projected regional changes in air temperature and precipitation.
• On average, the landforms are projected to lose more than one-fifth of their suitable environments under a moderate climate scenario (RCP4.5) and around one-third under a very high baseline emission scenario (RCP8.5) by midcentury.
• The sensitivity of pingos and ice-wedge polygons to the amount of rainfall was unambiguous, albeit TDD was almost as important for ice-wedge polygons.
• Alongside air temperatures, precipitation regime is suggested to be a key constraining factor of geomorphic development outside mountainous environments.
• For rock glaciers, air temperature and topography primarily controlled the occurrence potential. However, the approach had a limited ability to separate suitable environments for rock glaciers from those for other glacial or periglacial slope processes.

The predicted changes in environmental spaces suitable for the studied permafrost landforms will affect Arctic landscapes and their geomorphology across large areas. The potential implications should be considered when assessing impacts of climate change on natural and human systems.

This study was funded by the Academy of Finland (grants 315519, 307761). BMJ was supported by awards from the US National Science Foundation (OPP-1806213 and ICER-1927872) and GG was supported by ERC PETA-CARB (#338335), the HGF Impulse and Networking Fund (ERC-0013), and NASA NNX08AJ37G. The contribution by BE and KSL was funded by the University of Oslo, Norway. We wish to thank two anonymous reviewers for their highly valuable comments.

The data that support the findings of this study are available upon reasonable request from the authors.

OK, ML and JH developed the original idea. OK led the compilation of observational data with contributions from BE, GG, BMJ and KSL, and geospatial data processing with JH and JA. OK performed the statistical analyses with contributions from JA, JH and ML. OK wrote the manuscript with contributions from all the authors.

The authors declare that they have no competing interests.