Permafrost extent sets drainage density in the Arctic

Update 10/09/2023

A new release coincides with submitting revisions to PNAS! The repo now features a "Step 5" that calculates carbon stocks vulnerable to erosion with climate change.

Find the (un-revised) preprint here.

Summary

This coding project is the culmination of many months of learning Python geospatial workflows, Google Earth Engine, and parallel computing as part of my postdoctoral fellowship at the Neukom Institute for Computational Sciences at Dartmouth College. I am pleased to share the workflow with you here! Note that while all the satellite and topographic datasets used here are freely available, they are most easily downloaded/computed with Google Earth Engine. The conversion from Earth Enging server-side function and local client-side Geopandas/Pandas conversion is only possible through the amazing geemap package by Qiusheng Wu. The analysis of 10 m digital elevation models is only possible through the incredible LSDTopoTools software project.

I recommend two different environments to run these notebooks - I have a gee environment that is pretty much the included environment.yml. The second is a lsdtopy environment based on the lsdtopytools install instructions.

tl;dr you can download all 69,592 watershed results (nice) from watershed_export.csv to see what this is all about.

Code outline

1. Step 1 (which requires a Google Earth Engine account) uses the WWF HydroSHEDS Basins Level 10 dataset to find headwater watersheds from 25-90 degrees latitude broken up into 1 longitude slices. These geometries are then used to extract climate data from WorldClim BIO Variables V1, NDVI data from MOD13Q1.061 Terra Vegetation Indices 16-Day Global 250m, and elevation data from either SRTM Digital Elevation Data Version 4 (for watersheds < 60 latitude) or ArcticDEM Mosaic. The code then counts the number of pixels of river segments from the Hydrography 90m dataset. The channel segments and NDIV data are masked by the MOD44W.006 Terra Land Water Mask Derived From MODIS and SRTM Yearly Global 250m. A giant table of data is created based on the HYBAS ID of each watershed and turned into two shapefile files, arctic_h90.shp and midlat_h90.shp.
2. Step 2 calculates drainage density using the pixel counts of channels versus the number of pixels in the watershed from the Hydrography90m dataset. The code then intersects the centroids of the watersheds with shapefiles of both permafrost extent and glacial boundaries at the LGM and MIS 6. We can then create Figures 2 and 3 in the main text based on (1) permafrost extent, (2) glacial history, (3) watershed relief, (4) mean annual precipitation (MAP), and (5) mean annual temperature versus drainage density. The drainage density "heat map" is created by binning watersheds by relief and MAP, calculating the mean drainage density for the bin, and performing a binwise Mann-Whitney U test on the distribution of drainage densities of the two populations. This also performs the OLS regression of MAT versus MAP to calculate residual MAP values in Figure 3B an C. This script also makes many supplementary figures. There is also the random subsetting of non-permafrost and continuous permafrost watersheds for steps 3 and 4.
3. Step 3 queries Google Earth Engine to download either the USGS 3DEP 10m National Map Seamless (1/3 Arc-Second) DEM or a bilinearly resampled 10 m DEM from ArcticDEM for the subset watersheds. Each watershed gets a subdirectory within a directory lsdtt named after its HYBAS ID and a .tif is placed in the folder.
4. Step 4 does a lot - it uses LSDTT's lsdviztools to first run gio.convert4lsdtt on a watershed DEM to preprocess for the LSDTT command-line tools. Each watershed has the same lsdtt_parameters applied for lsdmw.lsdtt_driver with command_line_tool = "lsdtt-channel-extraction" in which I do things like get curvature, d-infinity flow accumulation, and area threshold and DrEICH channels. BUT I wrap this function up in a parsl dataflow to take advantage of the many cores available on Dartmouth's shared memory computers. YMMV. I then iterate through the resulting files to extract, for each watershed, data on drainage area, curvature, and channels. I then compare the two different channel datasets' drainage densities and create the curvature-cumulative-area-drainage-area plot of Figure 4.

" If I have seen further it is by standing on the sholders [sic] of Giants." - Isaac Newton