Skip to content

mankoff/sankey

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 

History

95 Commits
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Repository files navigation

On Ice Sheet Mass Flow

Table of contents

Introduction

Sankey diagrams for mass flow in Greenland and Antarctica.

Abstract & Methods

See ./ms.tex

Results

Note

All figures are width-proportional within and between each other.

./fig_aq_gl.png

./fig_aq_parts.png

Implementation

Greenland

CodeTermValueI/OPeriodSourceComment
RFRainfall45I2000-2019fettweis_2020
CDCondensation5I2000-2019fettweis_2020guesstimate
DPDeposition10I2000-2019fettweis_2020guesstimate
SFSnowfall685I2000-2019fettweis_2019
RFZRefreezing195IO2000-2019fettweis_2020RFZ = ME + RF - RU
EVEvaporation10O2000-2019fettweis_2020guesstimate
RURunoff440O2000-2019fettweis_2020
BMBasal melting20Osteadykarlsson_2020
DYNDischarge490-2000-2019mankoff_2020Submarine melting + calving
SUBSubmarine melting245Oenderlin_201350 % of discharge
ICECalving245Oenderlin_201350 % of discharge
GZRETGrounding line retreat ?5OSee methodsguesstimate
FRLOSSFrontal retreat50O2000-2020kochtitzky_2023
FRGAINFrontal advance0ONone in GL
SUSublimation60O2000-2019fettweis_2020
DDDrawdown-Derivedsum(O) - sum(I)
MGMass gain-Derivedsum(I) - sum(O)

Reset

trash G_GL
[[ -e ./G_GL ]] || grass -e -c EPSG:3413 ./G_GL

Basal melt

GZ retreat

From Millan (2022) http://doi.org/10.5194/tc-16-3021-2022

  • Gz retreat is ~0.13 km/yr (Fig. 3a)
  • Ice velocity is ~1200 m/yr (Fig. 3b) (not needed)
  • 20 km wide

Rates are higher per Ciraci (2023) http://doi.org/10.1073/pnas.2220924120, but

  • Ice surface close to flotation near GZ, and shelf is ~500 m thick, so estimate 600 m ice.

Therefore, gz retreat in Gt/year is width * thick * retreat rate * density

frink "0.13 km/yr * 20 km * 600 m * 917 kg/m^3 -> Gt/yr"
1.43052

Assume similar from other ice shelves too, for a total of ~5 Gt/yr GZ retreat in Greenland.

SMB

g.mapset -c MAR

ncdump -v TIME dat/MARv3.12-GRD-15km-annual.nc4 # 20-39 = 2000-2019
ncra --overwrite -d TIME,20,39 dat/MARv3.12-GRD-15km-annual.nc4 tmp/MAR_GL.nc

ncdump -v X10_110 tmp/MAR_GL.nc # 101
ncdump -v Y20_200 tmp/MAR_GL.nc # 181
g.region w=$(( -645000 - 7500 )) e=$(( 855000 + 7500 )) s=$(( -3357928 - 7500 )) n=$((-657928 + 7500 )) res=15000 -p

var=SF # debug
for var in SF RF RU SU ME SMB EVA CON DEP SUB MSK AREA; do
  r.in.gdal -o input=NetCDF:tmp/MAR_GL.nc:${var} output=${var}
  r.region -c map=${var}
done

r.mapcalc "GL_ice_all = (MSK > 50) & ((x()-y()) > 520000)" # Limit to ice and remove Canada
r.clump input=GL_ice output=clumps --o
main_clump=$(r.stats -c -n clumps sort=desc | head -n2 | tail -n1 | cut -d" " -f1)
r.mapcalc "GL_ice = if(clumps == ${main_clump}, 1, null())"
r.mask raster=GL_ice --o

# scale
## units are mm.w.eq. per grid cell. Grid cell areas are in km^2
## + mm.w.eq. -> m w.eq.: /1E3
## + m w.eq -> kg: *1E3
## + area in km^2 -> m^2: *1E3*1E3
## + kg -> Gt: /1E12
# ds = ds/1E3 * 1E3 * ds['AREA']*1E3*1E3 / 1E12
for var in SF RF RU SU ME SMB EVA CON DEP SUB; do
  r.mapcalc "${var} = (${var}/1000) * 1000 * (AREA * 1000*1000) / exp(10,12)"
done
r.mask -r

r.mapcalc "RFZ = ME + RF - RU"
for var in SF RF RU ME SMB EVA CON DEP SUB RFZ; do
  echo ${var} $(r.univar -g ${var} | grep sum)
done
SF sum=686.768815213334
RF sum=45.5535346610575
RU sum=440.665680238757
ME sum=589.542715610605
SMB sum=235.536411205988
EVA sum=7.9188290228966
CON sum=2.15906279235185
DEP sum=12.2697684982692
SUB sum=61.8983408836194
RFZ sum=194.430570032905

Discharge

import pandas as pd
df = pd.read_csv('/home/kdm/data/Mankoff_2020/ice/GIS_D.csv', index_col=0, parse_dates=True)

df = df['2000-01-01':'2019-12-31']
df.resample('YS').mean().mean().round().astype(int).values[0]
487

Antarctica

CodeTermValueI/OPeriodSourceComment
SMBSurface mass balance2021IO1979-2008rignot_2019 table 12098AQ-77Islands=2021
SMBSurface mass balance2582IO2000-2019fettweis_2020WITH shelves
DYNDischarge2229O2009-2017rignot_2019 table 12306-77

Davison 2023

DISDischarge1838O??davison_2023SUM(C7:C168)
BMBasal melt902OSUM(E7:E168)
CCalving1348OSUM(G7:G168)
SMBSurface mass balance412ISUM(I7:I168) shelf only

902+1348-412 = 1838

CodeTermValueI/OPeriodSourceComment
RFRainfall5I2000-2019fettweis_2020
CDCondensation1I2000-2019fettweis_2020
DPDeposition80I2000-2019fettweis_2020
SFSnowfall2750I2000-2019fettweis_2020
RFZRefreezing105IO2000-2019fettweis_2020
EVEvaporation5O2000-2019fettweis_2020
RURunoff10O2000-2019fettweis_2020
BMBasal melting70O-van-liefferinge_2013
DYNDischarge1090 + 1545-
SUBSubmarine melting1090Odavison_2023
ICECalving1545Odavison_2023
GZRETGrounding line retreat ?1O45 in Amundsen from Davison email
FRLOSSFrontal retreat79+122+145-1Ogreene_2022
FRGAINFrontal advance181+1+103Ogreene_2022
SUSublimation230O2000-2019fettweis_2020
DDDrawdown-Derived.= sum(O) - sum(I)
MGMass gain-Derived.= sum(I) - sum(O)
Exported aq_baseline to ./dat/aq_baseline.csv

Reset

trash G_AQ
[[ -e ./G_AQ ]] || grass -e -c EPSG:3031 ./G_AQ

Masks: East, West, Peninsula, Islands, Grounded and Shelves

grass ./G_AQ/PERMANENT

v.in.ogr input=${DATADIR}/NSIDC/NSIDC-0709.002/1992.02.07/IceBoundaries_Antarctica_v02.shp output=basins

g.region vector=basins res=10000 -pas

v.db.select map=basins|head
v.db.select -c map=basins columns=Regions | sort | uniq # East West Peninsula Islands
v.db.select -c map=basins columns=TYPE | sort | uniq # FL GR IS (float, ground, island)

v.to.rast input=basins output=east_GR use=val val=1 where='((Regions == "East") AND (TYPE == "GR"))'
v.to.rast input=basins output=east_FL use=val val=11 where='((Regions == "East") AND (TYPE == "FL"))'
v.to.rast input=basins output=west_GR use=val val=2 where='((Regions == "West") AND (TYPE == "GR"))'
v.to.rast input=basins output=west_FL use=val val=12 where='((Regions == "West") AND (TYPE == "FL"))'
v.to.rast input=basins output=peninsula_GR use=val val=3 where='((Regions == "Peninsula") AND (TYPE == "GR"))'
v.to.rast input=basins output=peninsula_FL use=val val=13 where='((Regions == "Peninsula") AND (TYPE == "FL"))'
v.to.rast input=basins output=islands use=val val=4 where='(TYPE == "IS")'
r.patch input=east_GR,east_FL,west_GR,west_FL,peninsula_GR,peninsula_FL,islands output=basins
r.category basins separator=":" rules=- << EOF
1:East (Grounded)
11:East (Floating)       
2:West (Grounded)
12:West (Floating)
3:Peninsula (Grounded)
13:Peninsula (Floating)
4:Islands	     
EOF

r.colors map=basins color=viridis

SMB (MAR)

g.mapset -c MAR

ncdump -v TIME dat/MARv3.12-ANT-35km-annual.nc4 # 20-39 = 2000-2019
ncra --overwrite -d TIME,20,39 dat/MARv3.12-ANT-35km-annual.nc4 tmp/MAR_AQ.nc

ncdump -v X tmp/MAR_AQ.nc # 176
ncdump -v Y tmp/MAR_AQ.nc # 148
g.region w=$(( -3010000 - 17500 )) e=$(( 3115000 + 17500 )) s=$(( -2555000 - 17500 )) n=$(( 2590000 + 17500 )) res=35000 -p

var=SF # debug
for var in SF RF RU ME SMB EVA CON DEP SUB MSK AREA; do
  r.in.gdal -o input=NetCDF:tmp/MAR_AQ.nc:${var} output=${var}
  r.region -c map=${var}
done

# scale
## units are mm.w.eq. per grid cell. Grid cell areas are in km^2
## + mm.w.eq. -> m w.eq.: /1E3
## + m w.eq -> kg: *1E3
## + area in km^2 -> m^2: *1E3*1E3
## + kg -> Gt: /1E12
# ds = ds/1E3 * 1E3 * ds['AREA']*1E3*1E3 / 1E12
for var in SF RF RU ME SMB EVA CON DEP SUB; do
  r.mapcalc "${var} = (${var}/1000) * 1000 * (AREA * 1000*1000) / exp(10,12)"
done

r.mapcalc "RFZ = ME + RF - RU"

Stats

r.mask --o raster=basins@PERMANENT --q maskcats="1 thru 3 10 thru 20" # drop 0 and Islands
for var in SF RF RU ME SMB EVA CON DEP SUB RFZ; do
  echo -n "${var}"
  r.univar -gt map=${var} zones=basins@PERMANENT | cut -d"|" -f2,13 | column -s"|" -t | sed 's/label.*//'
  r.univar -g ${var} | grep sum
  echo "#"; echo "#"
done
r.mask -r --q
[Raster MASK present]

SF
East (Grounded)       1383.5170055779
West (Grounded)       688.480301166503
Peninsula (Grounded)  270.324169435779
East (Floating)       172.41137746281
West (Floating)       180.27593549343
Peninsula (Floating)  56.6841289993761
sum=2751.6929181358

RF
East (Grounded)       0.53173174128475
West (Grounded)       0.22197445329555
Peninsula (Grounded)  2.1110471117335
East (Floating)       0.842541426357001
West (Floating)       0.4498711186449
Peninsula (Floating)  2.307238481508
sum=6.46440433282369

RU
East (Grounded)       1.53022074184155
West (Grounded)       0.0059355454226
Peninsula (Grounded)  1.8834556710495
East (Floating)       1.5089940427256
West (Floating)       0.0304982132294
Peninsula (Floating)  4.35827769837335
sum=9.317381912642

ME
East (Grounded)       15.4802688824357
West (Grounded)       4.3355515366436
Peninsula (Grounded)  17.3666975625616
East (Floating)       26.407263870097
West (Floating)       9.2284017518
Peninsula (Floating)  34.440989714197
sum=107.259173317735

SMB
East (Grounded)       1268.09953548228
West (Grounded)       678.756441138014
Peninsula (Grounded)  262.863634061747
East (Floating)       153.249402230902
West (Floating)       177.921656614902
Peninsula (Floating)  51.4267222532681
sum=2592.31739178111

EVA
East (Grounded)       0.60507218584105
West (Grounded)       0.195446970963
Peninsula (Grounded)  0.68838986544795
East (Floating)       0.70256713317005
West (Floating)       0.2422464141299
Peninsula (Floating)  0.70164322473235
sum=3.1353657942843

CON
East (Grounded)       0.00144321189675
West (Grounded)       0.00237202472355
Peninsula (Grounded)  0.01264862146645
East (Floating)       0.0031724865901
West (Floating)       0.0019547481581
Peninsula (Floating)  0.03522553443475
sum=0.0568166272697001

DEP
East (Grounded)       36.3995731586915
West (Grounded)       22.3853823370692
Peninsula (Grounded)  5.3339705229687
East (Floating)       5.70103389655939
West (Floating)       6.05853236904585
Peninsula (Floating)  1.5062480433876
sum=77.384740327722

SUB
East (Grounded)       150.624482553971
West (Grounded)       32.2385823068849
Peninsula (Grounded)  12.2348077688132
East (Floating)       23.4661462650309
West (Floating)       8.5418917438099
Peninsula (Floating)  3.94097993607855
sum=231.046890574587

RFZ
East (Grounded)       14.4817798818789
West (Grounded)       4.55159044451654
Peninsula (Grounded)  17.5942890032456
East (Floating)       25.7408112537284
West (Floating)       9.64777465721551
Peninsula (Floating)  32.3899504973317
sum=104.406195737917

[Raster MASK present]

Basal melt

Van Liefferinge (2013) http://doi.org/10.5194/cp-9-2335-2013

Convert MAT file to XYZ for importing into GRASS

import scipy as sp
import numpy as np
import pandas as pd

mat = sp.io.loadmat('/home/kdm/data/Van_Liefferinge_2023/Melt_Mean_Std_15exp.mat')
X = mat['X'].flatten() * 1E3 # convert from km to m
Y = mat['Y'].flatten() * 1E3
m = mat['MeanMelt'].flatten() / 10 # cm to mm

melt = pd.DataFrame(np.array([X,Y,m]).T, columns=['x','y','melt'])\
         .dropna()
melt.to_csv('./tmp/melt.csv', header=False, index=False)
melt.head()
xymelt
1487411.045e+06-2.14e+061e-09
1498591.03e+06-2.135e+060.00146608
1498601.035e+06-2.135e+060.000266042
1498611.04e+06-2.135e+061e-09
1498621.045e+06-2.135e+060.00045698
grass ./G_AQ/PERMANENT
g.mapset -c liefferinge_2023
r.in.xyz input=./tmp/melt.csv output=melt sep=, --o
echo "All: " $(r.univar -g map=melt | grep sum)
echo ""
r.univar -gt map=melt zones=basins | cut -d"|" -f2,13 | column -s"|" -t
All:  sum=69.3982306335468


label                 sum
East (Grounded)       45.7178033424208
West (Grounded)       18.0714170862276
Peninsula (Grounded)  2.93302497694997
Islands               0.279139711405429
East (Floating)       1.03624592705523
West (Floating)       0.781445329564939
Peninsula (Floating)  0.254017664974735

Antarctic Ice shelves

Submarine melt
import pandas as pd

fname = '~/data/Davison_2023/adi0186_table_s2.xlsx'

loc = pd.read_excel(fname, sheet_name='Total mass changes', index_col = 0, usecols = 'B,C,D', skiprows = 4)
loc = loc.drop('Antarctic Ice Shelves')


df = pd.read_excel(fname, sheet_name='Melt',
                   index_col = 1, skiprows = 3)

# drop uncertainty columns
unc = []
for c in df.columns:
     if type(c) == str:
          if c[0:8] == 'Unnamed:':
               unc.append(c)
df = df.drop(columns = unc)
df = df[df.columns[3:]]
df = df.iloc[1:]

df = pd.DataFrame(df.mean(axis='columns'))
df.columns = ['Mass']

df = loc.join(df)

import geopandas as gpd
fname = '~/data/NSIDC/NSIDC-0709.002/1992.02.07/IceBoundaries_Antarctica_v02.shp'
ew = gpd.read_file(fname)

df = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df['longitude'],df['latitude']), crs="EPSG:4326")
df = df.to_crs('epsg:3031')
e = ew.to_crs('epsg:3031')

idx = ew.sindex.nearest(df['geometry'], return_all=False)
df['Region'] = ''
for dfidx,ewidx in idx.T:
    arr = df.iloc[dfidx].copy(deep=True)
    arr['Region'] = ew.iloc[ewidx]['Regions']
    df.iloc[dfidx] = arr
df = df.drop(columns=['latitude','longitude'])
    
df.loc['Total'] = [df['Mass'].sum(), None, 'All']

df[['Mass','Region']].groupby('Region').sum().drop('Islands').round()
/tmp/ipykernel_1220053/4002790239.py:43: FutureWarning: The behavior of DataFrame concatenation with empty or all-NA entries is deprecated. In a future version, this will no longer exclude empty or all-NA columns when determining the result dtypes. To retain the old behavior, exclude the relevant entries before the concat operation.
  df.loc['Total'] = [df['Mass'].sum(), None, 'All']
RegionMass
All1087.2
East321.123
Peninsula172.74
West593.128
Calving

Same as above, different sheet. Reuses variables from above, run that first.

fname = '~/data/Davison_2023/adi0186_table_s2.xlsx'
df = pd.read_excel(fname, sheet_name='Calving',
                   index_col = 1, skiprows = 3)

# drop uncertainty columns
unc = []
for c in df.columns:
     if type(c) == str:
          if c[0:8] == 'Unnamed:':
               unc.append(c)
df = df.drop(columns = unc)
df = df[df.columns[3:]]
df = df.iloc[1:]

df = pd.DataFrame(df.mean(axis='columns'))
df.columns = ['Mass']

df = loc.join(df)

df = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df['longitude'],df['latitude']), crs="EPSG:4326")
df = df.to_crs('epsg:3031')
e = ew.to_crs('epsg:3031')

idx = ew.sindex.nearest(df['geometry'], return_all=False)
df['Region'] = ''
for dfidx,ewidx in idx.T:
    arr = df.iloc[dfidx].copy(deep=True)
    arr['Region'] = ew.iloc[ewidx]['Regions']
    df.iloc[dfidx] = arr
df = df.drop(columns=['latitude','longitude'])
    
df.loc['Total'] = [df['Mass'].sum(), None, 'All']

df[['Mass','Region']].groupby('Region').sum().drop('Islands').round()
/tmp/ipykernel_1220053/519284814.py:32: FutureWarning: The behavior of DataFrame concatenation with empty or all-NA entries is deprecated. In a future version, this will no longer exclude empty or all-NA columns when determining the result dtypes. To retain the old behavior, exclude the relevant entries before the concat operation.
  df.loc['Total'] = [df['Mass'].sum(), None, 'All']
RegionMass
All1543.55
East560.3
Peninsula236.416
West745.595
GZ retreat

Email from Davison

Ice ShelfMass change due to grounding line migration from 1997 to 2021 (Gt)Error (Gt)
Pine Island22040
Thwaites23025
Crosson20025
Dotson42080

(220+230+200+420)/(2021-1997) = 44.5833333333

Frontal Retreat from Greene 2022

[greene_Supplementary_Table_1.xlsx](https://github.com/user-attachments/files/15598602/greene_Supplementary_Table_1.xlsx)

I think the data in the attached spreadsheet from [Greene et al., 2022 ](https://doi.org/10.1038/s41586-022-05037-w) is everything needed for ice-shelf mass-change resulting from frontal advance/retreat, so in Excel `=BI189-O189` gives Antarctica’s net retreat from 1997 to 2021. Change the column to adjust the time period.

BI189 = 24596304.0 BI189 = 2021.2 Q189 = 24597630.0 Q189 = 2000.2

(24596304.0 - 24597630.0) / (2021.2-2000.2) = -63.1428571429

But we need to recreate this in code so we can split by east/west/peninsula

import pandas as pd
import geopandas as gpd
fname = "~/data/Greene_2022/data/greene_Supplementary_Table_1.xlsx"

df = pd.read_excel(fname, sheet_name='greene_iceshelf_area_and_mass',
                    index_col = 1, skiprows = 4)
df = df.rename(columns={'Unnamed: 2':'lat',
                        'Unnamed: 3':'lon'})

# drop uncertainty columns
unc = []
for c in df.columns:
    if type(c) == str:
        if c[0:8] == 'Unnamed:':
            unc.append(c)
df = df.drop(columns = unc)
df = df[['lat','lon',2000.2,2021.2]]
df = df.iloc[1:]

# Remove last two rows
aq = df.loc['Antarctica']
other = df.loc['Other']
df = df.iloc[:-2]
print(df.sum())
print("")
print(aq)
print("")
print(other)
lat       -12882.373098
lon         6279.268331
2000.2    682491.281291
2021.2    681213.775349
dtype: object

lat            -90
lon          every
2000.2    24597630
2021.2    24596304
Name: Antarctica, dtype: object

lat            NaN
lon            NaN
2000.2    23915136
2021.2    23915090
Name: Other, dtype: object
shelf = df.sum()
print("All AQ loss: ", (aq[2021.2] - aq[2000.2]) / (2021-2000))
print("Named shelf loss: ", (shelf[2021.2] - shelf[2000.2]) / (2021-2000))
print("Other loss: ", (other[2021.2] - other[2000.2]) / (2021-2000))
print("Named + Other: ", (((other + shelf)[2021.2] - (other + shelf)[2000.2]) / (2021-2000)))
print("Named %: ", 2.19/63.02*100)
All AQ loss:  -63.142857142857146
Named shelf loss:  -60.83361628651619
Other loss:  -2.1904761904761907
Named + Other:  -63.02409247699238
Named %:  3.4750872738813077
import geopandas as gpd
fname = '~/data/NSIDC/NSIDC-0709.002/1992.02.07/IceBoundaries_Antarctica_v02.shp'
ew = gpd.read_file(fname)
ew.drop(columns=['geometry']).head()
NAMERegionsSubregionsTYPEAsso_Shelf
0LarsenEPeninsulaIpp-JGRLarsenE
1Dawson_LambtonEastnanFLnan
2AcademyEastJpp-KGRFilchner
3Brunt_StancombEastK-AGRBrunt_Stancomb
4Riiser-LarsenEastK-AGRRiiser-Larsen
gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df['lon'],df['lat']), crs="EPSG:4326")

gdf = gdf.to_crs('epsg:3031')
ew = ew.to_crs('epsg:3031')

idx = ew.sindex.nearest(gdf['geometry'], return_all=False)
gdf['Region'] = ''
for gdfidx,ewidx in idx.T:
     arr = gdf.iloc[gdfidx].copy(deep=True)
     arr['Region'] = ew.iloc[ewidx]['Regions']
     gdf.iloc[gdfidx] = arr

gdf.head()

gdf.loc['Total'] = gdf.sum(axis='rows')
gdf.loc['Total', 'Region'] = 'All'

gdf['frontal change'] = (gdf[2021.2] - gdf[2000.2]) / (2021.2-2000.2)
pos = gdf[gdf['frontal change'] > 0]
neg = gdf[gdf['frontal change'] <= 0]
# gdf

print('neg', neg[['Region','frontal change']].groupby('Region').sum().round().abs())
print('')
print('pos', pos[['Region','frontal change']].groupby('Region').sum().round().abs())
print('')
print('all', gdf[['Region','frontal change']].groupby('Region').sum().round().abs())
neg            frontal change
Region                   
All                  61.0
East                 79.0
Peninsula           122.0
West                145.0

pos            frontal change
Region                   
East                181.0
Peninsula             1.0
West                103.0

all            frontal change
Region                   
All                  61.0
East                102.0
Peninsula           121.0
West                 42.0

Misc

Export tables to CSVs

(save-excursion
  (goto-char (point-min))
  (re-search-forward (concat "^#\\+name: " tbl) nil t)
  (next-line)
  (org-table-export (concat "./dat/" tbl ".csv") "orgtbl-to-csv")
  ;;(shell-command-to-string (concat "head " tbl ".csv"))
  (message (concat "Exported " tbl " to " (concat "./dat/" tbl ".csv")))
  )

Convert PDFs to PNG

About

Sankey diagram for ice sheet mass flow

Topics

Resources

Stars

Watchers

Forks

Packages

 
 
 

Contributors