Numerical approach of temperature distribution in a free boundary model for polythermal ice sheets

A numerical method for solving the thermal subproblem appearing in the modelization of polythermal ice sheets is described. This thermal problem mainly involves three nonlinearities: a reaction term due to the viscous dissipation, a Signorini boundary condition associated to the geothermic flux and an enthalpy term issued from the two phase Stefan formulation of the polythermal regime. The stationary temperature is obtained as the limit of an evolutive problem which is discretized in time with an upwind characteristics scheme and in space with finite elements. The nonlinearities are solved either by Newton-Raphson method or by duality techniques applied to maximal monotone operators. The application of the algorithms provides the dimensionless temperature distribution approximation and allows to identify the cold and temperate ice regions.