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 !

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

    2 months later

    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:

    1. 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?

    2. 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?

    3. 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!

        wudipao88 Hello !

        1. 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).
        2. '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?
        3. 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

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

            1. 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?

            2. 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?

            3. 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.

            1. I would like to know how the englacial strain heating term is calculated in the model? Where can I see this?

            2. If I want to ignore the advection term or the englacial strain heating term, how should I set up the model?

            Best wishes!

              Hi @wudipao88

              I highly recommend that you read:

              (among others). To answer your questions:

              1. 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).

              2. 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)

              3. 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)

              4. Look at src/c/classes/Elements/Penta.cpp function Penta::ViscousHeating

              5. 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

                8 days later

                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:

                1. Using the initial surface velocity from Rignot (2012) to invert the basal friction coefficient through data assimilation
                2. 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.
                3. Update the ice hardness based on the output steady-state temperature field.
                4. Re-invert the basal friction coefficient based on the output steady-state velocity field and the updated ice hardness.
                5. Repeat the above steps until convergence is achieved.
                6. Output the temperature field.

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

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