#!/usr/bin/env python
# coding: utf-8

# In[1]:


import pandas as pd
import numpy as np


# ## Read the data
# First, we read in a data file containing an estimate of Earth's contemporary glacier volume, separated in above and below sea-level:

# In[2]:


data = pd.read_hdf('rgi62_era5_itmix_df.h5', 'df')


# In[3]:


pd.options.display.float_format = '{:,.1f}'.format


# ## Farinotti et al. 2019
# In the article "A consensus estimate for the ice thickness distribution of all glaciers on Earth" by [Farinotti et al. (2019)](https://www.nature.com/articles/s41561-019-0300-3) the potential sea-level rise ($h_{sl}$) by melting all the ice ($V_{tot}$) instantaneously was estimated as follows:

# $$ h_{sl} = \frac{V_{tot} - V_{bsl}}{A_{ocean}} \frac{\rho_{ice}}{\rho_{ocean}}, \text{with} $$
# 
# $\rho_{ice}$ : density of ice
# 
# $\rho_{ocean}$ : density of ocean water
# 
# $V_{bsl}$: grounded volume below sea-level
# 
# $A_{ocean}$: area of Earth's oceans

# Let us calculate the total sea-level rise this would give, and the potential sea-level rise, if we neglected the fact that part of the volume is below sea-level, which has to be replaced by meltwater before sea-level will actually rise:

# In[4]:


def compute_slr(ice_vol_m3):
    """After Farinotti 2019"""
    rho = 900
    rho_oc = 1028
    A_oc = 362.5 * 1e9
    return ice_vol_m3 * rho / (A_oc * rho_oc)


v_asl = data.vol_itmix_m3 - data.vol_bsl_itmix_m3

total_slr = compute_slr(v_asl.sum())
total_slr_bsl = compute_slr(data.vol_bsl_itmix_m3.sum())
total_slr_wo_bsl = compute_slr(data.vol_itmix_m3.sum())

odf = pd.DataFrame(index=['F2019'])
odf['SLR (incl. V$_{bsl}$)'] = total_slr
odf['SLR (excl. V$_{bsl}$)'] = total_slr_wo_bsl
odf['% difference'] = 100*data.vol_bsl_itmix_m3.sum()/data.vol_itmix_m3.sum()
odf.T


# It can be seen from this exercise that roughly 15 % of the global ice volume is situated below sea-level and neglecting that fact would therefore lead to a considerable overestimation of potential sea-level rise of about 58 mm!

# We can do the same exercise for each region:

# In[5]:


def compute_slr_region(region):
    """slr2 regional"""
    vol_region = np.nansum(data.where(data["O1Region"]==region).iloc[:,-1])
    vol_region_bsl = np.nansum(data.where(data["O1Region"]==region).iloc[:,-2])
    slr = compute_slr(vol_region- vol_region_bsl)
    slr_bsl = compute_slr(vol_region_bsl)
    slr_wo_bsl = compute_slr(vol_region)
    return slr, slr_bsl, slr_wo_bsl

odf_reg = pd.DataFrame()
for region in data.O1Region.unique():
    slr_region = compute_slr_region(region)
    
    odf_reg.loc[region, 'SLR (incl. V$_{bsl}$)'] = slr_region[0]
    odf_reg.loc[region, 'SLR (excl. V$_{bsl}$)'] = slr_region[2]
    odf_reg.loc[region, '% difference'] = 100*slr_region[1]/slr_region[2]
    
odf_reg


# ## Potential issues with the SLR equation in F2019 

# Now we have an estimate of the glacier ice volume (in mm SLE) above and below sea-level. But concerning the estimation of potential sea-level rise melting that ice would cause, the equation above has two potential flaws:
# 
# 1. It removes the grounded volume below sea-level ($V_{bsl}$) in terms of the resulting ocean(!) water volume, which is $V_{bsl} \frac{\rho_{ice}}{\rho_{ocean}}$ , while in reality it should be removed as $V_{bsl} \frac{\rho_{freshwater}}{\rho_{ice}}$, because that is the volume of melted ice needed to replace the ocean water that was displaced by the ice below sea-level. The rest ($V_{tot} - V_{bsl} \frac{\rho_{freshwater}}{\rho_{ice}}$) will add volume as freshwater and thus has to be multiplied with $\frac{\rho_{ice}}{\rho_{freshwater}}$, resulting in the equation below. This is also equivalent to substracting $V_{bsl}$ from the total meltwater volume ($V_{tot} \frac{\rho_{ice}}{\rho_{freshwater}}$), since $V_{bsl}$ needs to be replaced by meltwater before the melt volume can actually increase the sea-level.
# 
# 
# 2. It divides by the density of ocean water, although melted ice will contain (nearly) no salt and therefore have a lower density, i.e. higher volume, than ocean water.

# ## A different approach
# With the reasoning above, we get following equation, which calculates the total sea-level height increase when all the ice volume is melted instantaneously, neglecting calving and submarine melt as well as the halosteric and thermosteric effects of melt occuring in the oceans:

# $$ h_{sl} = \frac{V_{tot} \cdot f_1 - V_{bsl}}{A_{ocean}}$$<br>
# 
# where $$f_1 = \frac{\rho_{ice}}{\rho_{freshwater}} $$

# Now we calculate the regional potential sea-level rise, and the difference to the Farinotti et al. (2019) equation as well as the percentage difference:

# In[6]:


def compute_slr2(v_tot, v_bsl):
    """slr2"""
    rho = 900
    rho_w = 1000
    A_oc = 362.5 * 1e9
    return (v_tot * (rho / rho_w) - v_bsl) / A_oc

total_slr2 = compute_slr2(data.vol_itmix_m3.sum(), data.vol_bsl_itmix_m3.sum())

def compute_slr_region2(region):
    """slr2 regional"""
    vol_region = np.nansum(data.where(data["O1Region"]==region).iloc[:,-1])
    vol_region_bsl = np.nansum(data.where(data["O1Region"]==region).iloc[:,-2])
    slr = compute_slr2(vol_region, vol_region_bsl)
    slr_bsl = compute_slr2(vol_region_bsl, vol_region_bsl)
    slr_wo_bsl = compute_slr2(vol_region, 0)
    return slr, slr_bsl, slr_wo_bsl

odf_reg_comp = odf_reg[[odf_reg.columns[0]]].copy()
odf_reg_comp.columns = ['SLR (F2019)']

for region in data.O1Region.unique():
    slr_region2 = compute_slr_region2(region)
    slr_region = compute_slr_region(region)
    
    odf_reg_comp.loc[region, 'SLR (New)'] = slr_region2[0]
    
total = odf_reg_comp.sum()
total.name = 'Global'
odf_reg_comp = odf_reg_comp.append(total)

odf_reg_comp['Diff'] = odf_reg_comp['SLR (New)'] - odf_reg_comp['SLR (F2019)']
odf_reg_comp['Diff (%)'] = 100 * odf_reg_comp['Diff'] / odf_reg_comp['SLR (F2019)']

odf_reg_comp


# It becomes evident that the global difference in estimated potential sea-level rise by melting all of Earth's glacier ice instantaneously is not very large (~0.8 %), but that there are considerable differences regionally. 
# 
# In regions without glacier ice volume below sea-level the estimate with the second equation is 2.8 % larger, because it divides by freshwater density, which is ~2.8 % lower than that of ocean water. In regions with much glacier ice below sea-level this is compensated for by the fact that the second equation removes the grounded volume below sea-level in terms of the meltwater needed to replace it, which is more than the ice that is situated below sea-level itself. This results in the second equation estimating less potential sea-level rise from the Antarctic periphery (region 19) where ~38 % of the glacier ice is located below sea-level.

# ## An expanded equation for sea-level rise due to glacier mass change
# 
# As already indicated in the blog post, there are more processes at play; more than we can cover in a blog post. But in this formula, we try to take into account calving, and halosteric/thermosteric effects of oceanic ice melt on global mean sea-level rise. If all the parameters were known, we could calculate the direct and delayed sea-level rise from glacier ice melt as follows:

# $$ h_{sl} = \frac{V_{tot} \cdot (f_c \cdot f_2 + (1 - f_c - f_{sm}) \cdot f_1 + f_{sm} \cdot (f_1 - f_T) + f_c \cdot f_2 \cdot f_h - f_c \cdot f_T) - V_{bsl}} {A_{ocean}}, \text{with} $$
# 
# $f_c$ : fraction of ice volume lost by calving
# 
# $f_2 = \frac{\rho_{ice}}{\rho_{saltwater}}$
# 
# $f_{sm}$: fraction of ice volume lost by oceanic/submarine melt
# 
# $f_h = \frac{\rho_{saltwater} - \rho_{freshwater}}{\rho_{freshwater}}$
# 
# $f_T$: fraction of volume contraction caused in the ocean due to
# cooling by the extraction of latent heat required for melting the ice volume<br>
# 
# In this equation the last two terms in the large parentheses refer to the
# delayed influence of glacier mass loss on sea-level rise induced by the
# melting of floating ice (i.e. icebergs). This can be positive, zero, or negative,
# depending on halosteric ($f_h$) vs. thermosteric ($f_T$) effect.
# 
# One parameter in this equation that is not easily calculated is the one of the thermosteric process. Here, we just set it to be slightly lower than the halosteric effect (~2 %). We take the fraction of calved and submarine melt volume to be 10 % and 5 %, respectively, and calculate the total sea-level rise this would lead to, the difference to the Farinotti et al. (2019) equation, the difference to the previous equation, the percentage difference to the Farinotti et al. (2019) equation, and the sea-level change that would occur delayed.

# In[7]:


def compute_slr3(v_tot, v_bsl):
    """slr3"""
    rho = 900
    rho_w = 1000
    rho_s = 1028
    f2 = rho/rho_s
    f1 = rho/rho_w
    A_oc = 362.5 * 1e9
    f_c = 0.10
    f_sm = 0.05
    f_t = 0.02
    f_h = (rho_s - rho_w) / rho_w
    total = (v_tot * (f_c * f2 + f1 * (1-f_c-f_sm) + f_sm * (f1 - f_t) + f_c * f2 * f_h - f_c * f_t) - v_bsl) / A_oc
    immediate = (v_tot * (f_c * f2 + f1 * (1-f_c-f_sm) + f_sm * (f1 - f_t)) - v_bsl) / A_oc
    delayed = (v_tot * (f_c * f2 * f_h - f_c * f_t)) / A_oc
    return total, immediate, delayed

total_slr3 = compute_slr3(data.vol_itmix_m3.sum(), data.vol_bsl_itmix_m3.sum())

odf = pd.DataFrame(index=['SLR'])
odf['Total'] = total_slr3[0]
odf['Immediate'] = total_slr3[1]
odf['Delayed'] = total_slr3[2]
odf['Diff. to F19'] = total_slr3[0] - total_slr
odf['Diff. to Eq. 2'] = total_slr3[0] - total_slr2
odf.T


# We see that this results in a slightly lower sea-level rise estimate than with the previous equation due to the thermosteric effect of submarine (and iceberg) melt. The delayed sea-level rise is very small, partly because the halosteric increase of melting an iceberg only affects the volume of the iceberg that is below sea-level, while the thermosteric effect affects the whole iceberg (neglecting atmospheric melt of the iceberg's part above sea-level). You can play around with the values of the parameters and look how that changes the results; e.g. put f_c = 1 and f_sm = 0 to see what would change if all the Earth's glacier ice would somehow calv at once.

# ## A simpler equation for (delayed) sea-level rise due to glacier mass change
# 
# To make the previous equation a bit slimmer, we neglect submarine melt and the thermosteric ($f_T$)
# effect. The equation then becomes:
# 
# $$ h_{sl} = \frac{V_{tot} \cdot (f_c \cdot f_2 + (1 - f_c) \cdot f_1 + f_c \cdot f_2 \cdot f_h) - V_{bsl}} {A_{ocean}} $$
# 
# This equation accounts for part of the ice mass being added to the oceans in the
# form of icebergs. The last term in the large parentheses is the delayed increase
# in ocean volume due to melting those via solely the halosteric effect. In this case we take the calved volume to be 15 %.

# In[9]:


def compute_slr4(v_tot, v_bsl):
    """slr4"""
    rho = 900
    rho_w = 1000
    rho_s = 1028
    f2 = rho/rho_s
    f1 = rho/rho_w
    A_oc = 362.5 * 1e9
    f_c = 0.15
    f_h = (rho_s - rho_w) / rho_w
    total = (v_tot * (f_c * f2 + f1 * (1-f_c) + f_c * f2 * f_h) - v_bsl) / A_oc
    immediate = (v_tot * (f_c * f2 + f1 * (1-f_c)) - v_bsl) / A_oc
    delayed = (v_tot * f_c * f2 * f_h) / A_oc
    return total, immediate, delayed

total_slr4 = compute_slr4(data.vol_itmix_m3.sum(), data.vol_bsl_itmix_m3.sum())

odf = pd.DataFrame(index=['SLR'])
odf['Total'] = total_slr4[0]
odf['Immediate'] = total_slr4[1]
odf['Delayed'] = total_slr4[2]
odf['Diff. to F19'] = total_slr4[0] - total_slr
odf['Diff. to Eq. 2'] = total_slr4[0] - total_slr2
odf.T


# As we do not take the extraction of ocean heat for ice melt into account with this formula, there is no difference in total sea-level rise compared to the second equation in this Notebook, but we get a simple estimate of delayed sea-level rise due to the halosteric effect.
# 
# ## Ultimately, we see that the differences to Farinotti et al. (2019) are below 1 %, but up to 4 % regionally. They become even smaller when we account for thermosteric effects. But I hope it was fun to think about this topic and play around with some equations!