Permafrost extent sets drainage density in the Arctic

We used the WWF HydroSHEDS Level 10 Basins dataset (59), a collection of watershed geometries hosted on Google Earth Engine to locate all headwater basins (defined as having no upstream contributing basin) between 23.5° and 90 °N latitude. We further filtered this dataset to find soil-mantled watersheds by only selecting watersheds with average annual normalized difference vegetation index (NDVI) $>$ 0.60, using each pixel’s annual maximum NDVI for the year 2021 as measured by MODIS’s 250 m Terra NDVI product (MOD13A1 v061). We masked the NDVI product by the MODIS water mask product (MOD44W v006). For remaining basins, we extracted mean annual temperature (MAT) and mean annual precipitation (MAP) as reported by WorldClim BIO variables v1 (60), covering data between 1960 and 1991. Because most Arctic watersheds had MAP $<$ 1,000 mm/y, and MAP $<$1,050 mm/y was found to be the boundary of the inflection point in Sangireddy et al. (36), we eliminated watersheds with $>$1,000 mm/y precipitation. We calculated basin relief by subtracting the highest and lowest elevations within the watershed of digital elevation models created by NASA’s 90-m resolution Shuttle Ray Topography Mission (SRTM) for 25${}^{°}$ to 60${}^{°}$ (its upper latitude limit) and 2-m resolution ArcticDEM for 60${}^{°}$ to 90${}^{°}$ (61). We calculated drainage density using the Hydrography90m dataset (62), which was created using the MERIT Hydro hydrologically conditioned digital elevation model which uses the SRTM dataset with a resolution of $\sim$90 m at the equator below 60${}^{°}$ latitude and the JAXA ALOS 30 m DEM from 60 to 90${}^{°}$ latitude (63) to derive a global network of 726 million stream segments with a minimum upstream contributing area of 0.05 km${}^2$. We masked this dataset with the MODIS water mask and calculated drainage density as the proportion of pixels associated with stream segments to the total number of flow accumulation pixels in the watershed boundary. This method of drainage density calculation is used to report data from Figs. 2 and 3. Once watershed attributes were determined from climate and satellite data, we intersected the centroids of watersheds with the Circum-Arctic Map of Permafrost and Ground Ice Conditions (64) to classify watersheds into categories of continuous, discontinuous, sporadic, isolated, and no permafrost.

We divided our data into “permafrost” (n = 31,974) and “no permafrost” (n = 37,618) based on the ground ice extent map. We binned each watershed of each group into one of 20 bins for MAP (range 0 to 1,000 mm/y) and relief (range 0 to 1,200 m) for a total of 400 bins. For bins that contained at least 10 watersheds in each of the permafrost and no permafrost groups, we calculated the ratio of the mean drainage density of permafrost versus non-permafrost datasets for each data bin (lower drainage density in permafrost watersheds results in a ratio $<$1.0) (Fig. 2A). To test the significance of the difference in ratios for each bin, we performed a Mann-Whitney U test using the Python package scipy (65) on the distributions of drainage density for permafrost and non-permafrost watersheds. The Mann-Whitney U test is a non-parametric test that does not presume a normal distribution between two populations and tests whether any sample of one population will be larger than any sample from the other population. We chose a threshold of $P<{10}^{-4}$ and differentiate bins that did not pass this significance test with hatches.

To confirm that glacial history is not a confounding variable in our study, we intersected the centroids with a map of the extent of the Last Glacial Maximum (LGM) at $\sim$24 ka at MIS 6 at $\sim$190 ka (66) to determine whether watersheds of different permafrost extents fell within the LGM boundary (which we assumed also included the MIS 6 extent), the MIS 6 boundary only, or unglaciated. We also performed Mann-Whitney U tests on boxplot pairs between permafrost extent categories for each of the glacial histories; all distribution differences are significantly different below a threshold of $p<{10}^{-4}$ (Fig. 2B).

We selected unglaciated watersheds in our dataset and separated those unglaciated sites into watersheds with (n = 16,002) and without (n = 23,220) permafrost. We eliminated permafrost watersheds with MAT $>$ 2.5 °C and non-permafrost watersheds with MAT $<$ $-$2.5 °C. We performed an ordinary least squares (OLS) regression between MAT and drainage density for each group of watersheds and found that generally, drainage density increases with MAT (Fig. 3A). To determine the role that MAP might play in the MAT-drainage density relationship (as temperature and precipitation covary), we performed an OLS regression between MAT and MAP for both datasets using the Python package statsmodel and calculated the difference between a site’s actual MAP and the regression fit prediction (residual) for each watershed. This procedure highlighted watersheds with particularly high (“wet”) or low (“dry”) MAP given their MAT. We then regressed these precipitation residuals against drainage density (Fig. 3 B and C).

Both the coarse resolution of the underlying MERIT Hydro digital elevation model (90 m) and the simplified approach to delineating channels by area thresholds are drawbacks to the use of Hydrography90m to calculate drainage density (62). To check whether our results are influenced by either of these factors, we randomly sampled all continuous permafrost and non-permafrost watersheds to generate a subset of unique non-permafrost (n = 476) and continuous permafrost (n = 460) watersheds that spanned the range of relief and MAP represented by the larger dataset (SI Appendix, Fig. S8). We limited the non-permafrost watersheds to the conterminous United States in order to use the 10-m USGS 3DEP National Map Seamless (1/3 Arc-Second) dataset. Of the 936 subset watersheds from continuous permafrost and non-permafrost landscapes, 504 had significant overlap in MAP-relief space and were therefore binned by MAP and relief as in Fig. 2A. We downloaded topographic data for each watershed boundary from the USGS dataset (non-permafrost) and used bilinear resampling to downscale the 2-m ArcticDEM dataset to 10 m to match the National Map dataset for continuous permafrost watershed boundaries. We then used the LSDTopoTools software suite (67) to preprocess the DEMs (remove invalid data, fill sinks, and remove dams) and calculate basic topographic metrics (see SI Appendix for algorithm parameters). We created flow accumulation rasters with the d-infinity algorithm (68) and calculated tangential curvature with a window of 100 m. We then used the DrEICH algorithm (69), which uses curvature and channel steepness thresholds to select channel heads, to create channel networks using constants $A_0$ = 1.0 and $m/n$ = 0.5. We calculated drainage density with the same pixel-counting method as the Hydrography90m dataset and found that while watersheds calculated with the Hydrography90m dataset contained 7$\times$ the proportion of channel pixels as the DrEICH algorithm channels (likely due to the difference in pixel resolution of the two datasets), the correspondence between the two drainage density metrics was generally good ($r^2$ from OLS regression=0.80; see SI Appendix, Fig. S7). The mean of the residual values for this regression is $<0$ in permafrost and $>0$ in non-permafrost, meaning that on average, Hydrography90m overestimates permafrost drainage density compared to the DrEICH algorithm and underestimates drainage density in non-permafrost watersheds. We then performed the same comparison between non-permafrost and permafrost drainage densities binned by relief and MAP as described in the text and shown in Fig. 2A and found all bins to exhibit lower drainage densities in continuous permafrost watersheds, with most bins exhibiting P values $<$ 0.005 (SI Appendix, Fig. S8). Despite shortcomings in both channel network delineations given poor remote sensing coverage in polar regions, these independent methods yield the same result.

We used the topographic metrics calculated from the high-resolution subset to asses the dominance of concavity (created by advection) and convexity (created by diffusion) across the landscape independent of channel delineation (Fig. 4 A and B). For each watershed in the subsample, we calculated a cumulative area distribution for flow accumulation (drainage area) binned by drainage areas of ${10}^2$ to ${10}^5$ m${}^2$ to determine the relative proportion of flow accumulation values in a watershed. We also found the median tangential curvature for each of these flow accumulation bins to determine the drainage area associated with a switch from convex hillslopes (negative tangential curvature) to concave river valleys (positive tangential curvature). For each type of permafrost extent, we calculated the median cumulative area distribution and tangential curvature for each drainage area bin.

For permafrost watersheds, we regressed MAT against drainage density (as done in Fig. 3A) to calculate that for every degree of higher MAT, a watershed’s channelized pixel percentage (the unit of drainage density used in this study) grows $\sim$0.15%. We then calculated this area in m${}^2$ for each watershed (dataset mean= 44,394 m${}^2$). We assigned each watershed an average carbon mass in the upper 100 cm by intersecting the watershed centroids with the Northern Circumpolar Soil Carbon Database version 2 (55). We multiplied this carbon mass by the per-degree channel expansion area for each watershed to calculate a summed total of 4.35 $\times$${10}^{10}$ kg C. The equivalent mass of CO${}_2$ equals about 34.6 million passenger cars driving for a year, according to the US EPA (70).

J.D. designed research; J.D. performed research; J.D. analyzed data; and J.D., M.C.P., and C.R.M. wrote the paper.

The authors declare no competing interest.