Skip to content

Latest commit

 

History

History
2389 lines (1913 loc) · 78.7 KB

workbook.org

File metadata and controls

2389 lines (1913 loc) · 78.7 KB
Raw

Freshwater Anomalies

Table of contents

Introduction

Summary

This workbook computes freshwater anomaly from Greenland and Antarctica for the NASA GISS ModelE.

  • Anomaly is defined as mass loss since 1990
  • Spatially, ModelE distributes a single value according to a ‘distribution field’ - a floating point mask that sums to 1.

Numbers and units

  • 1 Gt = 1 km3 (freshwater)
  • 360 Gt \(≈\) 1 mm SLR
  • 1 Sv = 31556 km3 yr-1 (\(∼\) π 104 km3 yr-1)

Overview numbers

import xarray as xr
ds = xr.open_dataset('/home/kdm/data/Mankoff_2021/461/MB_region.nc')\
    .sum(dim='region')\
    .resample({'time':'YS'})\
    .sum()['MB']

dd = ds.sel({'time':slice('1990-01-01','1999-12-31')}).mean(dim='time')
print('1990-1999: ', dd.data, dd.data/-362)

dd = ds.sel({'time':slice('1990-01-01','2019-12-31')}).mean(dim='time')
print('1990-2019: ', dd.data, dd.data/-362)

dd = ds.sel({'time':slice('2000-01-01','2019-12-31')}).mean(dim='time')
print('2000-2019: ', dd.data, dd.data/-362)

dd = ds.sel({'time':slice('2010-01-01','2019-12-31')}).mean(dim='time')
print('2010-2019: ', dd.data, dd.data/-362)

dd = ds.sel({'time':slice('2010-01-01','2019-12-31')}).resample({'time':'YS'}).mean()
dd = dd.where((dd.time.dt.year != 2019) & (dd.time.dt.year != 2012)).mean()
print('2010-2019 (excl 2012 & 2019): ', dd.data, dd.data/-362)
  • Greenland net ice mass loss citep:mankoff_2021
    • \(∼\)160 Gt yr-1 (0.44 mm yr-1) since 1990
    • \(∼\)200 Gt yr-1 (0.55 mm yr-1) since 2000
    • \(∼\)245 Gt yr-1 (0.68 mm yr-1) since 2010

{{{SKIPLINE}}}

  • \(∼\)500 Gt yr-1 is dynamic (outlet glacier discharge; citet:mankoff_2020_solid)
    • Distributed \(∼\)50 % ±40 % between submarine melt & icebergs citep:enderlin_2013
    • Fairly steady intra- and inter- annual
  • \(∼\)300–500 Gt yr-1 is surface melt and runoff citep:mankoff_2017_VHD,mankoff_2020_liquid
    • Some surface rivers, most submarine
    • Highly seasonal (most in summer)
    • Highly variable interannual

{{{SKIPLINE}}}

  • From the above, ice sheet mass gain (snowfall) should be \(∼\)600–800 Gt yr-1
  • Another \(∼\)120 Gt yr-1 land runoff (rain/snow on non-ice-sheet land) citep:mankoff_2020_liquid

Global ice mass change

citet:slater_2021 provides global ice mass change. This includes ice shelf mass loss, something not detected by GRACE and excluded from most mass change products. We use citet:slater_2021 in Antarctica.

{{{SKIPLINE}}}

In Greenland, we use citet:mankoff_2021 because it begins earlier, is continually updating, and ice shelf melting and thinning is negligible in Greenland relative to the other freshwater terms.

Slater 2021 Fig. 4 recreation

Original

Original

./fig/slater_2021_fig4_orig.png

Recreation

Recreation

S2021_kw = {'parse_dates':True, 'index_col':0}
SH_ice = pd.read_csv('~/data/Slater_2021/SH_seaice_cumul_1994_2017_annual.csv', **S2021_kw)\
           .rename(columns={"Cumulative mass change (Gt)":"Antarctic Sea Ice"})
NH_ice = pd.read_csv('~/data/Slater_2021/NH_seaice_cumul_1994_2017_annual.csv', **S2021_kw)\
           .rename(columns={"Cumulative mass change (Gt)":"Arctic Sea Ice"})
shelf_calving = pd.read_csv('~/data/Slater_2021/iceshelves_calving_cumul_1994_2017_annual.csv', **S2021_kw)\
                  .rename(columns={"Cumulative mass change (Gt)":"Ice Shelf Calving"})
shelf_thinning = pd.read_csv('~/data/Slater_2021/iceshelves_thinning_cumul_1994_2017_annual.csv', **S2021_kw)\
                      .rename(columns={"Cumulative mass change (Gt)":"Ice Shelf Thinning"})
AQ = pd.read_csv('~/data/Slater_2021/AIS_cumul_1994_2017_annual.csv', **S2021_kw)\
                      .rename(columns={"Cumulative mass change (Gt)":"Antarctica"})
GL = pd.read_csv('~/data/Slater_2021/GrIS_cumul_1994_2017_annual.csv', **S2021_kw)\
                      .rename(columns={"Cumulative mass change (Gt)":"Greenland"})
glacier = pd.read_csv('~/data/Slater_2021/Glacier_cumul_1994_2017_annual.csv', **S2021_kw)\
                      .rename(columns={"Cumulative mass change (Gt)":"Glaciers"})

S2021_kw = {'left_index':True, 'right_index':True, 'how':'outer'}
S2021 = SH_ice.merge(NH_ice, **S2021_kw)\
           .merge(shelf_calving, **S2021_kw)\
           .merge(shelf_thinning, **S2021_kw)\
           .merge(AQ, **S2021_kw)\
           .merge(GL, **S2021_kw)\
           .merge(glacier, **S2021_kw)
S2021.index.name = 'Date'

./fig/slater_2021_fig4.png

Issues

  • Although citet:slater_2021 provide global ice mass change, the data is only from 1994 through 2017
  • We will use other products with longer records where possible, but use citet:slater_2021 for QC
  • We will extrapolate citet:slater_2021 forward and backward if no other product exists
    • i.e., Antarctic ice shelves

Greenland

Mankoff 2021

citet:mankoff_2021 provide main ice sheet mass loss from 1840 through next week.

{{{SKIPLINE}}}

Most of the anomaly is SMB

./fig/mankoff_2021_fig2.png

citet:mankoff_2021

Fetch historical observations

mkdir input
wget -nc  https://dataverse.geus.dk/api/access/datafile/:persistentId?persistentId=doi:10.22008/FK2/OHI23Z/NBMCEK -O ./input/mankoff_2021.csv

wget -nc https://dataverse.geus.dk/api/access/datafile/:persistentId?persistentId=doi:10.22008/FK2/OHI23Z/XQHQOB  -O ./input/mankoff_2021.nc

Compare Mankoff 2021 & Slater 2021

Figure

./fig/mankoff_v_slater.png

Text

  • Upper panel: annual change (and difference between two products [dashed gray])
  • Lower panel: cumulative

Greenlandic Peripheral glaciers

  • Need to determine where Greenlandic peripheral glaciers are in the citet:slater_2021 data
    • [ ] Greenland?
    • [ ] Glaciers?
  • Presumably they’re included in the Greenland ice mass, because citet:slater_2021 uses GRACE, which cannot distinguish main ice from peripheral
  • citet:mankoff_2021 do not include peripheral glaciers, but can via a scaling factor
  • Note on previous slide - differences due to peripheral glacier treatment not detectable

Greenlandic icebergs vs. (runoff & submarine melt)

Text

Figure

./fig/Greenland-Iceberg-Chart.jpg

Greenland mass loss by ROI [Gt yr-1]

#+RESULTS[(2023-02-21 12:13:20) 3fac654a4b73926e98be6a010a19f62289822a8c]: mb_roi

regionMass
NE24.8409
CE3.15127
SE11.7459
SW23.1003
CW50.554
NW60.4492
NO36.0279

#+RESULTS[(2023-01-27 07:27:02) 67c3d0752fd16cd80ef00bbf168408a6e595a5b9]: iceberg_roi

regionDischarge
CE74.5944
CW79.8808
NE25.0692
NO24.1457
NW99.2541
SE148.667
SW20.5418

#+RESULTS[(2023-01-27 07:27:27) af0b8342115f1512478aac53eac87a5db1f49271]: stream_roi

regionRunoff
CE56.8306
CW43.505
NE49.1332
NO40.4823
NW50.519
SE74.169
SW118.234
regionMassDischargeRunoff
NE252549
CE17557
SE1114974
SW2321118
CW508044
NW609951
NO362440
TOTAL206472433
  • Recall: \(∼\)50 % discharge is submarine melt (runoff)
  • Distance from coast is all \(∼\)0 (from previous slide)

Greenland mass loss by ROI [%]

regionMassDischargeRunoff
CE11613
CW241710
NE12511
NO1759
NW292112
SE53117
SW11427
TOTAL999999

Greenland mass balance (annual)

timeMB
1990-137.596
1991-76.7293
199287.1372
1993-90.7483
1994-113.819
1995-211.913
1996132.208
19977.60407
1998-241.779
1999-47.0087
2000-77.1367
2001-26.1103
2002-142.517
2003-167.223
2004-165.81
2005-168.432
2006-239.821
2007-257.306
2008-201.225
2009-242.975
2010-376.771
2011-336.25
2012-429.338
2013-107.915
2014-184.639
2015-213.906
2016-255.964
2017-102.566
2018-75.8065
2019-425.994
2020-188.184
2021-211.337
2022-65.4334
20239.32125

Greenland mass balance (monthly)

All mass balance

MB
1990-014.72
1990-02-5.66418
1990-03-6.43677
1990-044.86196
1990-0520.2296
1990-06-78.8659
1990-07-155.311
1990-08-70.179
1990-0928.6134
1990-1029.0614
1990-1145.0175
1990-1246.3573
1991-0130.9776
1991-0237.0555
1991-038.78264
1991-04-10.715
1991-0530.2405
1991-06-70.9794
1991-07-170.2
1991-08-41.12
1991-099.11271
1991-1050.9149
1991-116.70501
1991-1242.4959
1992-0144.4852

Negative mass balance

MB
1990-010
1990-02-5.66418
1990-03-6.43677
1990-040
1990-050
1990-06-78.8659
1990-07-155.311
1990-08-70.179
1990-090
1990-100
1990-110
1990-120
1991-010
1991-020
1991-030
1991-04-10.715
1991-050
1991-06-70.9794
1991-07-170.2
1991-08-41.12
1991-090
1991-100
1991-110
1991-120
1992-010

Comments on previous slide

  • Doesn’t make sense to estimate from total mass balance
    • Feb and March with mass loss (probably) comes from steady ice discharge and low snowfall
  • Solid discharge roughly steady intra-annual
  • Better to estimate FW anomaly from stream discharge product
    • Need to define some monthly baseline
    • Then consider anomaly from baseline

Stream discharge: baseline [Gt/month] (1960–1989)

MonthCECWNENONWSESWTOTAL
100000000
200000000
300000000
400000000
511000125
697435111654
720152015182442155
8129759182686
9210015512
1000000001
1100000000
1200000000
TOTAL43333124336091314

Stream discharge: anomaly [Gt/month] ((1990–2019)-base)

MonthCECWNENONWSESWTOTAL
1-000000-0-0
2-0-0000000
30-00000-00
400000000
500-000112
6345453833
754101083748
8323234825
9200-00215
1000000001
1100000000
120000-0000
TOTAL13111816171326114

Stream discharge: anomaly [%] [(anom - base)/base]

MonthCECWNENONWSESW
1nannannannannannannan
2nannannannannannannan
3nannannannannannannan
4nannannannannannannan
5nan23nannannan5349
63861114126893154
724285164471417
821224847372031
91155nannannan4223
10nannannannannannannan
11nannannannannannannan
12nannannannannannannan

Summary

  • Some increase in solid ice discharge (10 %)
    • …of which 50 % is submarine melt
  • Most increase in liquid runoff
    • Only summer in JJA

{{{SKIPLINE}}}

  • Can build baseline distribution maps for solid & liquid
  • Can build anomaly distribution maps for solid & liquid
  • Can provide annual or monthly solid and liquid forcing

Output

# Contact: Ken Mankoff <ken.mankoff@nasa.gov>
# Code: https://github.com/NASA-GISS/freshwater-anomaly
# Data: Mankoff /et al./ (2021)
# Summary: Annual mass change from Greenland.
# Year	Gt/yr
1990	137.6
1991	76.7
1992	-87.1
1993	90.7
1994	113.8
1995	211.9
1996	-132.2
1997	-7.6
1998	241.8
1999	47.0
2000	77.1
2001	26.1
2002	142.5
2003	167.2
2004	165.8
2005	168.4
2006	239.8
2007	257.3
2008	201.2
2009	243.0
2010	376.8
2011	336.2
2012	429.3
2013	107.9
2014	184.6
2015	213.9
2016	256.0
2017	102.6
2018	75.8
2019	426.0

Output only positive values: notes

  • Use uncertainty from citet:mankoff_2021 published with data
  • Use published (1 σ) uncertainty because that uncertainty treatment is already conservative and if using 2 σ, mass gain occurs most year for lo estimates.

Output only positive values: results

timeMB_loMBMB_hiMB_err
199053.6137.6221.684
1991076.7164.787.9
199200077.6
199303.616785.8
19940113.8187.273.4
19954.8211.9288.576.5
199600084.3
1997002580.5
19980102323.681.9
1999047132.385.3
2000077.1154.977.8
2001026.111184.9
20020142.5225.783.2
20030167.2261.994.7
20040165.8249.683.8
200544.9168.4263.495
2006163.1239.8316.576.7
2007161.9257.3352.795.4
2008104201.2298.497.2
2009157.7243328.385.3
2010283.2376.8470.393.6
2011242.2336.2430.394.1
2012319.5429.3539.2109.8
201319.4107.9196.488.5
201491184.6278.393.7
2015118.5213.9309.395.4
2016154.2256357.8101.8
201711.7102.6193.490.9
2018075.8158.482.6
2019319.5426525.799.7
TOTAL2249.24890.17531.42641.3

Antarctica

IMBIE 2018

{{{SKIPLINE}}}

However, IMBIE does not have ice shelf thinning or calving

IMBIE 2018 vs Slater 2021

Text

Bottom panel: Blue line is repeated from middle panel (Slater 2021)

Figure

md5sum ~/data/IMBIE/*

./fig/slater_v_imbie.png

Antarctic icebergs

./fig/budge_2018.png

Iceberg distribution

import pandas as pd
import glob

root = '/home/kdm/data/Budge_2018/consol'
csvs = glob.glob(root + '/*.csv')

lon,lat = [],[]
for csv in csvs:
    df = pd.read_csv(csv, parse_dates=True, index_col='date')
    for i,c in enumerate(df.columns):
        if df[c].min() == df[c].max():
            continue
        if np.all(df[c] < 0):
            alat = df[c]
            alon = df[df.columns[i+1]]
            break

        lon = lon + list(alon)
        lat = lat + list(alat)

# plt.scatter(lon,lat)
pd.DataFrame(np.vstack((lon,lat)).T).to_csv('./tmp/lonlat.csv', index=False, header=None)

EPSG 3031

grass -c epsg:3031 ./G_AQ

eval $(m.proj -i input=tmp/lonlat.csv separator=comma,pipe |r.in.xyz input=- -sg)
g.region -pa n=$n s=-$n e=$(echo -1*$w|bc) w=$w res=100000 -s

m.proj -i input=tmp/lonlat.csv sep=comma,pipe | v.in.ascii input=- output=pts

# d.mon wx0
# d.vect pts

m.proj -i input=tmp/lonlat.csv sep=comma,pipe | r.in.xyz input=- output=bin method=n
r.mapcalc "bin10 = log(10,bin)"

# d.rast bin

EPSG 3031

grass -c epsg:4326 ./G_AQ

# g.region -pas res=0.125 s=-90 n=-45 e=-180 w=180
g.region -pas res=0.5 s=-90 n=-45 e=-180 w=180
v.in.ascii -n input=tmp/lonlat.csv sep=comma output=pts
cat tmp/lonlat.csv | sed 's/$/,1/' | r.in.xyz input=- separator=comma output=bin method=n
r.mapcalc "bin10 = log(10,bin)"
# d.rast bin

FW forcing distribution

Text

  • AQ solids
    • Spatial distribution of solid FW forcing may be unnecessary detail
    • Where icebergs are seen may be where they are not melting (survivorship bias)
    • Can distribute evenly around continent, or evenly where observed, or based on observation density
  • AQ liquids
    • Should be pegged to ice shelf edge
    • Do not have by shelf, only total for the continent

Fig

./fig/Survivorship-bias.svg.png

AQ by ROI

./fig/rignot_2019_table2.png

\(dM\) is anomaly including gains; does not distinguish calving vs. melt Baseline period is 14*360 = 5040 Gt or 5040/(2017-1979+1) = 130 Gt yr-1

citet:rignot_2019

AQ by ROI

./fig/rignot_2019_fig2.png

citet:rignot_2019

Output: notes

From citet:slater_2021 Table 1, uncertainties for Antarctica are:

ComponentUncertainty [Gt yr-1
Ice shelf calving36
Ice shelf thinning39
Antarctic land ice24
MEAN33
  • Assume these are 1 σ uncertainty
  • hi/lo estimates use the mean (33), max (39) or sum(99)

Output: code

# Contact: Ken Mankoff <ken.mankoff@nasa.gov>
# Code: https://github.com/NASA-GISS/freshwater-anomaly
# Data: Slater /et al./ (2021)
# Summary: Annual mass loss from Antarctica.
# Year	Gt/yr
1990	115.8
1991	115.8
1992	115.8
1993	115.8
1994	115.8
1995	110.2
1996	82.9
1997	154.4
1998	187.8
1999	499.9
2000	428.1
2001	430.4
2002	276.9
2003	424.6
2004	488.2
2005	453.6
2006	547.1
2007	679.5
2008	600.6
2009	534.7
2010	585.4
2011	416.6
2012	443.8
2013	499.6
2014	456.9
2015	381.3
2016	157.2
2017	331.8
2018	331.8
2019	331.8

Freshwater distribution mask

Introduction

  • ModelE distributes Greenlandic melt via a fractional mask (sums to 1)
  • We can use the same mask to distribute the freshwater anomaly
  • However, the anomaly is distributed differently than the baseline
  • For example, Fig. 1 of citet:mankoff_2021 shows the SE (included below) has no net mass loss, and no change from <1990 baseline. However, SE should still have an annual baseline meltwater volume flow rate, because winter snowfall is offset by summer melt.

./fig/mankoff_2021_fig1_orig.png

Existing mask

{{{SKIPLINE}}}

Land classification mask (could be useful when creating new mask)

  • 2.5x2 : Z2HX2fromZ1QX1N.BS1.nc
  • 1x1 : OZ1QX1N.BS1.nc
  • 1/60 : etopo_ice_g1m.nc

Contents of GLMELT_144X90_gas.OCN.nc

netcdf GLMELT_144X90_gas.OCN {
dimensions:
	lon = 144 ;
	lat = 90 ;
variables:
	float lon(lon) ;
		lon:units = "degrees_east" ;
	float lat(lat) ;
		lat:units = "degrees_north" ;
	float mask(lat, lon) ;
data:

 lon = -178.75, -176.25, -173.75, -171.25, -168.75, -166.25, -163.75, 
    -161.25, -158.75, -156.25, -153.75, -151.25, -148.75, -146.25, -143.75, 
    -141.25, -138.75, -136.25, -133.75, -131.25, -128.75, -126.25, -123.75, 
    -121.25, -118.75, -116.25, -113.75, -111.25, -108.75, -106.25, -103.75, 
    -101.25, -98.75, -96.25, -93.75, -91.25, -88.75, -86.25, -83.75, -81.25, 
    -78.75, -76.25, -73.75, -71.25, -68.75, -66.25, -63.75, -61.25, -58.75, 
    -56.25, -53.75, -51.25, -48.75, -46.25, -43.75, -41.25, -38.75, -36.25, 
    -33.75, -31.25, -28.75, -26.25, -23.75, -21.25, -18.75, -16.25, -13.75, 
    -11.25, -8.75, -6.25, -3.75, -1.25, 1.25, 3.75, 6.25, 8.75, 11.25, 13.75, 
    16.25, 18.75, 21.25, 23.75, 26.25, 28.75, 31.25, 33.75, 36.25, 38.75, 
    41.25, 43.75, 46.25, 48.75, 51.25, 53.75, 56.25, 58.75, 61.25, 63.75, 
    66.25, 68.75, 71.25, 73.75, 76.25, 78.75, 81.25, 83.75, 86.25, 88.75, 
    91.25, 93.75, 96.25, 98.75, 101.25, 103.75, 106.25, 108.75, 111.25, 
    113.75, 116.25, 118.75, 121.25, 123.75, 126.25, 128.75, 131.25, 133.75, 
    136.25, 138.75, 141.25, 143.75, 146.25, 148.75, 151.25, 153.75, 156.25, 
    158.75, 161.25, 163.75, 166.25, 168.75, 171.25, 173.75, 176.25, 178.75 ;

 lat = -89, -87, -85, -83, -81, -79, -77, -75, -73, -71, -69, -67, -65, -63, 
    -61, -59, -57, -55, -53, -51, -49, -47, -45, -43, -41, -39, -37, -35, 
    -33, -31, -29, -27, -25, -23, -21, -19, -17, -15, -13, -11, -9, -7, -5, 
    -3, -1, 1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 
    35, 37, 39, 41, 43, 45, 47, 49, 51, 53, 55, 57, 59, 61, 63, 65, 67, 69, 
    71, 73, 75, 77, 79, 81, 83, 85, 87, 89 ;

 mask =
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
[...]

Visualization of GLMELT_144X90_gas.OCN.nc

./.ob-jupyter/4ae2299e0176cc83b7260b90276f2c9a1f22f10f.png

Mask editing

Make mask QGIS compatible

import xarray as xr

fname = "GLMELT_144X90_gas.OCN.nc"
ds = xr.open_dataset("./dat/" + fname)

ds["crs"] = True
ds["crs"].attrs["grid_mapping_name"] = "latitude_longitude"

ds["mask"].attrs["grid_mapping"] = "crs"
ds["mask"].attrs["_FillValue"] = 0 # optional

ds.to_netcdf("./dat/" + fname[:-3] + ".crs.nc")

ncap2 -h -O -s 'crs=1B' ./dat/GLMELT_144X90_gas.OCN.nc ./dat/GLMELT_144X90_gas.OCN.v2.nc

ncatted -h -O \
    -a grid_mapping_name,crs,c,c,'latitude_longitude' \
    -a grid_mapping,mask,c,c,'crs' \
    -a _FillValue,mask,c,c,0 \
     ./dat/GLMELT_144X90_gas.OCN.v2.nc
ncdump -chs ./dat/GLMELT_144X90_gas.OCN.crs.nc

Print locations (NH; 4x5)

import xarray as xr
ds = xr.open_dataset('./dat/GLMELT_4X5.OCN.nc')

ds = ds.where((ds['mask'] == 1) & (ds['lat'] > 0))
da_stacked = ds['mask'].stack(notnull=['lat','lon'])
dd = da_stacked[da_stacked.notnull()]
# print(dd)
df = pd.DataFrame([_ for _ in zip(dd['lat'].values,dd['lon'].values,dd.values)],
                  columns=['lat','lon','mask'])

# df.set_index(['lat','lon'])
df
latlonmask
062-52.51
162-37.51
266-57.51
366-37.51
466-32.51
566-27.51
670-57.51
770-22.51
874-57.51
974-17.51

Print locations (NH; 2.5x2)

import xarray as xr
melt = xr.open_dataset('./dat/GLMELT_4X5.OCN.nc')
cl = xr.open_dataset('./dat/Z2HX2fromZ1QX1N.BS1.nc')
ds = melt.merge(cl)

ds = ds.where((ds['mask'] == 1) & (ds['lat'] > 0))

st = ds.stack(notnull=['lat','lon'])
da_stacked = ds['mask'].stack(notnull=['lat','lon'])
dd = da_stacked[da_stacked.notnull()]
# print(dd)
df = pd.DataFrame([_ for _ in zip(dd['lat'].values,dd['lon'].values,dd.values)],
                  columns=['lat','lon','mask'])

# df.set_index(['lat','lon'])
df
latlonmask
062-52.51
162-37.51
266-57.51
366-37.51
466-32.51
566-27.51
670-57.51
770-22.51
874-57.51
974-17.51

New mask creation

50, 100, 250, 500 km buffer around Greenland (raw)

  • Generate on 1/8 ° (lat,lon) grid grid
  • Can then be resampled as needed
grass -c epsg:4326 ./G
g.region n=90 s=-90 w=-180 e=180 res=0.125 -pas # 1/8 = 0.125

mkdir -p tmp
ogr2ogr -f KML tmp/M2019.kml ${DATADIR}/Mouginot_2019/Greenland_Basins_PS_v1.4.2.shp


# d.mon wx0
# d.vect M2019

v.in.ogr input=tmp/M2019.kml output=M2019
# db.select table=M2019 | head

v.db.addcolumn map=M2019 columns="REGION, INT"
# Encode using clock face numerics
db.execute sql='UPDATE M2019 SET REGION=1 where SUBREGION1 = "NE"'
db.execute sql='UPDATE M2019 SET REGION=12 where SUBREGION1 = "NO"'
db.execute sql='UPDATE M2019 SET REGION=11 where SUBREGION1 = "NW"'
db.execute sql='UPDATE M2019 SET REGION=3 where SUBREGION1 = "CE"'
db.execute sql='UPDATE M2019 SET REGION=9 where SUBREGION1 = "CW"'
db.execute sql='UPDATE M2019 SET REGION=5 where SUBREGION1 = "SE"'
db.execute sql='UPDATE M2019 SET REGION=7 where SUBREGION1 = "SW"'
v.to.rast input=M2019 output=M2019 type=area use=attr attribute_column=REGION

# d.mon wx0
# d.rast M2019

r.grow.distance -m input=M2019 value=M2019_grow distance=distance metric=geodesic

# r.mapcalc "ocean5km = if((distance > 0) & (distance < 50000), M2019_grow, null())"
r.mapcalc "ocean_gl_500 = if((distance > 0) & (distance < 500000), M2019_grow, null())"
r.mapcalc "ocean_gl_250 = if((distance > 0) & (distance < 250000), M2019_grow, null())"
r.mapcalc "ocean_gl_100 = if((distance > 0) & (distance < 100000), M2019_grow, null())"
r.mapcalc "ocean_gl_050 = if((distance > 0) & (distance < 50000), M2019_grow, null())"

# r.out.gdal format=netCDF input=ocean5km output=tmp/ocean5km.nc
r.out.gdal format=netCDF input=ocean_gl_500 output=tmp/ocean_gl_500.nc
r.out.gdal format=netCDF input=ocean_gl_250 output=tmp/ocean_gl_250.nc
r.out.gdal format=netCDF input=ocean_gl_100 output=tmp/ocean_gl_100.nc
r.out.gdal format=netCDF input=ocean_gl_050 output=tmp/ocean_gl_050.nc

50, 100, 250, 500 km around Antarctica

# grass ./G/PERMANENT

ogr2ogr -f KML tmp/aq.kml ${DATADIR}/NSIDC/NSIDC-0709.002/1992.02.07/IceBoundaries_Antarctica_v02.shp
v.in.ogr input=tmp/aq.kml output=boundary

ogr2ogr -f KML tmp/shelf.kml ${DATADIR}/NSIDC/NSIDC-0709.002/1992.02.07/IceShelf_Antarctica_v02.shp
# v.in.ogr input=tmp/shelf.kml output=shelf_all
# v.db.droprow input=shelf_all where='Name == "Filchner" OR Name == "Ross_West" OR Name == "Ross_East"' output=shelf
v.in.ogr input=tmp/shelf.kml output=shelf

v.to.rast input=shelf output=shelf type=area use=val val=1 # attr attribute_column=REGION
v.to.rast input=boundary output=boundary type=area use=val val=1

r.grow.distance -m input=shelf value=shelf_grow distance=distance_aq metric=geodesic

r.mapcalc "ocean_aq_500 = if((distance_aq > 0) & (distance_aq < 500000) & isnull(boundary) & (y() > -80), shelf_grow, null())"
r.mapcalc "ocean_aq_250 = if((distance_aq > 0) & (distance_aq < 250000) & isnull(boundary) & (y() > -80), shelf_grow, null())"
r.mapcalc "ocean_aq_100 = if((distance_aq > 0) & (distance_aq < 100000) & isnull(boundary) & (y() > -80), shelf_grow, null())"
r.mapcalc "ocean_aq_050 = if((distance_aq > 0) & (distance_aq < 50000) & isnull(boundary) & (y() > -80), shelf_grow, null())"

r.out.gdal format=netCDF input=ocean_aq_500 output=tmp/ocean_aq_500.nc
r.out.gdal format=netCDF input=ocean_aq_250 output=tmp/ocean_aq_250.nc
r.out.gdal format=netCDF input=ocean_aq_100 output=tmp/ocean_aq_100.nc
r.out.gdal format=netCDF input=ocean_aq_050 output=tmp/ocean_aq_050.nc

Realistic distribution

Antarctica

From citet:rignot_2019 Table 2, 2009 to 2017

LocMass%
W dM + dD-159-132 = -291291/466 = 0.624463519313
P dM + dD-42 - 45 = -8787/466 = 0.18669527897
E dM + dD-51 - 37 = -8888/466 = 0.188841201717
Total-291 - 87 - 88 = -46662+19+19 = 100

Experiment

LocationAmount
West70
East15
Peninsula15

Greenland

Solid and liquid discharge combined by sector

From [[id:20230127T114208][Greenland mass loss by ROI [%]​]]

regionDischargeRunoff
CE1613
CW1710
NE511
NO59
NW2112
SE3117
SW427
TOTAL9999
 
import numpy as np
import pandas as pd

dr = np.array(dr)
df = pd.DataFrame(data=dr[1:-1,1:].astype(int), columns=dr[0,1:], index=dr[1:-1,0])

df['Sum'] = df['Discharge'] + df['Runoff']
df['Sum'] = (df['Sum']/df['Sum'].sum()*100).round()
df.loc['Total'] = df.sum()
print(df)
       Discharge  Runoff    Sum
CE          16.0    13.0   15.0
CW          17.0    10.0   14.0
NE           5.0    11.0    8.0
NO           5.0     9.0    7.0
NW          21.0    12.0   17.0
SE          31.0    17.0   24.0
SW           4.0    27.0   16.0
Total       99.0    99.0  101.0

Code

# grass ./G/PERMANENT

ogr2ogr -f KML tmp/aq.kml ${DATADIR}/NSIDC/NSIDC-0709.002/1992.02.07/IceBoundaries_Antarctica_v02.shp
v.in.ogr input=tmp/aq.kml output=boundary

r.mapcalc "peninsula = if((x() > -80) & (x() < -62) & (ocean_aq_100 == 1), 1.0, 0)"
eval $(r.univar -g peninsula)
P=$(echo 15/${sum} | bc -l) # peninsula percent, spread over number of peninsula cells
r.mapcalc "peninsula = if(peninsula == 1, ${P})"

r.mapcalc "east = if((x() > 20) & (x() < 180) & (ocean_aq_100 == 1), 1.0, 0)"
eval $(r.univar -g east)
E=$(echo 15/${sum} | bc -l) # east percent, spread over number of east cells
r.mapcalc "east = if(east == 1, ${E})"

r.mapcalc "west = if((x() > -150) & (x() < -90) & (ocean_aq_100 == 1), 1.0, 0)"
eval $(r.univar -g west)
W=$(echo 70/${sum} | bc -l) # west percent, spread over number of west cells
r.mapcalc "west = if(west == 1, ${W})"



r.mapcalc "CE = if(ocean_gl_100 == 3, 1.0, 0)"
eval $(r.univar -g CE)
x=$(echo 15/${sum} | bc -l)
r.mapcalc "CE = if(CE == 1, ${x})"

r.mapcalc "CW = if(ocean_gl_100 == 9, 1.0, 0)"
eval $(r.univar -g CW)
x=$(echo 14/${sum} | bc -l)
r.mapcalc "CW = if(CW == 1, ${x})"

r.mapcalc "NE = if(ocean_gl_100 == 1, 1.0, 0)"
eval $(r.univar -g NE)
x=$(echo 8/${sum} | bc -l)
r.mapcalc "NE = if(NE == 1, ${x})"

r.mapcalc "NO = if(ocean_gl_100 == 12, 1.0, 0)"
eval $(r.univar -g NO)
x=$(echo 9/${sum} | bc -l)
r.mapcalc "NO = if(NO == 1, ${x})"

r.mapcalc "NW = if(ocean_gl_100 == 11, 1.0, 0)"
eval $(r.univar -g NW)
x=$(echo 17/${sum} | bc -l)
r.mapcalc "NW = if(NW == 1, ${x})"

r.mapcalc "SE = if(ocean_gl_100 == 5, 1.0, 0)"
eval $(r.univar -g SE)
x=$(echo 24/${sum} | bc -l)
r.mapcalc "SE = if(SE == 1, ${x})"

r.mapcalc "SW = if(ocean_gl_100 == 7, 1.0, 0)"
eval $(r.univar -g SW)
x=$(echo 16/${sum} | bc -l)
r.mapcalc "SW = if(SW == 1, ${x})"


r.mapcalc "ocean_aq_rignot = (east + west + peninsula)"
r.mapcalc "ocean_gl_mankoff = (CE +  CW + SE + SW + NE + NW + NO)"

r.null map=ocean_aq_rignot null=0
r.null map=ocean_gl_mankoff null=0
r.mapcalc "ocean_aq_gl = ocean_aq_rignot + ocean_gl_mankoff"
r.null map=ocean_aq_gl setnull=0

r.report -n ocean_aq_gl units=k
r.univar ocean_aq_gl # sum = 200

r.out.gdal format=netCDF input=ocean_aq_gl output=tmp/ocean_aq_gl_nonuniform.nc

ncrename -v Band1,mask tmp/ocean_aq_gl_nonuniform.nc tmp.nc
mv tmp.nc out/ocean_aq_gl_nonuniform.nc

Scale to 288x144

import numpy as np
import xarray as xr

ds = xr.open_dataset('./out/ocean_aq_gl_nonuniform.nc')
ds['lon'].attrs['units'] = 'degrees_east'
ds['lat'].attrs['units'] = 'degrees_north'

ds['crs'] = True
ds['crs'].attrs['grid_mapping_name'] = 'latitude_longitude'

ds['mask'].attrs['grid_mapping'] = 'crs'
ds['mask'].attrs['_FillValue'] = 0

ds.attrs = {}        
ds.attrs['Creator'] = 'Ken Mankoff'
ds.attrs['email'] = 'ken.mankoff@nasa.gov'
ds.attrs['source'] = 'https://github.com/NASA-GISS/freshwater-anomaly'

latbins = np.arange(-90,90+1)
lonbins = np.arange(-180+0.625, 180+0.625*2, 1.25)
ds2 = ds.groupby_bins('lat', latbins, labels=latbins[:-1]+0.5)\
        .max()\
        .groupby_bins('lon', lonbins, labels=lonbins[:-1])\
        .max()\
        .rename({'lat_bins':'lat','lon_bins':'lon'})\
        .transpose() # Model expects variables to be (lat,lon)
# ds2 = ds2.rename({'ones':'mask'})
ds2['lon'].attrs['units'] = 'degrees_east'
ds2['lat'].attrs['units'] = 'degrees_north'
ds2['mask'].attrs['_FillValue'] = 0
ds2 = ds2.fillna(0)
ds2.to_netcdf('./out/fw_nonuniform_288x180.nc', format='NETCDF3_CLASSIC')

# Trickier: non-uniform lat bins
# ds['ones'].coarsen(lat=int(4/(1/8))).mean().coarsen(lon=int(5/(1/8))).mean().to_netcdf('./dat/fw_4x5.nc')
print(ds2)

Build mask with scaled distribution

  • See slides above for [[id:20230127T114208][Greenland mass loss by ROI [%]​]]

WARNING: This spreads the distribution per sector over the number of grid cells. Because I’m working on my own map of Greenland at high res, some of these grid cells will be land-cells in ModelE. The relative number of land cells per sector should be similar between sectors (large sectors means long coastline meaning many land cells, vs. small sector short coast few cells), meaning this error may not be significant to final results.

import numpy as np
import xarray as xr

ds = xr.open_dataset('./tmp/ocean_gl_050.nc')
ds['lon'].attrs['units'] = 'degrees_east'
ds['lat'].attrs['units'] = 'degrees_north'

# record the RIO for each cell
ROI = np.empty(ds['Band1'].data.shape, dtype='U2')
b1 = ds['Band1'].values
ROI[b1 == 12] = 'NO' # cloc face
ROI[b1 == 1] = 'NE'
ROI[b1 == 3] = 'CE'
ROI[b1 == 5] = 'SE'
ROI[b1 == 7] = 'SW'
ROI[b1 == 9] = 'CW'
ROI[b1 == 11] = 'NW'
ROI[b1 == b1[0,0]] = ''
ds['ROI'] = (('lat','lon'), ROI)
ds['ROI'].attrs['ROI_source'] = 'Mouginot /et al./ (2019); https://doi.org/10.7280/D1WT11'
ds = ds.drop_vars('Band1')

for dist in ['050', '100', '250', '500']:
    gl = xr.open_dataset('./tmp/ocean_gl_'+dist+'.nc')
    aq = xr.open_dataset('./tmp/ocean_aq_'+dist+'.nc')
    ds['ones_'+dist] = (gl['Band1'].notnull()+aq['Band1'].notnull()).astype(np.int32)
    ds['ones_'+dist].attrs['description'] = 'Grid cells w/in '+dist+' km of coast with freshwater runoff, submarine melt, or iceberg forcing'

# ds['base'] = ds['ones']
# ds['base'].attrs['description'] = 'Baseline freshwater runoff, submarine melt, or iceberg forcing, scaled by sector contribution'

# ds['base_solid'] = ds['ones']
# ds['base_solid'].attrs['description'] = 'Baseline iceberg forcing, scaled by sector contribution'

# ds['base_liquid'] = ds['ones']
# ds['base_liquid'].attrs['description'] = 'Baseline freshwater runoff and submarine melt scaled by sector contribution'

# ds['anom'] = ds['ones'].astype(np.float32)
# ds['anom'].attrs['description'] = 'Anomaly freshwater runoff, submarine melt, or iceberg forcing, scaled by sector contribution'
# contrib = {'CE':1, 'CW':24, 'NO':12, 'NE':17, 'NW':29, 'SE':5, 'SW':11}
# ds['anom'].attrs['distribution_pct'] = str(contrib)
# # ds['anom'] = ds['anom'].where(ds['ROI'] != 'CE', other=(1  / ds['ones'].where(ds['ROI'] == 'CE').sum().data))
# # ds['anom'] = ds['anom'].where(ds['ROI'] != 'CW', other=(24 / ds['ones'].where(ds['ROI'] == 'CW').sum().data))
# # ds['anom'] = ds['anom'].where(ds['ROI'] != 'NO', other=(12 / ds['ones'].where(ds['ROI'] == 'NO').sum().data))
# # ds['anom'] = ds['anom'].where(ds['ROI'] != 'NE', other=(17 / ds['ones'].where(ds['ROI'] == 'NE').sum().data))
# # ds['anom'] = ds['anom'].where(ds['ROI'] != 'NW', other=(29 / ds['ones'].where(ds['ROI'] == 'NW').sum().data))
# # ds['anom'] = ds['anom'].where(ds['ROI'] != 'SE', other=(5  / ds['ones'].where(ds['ROI'] == 'SE').sum().data))
# # ds['anom'] = ds['anom'].where(ds['ROI'] != 'SW', other=(11 / ds['ones'].where(ds['ROI'] == 'SW').sum().data))
# for k,v in contrib.items():
#     ds['anom'] = ds['anom'].where(ds['ROI'] != k, other=(v  / ds['ones'].where(ds['ROI'] == k).sum().data))

# ds['anom_solid'] = ds['ones'].astype(np.float32)
# ds['anom_solid'].attrs['description'] = 'Anomaly iceberg forcing, scaled by sector contribution'
# contrib = {'CE':16, 'CW':17, 'NO':5, 'NE':5, 'NW':21, 'SE':31, 'SW':4}
# ds['anom_solid'].attrs['distribution_pct'] = str(contrib)
# for k,v in contrib.items():
#     ds['anom_solid'] = ds['anom_solid'].where(ds['ROI'] != k, other=(v  / ds['ones'].where(ds['ROI'] == k).sum().data))

# ds['anom_liquid'] = ds['ones'].astype(np.float32)
# ds['anom_liquid'].attrs['description'] = 'Anomaly freshwater runoff and submarine melt scaled by sector contribution'
# contrib = {'CE':13, 'CW':10, 'NO':11, 'NE':9, 'NW':12, 'SE':17, 'SW':27}
# ds['anom_liquid'].attrs['distribution_pct'] = str(contrib)
# for k,v in contrib.items():
#     ds['anom_liquid'] = ds['anom_liquid'].where(ds['ROI'] != k, other=(v  / ds['ones'].where(ds['ROI'] == k).sum().data))

ds['crs'] = True
ds['crs'].attrs['grid_mapping_name'] = 'latitude_longitude'

for v in ds.data_vars:
    ds[v].attrs['grid_mapping'] = 'crs'
    if (v != 'crs') and (v != 'ROI'):
        ds[v].attrs['_FillValue'] = 0

ds.attrs = {}        
ds.attrs['Creator'] = 'Ken Mankoff'
ds.attrs['email'] = 'ken.mankoff@nasa.gov'
ds.attrs['source'] = 'https://github.com/NASA-GISS/freshwater-anomaly'

ds.to_netcdf('./out/fw.nc')

# current resolution: 1/8

latbins = np.arange(-90,90+1)
lonbins = np.arange(-180+0.625, 180+0.625*2, 1.25)
ds2 = ds.groupby_bins('lat', latbins, labels=latbins[:-1]+0.5)\
        .max()\
        .groupby_bins('lon', lonbins, labels=lonbins[:-1])\
        .max()\
        .rename({'lat_bins':'lat','lon_bins':'lon'})\
        .transpose() # Model expects variables to be (lat,lon)
# ds2 = ds2.rename({'ones':'mask'})
ds2['lon'].attrs['units'] = 'degrees_east'
ds2['lat'].attrs['units'] = 'degrees_north'
ds2.to_netcdf('./out/fw_288x180.nc', format='NETCDF3_CLASSIC')

# Trickier: non-uniform lat bins
# ds['ones'].coarsen(lat=int(4/(1/8))).mean().coarsen(lon=int(5/(1/8))).mean().to_netcdf('./dat/fw_4x5.nc')
print(ds2)

Variable to mask

cdo chname,ones_100,mask fw.nc fw_100.nc
ncrename -v ones_100,mask fw.nc fw_100.nc

NOT USED

ds = xr.Dataset()

step = 1/60
ds['lon'] = np.arange(-180 + (step/2),180, step, dtype=np.float64)
ds['lon'].attrs['units'] = 'degrees_north'

ds['lat'] = np.arange(-90 + (step/2), 90, step, dtype=np.float64)
ds['lat'].attrs['units'] = 'degrees_east'

ds['crs'] = True
ds['crs'].attrs['grid_mapping_name'] = 'latitude_longitude'

ds['ROI'] = (('lon','lat'), np.zeros((ds['lon'].size, ds['lat'].size), dtype='U2'))
# ds['ROI'].attrs['_FillValue'] = ''
ds['ROI'].attrs['grid_mapping'] = 'crs'
ds['ROI'].attrs['ROI_source'] = 'Mouginot /et al./ (2019); https://doi.org/10.7280/D1WT11'

## Example: Can set fields with:
# ds['ROI'][:,np.where(ds['lat'] == 80)[0][0]:np.where(ds['lat'] == 90)[0][0]] = 'NO'
## Let's create a small function to help with that.

def lati(val): return np.where(ds['lat'] == val)[0][0] # lat index
def loni(val): return np.where(ds['lon'] == val)[0][0]
# use like: ds['ROI'][:,lati(80):lati(90)+1] = 'NO'

ds['ROI'][loni(ds.isel(27,38] = 'SW'
# ds['ROI'][ds.isel(27,38] = 'SW'
# ds['region'][26,39] = 'SW'
# ds['region'][25,40] = 'SW'

# ds['region'][27,38] = 'SE'
# ds['region'][27,38] = 'SE'
# ds['region'][27,38] = 'SE'


ds.attrs['Creator'] = 'Ken Mankoff'
ds.attrs['email'] = 'ken.mankoff@nasa.gov'
ds.to_netcdf('icesheet_runoff_mask.nc')
print(ds)

Greenland mask

  • ones : All 1 within 100 km of Greenland coast
  • anom : MB anomaly by ROI
  • anom_liquid SMB anomaly by ROI
  • anom_solid Discharge anomaly by ROI

Antarctic mask ???

Even distribution w/in 500 km of coast

  • Solid discharge at ice shelf fronts
    • Scaled by ice shelf area??
  • Ice shelf melting at ice shelf fronts
    • Scaled by ice shelf area??

{{{SKIPLINE}}}

Doesn’t reflect reality - most of the mass loss from small shelves, not Ross and Filchner-Ronne (mass gain).

  • Exclude Ross and Filchner-Ronne?

{{{SKIPLINE}}}

citet:NSIDC_0709,rignot_2013 See https://nsidc.org/data/nsidc-0709

New mask contents

./fig/panoply_fw.nc.png

New mask visualization (Greenland)

./.ob-jupyter/04adf306af466580479853e74b920e0d45dbc623.png

New mask visualization (Antarctica)

./.ob-jupyter/f6294105a78024b0fca36d4b9a3369064688cc88.png

Summary

Questions & To Do

  • [ ] Updated baseline mask
    • [ ] Include N. Greenland
    • [ ] Treat spatial variability
    • [ ] Separate (liquid runoff & submarine melt) v. solid discharge
  • [ ] Anomaly: determine regional distribution of
    • [ ] Liquid runoff & submarine melt
    • [ ] Solid ice discharge
  • [ ] Build anomaly distribution masks for each of these

Appendix

References

\printbibliography[heading=none]

About This Document

This document is an Emacs Org Mode plain-text file with code and text embedded. If you are viewing:

  • A PDF, HTML, or DOC file, then it was generated by exporting from Org. Not all of the Org parts (code, results, comments, etc.) were exported. The Org source file is available upon request, and may be embedded in the PDF. You can access files embedded in PDF files with from within your PDF viewer.
  • A file with an org extension in something other than Emacs, then you are seeing the canonical version and the full source, but without any syntax highlighting, document structure, or the ability to execute the code blocks.
  • An Org file within Emacs, then this is the canonical version. You should be able to fully interact and reproduce the contents of this document, although it may require 3rd-party applications (Python, etc.) and a similar Emacs configuration. This is available upon request.

LaTeX Header

Beamer

Title Page

Theme

Other

References

Hyperref

Tweak References

Background Block

Code export

Local Variables