Description:

For the periodic simulations, the air temperature oscillates around the true value during surging (see figure). The oscillation is in either direction (above/below reference) and is nearly symmetric in time. This is obviously nonphysical as the air temperature / elevation relationship shouldn't be affected by surface speed or the timestep.

Here's an alternate view of the problem. This is the surface air temperature at the ELA, which is fixed elevation so it should have a constant air temperature. You can see the asymmetric but equal amplitude oscillations in either direction, which occur at the same interval as the forced surges (30 years). The 1 degree amplitude is huge in the context of our modeling, and too large to really ignore.

Diagnosis:

I'm fairly confident this problem has to do with variable timesteps used for the periodic simulation. To confirm I've run a test case with the a constant during quiescence and surging, for both the dynamic and thermodynamic solvers, of dt=0.05. Oscillation go away, though there is some noise in the mean value it's relatively minor and acceptable for where we are in the timeline.

Proposed solution:

My initial inclination is that this stems from a problem in converting time into day of year to evaluate the air temperature and melt functions. In order to fix this, I will re-parameterize all the code to be a function of time in years, instead of day of year.

Implementing the Vectorized (and threaded) Stokes Solver (Stokes-Vec) would result in significant simulation speed ups. There have be numerous reasons barring us from its implementation, some resolved some still unresolved.

Resolved Reasons:

• .result file parsing:
  • An early attempt use Stokes-Vec was abandoned because it wrote it's solution vectors to .result files differently, than the legacy Stokes solver. elmer2nc was not built out to handle these differences, so Stokes-Vec was abandoned.
  • With the adoption of NetcdfUGRIDOutputSolver Solver, this no longer an obstacle. NetCDF files with solution vectors are directly written by the .sif file so don't need to parse .result files anymore.

Unresolved Reasons:

• Structured Mesh Mapper:
  • It appears as though the Stokes-Vec solver does not work well with the Structured Mesh Mapper.

In the surface boundary process function (Surface_Processes in src/elmer_UDF/SurfaceBoundary.f90) I calculate the volumetric heat source (Q_lat) based on Eqn. (9) from Wilson and Flowers (2013):

My initial attempts to prescribe Q_lat as prescribing a volumetric source (J/m^3) along the surface boundary conflicted with the Dirichlet B.C. based on air temp and in-effect Q_lat was ignored by Elmer. Previously as solution this problem, I prescribed Q_lat as a volumetric heat source one node below the surface, which worked but was nonphysical.

Instead I've edited the code to convert Q_lat to a mass heat source (J/kg) which can be used to augment the Dirichlet B.C. for enthalpy along the free surface. See Eqn. (6) from Licciulli et al. (2019) for an example of this approach.

I'm unsure what density should I use to convert Q_lat from a volumetric (J/m^3) to mass (J/kg) heat source. I'm currently just using the density of water:

Per communication with Adrien Gilbert, the enthalpy solver in Elmer can only achieve steady-state convergence with a tolerance of 1e-3 or 1e-4. Adrien has suggested the lack of convergence is due to the discontinuity in the enthalpy diffusivity between cold and temperate ice.

The discontinuous enthalpy diffusivity has caused problems in other ice-sheet model, ISSM dealt with the problem by implementing a solver specific stabilization method (Ruckamp et al, 2020). This is out of the scope of what we can do in Elmer on our timescale.

A simpler approach might be fitting a smooth function between the cold and temperate diffusivity as a form of regularization (see Wilson and Flowers, 2013).

What we want to do:

Periodic surging experiments would greatly benefit from more efficient timesteping. Using the minimum timestep, which is typically the Stokes solver timestep during surging, for the entire 2000+ year simulation is computational prohibitive.

During the quiescence the Stokes solver probably only needs dt=1.0, whereas it's probably run with dt=0.05 during surging. Furthermore, the Enthalpy solver probably only needs to be run with dt=0.1 during quiescence, whereas it's also run with dt=0.1 during surging.

Because of the regular return intervals, this should be very doable in Elmer/Ice. Something like:

Simuluation
.
.
.
  Timestep Intervals(6)          = Integer 100 450 100 450 100 450 ![a]
  Timestep Sizes(6)              = Real    0.05 0.1 0.05 0.1 0.05 0.1  ![a]
.
.
.
End 

...

Solver X
  Exec Interval(6)            = Integer 1 10 1 10 1 10
  Timestep Scale(6)           = Real    1 10 1 10 1 10
 .
 . 
  .
End Solver

where Solver X is a dynamics solver which only needs to be run w/ dt=1.0 during quiescence. Thermodynamic solver would be executed with dt=0.1. This simulation would be 3-surge cycle with a 5-year active phase and a 45-year quiescent phase.

The Problem:

The Timestep Scale keyword from the solver section only excepts constants. See the "elmer/fem/src/MainUtils.f90" file:

     IF(.NOT. ListGetLogical( Params,'Auxiliary Solver',Found)) THEN
       DTScal = ListGetConstReal( Params, 'Timestep Scale', Found )
       IF ( .NOT. Found ) DTScal = 1.0_dp

This a problem, given that Exec Intervals needs an array the same length as Timestep Sizes and Timetep Intervals. Therefore, we can execute solvers at different intervals but the timestep will not be correctly scaled.

How to fix:

This is a problem with the "elmer/fem/src/MainUtils.f90" file which is beyond the scope of what I'm comfortable editing

I've posted on the Elmer forum requesting help to see if anyone can update the code to fix this. Elmer Forum post can be found here:

http://elmerfem.org/forum/viewtopic.php?t=7836

How do we determine where the maximum water content should be the value for glacier ice (~3%) versus firn (~10%)?

For now I've set the maximum water content based on density:

Gwenn made a good point that maybe firn values of maximum water content should be confined to the "firn aquifer". The firn aquifer thickness is a free parameter which we need to set and have been using a reference value of ~3 m for experiments . So in practice, because of our vertical grid cell spacing, this means only the surface nodes would have firn water content values.

My thought process was given the higher porosity of firn we should allow for more water content. That being said the whole idea of setting firn aquifer thickness is that's the depth over which melt percolates, below which some unmodeled physical processes initiate instantaneous run-off. From the tests I was running with sample data using the firn model, rarely did melt actually reach an impermeable layer (i.e. the ice base) more often in my tests all the melt was consumed to heat the snow pack. I didn't allow for ice lenses when testing the firn model, but that seems beyond our needed level of complexity.

Not really sure about this one, guess some test to see how important this actually is would be useful.

Given the higher mesh densities at the boundaries (and therefore finer mesh resolutions), the simple approach for filtering the fictitious ice thickness added by Elmer ( $h_{\rm min}$ ) no longer works correctly:

If the vertical mesh resolution is fine enough, the ice thickness in the second (or even third) vertical nodes from the bottom will still be less then $h_{\rm min}$ and therefore forced to zero.

Proposed Solution:

• Can we use the calculate_length function to aid in the filtering? Since it only used the ice-thickness along the free surface it shouldn't be subject to this same bug.

The variable name used to set the initial condition for the Surface_Enthalpy field was incorrectly spelled. No underscore included (i.e. Surface Enthalpy).

https://github.com/andrewdnolan/thermal-structure/blob/0fa70ed5d1041d26691d467032663b005d93f9b2/expr/03_PeriodicSurge/sifs/periodic_surge.sif#LL142C30-L142C30

This meant the field was never updated in time, and the value from the restart simulation was used for the entire diagnostic simulation, which in effect is a position dependent air temperature. This is only a problem as of the "new" heat pump parameterization, where the restart need to be specified in order for the fields to be updated in time.