Multi-Scale Sea Ice Kinematics Modeling with a Grid Hierarchy in Community Earth System Model (version 1.2.1)

High-resolution sea ice modeling is becoming widely available for both operational forecasts and climate studies. Sea ice kinematics is the most prominent feature of high-resolution simulations, and with rheology models such as ViscousPlastic, current models are able to reproduce multi-fractality and linear kinematic features in satellite observations. In this study, we carry out multi-scale sea ice modeling with Community Earth System Model (CESM) by using a grid hierarchy (22 km, 7.3 km, and 2.5 km grid stepping in the Arctic). By using atmospherically forced experiments, we simulate consistent sea 5 ice climatology across the 3 resolutions. Furthermore, the model reproduces reasonable sea ice kinematics, including multifractal deformation and scaling properties that are temporally changing and dependent on circulation patterns and forcings (e.g., Arctic Oscillation). With the grid hierarchy, we are able to evaluate the model’s effective spatial resolution regarding the statistics of kinematics, which is estimated to be about 6 to 7 times that of the grid’s native resolution. Besides, we show that in our model, the convergence of the Elastic-Viscous-Plastic (EVP) rheology scheme plays an important role in reproducing 10 reasonable kinematics statistics, and more strikingly, simulates systematically thinner sea ice than the standard, non-convergent experiments in landfast ice regions of Canadian Arctic Archipelago. Given the wide adoption of EVP and subcycling settings in current models, it highlights the importance of EVP convergence especially for climate studies and projections. The new grids and the model integration in CESM are openly provided for public use.


Introduction
Sea ice is the interface between the polar atmosphere and ocean, and therefore an important modulating factor of the polar air-sea interactions. The momentum input into the sea ice from the atmosphere and the ocean causes sea ice drift, as well as failures and deformations at a wide range of temporal and spatial scales. As revealed by high-resolution, kilometer-scale satellite observations with Synthetic Aperture Radars (SAR), the deformation of the sea ice is shown to be multi-fractal, with 2 Grid Generation and Model Integration

TS grids -a tripolar grid hierarchy
We design a new tripolar grid generation method for global ocean-sea ice modeling. The generated grid is orthogonal and compatible with many existing ocean and sea ice models, including POP, CICE, and MOM. As shown in Figure 2, it consists 10 of two patches, the southern patch (SP) and the northern patch (NP), divided by a certain latitudinal circle (φ) in the northern hemisphere. For SP, the grid lines are purely zonal or meridional. For NP, there are two north poles for the grid, which are placed on the Eurasia and North America land masses. However, unlike many tripolar grids such as TX0.1V2 in CESM, this new grid features smooth grid scale transition on the boundary of SP and NP. In order to achieve this, a non-trivial grid generation process is carried out on the stereographic projection of NP. Specifically, a series of embedded ellipses are constructed based 15 on resolution requirements, with: (1) the outermost ellipse as the projection of the latitudinal circle at φ (i.e., a circle), and (2) the innermost ellipse as the line linking the two grid poles. The foci of the outermost ellipse are both at the North Pole, and with the progression to inner ellipses, they gradually move towards the two grid poles, respectively. The ellipses form the "zonal" grid lines in NP. After the ellipses are constructed, the "meridional" grid lines in NP are constructed, starting from the southern boundary of NP, down to the innermost ellipse. During this process, it is ensured that: (1) the grid lines are 20 constructed consecutively between adjacent ellipses, with starting points on the outer ellipse; and (2) they are linked with the Sub-mesoscale capable meridional grid lines in SP to ensure overall continuity of the global grid. The details of the construction of ellipses and the smooth transition of meridional grid scales are further described in Appendix A. Furthermore, the meridional grid scales in both SP and NP are constructed to alleviate the grid aspect ratio, reducing meridional grid sizes in higher latitudes.
By using the grid generation method, we generate a suite of tripolar grids: TS045, TS015 and TS005. These grids all have the the boundary between SP and NP (φ ) at 10 • N , and grid poles at: (63 • N, 104 • W ) and (59.5 • N, 76 • E). Table 1 shows the 5 detailed configuration of these grids. TS045 is the coarsest grid with the nominal resolution of 0.45 • , and it targets at climate modeling and the typical resolution range for ocean component models in Coupled Model Intercomparison Projects (CMI).
TS015 and TS005 are about 3 times (or nominal 0.15 • ) and 9 times (or nominal 0.05 • ) the resolution of TS045, respectively. Figure 3 shows the horizontal grid scale (s = (dx 2 + dy 2 )/2) in polar regions for TS015 and TS005. The average grid scale in TS005 is 2.45 km in the Arctic oceanic regions (north of 65 • N ), and the scales of TS015 and TS045 are 7.34 km 10 and 22.01 km, respectively. From the ocean modeling perspective, we use the Rossby deformation radii (R) as the proxy for mesoscale (Chelton et al., 1998), and investigate the capabilities of each grid. In specific, the criterion in Hallberg (2013) is adopted: the mesoscale resolving is attained when s is smaller than the half of the local value of R. Based on the annual mean WOA13 climatology of salinity and temperature Zweng et al., 2013), we construct the global distribution of R. In Table 1 we note that TS015 is "almost eddying", because this grid attains mesoscale-resolving for 65% of remapping based advection (Dukowicz and Baumgardner, 2000), and (4) EVP rheology model (Hunke and Dukowicz, 1997), the ridging/rafting scheme and ice strength model (Lipscomb et al., 2007;Rothrock, 1975). For the new grids, we choose shorter thermodynamics and dynamics time steps for grids with higher resolution (Tab. 2). Since sea ice kinematics is the focus of this study, a series of EVP subcycle numbers are chosen for each grid. While 120 subcycles per hour is adopted for some 1 • resolution CMIP simulations (Jahn et al., 2012;Xu et al., 2013), we experiment with larger values of cycle counts (2880 5 subcycles per hour or higher), as shown in Table 2.
In our study, CICE is coupled to the slab ocean model (SOM) in CESM, which provides a climatological seasonal cycle of ocean mixed layer depth and heat potential. The same configuration for SOM is used for experiments with TS045, TS015 and TS005. The reason why we use SOM instead of a fully dynamic-thermodynamic ocean model (such as POP) is three fold.
First, in this study we focus on the simulation of sea ice kinematics and the inter-comparison across different resolution settings.

10
Therefore, by using a single-column model for the ocean, we eliminate the factors that may compromise the comparability, including the inconsistency in modeled ocean status across the resolution range, ocean and coupled turbulence, etc. Second, since atmospheric forcing is the major driver of sea ice drift and kinematics, we consider SOM eligible for the purpose of  6 https://doi.org/10.5194/gmd-2020-160 Preprint. Discussion started: 2 June 2020 c Author(s) 2020. CC BY 4.0 License. Ocean Model Intercomparison Project (Griffies et al., 2016). Specifically, CORE dataset is a Normal-Year Forcing (NYF) with the climatological annual cycle based on NCEP atmospheric reanalysis, and it has a spatial resolution of about 2 • (T62) with four-times daily wind stress fields. Following the common practice in CESM, for the coupling of atmospheric state variables (such as air temperature, humidity, etc.), a bilinear interpolator is used between T62 and each TS grid. For wind stress, the patch-recovery algorithm is adopted to ensure good structure of wind fields on the ocean-sea ice grid. For fluxes, a first-order 10 conservative interpolator is utilized. All these interpolators (in total 6) are generated through CESM mapping toolkit and ESMF regridding toolkit (https://www.earthsystemcog.org/projects/esmf/).
All the model integration for TS grids in CESM, including grid files in the format of POP, and interpolation files, are openly available (details in Code and data availability). The time stepping are based on 240 EVP subcycles per dynamics time step (NDTE) for TS045 and TS015, following the settings in Table 2. The experiments are outlined in Figure 4.a. Specifically, both experiments with TS045 and TS015 starts on Jan. 1st with no sea ice, and the model is gradually reaching an equilibrium state for both sea ice coverage and volume. We As shown in Figure 4.b, for both TS045 and TS015, the Arctic sea ice extent and volume approach equilibrium at about the 5th and 30th year, respectively. In their equilibrium states ( With good consistency of modeled ice climatology between TS045 and TS015 experiments, we adopt a warm start-up 10 strategy for TS005 ( Fig. 4.a). We map the CICE restarting file from TS015 in year 36 (0:00 on Jan. 1st) onto TS005, and continue the experiment with TS005 to year 42. Therefore, the timeline of TS005 simulation is aligned to year 36 of TS015 in Figure 4.b. In Figure 5.a, b, e and h, we also show the results with TS005. With the warm start-up, the experiments with TS005 approaches equilibrium towards year 42. Similar to TS045 and TS015, the experiment with TS005 produces reasonable sea ice climatologies, but with a minor decrease with respect to TS045 and TS015. The major difference is the reduced thickness for 15 all the seasons, and the smaller sea ice coverage during summer. The overall sea ice coverage and volume of TS005 is also in good agreement with satellite observations and PIOMAS dataset. Since we adopt the same configurations for parameterization schemes across all the three grids, we consider the results reasonable and further carry out the analysis of sea ice kinematics with high-frequency outputs of these experiments.

Experiments and Analysis
One major obstacle of running high-resolution models is the huge amount of computational overhead and the long durations 20 for climate simulations. In this study, we utilize an Intel processor-based cluster with 24 cores/node for all the experiments.
For TS045 and TS015, more than 5 simulated year per day (SYPD) can be attained with less than 1000 cores for all the experiments. For TS005, the simulation speed is about 1 SYPD with 1920 cores and NDTE=240, which is reasonable given the high resolution of the grid. The utilization of larger computational facilities for TS015 and TS005 remains an important direction of future work.

Sea ice kinematics and scaling on representative days
In order to study the sea ice kinematics, we output two years' daily mean sea ice fields for all three TS grids, covering year 41 to 42 as in Figure 4.a. The deformation rates in their invariant form are defined in Equation 1 through Equation 3, including: shearing rate (ε shear ), divergence rate (ε div ) and total deformation rate (ε total ). They are computed from the daily mean prognostic sea ice drift speeds (u's and v's) defined on Eulerian grid locations (with Arakawa-B staggering in CICE).

TS015 TS005
Year 0 We have manually chosen two typical days that are representative of Arctic sea ice drift patterns: Dec. 20th and Feb. 6th % of the total variance, with the normalized spatial pattern (unitless) and the principle component (time series) shown in Figure   S1. The wintertime AO index of the NYF dataset is shown in Figure 6.c. The winter of the NYF forcing dataset corresponds to 10 https://doi.org/10.5194/gmd-2020-160 Preprint. Discussion started: 2 June 2020 c Author(s) 2020. CC BY 4.0 License.
an overall neutral AO status, and the variability of wintertime weekly AO index (25 hPa as in Fig. 6.c) is also on par with the average intra-seasonal variability of the 50-year NCEP reanalysis data (22 hPa as in Fig. S1.b).
The two representative days are: (1) Dec. 20th, on which date the high-pressure center resides in Beaufort Sea, and a negative AO index is witnessed, and (2) Feb. 6th, on which date the high-pressure center shifts towards the eastern hemisphere, the lowpressure system in the Atlantic sector extends further into the Arctic, and a slightly positive AO index is present. Furthermore, 5 because asymptotic convergence of sea ice kinematic to increase in the EVP subcycle count is witnessed in existing studies (Lemieux et al., 2012;Koldunov et al., 2019), we limit the analysis of kinematics and scaling to the experiments with the largest EVP sybcycle count (NDTE=960) for each grid.  For comparison, the run with TS015 produces much fewer features, while that with TS045 only produces the major 1 or 2 deformation features. However, on both days, there is also evidence that, even at TS005, the model is limited for resolving fine-scale kinematic structures. Some large kinematic features, especially in the central Arctic, have 30 ends that are not well defined, which is clearly not bounded by the grid sizes. This applies to both shearing and divergence, and it indicates that the effective resolution of the model needs to be evaluated against the grid's native resolution, in terms of resolving kinematic features.
We examine whether the distribution of total deformation rates follow power-law distributions, as previously shown with RGPS data in Marsan et al. (2004). Specifically, the region outlined in Figure 7 is used for study, and we show the daily and 35 11 https://doi.org/10.5194/gmd-2020-160 Preprint. Discussion started: 2 June 2020 c Author(s) 2020. CC BY 4.0 License.
3-day complementary cumulative distribution function (C-CDF) for total deformation rates in Figure 9. The 3-day deformation data is computed from the 3-day mean velocity fields for the period of Dec. 19th to Dec. 21th, and Feb. 5th to Feb. 7th, respectively. For both daily and 3-day cumulative distribution, we attain multi-fractality across all the three resolutions. With higher resolution, the model consistently simulates more extreme deformation events, and a flatter C-CDF for the deformation rates.

5
Since the three grids have nearly exact grid stepping ratio of 1:3:9, we carry out: (1) the spatial scaling of TS005 onto TS015 and TS045, and (2) that of TS015 onto TS045. At larger spatial scales, the slopes become steeper for both TS005 and TS015 on Dec. 20th, but remain unchanged on Feb. 6th. The C-CDFs for each resolution, including modeled with the native resolution as well as scaled from higher resolutions, are shown in each panel in Figure 9 with same symbol. The slopes of C-CDFs from scaled rates are shallower than non-scaled rates, indicating that for each grid, the effective resolution is lower 10 than the grid's native resolution. For example, at 0.45 • , the slopes from TS005 are shallower than those from TS015 and TS045 (-1.6, -2.1 and -2.7 for Dec. 19th to Dec. 21th). This indicates that at the scale of 22 km, more extreme deformation events are present in the run with TS005 than TS015, and in turn, TS045. The main cause may be due to that the effective resolution of TS015 is lower than 0.45 • (i.e., TS045), in terms of simulating realistic shape of the power-law distribution for deformation rates. We further scale down the results of TS005 and TS015 to attain the slopes at spatial scales coarser than 0.45 • . By using 15 the results of TS005 as the standard, we attain same slope with TS015 as TS005 at the scale of 42 km and 50 km for Dec.
20th and Feb. 6th, respectively. This result gives certain hint of the effective resolution of TS015 (with grid resolution of about 7.3 km). Regarding the CDF of sea ice kinematics, the effective resolution is variable among days with different kinematic features, and about 6 to 7 times the native resolution of TS015. Also, since we witness flatter C-CDF slopes in scaled datasets, we expect the "real" tail of C-CDF at 0.05 • is flatter than the modeling result from TS005, which have the slope of -0.9 for 20 3-day fields around Dec. 20th, and -0.6 around Feb. 6th.
As reported in Marsan et al. (2004), the slopes for the tails of the C-CDF for 3-day revisit interval are -2.5 and -3.6 for the nominal spatial scale of 16 km and 280 km, respectively. The data analyzed in Marsan et al. (2004) are around Nov. 6th, 1997, which corresponds to a negative AO index of similar magnitude to Dec. 20th in our experiments. Since 0.45 • in our experiment corresponds to the spatial scale of 22 km, there are slightly flatter C-CDF for the 3-day deformation fields as modeled with 25 TS045 (-2.7), TS015 (-2.1) and TS005 (-1.6) than Marsan et al. (2004).
We further carry out systematic analysis for the spatial scaling for these two representative days (Fig. 10). Specifically, the moments (q) that are adopted to evaluate structure function of scaling and the multi-fractality are: 0.5, 1, 2 and 3. The polynomial for the least-square fitting of the structure function is: β(q) = a · q 2 + b · q. At q = 1, the three resolutions show slight differences for β: for both days, and higher resolution starts with a slightly higher mean deformation rate (0.15 and 30 0.12 for TS005 on Dec. 20th and Feb. 6th). At q = 3, the structure function is in the range of 0.9 (TS045) and 1.4 (TS005) on Dec. 20th, and about 1.3 on Feb. 6th for all three grids. At higher orders (e.g., q = 2 or q = 3), high-resolution produces faster decay of deformation rates with spatial scales (larger β) on Dec. 20th, but on Feb. 6th, the differences of decay rates are less pronounced between the grids. At 3-day scale, the average deformation rate is generally lower than at 1-day scale, and the difference on higher orders is also evident between Dec. 20th and Feb. 6th. Furthermore, no positive value of β is detected at q = 0.5, which is consistent with Marsan et al. (2004) (Fig. 4 of the reference).
The structure function of β(q) shows a strong curvature, and this applies to all the grids on both days, as well as both temporal scales (1-day and 3-day). This further indicates that the model simulates multi-fractal sea ice kinematics. On Dec.
20th, the curvature level (a) is higher in TS005 (0.15) than TS015 (0.1) and TS045 (0.08). But on Feb. 6th, the curvature 5 levels are more consistent across the 3 resolutions (0.18 for TS045, 0.2 for TS015, and 0.21 for TS005). The differences in the statistics of scaling is the result of both the localization of the specific deformation field, as well as the resolution of the grid. Besides, there is slight drop in β when evaluating 3-day deformation fields with 1-day counterparts. Temporal scaling is indicated by the decrease in β for the daily field and 3-day field for Dec. 20th, and not evident for Feb. 6th.
Based on the numerical experiments with NYF dataset, we show that the sea ice kinematics, including the deformation 10 and its scaling, are distinctive on different days. On the two representative days, the deformation fields show multi-fractality, which is consistent with existing studies with observational datasets and modeling results. On Feb. 6th, there are more larger deformation events than Dec. 20th: higher 95-th percentile is present for all 3 resolutions, including scaled results (Fig. 9).
Also the structure function for spatial scaling is more convex on Feb. 6th and Dec. 20th, indicating less dominant large-scale features on Feb. 6th. Furthermore, there is more effective temporal scaling on Dec. 20th than Feb. 6th, as shown for C-CDFs in 15 Figure 9 and structure functions in Figure 10. Large-scale sea ice drift patterns were found to be greatly accurately determined by geostropic winds and the associated SLP field and AO indices (Rigor et al., 2002). The associated sea ice deformations at smaller spatial and temporal scales, due to the multi-fractality and large scale-small scale linkage, is highly dependent on the atmospheric forcings which contain inherent variability at different scales. Furthermore, the sea ice deformation is also greatly dependent on sea ice status, such as ice thickness and strength. Since the experiments in this study target at climatologlical sea 20 ice states, the ice thickness and multi-year ice coverage is higher than current satellite observations such as RGPS (Lindsay and Stern, 2003). Historical simulations which are driven by high-resolution, inter-annually changing atmospheric forcings, as well as coupling to dynamical atmospheric and oceanic models, are needed to compare with coinciding observations such as RGPS.

25
In this section we evaluate the modeled kinematics to the EVP subcycling and convergence. In Figure 11 we show the probability density (or PDF) of daily deformation rates during wintertime (Dec. to Feb.) for the three grids. All the simulations attain a good shape for the tail of of the PDF, approaching the slope of -3. For total deformation rate, there exists a well defined mode at 0.1%/d to 0.2%/d for runs with NDTE=960, and at NDTE=120, a slight shift of mode to higher values (between 0.5%/d to 1%/d for TS015 and TS005). At different resolutions, the EVP subcycling count plays a different role in the shape of PDF.

30
For TS045, the PDF is consistent between different subcycle counts. But for TS015 and TS005, there are much more evident differences: (1) with larger NDTE, there are more regions with smaller shearing and divergence (less than 0.1 %/d); (2) there is general convergent behavior of the shape of the PDF when NDTE is large (e.g., NDTE>480 for TS015). This behavior also 13 https://doi.org/10.5194/gmd-2020-160 Preprint. Discussion started: 2 June 2020 c Author(s) 2020. CC BY 4.0 License.

Dec-20
Feb-6 (a) (b) (c) applies for specific days (see Fig. S2 and Fig. S3 for the PDF on Dec. 20th and Feb. 6th, respectively). This indicates that insufficient EVP subcycle count leads to overall non-convergent deformation rate distributions.
A visual inspection of the deformation fields reveals the loss of the kinematic features when the convergence is not attained. Figure 7 shows the daily total deformation fields on Dec. 20th as simulated with different NDTE values for each grid. The more remarkable difference is between the runs with NDTE=120 and NDTE=960 with the highest resolution (i.e., TS005).

5
Although the linear kinematic features are well defined with NDTE=960, the deformation field is much noisier for the run with NDTE=120, with only the larger features detectable. The noise level is at about 1 %/d or lower, which corresponds to the mode of PDF in Figure 12. With larger values of NDTE, the noise level decreases and drops lower than the physical deformation rate, and a convergent PDF and linear feature map is attained. The results show that with the traditional EVP implementation, with higher resolution, even more subcycles are needed to reach convergence for the kinematics. This causes further deterioration 10 of simulation speed, given that the dynamics time step is already decreased in high-resolution runs, due to CFL conditions in processes such as advection (Tab. 2).
For TS045, although the PDF of the deformation rates are similar among runs with different NDTE values, the deformation field is also slightly more noisy in runs with NDTE=120 and NDTE=240 than that with NDTE=960. The region with the biggest difference is within the Canadian Arctic Archipelago. We further evaluate the effect of subcycling by carrying out 15 the experiments with NDTE=240 and NDTE=960 for TS045 further to year 50 (Fig. 4). We show the difference of March (September) sea ice thickness between NDTE=960 run and NDTE=240 run on year 50 in Figure 13.a (13.b). Although there is only less than 1% difference in the basin-scale ice volume between the two experiments, the sea ice is remarkably thinner   Figure 9. Cumulative distribution functions for the total deformation rates on Dec. 20th and Feb. 6th. Region of study is outlined in Fig. 7 and Fig. 8. Daily rates are shown in the first row and 3-day rates on the second row. Slopes are computed for the range of the cumulated probability between 0.05 and 0.25 for all C-CDF's, which correspond to 95th and 75th percentile of the deformation fields, respectively.
in Canadian Arctic Archipelago (CAA) in the run with NDTE=960. The difference in ice thickness is all-year around (i.e., in both September and March). Specifically, in the run of NDTE=960, ice arch forms by the eastern end of the Lancaster Sound and near Amund Ringnes Island, resulting in thicker ice in these regions. In other parts of CAA, the ice is thinner by about 15 to 20 cm on average, and up to 1 meters thinner in certain areas. As shown in Figure 13.c, the run with NDTE=240, which started from year 1, has already attained equilibrium in sea ice thickness and volume in CAA well before year 41. 5 However, in the run with NDTE=960 which starts from year 41, gradually shifts to another equilibrium status for ice thickness towards year 50. The overall volume difference is uniformly 100 km 3 in March and September, which consists of about 5% of the total volume in September. We further attribute the difference of winter ice thickness to the thermodynamic (Fig. 13.d), advection ( Fig. 13.e), and dynamical ridging/rafting (Fig. 13.f) contributions during the freeze-up seasons, computed as multi- year mean fields for December, January and February (DJF) of year 41 and 50. As compared with the run with NDTE=960, 10 there exist small deformation events in the CAA in the run with NDTE=240, which are arguably due to non-convergent EVP solutions (see also Fig. 12 for comparisons of deformation fields). Therefore, there is more ice thickness increase in the run with NDTE=240 due to these events (Fig. 13.f). Also, since September ice is thinner in the run with NDTE=960, it is intuitive that the winter ice growth should be higher. On the contrary, the thermodynamic ice thickness growth is in general lower in  Dec-19 ~ Dec-21 Feb-5 ~ Feb-7 Figure 10. Spatial scaling of total deformation rate on 1-day and 3-day scale. The daily mean velocity fields on Dec. 20th and Feb. 6th, as well as the 3-day mean velocity fields centering on Dec. 20th and Feb. 6th are used to compute the deformation rates and scaling curves.
Each panel contains the scaling curves and the structure function relating the scaling coefficients (β's) to the moment of order (q).
the run with NDTE=960, except for regions with the biggest thickness decrease in September as compared with the run with NDTE=240 ( Fig. 13.d and Fig. 13.b). This is mainly due to that the thermodynamic ice growth is also closely tied to the kinematics processes. With more noisy sea ice movement fields in the run with NDTE=240, the sea ice formation and growth is also promoted with very small deformation events, resulting in more thermodynamic growth in the run with NDTE=240. As compared with ridging and rafting, the sea ice advection is mainly responsible for redistributing the ice mass and thus plays a 5 minor role for the overall ice volume in CAA in our experiments (Fig. 13.e). Since we adopt a climatological forcing for our experiments, the sea ice in CAA is mostly fast ice during winter.
The dependence of the modeled sea ice thickness in CAA on EVP convergence in our study is purely numerical, but calls for further attention from both modelers and modeling data users. It highlights the importance of numerical convergence of EVP (or any other candidate solvers to VP) on the modeled sea ice climatology in fast ice regions. Since the experiments are 10 idealized in this study, the effects in realistic simulations may be also subjected to grid resolution and staggering, as well as coupling and feedback processes.

TS045
TS015 TS005 NDTE=120 NDTE=240 NDTE=480 NDTE=960 Figure 11. Probability density of modeled daily deformation rates during winter months (Dec. to Feb.) with TS045 (first row), TS015 (second row) and TS005 (third row). The region of study is outlined in Fig. 7. Shearing rate (left column), divergence rate (central column) and total deformation rate (right column) are shown. Runs with different EVP configurations (NDTE from 120 to 960) are marked by the same color as in Fig. 4.a.
In this paper we introduced sea ice modeling with a multi-resolution framework in Community Earth System Model. A grid hierarchy is constructed with the resolution range spanning climate simulations (0.45 • ) to sub-mesoscale modeling (0.05 • ).
At 0.05 • , the grid stepping in the Arctic region is approximately 2.45 km. The grid hierarchy is incorporated into CESM, and by using atmospherically forced experiments, we simulate and evaluate multi-scale sea ice kinematics and scaling properties. 5 We have found good consistency of the Arctic sea ice climatology and kinematics across the resolution range. Multi-fractal sea ice deformation is accurately modeled by all three resolutions, with good agreement with observational works in terms of scaling properties. In our study, high-resolution (0.05 • ) runs yield the most trustworthy kinematic features, and the multiresolution simulations provide a unique approach for evaluating sea ice kinematics in lower resolutions. Furthermore, the convergence of Elastic-Viscous-Plastic rheology model is evaluated, which show significant impact of EVP convergence on 10 kinematic statistics, as well as landfast ice in Canadian Arctic Archipelago. The framework of utilizing the grid hierarchy, including TS045, TS015 and TS005, provides an infrastructure for multi-scale modeling studies for both ocean and sea ice in the future. The three grids and the model integration are openly provided for public use for the version 1.2.1 of CESM.
Sea ice kinematics has been the focus of the high-resolution sea ice modeling community in recent years. Based on SAR remote sensing such as ASAR and RADARSAT, kilometer-level or finer observations of sea ice drift and deformation are 15 made possible, and serve as the backbone for the validation of high-resolution sea ice simulations. Modeling studies with structured grids such as MITGcm (Spreen et al., 2017;Hutter et al., 2018) are similar to our study, although model specifics are different. The sea ice deformation and scaling can be easily derived with model output for spatial scaling, but Lagrangian based diagnostics may be needed for the full analysis of temporal scaling. In our study, we mainly utilized model output for the study of spatial scaling properties in Section 3.2, with initial study with temporal scaling analysis with 3-day mean drift 20 fields. Specifically, in our study, the process dependent scaling properties were found in our study for representative days for sea ice drift and AO index. As future work, a full exploration of AO and ice condition dependent scaling analysis is planned with the grid hierarchy and CESM. Besides, the scaling analysis remains a traditional tool for evaluating sea ice kinematics, but it maybe insufficient for fully evaluating the sea ice deformation properties. Novel statistics based on linear kinematics features (LKF) are proposed in recent studies, such as Hutter and Losch (2020) and Ringeisen et al. (2019). The utilization of 25 a full suite of diagnostics in our modeling framework, including temporal-spatial scaling and LKF based approaches, serve as an important direction for future work.
Unstructured grid based models and Lagrangian models have unique capabilities in modeling sea ice. Regionally focused studies can be easily carried out with variable grid sizes, such as the Arctic simulation with FESOM (Wang et al., 2018;Koldunov et al., 2019). Purely Lagrangian models such as neXtSIM (Rampal et al., 2019) are potentially free of the resolution 30 issues for resolving small deformation features. Specifically, neXtSIM utilizes elasto-brittle rheology (Girard et al., 2011), which are shown to simulate reasonable sea ice kinematics and scaling properties even at moderate resolution settings (Rampal et al., 2019). Although traditional VP rheology can produce deformation features that are consistent with observations at resolution of 1 km (Hutter and Losch, 2020), there exists efforts to improve certain aspects such as reproducing observed anisotropy in deformation fields as well repeating deformation features (Tsamados et al., 2013). Regarding rheology models, we plan to carry out simulations and comparative studies with anisotropic rheology models such as EAP, with updated versions of the component models or the coupled model.
While EVP provides a numerically stable and easy-to-implement solver for the traditional VP model, the convergence of EVP solutions to VP model is the focus of many recent efforts. Based on a traditional implementation of EVP in our model, 5 we have witnessed asymptotic and converging behavior of modeled kinematics by the increasing EVP subcycle count. But this comes at a large computational overhead: at 960 subcycles per step for TS005, the simulation speed has halved compared with 240 subcycles per step, and over 50% of the computational time is consumed in EVP, with less than 0.5 SYPD with 1920 processor cores. Ideally, an efficient and scalable implicit solver promises a more elegant and numerically sound solution to the solving problem of VP model (Lemieux et al., 2010). Also, adaptive methods that complements the convergence and efficiency 10 problems of traditional EVP solver, such as Kimmritz et al. (2015) and Koldunov et al. (2019), are considered in our future work for the integration of the upgraded version of the coupled or sea ice model.
In Sec. 3.3, we have shown that the mean states for ice thickness and volume can be systematically shifted due to the numerical behavior of EVP solver. Since this issue is purely numerical, the uncertainty caused by it is different from other factors, such as the choice of ice strength parameterization scheme. In our study, the region which is most sensitive is landfast 15 ice in CAA, where significant decrease of ice volume is witnessed when solver convergence is attained with enough subcycle counts. Furthermore, more subcycling is required for higher resolution to reach convergence. Given the wide adoption of VP and EVP in climate models, using the model outputs for climate research and applications (such as the projection of shipping routes) could face potential compromises if the convergence issue is overlooked. Especially, careful choices should be made for configuring EVP at different resolutions (Kiss et al., 2020). In Koldunov et al. (2019) the authors also discovered sensitivity of 20 modeled ice thickness to EVP subcycling, but an increase in ice thickness with more EVP subcycles (Fig. 2 of the reference).
Also the region with most significant thickness change is different in Koldunov et al. (2019), covering both regions in CAA and north of CAA and Greenland. The different behavior to EVP subcycling may be due to the differences in the numerical experiments, as well as model physics, including grid resolutions, ice thickness distribution settings, etc. Although both show relationship of ice thickness to rheology model, more analysis is needed for the attribution and explanation of aforementioned 25 differences. Besides, tidal processes and interactions with sea ice are potentially important for the simulation of landfast ice (Lemieux et al., 2018). Since these processes are absent in our current model, we plan to include parameterization schemes that account for the their influence on sea ice kinematic in the future.
With the wide availability of high-performance computing utilities and the progress in model developments, multi-scale simulations is becoming more common for climate modeling community. While 1 • models are still dominant in Coupled resolutions are adopted, including OM4p5 (0.5 • ) and OM4p25 (0.25 • ). Parameterization schemes in the ocean and the sea ice models are chosen and tuned to each specific resolution. For example, mesoscale eddy induced mixing parameterization is 35 usually adopted for low-resolution ocean models (0.5 • or coarser), but inactive for higher resolutions. In our study, we use three resolutions for the study of Arctic sea ice kinematics, including: 0.45 • , 0.15 • and 0.05 • . As is shown in the scaling analysis, this multi-resolution framework enables comparative analysis across the resolution. Given that the modeled sea ice climatology is reasonable and consistent among the three resolutions, we consider the results adequate for the analysis of kinematics and scaling based on the coupling to the Slab Ocean Model. Beyond the lower computational overhead of this approach, we can 5 also attain an equilibrium status for the sea ice with fewer model years. Especially for 0.05 • grid (TS005), the computational cost and the duration is prohibitively high to fully spin-up the ocean-sea ice coupled model. Furthermore, it reduces the uncertainty of ocean's modeling on the sea ice, including: (1) ocean model's parameterization schemes that are not aligned and potentially not well-tuned between the different resolutions; and (2) the avoidance of ocean and ocean-sea ice coupled internal variability that may potentially compromise the comparability across the resolutions. In the future, we plan to carry out spin-10 up and long-term simulation of the ocean-sea ice coupled system using the proposed multi-resolution framework, which will introduce variability to the modeled sea ice states from the oceanic dynamic processes and coupled processes. Specifically, the migration of the status from cheap, low-resolution runs (TS045) to TS015 and TS005 is planned, in order to facilitate faster spin-up for high-resolution runs. This is accompanied by historical, inter-annual forcing datasets, for the production of realistic simulations and comparison with coinciding satellite observations of the sea ice kinematics. innermost ellipse is achieved, by controlling the semi-major axes (a's) and semi-minor axes (b's) of the ellipses. Without the loss of generality, we rescale the projected NP, so that the outermost ellipse is the unit circle and the innermost ellipse resides on the x-axis, and the layout is shown in Figure A1.a. Therefore, we have: (1) for the outermost ellipse, a = b = 1; and (2) for the innermost ellipse, a = α and b = 0, where α is half of the distance between the two grid poles. To ensure a continued change of meridional grid scales on NP-SP boundary, the change in b should be equal to that in a on the outermost ellipse. In Second, with the embedded ellipses, the orthogonal grid can be constructed from the outermost circle. Since NP is directly linked to SP at latitude of φ, we specify the grid points on this boundary (i.e., the outermost circle), and extend into NP to form the meridional grid lines. The construction process is iterative, starting from these points on the outermost ellipse, and for each step constructing lines between two adjacent ellipses. In each step, the ending locations of the lines are located on the inner 20 ellipse, under orthogonality constraints. Fig. A1.a shows a specific example between two adjacent ellipses (marked by blue and dash-blue). This whole approach is similar to the grid generation methods that are based on numerical integration processes, as in Madec and Imbard (1996) and Xu et al. (2015). Figure A2 show the meridional and zonal grid scales (dx and dy) along typical meridional grid lines for TS005. As is shown, there is smooth transition of the meridional grid scale on the boundary of SP and NP (at about J=2520 for TS005). Also the 25 overall grid scale anisotropy is kept lower than 1.5 in the oceanic areas.  Figure 2 with the same color coding. Note that when I=1, dy (meridional grid size) approaches zero when J is close to the northern end (J=5400), which does not affect the model's time stepping since it is near the grid pole which resides on the land.