Implementation of the Depth Integrated Viscosity Approximation in SICOPOLIS

Implementation of the Depth Integrated Viscosity Approximation in SICOPOLIS First tests and results 17 September 2024 Internship, M2 SOAC intern: Félix Grandadam supervisor: Prof. Ralf Greve (ILTS) university supervisor: Dr. Louis Gostiaux Abstract In this report, we describe the first iteration of simulations following the implementation of the Depth Integrated Viscosity Approximation (DIVA) into SICOPOLIS, an open-source standalone 3D coupled thermomechanical ice sheet model. The approximation regards the dynamical part of the model, mainly the momentum balance equation, it was first developed by Goldberg (2011). Most of our internship consisted in properly coding this new dynamic, and the set of experiment performed in this report aims to identify any flaw in the implementation. Those experiments come from the European Ice Sheet Modeling Initiative (EISMINT) Payne et al. (2000) and creates a circular ice sheet, of which we studied the final steady state. DIVA was found to be behaving on par with the other dynamics present in the model, showing great adaptability in the different limit case it was subjected to. Some marginal concerns were raise as to how DIVA handles its ice front, creating seemingly more instabilities than the other dynamics, but with little repercussion on the ice-sheet global state. Hence, we advice to take the DIVA dynamics a step further in more complex simulations together with the Hybrid dynamics, like the Greenland spin-up (figure 19), to better discriminate their respective performances. Cite as: Grandadam, F. 2024. Implementation of the Depth Integrated Viscosity Approximation in SICOPOLIS. Internship Report, Claude Bernard University Lyon 1, France, and Hokkaido University, Sapporo, Japan (doi: 10.5281/zenodo.14732938). ii Contents 1 Introduction 1 2 Theory & model 1 2.1 Icedynamicstheory.................................... 1 2.1.1 incompressible thermovisquous fluids . . . . . . . . . . . . . . . . . . . . . 2 2.1.2 Isotropic, non-linear viscous fluid . . . . . . . . . . . . . . . . . . . . . . . . 2 2.1.3 Navier Stokes and hydrostatic approximation . . . . . . . . . . . . . . . . . 2 2.1.4 Depth Integrated Viscosity Approximation (DIVA) . . . . . . . . . . . . . . 5 2.1.5 Icethicknessequation .............................. 6 2.2 SICOPOLIS......................................... 7 3 Simulations details 7 4 Results & discussion 8 4.1 ExperimentA ....................................... 10 4.2 ExperimentG ....................................... 12 4.2.1 ExperimentG2 .................................. 15 4.3 ExperimentH ....................................... 16 4.4 Discussion ......................................... 19 5 Conclusion 19 6 Acknowledgment 21 7 Bibliography 22 7.1 Sigmatransformation................................... 24 7.2 Code ............................................ 24 7.3 expA higher resolution figure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 iii Félix Grandadam Internship, M2 SOAC 1 Introduction The accurate modeling of ice sheet dynamics is crucial for understanding their role in global climate regulation and predicting future sea-level rise. Given the complexities of ice sheet evolution, especially at a continental scale, direct measurements alone are insufficient for developing a comprehensive 3D mapping of these vast ice masses. This limitation is particularly pressing in the current climatic context, where ice sheets represent one of the most significant uncertainties in climate projections Bamber et al. (2019), Alley et al. (2005). The implementation of the Depth Integrated Viscosity Approximation (DIVA) particularly was motivated as it was designed to balance the need for detailed physical representation with computational efficiency, ensuring that it can be applied to large-scale ice sheets without prohibitive computational costs. Indeed, it achieves higher-order precision while needing a similar 2D partial differential equation solving, as lower orders dynamics already implemented in SICOPOLIS. By improving the accuracy of the dynamics in the model, we aim to provide more reliable projections of ice sheet behavior, which is essential for anticipating their impact on global sea levels and climate systems Alley et al. (2005). During the last decade, implementation of higher-order dynamics represented a significant step forward in ice sheet modeling, offering tools that can better capture the complex processes governing ice sheet evolution Goldberg (2011), Lipscomb et al. (2019). 2 Theory & model 2.1 Ice dynamics theory To further introduce the work realised during this internship, we first present key part of the stoke flow problem and common associated approximation in ice dynamics; This section is largely based on Greve and Blatter (2009). The stresses under which the ice is found plays a predominant role in the general flow of an ice sheet. The Cauchy stress tensor tdescribes such stresses at a particular point (or rather, at the surfaces of a small volume) and is defined as such : t=  txx txy txz tyx tyy tyz tzx tzy tzz (1) This tensor is symmetric, so for example, we wont distinguish txz from tzx. if n is the normal to a surface, then the stresses under which this surface is found is :  tn=t:n . For an ice sheet, a stress free upper boundary condition is inferred, so at the free surface h(x, y):  tn=t:n = 0(2) at the base of the ice sheet b(x, y), the dynamic boundary condition relates the basal drag τb(the shear component of the full basal stress  tn) to the basal velocity thanks to a sliding law, which can be written as : τb=β(vb, P)vb(3) with βa drag coefficient depending on the type of sliding law chosen. In SICOPOLIS, it is an empirical Weertman-Budd-type law, see J.Weertman (1957), Greve and Blatter (2009). 1 Félix Grandadam Internship, M2 SOAC We now introduce the strain-rate tensor D, its element are defined as such : Di,j =1 2∂vi ∂j +∂vj ∂i (i, j)∈ {x, y, z}(4) With the diagonal element associated with the stretching rates, and the off diagonal the shear rates (along the coordinates axes). 2.1.1 incompressible thermovisquous fluids To solve the mass, momemtum (and energy) balance, a closure relation is needed. Thus, we confine ourselves to incompressible thermovisquous fluids. Hence, the stress tensor tis only a function of (D, T, −−→ grad(T), P)with Tthe temperature and P the pressure. The incompressibility condition makes it that the pressure does not contribute to the deformation of the body; So, we define the deviatoric stress tD, the stress responsible for the deformation as : tD(D, T, −−→ grad(T), P) = t+PI(5) (with Pwhich can be defined as the negative average of the normal stresses : P=−1 3tr(t).) 2.1.2 Isotropic, non-linear viscous fluid On the anisotropy of ice : we suppose that the orientation distribution of the crytals composing a macroscopic compound is random, then, the anisotropy of the crytals balances out and ice behaves as an isotropic material. To then approach the behavior of ice, and particularly it’s tendency to undergo slow deformation while under persistent stress, a phenomenon known as creep; We use Glen’s flow law (Glen and Perutz (1955)), linking the deformation of ice to its viscosity and the strain-rate it undergoes : tD= 2ηD(6) with ηthe viscosity : η=1 2A(T′) −1 nd−(1−1 n) e(7) With A(T′)the rate factor, an Arrhenius type law function of T′the melting point of the ice. deis known as the effective strain rate : de=q1 2tr(D.D). And with n the stress exponent, usually equal to 3, to fit the observed deformation of ice under continuous stress (creep). The component form of the effective strain rate is as follow (using the mass conservation eq 9): de=qD2 xx +D2 yy +DxxDyy +D2 xy +D2 xz +D2 yz (8) 2.1.3 Navier Stokes and hydrostatic approximation Under these conditions, the full Stokes flow problem is then : Mass balance : div(−→ v) = 0 (9) Momentum balance : ρD−→ v Dt =div(t) + −→ f(10) 2 Félix Grandadam Internship, M2 SOAC which then leads in component form, using the definition of the deviatoric stress eq.5, to : ρDvi Dt =−∂P ∂i +∂tD ii ∂i +∂tD ij ∂j +∂tD ik ∂k +fi(i, j, k)∈ {x, y, z}(11) Then, characteristic values for an ice sheet horizontal and vertical extent is established table 1 : horizontal extent [L] 1000 km vertical extent [H] 1 km horizontal velocity [U] 100 m.a−1 vertical velocity [W] 0.1 m.a−1 Pressure [P] ρg[H]≈10 MPa time-scale [t] [L]/[U] = [H]/[W] = 10 000 a Table 1: Characteristic values of an ice sheet, from Greve and Blatter (2009) From this values, one can realise that : the Froude number ( F r =q[U] g[H]the ratio of the flow inertia to the gravitational force field ) is very small (≈10−15 at most). Hence, the acceleration term in the momentum balance is negligible for the flow of ice sheets. Also, the Coriolis term, when compared to the pressure gradient appears to be negligible : 2ρΩ[U] [P]/[L]≈10−8, so the external force field can be approximated to only the gravitational field. Moreover, with the hydrostatic condition, one can show that the vertical normal stress is tzz ≈ [P]≈10MPa, while the shear vertical stresses txz, tyz ≤100kPa, so tzz ≫txz, tyz This leads to a simplified momentum balance from eq 10, known as the hydrostatic full Stoke force balance : ∂txx ∂x +∂txy ∂y +∂txz ∂z = 0 ∂tyy ∂y +∂txy ∂x +∂tyz ∂z = 0 ∂tzz ∂z =ρg With, in green the gradients of the longitudinal and lateral shear stresses respectively, also known together as membrane stresses, and in red the gradient of the vertical shear stress. terms associated with either of these 2 stresses will conserve this color code from now on. Then, by integrating the vertical component we directly have, using the stress free boundary condition eq2 : tzz =ρg(z−h), by then using the definition of the deviatoric stress eq5 a relationship between tD xx,tD yy and Pcan be found. Furthermore, by comparing the velocity gradients (composing the deviatoric stresses), one can observe that, due to the low aspect ratio ϵof an ice sheet, the horizontal gradient of the vertical velocity is negligible in front of the vertical gradient of the horizontal velocity : ∂hvz ∂zvh ≈[H][W] [U][L]=ϵ2= 10−6(12) 3 Félix Grandadam Internship, M2 SOAC These approximations lead to, using the deviatoric stress definition and the Glen’s flow law eq 5, 6 : ∂ ∂x 2η2∂vx ∂x +∂vy ∂y +∂ ∂y η∂vx ∂y +∂vy ∂x +∂ ∂z η∂vx ∂z =ρg ∂h ∂x ∂ ∂y 2η2∂vy ∂y +∂vx ∂x +∂ ∂x η∂vx ∂y +∂vy ∂x +∂ ∂z η∂vy ∂z =ρg ∂h ∂y (13) This set of approximation is known as the Blatter-Pattyn approximation ( Blatter (1995), Pattyn (2002) ). ( The vertical velocity can then be found with the mass conservation as it is fully decoupled from the horizontal mass balance, see Greve and Blatter (2009) for the full resolution.) The effective strain-rate in this case is as follow : dBP e=r(∂xvx)2+ (∂yvy)2+∂xvx.∂yvy+1 4(∂yvx+∂xvy)2+1 4(∂zvx)2+1 4(∂zvy)2(14) The Blatter-Pattyn approximation needs the resolution of a system of Partial Differential Equations (PDE) in 3 dimensions which is costly in term of processing power, so further approximation have been made to avoid solving such a system. 2 approximations are especially worth mentioning : The Shallow Ice Approximation (SIA) and the Shallow Shelf Approximation (SSA). The Shallow Ice Approximation (SIA) is the oldest, simplest and most widely used approximation of Blatter-Pattyn, where the membrane stress (terms in green) are neglected from equation 13. This leads to a balance between the vertical shear stress and the gravitational driving stress only : ∂ ∂z η∂vx ∂z =ρg ∂h ∂x ∂ ∂z η∂vy ∂z =ρg ∂h ∂y (15) with the effective strain-rate composed of only the shear terms : dSIA e=r1 4(∂zvx)2+1 4(∂zvy)2(16) The SIA is representing the typical flow of a cold, slow flowing ice with no (or little) sliding where creeping is predominant. On the opposite end, there is the Shallow Shelf Approximation (SSA) MacAyeal (1989), a dynamic originating from the plug flow behaviour of an ice shelf (floating ice) where the horizontal velocities have no vertical dependency due to negligible basal drag. In SSA, by dropping the vertical shear stress on equation 13 ( terms in red), we allow the resolution to happen for a depth invariant velocity field, leading to a system of PDE in 2D for the vertical averaged of the velocity field simply by integrating over the depth the modified equation 13. This leads to : 4 Félix Grandadam Internship, M2 SOAC ∂ ∂x 2H¯η2∂¯vx ∂x +∂¯vy ∂y +∂ ∂y H¯η∂¯vx ∂y +∂¯vy ∂x =ρgH ∂h ∂x ∂ ∂y 2H¯η2∂¯vy ∂y +∂¯vx ∂x +∂ ∂x H¯η∂¯vx ∂y +∂¯vy ∂x =ρgH ∂h ∂y (17) With H the ice thickness (defined as H(x, y) = b(x, y)−h(x, y)), and ¯vxthe mean value of vxover the depth : ¯vx=1 HRh bvx(z)dz and similarly for ¯vy,¯η. Similarly, the effective strain rate uses the mean velocity and the vertical shear terms are discarded: dSSA e=r(∂x¯vx)2+ (∂y¯vy)2+∂x¯vx.∂y¯vy+1 4(∂y¯vx+∂x¯vy)2(18) It is to be noted that, we can add back a basal drag term over the mean velocity, so acting on the whole column, to add a tuning parameter and slow down the flow where needed. This modification is sometimes called the Shelfy Stream Approximation. (SSTA) The plug-flow we get from these approximations (SSA and SSTA) are considered valid for fastflowing ice due to high basal sliding. A third dynamic used a lot in more recent models, is an hybrid dynamics (Bernales et al. (2017)), which acknowledges that the 2 previously shown approximations (SIA and SSTA/SSA) describe 2 opposite and complementary flow regimes. Thus by merging numerically the output of the 2 approximations, we can derive a result adequate for both fast and slow flowing ice. This numerical solution is then simply called the hybrid dynamics. Still, to achieve a near-similar level of accuracy as the BP approximation with an higher order, theory-derived equation, while solving a 2D PDE system only, Goldberg and others, defined the Depth Integrated Viscosity Approximation (DIVA) Goldberg (2011). 2.1.4 Depth Integrated Viscosity Approximation (DIVA) The DIVA dynamics has been developed "by making approximations to the functional that yields the Euler–Lagrange equations, instead of to the equations themselves" Goldberg (2011). The work presented in this part is mostly based on Lipscomb et al. (2019). This development, beyond the scope of this report, allows one to write a momentum balance and an effective stress, similar to Blatter-Pattyn eq.13, with the membrane stress terms computed with the mean velocities, as in SSA : Hence the effective strain rate for DIVA : dDIV A e=r(∂x¯vx)2+ (∂y¯vy)2+∂x¯vx.∂y¯vy+1 4(∂y¯vx+∂x¯vy)2+1 4(∂zvx)2+1 4(∂zvy)2(19) And the momemtum balance equation : 1 H ∂ ∂x 2H¯η2∂¯vx ∂x +∂¯vy ∂y +1 H ∂ ∂y H¯η∂¯vx ∂y +∂¯vy ∂x +∂ ∂z η∂vx ∂z =ρg ∂h ∂x 1 H ∂ ∂y 2H¯η2∂¯vy ∂y +∂¯vx ∂x +1 H ∂ ∂x H¯η∂¯vx ∂y +∂¯vy ∂x +∂ ∂z η∂vy ∂z =ρg ∂h ∂y 5 Félix Grandadam Internship, M2 SOAC Then, the membrane stress is independent of depth, so we can integrate over the depth, and use the boundary conditions from equations 3, 2 to substitute the vertical shear term for the basal stress only : ∂ ∂x 2H¯η2∂¯vx ∂x +∂¯vy ∂y +∂ ∂y H¯η∂¯vx ∂y +∂¯vy ∂x −τb,x =ρgH ∂h ∂x ∂ ∂y 2H¯η2∂¯vy ∂y +∂¯vx ∂x +∂ ∂x H¯η∂¯vx ∂y +∂¯vy ∂x −τb,y =ρgH ∂h ∂y (20) With τb,i =βvb,i the basal shear stress in the i direction as defined by the boundary conditions equation 3. For this equation to be a 2D PDE, we must express the basal stress as a function of the mean velocity. Following Arthern et al. (2015) we define some useful integral Fnfor clarity : Fn=Zh b 1 ηh−z Hn dz (21) Then, by re-arranging a few terms, Goldberg (2011) showed that the vertical profile of viis related to the basal velocity by : vi=vb,i +βvb,i Zz b 1 ηh−z′ Hdz′(22) By integrating, Arthern et al. (2015) as demonstrated that we can relate the basal velocity vb,i to the mean velocity ¯vi: ¯vi=vb,i(1 + βF2)(23) We can then define an effective friction parameter βeff =β 1+βF2to be able to express, from equation 23 and 3 the basal stress as a function of the mean velocity : τb,i =βvb,i =βeff ¯vi(24) Which means the PDE to solve (eq 20) can be expressed with the mean velocity exclusively and is thus 2D. 2.1.5 Ice thickness equation The horizontal velocity is the key component to compute the ice volume flux terms in the ice thickness equation : ∂tH=−(∂x(Qx) + ∂y(Qy)) + as(25) Where Hdenotes the ice thickness (H=h−b), asthe accumulation rate, and Qithe icomponents of the volume flux (vertically integrated horizontal velocity). See Greve and Blatter (2009). 6 Félix Grandadam Internship, M2 SOAC Figure 7: top view elevation map of the 4 dynamics (diva, hyb, sia, ssta) for experiment G at steady state. Regarding the horizontal velocity, the simulations stay in very good agreement with each other at first glance : Both the mean velocity distribution figure 8 and the vertical profile of the horizontal velocity at y=0.0 figure 9 show very little difference : Figure 8: mean horizontal velocity distribution of the 4 dynamics (diva, hyb, sia, ssta) for experiment G at steady state. The mean velocity is at around 55m/a with a maximum not going over 80 m/a across all smulations. Seemingly, the vertical profile figure 9 shows almost no difference between the different dynamics. This can be explain by the fact that in this limit case, the velocity is not driven by internal deformation but mainly by basal sliding, setup with the same Weertman-type sliding law across all simulations. So the dynamics do not play a predominant role in this experiment. 13 Félix Grandadam Internship, M2 SOAC Figure 9: side view of the absolute horizontal velocity profile at y=0.0 of the 4 dynamics (diva, hyb, sia, ssta) for experiment G at steady state. Indeed, if we look at figure 10, where is represented the 2D maps of basal and surface horizontal velocity respectively. The aspect and values of all dynamics are of the same order of magnitude. Figure 10: top view maps of the basal and surface horizontal velocities respectively of the 4 dynamics (diva, hyb, sia, ssta) for experiment G at steady state. Going further, if we subtract the basal velocity to the surface one, we obtain the surface velocity map due to internal deformation. This is what is depicted figure 11. Then it becomes apparent than, for one : the hybrid simulation is dominated by the SSTA regime and almost display a complete plug-flow regime (as is SSTA). For two : the internal ice deformation is responsible, both in DIVA and SIA, for only 5m/a at most in the surface velocity, which is little, but not negligible when striving for precision. Hence, DIVA displays a good response to the experiment, and is, on the contrary to experiment A remarkably stable, even at the margins. 14 Félix Grandadam Internship, M2 SOAC Figure 11: top view map of the difference between surface and basal velocity of the 4 dynamics (diva, hyb, sia, ssta) for experiment G at steady state. Still, in order to better discriminate the response of the different dynamics in this setup, we ran an additional simulation where we tuned down the basal sliding coefficient, to still allow sliding everywhere, but have it of the same order of magnitude as the internal deformation. This is experiment G2 : 4.2.1 Experiment G2 This slower setup shows the limit of using SSTA dynamics : Figure 12, showing the elevation profiles of the different simulation, shows that the SSTA is drastically different, with a summit at more than 4500m whereas all the other simulation are lower than 3500m. Figure 12: sideview of the elevation, through the middle of the ice sheet (x=0) for experiment G2 at steady state. The dots represent the data points (25km resolution). This difference comes from the imposed plug flow regime of the SSTA dynamics which fails to take into account ice deformation. Indeed if one looks at the velocity profile figure 13, this issue becomes apparent with all simulation but SSTA showing a non negligible velocity variation depthwise, especially at the margin where internal deformation is responsible for the majority of the surface velocity: from 20m/a near the ice base to around 60m/a close to the surface. 15 Félix Grandadam Internship, M2 SOAC Figure 13: side view of the absolute horizontal velocity profile at y=0.0 of the 4 dynamics (diva, hyb, sia, ssta) for experiment G2 at steady state. On the global distribution of the mean velocity observed figure 14, all simulations but SSTA seem to deal well with the ice deformation being not negligible. DIVA appears to have higher maximum values as in experiment A but to a lesser degree (less than 15m/a difference). Figure 14: mean horizontal velocity distribution of the 4 dynamics (diva, hyb, sia, ssta) for experiment G2 at steady state. 4.3 Experiment H In this experiment, basal sliding is authorized only where the ice base is at its melting point. This hypothesis is deemed more realistic as basal melt-water likely act as a lubricant at the ice base. This setup induces a discontinuity between cold and temperate ice, which translate to more complex behavior in the ice sheet. The radial symmetry is broken and a spike patterning is observed across all simulation. This effect happened for almost all model originally tested in Payne et al. (2000). 16 Félix Grandadam Internship, M2 SOAC Figure 15: top view elevation map of the 3 dynamics (diva, hyb, sia) for experiment H at steady state. Experiment H, produces a stable steady state with ridges and valleys on the ice sheet, visible on figure 15. These ridges follow as one could expect the low velocity area, where the spike patterning is really visible : see figure 16 Figure 16: top view maps of the basal and surface horizontal velocities respectively of the 3 dynamics (diva, hyb, sia) for experiment H at steady state. 17 Félix Grandadam Internship, M2 SOAC The specific patterns are simulations dependent. SIA dynamics is the most different whereas DIVA and Hybrid seems to have a somewhat similar distribution. But, when looking at the surface velocity distribution due to ice deformation figure 17;The Hybrid simulation, due to its nature, seems to really be stretch between a slow, deformation-driven ice flow at the center and ridges of the ice sheet and a fast plug-flow regime in the valleys and margin, with almost no ice deformation. DIVA on the other hand, seems to have a smoother transition between slow and fast flowing ice, and a somewhat uniform ice deformation in all the ice sheet. Still, DIVA has rather chaotic margins when compared to the other simulations. Figure 17: top view map of the difference between surface and basal velocity of the 3 dynamics (diva, hyb, sia) for experiment H at steady state. Although the velocity and patterning differences, all the simulation remain clustered around the same volume/surface area within a 1.5% difference of one another. And regarding the computation time, DIVA appears to be more costly, taking 20% more time to complete than the Hybrid simulation. Simulation Total Volume (m3) Total Area (m2) computation time SIA 1.695e+15 1.031e+12 12min 03s Hybrid 1.717e+15 1.031e+12 21min 43s DIVA 1.719e+15 1.038e+12 27min 18s Average 1.710e+15 1.033e+12 / Table 5: Steady state values for Experiment H In the same fashion as the other experiments, when looking at the mean velocity distribution figure 18, on one hand the global mean horizontal velocity of DIVA is lower, being at around 50m/a while Hybrid and SIA are at around 60m/a. On the other hand the maximum velocity reach using DIVA is well above the one in the Hybrid simulation. Those maximum are happening near the margin and are probably not physical but more due to instabilities, which we will discuss in the next section. 18 Félix Grandadam Internship, M2 SOAC Figure 18: mean horizontal velocity distribution of the 3 dynamics (diva, hyb, sia) for experiment G at steady state. 4.4 Discussion The main issue DIVA is dealing with are the "messy" margin some setups can produce. Looking further into this, it seems that, in experiment A and H, DIVA does not reach a perfect steady state and has little elevation variation near the margins. This would explain a few of the velocity maximum observed near the margins in experiment A and H (figures 4, 16). These isolated maxima are also highlighted by the boxplots (figures 6, 18 ). All simulations seems to have a certain amount of instabilities at their margin, this could be partly linked to the steep slope created at the margin, which can render some of the approximation leading up to the dynamics used inadequate, especially at lower spatial resolution. It seems also linked to the discontinuity induce between temperate and cold ice base, as the pattern of some of the instabilities observed are closely correlated to the cold-temperate ice base distribution. Our best guess as to why this instability occurs more in DIVA is linked to the Arakawa grid system used in SICOPOLIS. Indeed, in DIVA, multiple variables need to be interpolated to the staggered or main grid. But interpolating linearly a non linear equation, like equation 21 for example, can lead to error, especially at lower resolution and when on the margin, as the 2 points used for interpolation could very well not both be glaciated points, or in the same thermodynamic state. This cases are hard to catch and can potentially lead to non negligible errors, especially as it can rely on variables being defined outside of their physical domain. To further support our reasoning, we tried running these experiments at a higher resolution (5km instead of 25km) for DIVA and observe a net amelioration in the symmetry and cleanness of the ice sheet near the margins (see Appendix). 5 Conclusion To conclude, during this internship, a new dynamics was implemented in the SICOPOLIS model, this new tool should theoretically allow 3D modelling with a better accuracy and limited additional computation time. In this report, we presented a first iteration of test from EISMINT Payne et al. (2000) to ensure DIVA is properly implemented in SICOPOLIS. These different experiment, showed 19 Félix Grandadam Internship, M2 SOAC that DIVA was on par with the other dynamics. Showcasing good adaptability to the different limit cases (fast and slow flowing ice) it was subjected to. This apparent versatility is a really good point for DIVA, as more complex ice sheet obviously share fast and slow flowing region. Having a single, theory derived dynamics being able to encapsulate both of these state should add stability. Likely observation of this better handling of transitional region from slow to fast flowing ice has been made in both experiment A and H. Up to now, SICOPOLIS used the hybrid dynamics to numerically merge SIA and SSTA dynamics into a single simulation. On the contrary, DIVA seems to struggle to maintain a full steady state at its margin, which creates an unclear ice front and velocity profile. This potential issue could be linked to how the approximation is coded in SICOPOLIS as discussed previously. If this turns out to be the case, then this issue is surely fixable in the future by trying to deal more carefully with interpolation. Still, although visually obvious in our set of experiment on any 2D top view map, and creating few maxima on the velocity field, there is little evidence as if this issue creates any durable changes to the final steady state of the ice sheet. The surface area and volume of the different simulation being very clustered for every experiment. Regarding the actual performance of DIVA, it is hard to quantify the gains of this new dynamics when compared to the others as the experiments were not meant for this and provided little discrimination between the simulations. A future set of more complex experiment and comparison should be perform to quantify the potential precision gain of the depth integrated viscosity approximation. Looking at the computation time, DIVA seems less efficient than hybrid, being around 5min or 20% slower in experiments G and H. Still it is hard to quantify how much more costly it actually is with only these few experiments, as, this added time could very well be not linear with more complex, longer simulations. The set of experiment performed in this report was meant to identify any flaw in the implementation of the new approximation. It did raise marginal concern as to how the ice front is handled, but, the overall results being on par with the other dynamics, DIVA looks ready to be tested for heavier and more complex simulation. Especially, we intend to use this new dynamics for a paleo spin-up simulation of Greenland, which aims to recreate the realistic, 3D present day Greenland ice sheet, as shown on the first page of this report figure 19 with the Hybrid dynamics. This kind of spin-up are capital to accurately predict future evolution as the history of an ice sheet can be felt during thousands of year through for example : ice temperature and isostasy adjustment. Hence the need for a dynamic able to accurately and efficiently model 3D ice sheet like DIVA. Figure 19: Observed (left) vs. simulated (right) surface velocity of the present-day GrIS. Hybrid dynamics, and two resolutions (low/40 km, high/4 km) were employed. From "A multi-phase spin-up method for the Greenland ice sheet, and its influence on future changes of the ice sheet" (2024.9.18 JSSI & JSSE Joint Conference, by R.Greve, C.J.Berends, J.Bernales, F.Grandadam) 20 Félix Grandadam Internship, M2 SOAC 6 Acknowledgment I would like to thanks the ILTS institute and my supervisor Prof. Ralf Greve for the opportunity they gave me with this internship and the caring feedback provided. Ms. Arakawa for all the work she did to prepare my stay in Japan. Dr. Combot, Dr. Kuhfeld, Dr. Do, Imazu, Yazawa and all the other students and P.H.Ds in my office for their support. Finally i would like to thanks the M2 ocean promotion for this last year, i would have not done it without them. 21 Félix Grandadam Internship, M2 SOAC References R. Alley, P. Clark, P. Huybrechts, and I. Joughin. Ice-sheet and sea-level changes. Science, 310:456 – 460, 2005. doi: 10.1126/SCIENCE.1114613. A. Arakawa and V. R. Lamb. Computational design of the basic dynamical processes of the ucla general circulation model. In J. CHANG, editor, General Circulation Models of the Atmosphere, volume 17 of Methods in Computational Physics: Advances in Research and Applications, pages 173–265. Elsevier, 1977. doi: https://doi.org/10.1016/B978-0-12-460817-7. 50009-4. URL https://www.sciencedirect.com/science/article/pii/ B9780124608177500094. R. J. Arthern, R. C. A. Hindmarsh, and C. R. Williams. Flow speed within the antarctic ice sheet and its controls inferred from satellite observations. Journal of Geophysical Research: Earth Surface, 120(7):1171–1188, 2015. doi: https://doi.org/10.1002/2014JF003239. URL https://agupubs. onlinelibrary.wiley.com/doi/abs/10.1002/2014JF003239. J. Bamber, M. Oppenheimer, R. Kopp, W. Aspinall, and R. Cooke. Ice sheet contributions to future sea-level rise from structured expert judgment. Proceedings of the National Academy of Sciences of the United States of America, 116:11195 – 11200, 2019. doi: 10.1073/pnas.1817205116. J. Bernales, I. Rogozhina, R. Greve, and M. Thomas. Comparison of hybrid schemes for the combination of shallow approximations in numerical simulations of the antarctic ice sheet. The Cryosphere, 11(1):247–265, 2017. doi: 10.5194/tc-11-247-2017. URL https://tc.copernicus.org/ articles/11/247/2017/. H. Blatter. Velocity and stress fields in grounded glaciers: a simple algorithm for including deviatoric stress gradients. Journal of Glaciology, 41(138):333–344, 1995. doi: 10.3189/S002214300001621X. J. W. Glen and M. F. Perutz. The creep of polycrystalline ice. Mathematical and Physical Sciences, (228), 1955. D. N. Goldberg. A variationally derived, depth-integrated approximation to a higher-order glaciological flow model. Journal of Glaciology, 57(201):157–170, 2011. doi: 10.3189/002214311795306763. R. Greve and H. Blatter. Dynamics of ice sheets and glaciers. 2009. URL https://api. semanticscholar.org/CorpusID:128734526. J.Weertman. On the sliding of glacier. Journal of Glaciology, 3(21), 1957. W. H. Lipscomb, S. F. Price, M. J. Hoffman, G. R. Leguy, A. R. Bennett, S. L. Bradley, K. J. Evans, J. G. Fyke, J. H. Kennedy, M. Perego, D. M. Ranken, W. J. Sacks, A. G. Salinger, L. J. Vargo, and P. H. Worley. Description and evaluation of the community ice sheet model (cism) v2.1. Geoscientific Model Development, 12(1):387–424, 2019. doi: 10.5194/gmd-12-387-2019. URL https://gmd. copernicus.org/articles/12/387/2019/. D. R. MacAyeal. Large-scale ice flow over a viscous basal sediment: Theory and application to ice stream b, antarctica. Journal of Geophysical Research: Solid Earth, 94(B4):4071–4087, 1989. doi: https://doi.org/10.1029/JB094iB04p04071. URL https://agupubs.onlinelibrary. wiley.com/doi/abs/10.1029/JB094iB04p04071. F. Pattyn. Transient glacier response with a higher-order numerical ice-flow model. Journal of Glaciology, 48(162):467–477, 2002. doi: 10.3189/172756502781831278. 22