Cumulative Elevation Calculation

[Jez-Carter]

I've noticed some odd behaviour with cumulative elevation. That is cumulative elevation during the reference spin-up phase does not flatten (e.g. I would expect eventually cumulative elevation would stop changing as the state reaches an equilibrium). To show this behaviour I've created the following example script, where I've used the example csvs for the Summit station and have extended the spinup time even further to emphasise the result. Note, since the average snowfall is 0.23 mIeq/yr it's expected the number of years required to completely refresh the column is 120/0.23=522 (120m is the domain depth).

import glob
import numpy as np
import pandas as pd
import json
import subprocess
### Extending SpinUP
for f in glob.glob('/CommunityFirnModel/CFM_main/CFMinput_example/example*'):
    df = pd.read_csv(f).T
    df_spin = df[df.index.astype('float32') < 1980]
    df_extended = pd.concat([df_spin,df])
    old_index = df_extended.index.astype('float32')
    step_size = old_index[1]-old_index[0]
    new_index = np.arange(0,len(old_index))*step_size + old_index.max() - len(old_index)*step_size
    df_extended = df_extended.set_index(new_index)
    new_filename = f.split('/')[-1].replace('example','extended')
    df_extended.T.to_csv(f'/CommunityFirnModel/CFM_main/CFMinput_example/{new_filename}',index=False, header=True)

### Loading and editing JSON file for CFM run
with open('/CommunityFirnModel/CFM_main/example.json', 'r') as read_file:
    adjusted_json = json.load(read_file)
adjusted_json['InputFileNameTemp'] = 'extended_TSKIN.csv'
adjusted_json['InputFileNamebdot'] = 'extended_BDOT.csv'
adjusted_json['InputFileNamemelt'] = 'extended_SMELT.csv'
adjusted_json['InputFileNameRain'] = 'extended_RAIN.csv'    
adjusted_json['input_type'] = 'csv'
adjusted_json['TWriteStart'] = 980.0
adjusted_json['residual_strain'] = 0.0
adjusted_json['doublegrid'] = False

with open('/CommunityFirnModel/CFM_main/adjusted.json', 'w') as json_file:
    json.dump(adjusted_json, json_file, indent=0)

### Running CFM with edited JSON
subprocess.call(['python','/CommunityFirnModel/CFM_main/main.py','/CommunityFirnModel/CFM_main/adjusted.json','-n'])

The result of running the above script is the following, where it can be seen no steady state is reached.

import xarray as xr
ds = xr.open_dataset('/CommunityFirnModel/CFM_main/CFMoutput_example/df/CFMresults.hdf5')
dimensions_dict = {ds.density.dims[0]:'time',ds.density.dims[1]:'cell'}
ds = ds.rename_dims(dimensions_dict)
year_data = ds.isel(cell=0).density.data
ds = ds.assign_coords(year = ('time',year_data))
ds = ds.isel(cell=slice(1,None))
ds['cumulative_elevation']=ds.DIP[:,6]
ds['cumulative_elevation'].plot.line(x='year')

The equation for elevation change is given in firn_density_no_spin.py as:

self.dHcorr = (self.sdz_new - self.sdz_old) + self.dh_acc + self.dh_melt - (iceout_corr*self.t[iii]) # iceout has units m ice/year, t is years per time step.

This doesn't make complete sense to me as it feels like melt is being counted twice since (a similar issue as mentioned in issue #9), since self.sdz_new is labelled as including the impact of melt:

self.sdz_new = np.sum(self.dz) #total column thickness after densification, melt, horizontal strain, before new snow added

I tried adjusting the code to leave out self.dh_melt from the above:

#self.dHcorr = (self.sdz_new - self.sdz_old) + self.dh_acc + self.dh_melt - (iceout_corr*self.t[iii]) # iceout has units m ice/year, t is years per time step. 
self.dHcorr = (self.sdz_new - self.sdz_old) + self.dh_acc - (iceout_corr*self.t[iii]) # iceout has units m ice/year, t is years per time step.

Doing this I get the following result, where minimal difference is made since Summit receives minimal melt:

I've also tried replacing the above with a commented out line in the original code, which seems like it should work:

#self.dHcorr = (self.sdz_new - self.sdz_old) + self.dh_acc + self.dh_melt - (iceout_corr*self.t[iii]) # iceout has units m ice/year, t is years per time step. 
self.dHcorr = self.z[-1] - self.z_old[-1]

This does seem to result in a steady state reached at around 500 years:

Although, this doesn't seem to work for high-melt sites, which I've simulated using the following example adjustment to the input data:

### Extending SpinUP
for f in glob.glob('/CommunityFirnModel/CFM_main/CFMinput_example/example*'):
    df = pd.read_csv(f).T
    df_spin = df[df.index.astype('float32') < 1980]
    ###############
    if 'SMELT' in f:
        df_spin = df_spin+2
        df = df+5
    elif 'BDOT' in f:
        df_spin = df_spin+1
        df = df+1
    ###############
    df_extended = pd.concat([df_spin,df])
    old_index = df_extended.index.astype('float32')
    step_size = old_index[1]-old_index[0]
    new_index = np.arange(0,len(old_index))*step_size + old_index.max() - len(old_index)*step_size
    df_extended = df_extended.set_index(new_index)
    new_filename = f.split('/')[-1].replace('example','extended')
    df_extended.T.to_csv(f'/CommunityFirnModel/CFM_main/CFMinput_example/{new_filename}',index=False, header=True)

The above adjustment produces the following, which while it reaches an equilibrium shows no step change at the year 1452 where melt goes from approximately 2 to 5 m iEq/yr.

Am I missing something in all the above or is there an issue with the calculation of cumulative elevation? Best, Jez

[maximusjstevens]

Hi @Jez-Carter, thanks for digging into this. Sorry, I have not had time in the last month to dig into your previous issue #9. I think you are not missing anything; that dh code was actually developed prior to the CFM getting melt physics (i.e., it is pretty old!). If you want, you can dive into the code and try to fix the issue and make a pull request, or hopefully I can find a moment later this week or next week to take a look.

[Jez-Carter]

Hi @maximusjstevens, thanks for offering to take a look, I think that would be best and to get your further thoughts. I'll also get in contact with Vincent as I believe he's worked with cumulative elevation before. I won't be able to spend much time digging into the code myself over the next couple of weeks as I have another deadline approaching for a different piece of work but will keep you updated if I am able to.

[vverjans]

Hello everyone, Jezz, I hard-coded my own variable for computing elevation changes. The CFM structure has changed since then, so I cannot simply push my (very) old code. However, I am happy to look into that in the coming days. I can probably work on this aspect in the coming days (or early next week). Jezz, I will email you to get all your input files.

[vverjans]

I think that I fixed the problem. See the email I just sent to Max and you.