2D and 3D Simulation Data for: Numerical Determination of the Nonlinear Convection Threshold in Nitrogen Ice: Application to the Margins of the Cellular Plains in Sputnik Planitia, Pluto

Numerical simulations of bottom-heated convection in a layer of nitrogen ice under temperature conditions expected on Sputnik Planitia, Pluto. Numerical calculations are performed in a 12x1 box in 2D and a 12x12x1 box in 3D, where the top and bottom boundaries are maintained at constant temperatures of T_s and T_b, respectively. The top boundary is free-slip, whereas the bottom boundary is no-slip. The vertical boundaries are free-slip and thermally insulated. In all of the calculations shown here except one, the viscosity follows the power-law creep function, i.e., it depends on both temperature and stress, and the experimentally determined flow-law parameters for power-law creep in nitrogen ice are used (Yamashita et al., 2010). The 2D calculations are conducted using the Citcom code (Moresi & Solomatov, 1995) and 3D calculations using the CitcomCU code (Moresi and Gurnis, 1996; Zhong, 2006). The numerical domain is divided into cube-shaped finite elements, and the number of elements in the vertical direction is 32. For a 12x12x1 box, this resolution translates into 4,718,592 finite elements. We investigate the stability of the localized plume over a range of conditions relevant to Sputnik Planitia. These include T_s = 32 – 40 K, T_b = 47 – 63 K, and a range of experimentally constrained values for the flow-law parameters n, E, and A for power-law creep in nitrogen ice, where n is the stress exponent (n = 2.2 +/- 0.2), E is the activation energy (E = 3.5 +/- 0.5 kJ/mol), and A is the scaling coefficient. In the 3D calculations, we focus on only three sets of permissible values of these parameters, which are: n = 2.2, E = 3.5 kJ/mol, and A = 3.74x10^(-12) s^(-1) Pa^(-n) (reference values); n = 2.0, E = 3.2 kJ/mol, and A = 2.18x10^(-11) s^(-1) Pa^(-n); and n = 2.4, E = 3.8 kJ/mol, and A = 6.41x10^(-13) s^(-1) Pa^(-n). For each. set of assumptions on T_s, T_b, n, and E, the stability boundaries [in the Rayleigh number (Ra)-space] of the localized convection planform are determined using the approach described in previous studies (e.g., Solomatov, 2024). The results of our 2D and 3D simulations are provided in the attached document. The document contains 5 sheets. In the first sheet (titled '2D-Ra_loc,c'), we provide the values of the highest Rayleigh number at which localized convection decays, denoted by Ra(-)_loc,min, the lowest and highest Rayleigh numbers - Ra(+)_loc,min and Ra(-)_loc,max, respectively - at which a localized plume is stable, and the lowest Rayleigh number - Ra(+)_loc,max - at which the localized convection planform is unstable and transitions to widespread convection for each set of assumptions on n, E, T_s, and T_b in 2D simulations. The stability boundaries for localized convection in 3D are given in the sheet title '3D-Ra_loc,c'. The nondimensional values of the surface observables for the localized convection planform output by our 2D simulations, namely the width of the localized cell, w*_h, the peak amplitude of the surface topographic anomaly, A*_h, the average surface heat flux, F*_avg, and the average surface convective velocity, U*_avg, are provided in the sheet titled '2D-nondimen-params'. The same parameters output by our 3D simulations are provided in the sheet '3D-nondimen-params'. For each surface observable, the parameter values at the lowest and highest Ra for stable localized convection are given. The formulae for deriving the corresponding dimensional values are explained in the sheet "dimensionalization".