basal temperature calculation

wudipao88

Hi everyone

I would like to know whether it is possible to calculate the basal temperature based on known information (including geothermal heat flux, geometry information, surface information, etc.) using the thermal model? There is no tutorial example in the Userguide, so I would like to know how to implement it? What are the specific model inputs required?

Thanks !

mathieumorlighem

Hello @wudipao88

yes it is possible to do that. We have not had time to add tutorials on the thermal model yet, but you can for example take a look at tests/NightlyRun/test327.m. Make sure to set boundary conditions properly. Since you generally have advection, you need to constrain the temperature at all boundary where ice is coming in.
All the best
Mathieu

wudipao88

mathieumorlighem Hi Mathieu

I recently tried the calculation of the basal temperature according to this thermal model example (i.e. tests/NightlyRun/test327.m), and I encountered some confusion, such as:

If my model domain is located at ice divides, in this case, the ice in the model domain flows very slowly, how should I constrain the temperature where the ice flows in?
I want to solve for the steady-state basal temperature, so I change the transient solution to a steady-state solution. However, there seem to be two ways to do a steady-state solution. The first is to run md = solve(md, 'Thermal') when md.timestepping. time_step = 0, and the second is to run md = solve(md, 'steadystate'). However, the steady-state temperature fields obtained by these two methods are different. Why is that? Which method solves the temperature field correctly?
When performing a steady-state solution (here md.thermal.isenthalpy = 1 and md.thermal.isdynamicbasalspc = 1), I think the steady-state basal temperature (i.e. the temperature of the bottom grid) calculated using different md.basalforcings.geothermalflux should be different. However, whether I use md = solve(md, 'Thermal') or md = solve(md, 'steadystate'), the steady-state basal temperature I calculated based on different md.basalforcings.geothermalflux is the same. Why is this? Note that the temperatures of the non-bottom-most grids are different, only the bottom-most grid is the same.
However, when I do not use the enthalpy formula (i.e., md.thermal.isenthalpy = 0), the steady-state basal temperature calculated using different md.basalforcings.geothermalflux is different. Therefore, I don't know the reason, and I don't know whether I should turn on the enthalpy formula when calculating the basal temperature.

Looking forward to your reply, thank you!

mathieumorlighem

wudipao88 Hello !

if you are looking at a divide, you don't have any "incoming" ice and so you don't need to worry about lateral boundary conditions (just the surface temperature matters).
'thermal' takes the velocity from md.initialization.vx/vy to compute the advection, whereas 'steadystate' solves iteratively a stress balance (new vx and vy) and then temperature, until they both converge. So it is normal if they look different because most likely md.results.SteadystateSolution.Vx/Vy is potentially very different from md.initialization. I hope this makes sense?
changing the geothermal heat flux should definitely change your basal temperature. Make sure you are using the right units and that sliding speed is not too high.

I hope this helps
Mathieu

wudipao88

mathieumorlighem Thank you for your explanation, it is helpful! I still have some confusion:

If using md = solve(md, 'Thermal'), does md.initialization.vx/vy in the advection term correspond to the surface velocity, basal sliding velocity, or the depth-averaged velocity of the entire ice column?
How is the basal Neumann boundary condition set in the thermal model? Since this boundary condition is not directly visible in the model, is it equal to the geothermal heat flux plus the basal frictional heat flux? If so, is the basal frictional heat flux calculated by the basal shear stress multiplied by the basal sliding velocity? If using md = solve(md, 'Thermal'), is the basal sliding velocity then replaced by md.initialization.vx/vy? If using md = solve(md, 'steadystate'), is the basal sliding velocity corresponding to the velocity at the bottommost mesh obtained from the stress balance solution?
How is md.results.ThermalSolution.BasalforcingsGroundediceMeltingRate calculated in the model? Is it calculated using the following formula? When calculating it, does it not consider whether the basal temperature is above the pressure melting point? I noticed that even if the basal temperature is not above the pressure melting point, there is still a corresponding basal melting rate.

I would like to know how the englacial strain heating term is calculated in the model? Where can I see this?
If I want to ignore the advection term or the englacial strain heating term, how should I set up the model?

Best wishes!

mathieumorlighem

Hi @wudipao88

I highly recommend that you read:

(among others). To answer your questions:

they are the 3d velocities and you decide what you want to prescribe. By default extrude will copy the values that they have over the entire column (plug flow).
it is done automatically, take a look at src/c/analyses/EnthalpyAnalysis.cpp and EnthalpyAnalysis::CreatePVectorSheet we impose the following heat flux: heatflux=(basalfriction+geothermalflux). The velocity used are all from md.initialization except if you use steady state in which case they will be calculated until you reach a thermo-mechanical steady state (i.e. the flow speed is consistent with the thermal state and ∂T/∂t = 0)
No it is not calculated using the formula you are showing because this formula assumes that all the heat is converted to latent heat, whereas some of the heat can be sensible heat (and warms up the ice without melting). You can look at EnthalpyAnalysis::ComputeBasalMeltingrate (which follows the papers I highlighted above)
Look at src/c/classes/Elements/Penta.cpp function Penta::ViscousHeating
Then you will need to hard code it yourself in EnthalpyAnalysis::CreatePVectorVolume(setphi` to 0 after it is calculated and recompile) or set the velocity to 0 to stop advection.

I hope this helps
Mathieu

wudipao88

mathieumorlighem Thanks a lot for your patient help! This information is valuable to me!

wudipao88

mathieumorlighem Hi Mathieu,
According to the information you provided last time, I have learned a lot about thermal modeling. Now I am trying to reproduce the steady-state temperature simulation results (i.e. Figure 1c) from your previous research (i.e. Dependence of century-scale projections of the Greenland ice sheet on its thermal regime), so as to verify my understanding and use of the thermal model. My understanding of how to obtain Figure 1c is to follow the following iterative steps:

Using the initial surface velocity from Rignot (2012) to invert the basal friction coefficient through data assimilation
Input the inverted basal friction coefficient back into the model, then perform a thermo-mechanical steady-state simulation (i.e., using md = solve(md, 'steadystate') in ISSM), and output the velocity and temperature fields.
Update the ice hardness based on the output steady-state temperature field.
Re-invert the basal friction coefficient based on the output steady-state velocity field and the updated ice hardness.
Repeat the above steps until convergence is achieved.
Output the temperature field.

I am unsure if my understanding is correct? I await your response, thank you!

mathieumorlighem

Hello, yes that's correct. We needed to do a few iterations to make sure the velocity and temperature were consistent.

wudipao88

mathieumorlighem Hi Mathieu,
I now want to turn off horizontal heat conduction in the thermal model (but keep vertical conduction). Where should I change the source code? In addition, if I want to turn off the strain heat term, can I just set phi to 0 in EnthalpyAnalysis::CreatePVectorVolume?

mathieumorlighem

to remove the horizontal conditions you can change EnthalpyAnalysis::CreateKMatrixVolume:

      /*Conduction: */
      factor = D_scalar*kappa/rho_ice;
      for(int i=0;i<numnodes;i++){
         for(int j=0;j<numnodes;j++){
            Ke->values[i*numnodes+j] += factor*(
                     dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]);
         }
      }

change to

      /*Conduction: */
      factor = D_scalar*kappa/rho_ice;
      for(int i=0;i<numnodes;i++){
         for(int j=0;j<numnodes;j++){
            Ke->values[i*numnodes+j] += factor*(
                     dbasis[2*numnodes+j]*dbasis[2*numnodes+i]);
         }
      }

Note that some papers only make this approximation for convenience (to speed up the calculation, but we use a different approach so it would not change anything in the case of ISSM in terms of computing time), I personally don't think it makes sense.

and yes, you can set phi to 0 to remove strain heating.
Cheers
Mathieu