### Data

Table of Contents

Sea ice draft data were obtained from upward looking sonar (ULS) moored in the East Greenland Current in the western Fram Strait. The dataset continuously covers the last three decades (1990–2019) with some short temporal gaps. Four ULSs were zonally aligned at approximately 79° N from 3° W to 6.5° W (Fig. 1a). The latitude of the mooring array changed from 79° N to 78.8° N in 2001. The zonal positions (names) of the moorings equipped with ULS are 3° W (F11), 4° W (F12), 5° W (F13) and 6.5° W (F14), respectively. There are three main temporal data gaps during the three decades of measurements, that is, 1996, 2002 and 2008. Except for these gaps, ULSs were in operation although their number varied from time to time. The ULS measures the travel time of the sound reflected at the bottom of the floating sea ice, from which we calculate the ice draft, the underwater fraction of sea ice^{54}. The raw data were processed to ice draft using procedures described in earlier literature^{55,56}. The accuracy of each draft measurement ranges from 0.1 m (ice profiling sonars (IPS) deployed after 2006) to 0.2 m (ES300 instruments before 2006), while the uncertainty of each individual measurement is not subject to bias errors and the summary error statistics of monthly values are less than 0.1 m^{57}.

The daily mean sea ice motion product provided by the National Snow and Ice Data Center (Polar Pathfinder Daily 25 km EASE-Grid Sea Ice Motion Vectors v.4.1, hereafter NSIDCv4 ice drift) was applied for backward trajectory calculation of sea ice floes, estimating residence time of sea ice in the Arctic and analysis of spatial average of ice drift speed. The product was derived from a combination of various ice motion estimates from remote sensing and tracks of ice-tethered buoys^{52}. We applied the motion vectors from 1984 to 2019 for the backward trajectory calculation. Sea ice concentration data were taken from the Global Sea Ice Concentration Climate Data Records of the EUMETSAT OSI SAF^{51} (OSI-409 v.1.2 for the period from January 1984 to April 2015 and OSI-430 v.1.2 for the period after April 2015). The ice concentration data were used for the sea ice trajectory calculations and to analyse spatial and temporal changes of ice concentration.

Sea surface temperatures obtained from the NOAA/NESDIS/NCEI Daily Optimum Interpolation Sea Surface Temperature v.2.1 (ref. ^{58}) were used to examine temporal variation in sea surface temperature in areas of ice formation (Extended Data Fig. 1a,b). A dataset of Arctic Ocean in situ hydrographic observations^{59} was also used to examine the change in ocean surface temperature between 1990 and 2006 and between 2007 and 2019. The dataset was revised using recent observations^{60} and gridded to 110 × 110 km cells covering the whole Arctic Ocean. The mean temperature in each period (1990–2006 and 2007–2019) in each cell was defined by an average of all available measurements for every 3-month period (January–March, April–June, July–September and October–December). If the number of available data in a cell was less than four, a missing value was assigned to the cell. The difference in summer sea surface temperature (July–September, 0–20 m) between two periods (1990–2006 and 2007–2019) was used for the analysis (Extended Data Fig. 1c). Atmospheric data were taken from the European Center for Medium-Range Weather Forecasts Reanalysis v.5 (ref. ^{61}) (ERA5). Daily and monthly mean sea-level pressure and 10 m wind were used in the analyses.

### Ice thickness distribution

Ice thickness distributions were derived on a monthly basis. All sea ice draft measurements in the Fram Strait from 1990 to 2019 were classified into draft thickness bins of 0.1 m, ranging from 0 to 8 m (80 bins in total). The number of data samples (ice draft measurements) used to derive the distributions varied from time to time. From 1990 to 2005, *O*(10^{4}) samples were used to derive the distribution functions on a monthly basis (measurement interval of 240 s in most cases); after 2006, *O*(10^{6}) samples were used (interval of 2 s). The number of samples used were sufficient to derive thickness distribution on a monthly basis^{57}. The number of data in each bin were divided by the total number of measurements to derive the distribution function. The open water fraction (that is, zero ice thickness) was excluded when deriving the function. Distribution functions, including the open water fraction, are shown in Extended Data Fig. 3. If the temporal coverage of the data samples was less than 15% of a monthly coverage, the distribution function was not defined and removed from the analysis. The draft thickness distributions were converted to ice thickness distributions by an average ratio of draft to thickness in the Fram Strait, 1.136 (ref. ^{62}). Although the ratio has some seasonal variability and might have slightly changed due to changes in ice density and snow load in different seasons and years, we assumed that the change was not considerable for the aim of the current analyses. A composite time series of ice thickness distribution (Fig. 2a and Extended Data Fig. 3) was obtained from a combination of available distribution functions from F11 to F14 . More specifically, one of the distribution functions from F11 to F14 was used with a priority order of F13, F14, F12, F11. The composite time series is shown in Fig. 2a, while the time series at each site are shown in Extended Data Fig. 4. The fractions of sea ice thicker than two thresholds, 4 and 5 m, were calculated as the cumulative function of all available distributions (that is, all distributions from F11 to F14) in each month and are shown in Fig. 2d.

The uncertainties on the estimates of monthly fractions of ice thicker than a threshold and position of the modal ice peak were assessed numerically using a moving block bootstrap approach^{63}. Bootstrapping is a family of resampling techniques used to derive uncertainties on various complex estimators for large datasets and uses random sampling with replacement^{64}. The presence of autocorrelation in ice draft series from IPS (ULSs deployed after 2006) suggested using the moving block bootstrap approach^{63}. The method splits the original monthly series of ice drafts *O*(10^{6}) samples of length *N* into *N* *−* *K* + 1 overlapping blocks of length *K* each. The block length was set to 30 samples that approximately corresponded to a distance of 10 m covered by ice travelling at 0.3 knots. It roughly corresponded to the lower limit of a horizontal spatial scale of ice ridges. For ES300 instruments (ULSs deployed until 2005) with a lower sampling rate of 240 s, ordinary bootstrapping was used.

At each of *M* steps of bootstrap sampling, *N*/*K* blocks were drawn at random, with replacement, from the constructed set of *N* *−* *K* + 1 blocks, making a new bootstrap data sample for the month. A Gaussian noise of *N*(0,1) was further added to the data to account for measurement uncertainty. The mean and s.d. of the fractions of thicker ice and modal ice peak position were then calculated directly from the *M* estimates derived at each step of the procedure. The results suggested that the monthly coefficient of variation or a ratio of the s.d. of the estimate to its mean, for the IPS data varies from 1% to 3% for both fractions, being on average lower (1–2%) for a fraction of ice greater than 4 m thick. For ES300, the coefficient of variation was slightly higher at about 4(6)% for the fractions of ice thicker than 4(5) m. The same applied to the position of the modal peak, which showed a coefficient of variation of 0–3% being typically closer to 0 for the IPS data. It suggested that the selected bin width was large enough to accommodate uncertainties related to the approach and data. For the ES300 data, the monthly s.d. of the modal peak position was higher, up to 30 cm, and the coefficient of variation was within 9%. Therefore, we postulate that the inferred uncertainties are far too low to have any noticeable influence on the results of the shift detection analysis.

### Modal peak and variance of ice thickness distributions

As statistical measures of ice thickness distributions, we examined modal thickness and variance. The monthly mean ice thickness distributions (740 samples in total) were fitted to log-normal functions:

$$F(x)=\frac{1}{\sqrt{2{\rm{\pi }}{\sigma }_{{\rm{f}}}^{2}}x}\exp \left(\frac{-{({\rm{l}}{\rm{n}}x-{\mu }_{{\rm{f}}})}^{2}}{2{\sigma }_{{\rm{f}}}^{2}}\right)$$

(4)

where *σ*_{f} and *μ*_{f} are the fitting parameters, *x* is the ice thickness bin and *F* is the distribution function. To detect the second peak of the distributions that represents multi-year ice travelled across the Arctic Basin, a cut-off threshold was introduced. The threshold was used to exclude thin sea ice fraction, which is supposedly formed in the vicinity of the Fram Strait and is not representing basin-wide changes of ice properties in the Arctic. We defined the threshold by the minimum between the first and second peak of each monthly mean distribution. A set of two consecutive negative gradients (towards thicker bins) followed by two consecutive positive gradients was used to detect the minimum (after applying 3-bin smoothing), while a threshold of 1.53 m (corresponding to 1.3 m of ice draft) was applied when an estimated threshold was thicker than 3 m. Function values ranging lower than the threshold were set to zero (zero case) or excluded from the fitting (NaN case). A least-square minimization was applied to fit a log-normal function to the distribution. In general, the fitted log-normal functions represent the distribution very well. The NaN case slightly underestimates the modal peak, while the zero case captures the peak very well. Examples of distribution functions, together with cut-off thresholds and fitted log-normal functions, are shown in Extended Data Fig. 5. The modal peak of the log-normal function roughly gives the thickness of thermodynamically grown sea ice, while the variance of the function quantifies the deformed fraction of sea ice (dynamically thickened thickness). Changes in modal peak height and variance of the fitted log-normal functions, var(*x*) = exp (2*μ*_{f} + *σ*_{f}^{2}) (exp (*σ*_{f}^{2}) − 1), for the last three decades are summarized in Fig. 2b,c. The time series of modal thickness and the fitting parameters *σ*_{f} and *μ*_{f} are summarized in Extended Data Fig. 6.

### Regime shifts detection

A sequential algorithm for regime shift detection^{23} was applied to all time series. The method identifies discontinuities in a time series using a data-driven approach that does not require an a priori assumption on the timing of the regime shifts. The method first identifies potential change points sequentially by checking if the anomaly of the data point is statistically significant from the mean value of the current regime. If it is significant, the following data points are sequentially used to assess the confidence of the shift, using a regime shift index (RSI). RSI represents a cumulative sum of normalized deviations from the hypothetical mean level for the new regime, for which the difference from the mean level of the current regime is statistically significant according to a Student’s *t*-test. If the RSI is positive for all points sequentially within the specified cut-off length, the null hypothesis of a constant mean is rejected. This led us to conclude that the regime shift might have occurred at that point in time^{65}. If multiple data are available at a certain point in time (that is, multiple sites from F11 to F14), the mean value is applied in the time series. Before testing, the temporal gaps of the time series were interpolated by the average of all available data (modal peak height (Fig. 2b), variance (Fig. 2c) and modal thickness (Extended Data Fig. 6a) of ice thickness distributions, fraction of thick sea ice (Fig. 2d) and residence time of sea ice (Fig. 3a and Extended Data Fig. 2)). The cut-off length was set to 7 years (84 months) to cover the advection timescale (travelling time across the Arctic) of sea ice, while at the same time, detecting shifts occurring at a timescale shorter than a decade. Other cut-off lengths (3, 4, 5, 6, 8 and 10 years) were also tested to see the sensitivity and robustness of the results. A summary of the test results is given in Extended Data Table 1; the timing of the detected shifts is shown in all time series except for Fig. 2. The timing of the detected shifts of modal peak height and variance of ice thickness distributions are shown in Extended Data Fig. 6b,c, while those of the fraction of thick sea ice are shown in Extended Data Fig. 8. RSIs, respective *P* values and the shift of the means are summarized in Extended Data Table 2.

### Sea ice trajectory analysis

To investigate changes of pathways and residence time of sea ice in the Arctic basins, sea ice trajectories were calculated for the last three decades. Eight pseudo-ice floes were settled in the western half of the Fram Strait section (from prime meridian to 10° W) at the same time and advected backwards in time. The calculations started on the 15th of every month from 1990 to 2019. Daily sea ice motion vectors from the NSIDCv4 were used to update daily position of ice floes backwards in time. Ice motion vectors at the respective floe positions were calculated by interpolation of surrounding points with Gaussian-type weighting (e-folding scale of 25 km). Each trajectory calculation was performed 6 years back in time, while it was terminated if no motion vector was available within a 25-km distance or sea ice concentration at the floe position was lower than 15%. The sea ice concentration at the ice floe position was obtained from OSI-409/OSI-430 with a Gaussian-type weighting (e-folding scale of 12.5 km). The position of each trajectory termination was used to define the location of ‘initial sea ice formation’. Trajectories shorter than three months were excluded from the analysis because they represent ice floes formed in the vicinity of the Fram Strait.

Uncertainty of the daily position of the pseudo-floes was assessed by comparisons with ice-tethered buoy tracks obtained from the International Arctic Buoy Program (IABP)^{53}. We used 83 buoy tracks that arrived in the Fram Strait from 2000 to 2018 and calculated the corresponding pseudo-buoy tracks backwards in time. The comparisons showed that the mean error of the daily pseudo-buoy positions can be reasonably approximated by a linear function of backtracking days^{19}, error = 50 + (backtracking days)/2 km. We applied this empirical formula as an error of the daily position of the backward trajectories from 0 to 500 backtracking days, which corresponds to a 200 (300) km error after 300 (500) backtracking days. Note that this error estimate may underestimate the uncertainty because IABP buoy tracks have been included in the NSIDCv4 ice motion product. However, comparisons between non-IABP buoys and pseudo-buoy tracks derived from the NSIDCv4 with error estimates by a bootstrap method showed that pseudo-tracks are largely parallel to the corresponding buoys and the error does not monotonically increase over time^{66}. The estimated error circles (approximately 300 km) of ice formation location in the present study are sufficiently small compared to the polygons in Fig. 3b (greater than 1,500-km width), which guarantees the robustness of the analysis.

The residence time of sea ice in the Arctic basins was defined by the period from the start to the termination date of each trajectory. We calculated an average of residence time of eight pseudo-ice floes that arrived in the Fram Strait at the same time and used it to define the mean residence time of ice floes for each month. The uncertainty of the residence time was defined by the s.d. of the residence time of the eight pseudo-ice floes.

### Stochastic model of dynamic ice thickening

The log-normal form of the ice thickness distribution can be obtained from a simple proportionate growth process, \({X}_{m}={X}_{0}{\Pi }_{i=0}^{i=m-1}{a}_{i}\). If ln *X*_{0} and ln *a*_{0}, ln *a*_{1}, .., ln *a*_{i} are uncorrelated; thus, the probability function of *X*_{m} for large *m* is given by:

$$f({X}_{m})=\frac{1}{\sqrt{2{\rm{\pi }}\mathop{\sigma }\limits^{ \acute{} }{(m)}^{2}}{X}_{m}}\exp \left(\frac{-{({\rm{l}}{\rm{n}}({X}_{m}/\mathop{X}\limits^{ \acute{} }(m)))}^{2}}{2\mathop{\sigma }\limits^{ \acute{} }{(m)}^{2}}\right)$$

(5)

where \(\acute{X}(m)={e}^{\nu m}\) and \(\acute{\sigma }(m)=\sigma {m}^{1/2}\) and *ν* and *σ*^{2} are the mean and variance of the population distribution of ln *a*_{m} (including *X*_{0}), respectively^{67}.

In this study, we provide a concept and description of a stochastic model that formulates sea ice thickening associated with dynamic ice deformation. The model formulates three features of dynamic sea ice thickening by ridging and/or rafting: (1) dynamic ice thickening is a stochastic process (areal and thickening stochasticity); (2) thicker ice has a larger potential to get thicker than thinner ice at a dynamic event (proportionate ice thickening); and (3) thinner ice has a higher probability of dynamic deformation due to its weaker ice strength (preferential deformation of thinner over thicker ice types). The first point consists of two types of stochasticity in the dynamic ice thickening process. One is areal stochasticity, corresponding to the fact that ice deformation only occurs for a small fraction of the pack ice while the rest of the ice is unchanged when a dynamic event occurs. The other is thickening stochasticity, representing the fact that the thickness gain by ridging/rafting varies in space and differs between events. The second point represents a sea ice characteristic that thicker ice is tolerant and can exert stronger compressive force on the ice forming ridges and/or rafts; hence, more energy is potentially available for the dynamic thickening^{30,31}. The third point represents the fact that the thinner part of the pack ice is preferentially ridged/rafted when a dynamic event occurs^{68,69}. This also takes into account the effect of ice thickness changes on the dynamic thickening process, for example, the thinner ice condition in recent years increases the likelihood of ice deformation.

The two first points can be formulated by a proportionate thickening of stochastic ice thickness:

$${X}_{i}={a}_{i-1}{X}_{i-1}$$

(6)

where *X* is the stochastic ice thickness at a certain location, *i* denotes the time index counting sporadic dynamic events (for example, passage of a storm) and *a*_{i−1} is the conditional proportionate thickening increment due to the event at *i*−1. The increment *a*_{i−1} is a stochastic variable representing both areal and thickening stochasticity. They are implemented as:

$${a}_{i}=\left\{\begin{array}{cc}1+b{r}_{i} & {\rm{w}}{\rm{i}}{\rm{t}}{\rm{h}}\,{\alpha }_{i}({\rm{ \% }})\,{\rm{p}}{\rm{r}}{\rm{o}}{\rm{b}}{\rm{a}}{\rm{b}}{\rm{i}}{\rm{l}}{\rm{i}}{\rm{t}}{\rm{y}}\\ 1 & {\rm{w}}{\rm{i}}{\rm{t}}{\rm{h}}\,1-{\alpha }_{i}({\rm{ \% }})\,{\rm{p}}{\rm{r}}{\rm{o}}{\rm{b}}{\rm{a}}{\rm{b}}{\rm{i}}{\rm{l}}{\rm{i}}{\rm{t}}{\rm{y}}\end{array}\right.$$

(7)

where *b* is a proportionate thickening constant, *r*_{i} is a stochastic thickening increment that represents the thickening stochasticity of *i*-th dynamic event and *α*_{i} is the areal probability of dynamic thickening that represents areal stochasticity and gives the probability of ice thickening occurrence. This formula indicates that when a dynamic event occurs, *α*(%) area of pack ice experiences dynamic thickening (ridging/rafting), while the rest (1 − *α*(%)) is unchanged (areal stochasticity). The thickness gain, *br*_{i}, in the dynamic thickening area, is also a stochastic variable: the possible maximum gain is *b* while the minimum is 0 (*r*_{i} is random, so 0 ≤ *r*_{i} < 1) (thickening stochasticity).

We applied the proportionate thickening constant as *b* = 0.4. The choice of *b* = 0.4 implies that sea ice in the ridging/rafting area gains 0.4 m thickness at maximum (0.2 m on average) when the ice is initially 1-m thick, while it gains a 1.2 m thickness at maximum (0.6 m on average) when initially the ice is 3-m thick. The value of the parameter *b* comes from a recent high-resolution observation (5 × 5 m resolution covering a 9-km^{2} area) of single ice deformation event north of Svalbard, which describes changes in the sea ice freeboard just before and right after a storm event^{70}. According to this study, the change in sea ice freeboard in a converging area is 0.07 m (from 0.36 m to 0.43 m) on average, corresponding to 0.58 m dynamic thickening of ice by ridging and/or rafting (assuming the freeboard to thickness ratio = 8.35)^{62}. The gain relative to the mean ice thickness is estimated by *b* = 0.58/1.45 = 0.4 (the mean ice thickness in the survey area = 1.45 m). We applied the value for the proportionate thickening constant, *b*, as the first approach to develop the model. Although the study captured detailed spatial change in the sea ice freeboard, the estimate of *b* comes from a one-time event and hence needs further assessments by future observations. It should be noted that *b* = 0.4 is an areal average estimated from the observations, whereas we applied *b* as the upper bound of the proportionate thickening. This is because the probability function of the stochastic thickening increment, *r*_{i}, is not known so far; hence, we assumed a constant probability of the increment between 0 and *b*, potentially causing excessive thickening near the upper bound. The effect of the choice of *b* is discussed below.

The areal probability of dynamic thickening, *α*, is included to take the third point into account, that is, thinner ice has more chance to be ridged and/or rafted than thicker ice when a dynamic event occurs. To implement this feature, *α* is given by a function of ice thickness: the areal probability is inversely proportional to the stochastic ice thickness *X*_{i}:

$${\alpha }_{i}({X}_{i})=\frac{8}{{X}_{i}+1}({\rm{ \% }}).$$

(8)

The formula indicates that 1-m thick ice experiences dynamic thickening at 4% areal probability, while 3-m ice experiences dynamic thickening at 2% areal probability. Our first implementation of this formula is based on an observational estimate of areal fraction of dynamic thickening^{70}. According to the high-resolution survey of a single dynamic event, thickening occurred in 4% of the survey area with a mean ice thickness of 1.45 m. This formula also needs further evaluation by comparing with future observations that address the relationship between areal probability of deformation and ice thickness.

Another parameter necessary for the model is the number of dynamic events, *m*, that is, external forcing that could cause mechanical fracturing of sea ice and consequent ridging and/or rafting. We used the number of Arctic cyclones that passes over the ice pack as a first-order indicator of the number of dynamic events. Typically 90–130 cyclones per year occur in the Arctic Ocean (40–60 cyclones in winter, 50–70 cyclones in summer)^{71}. A typical size of an Arctic cyclone is approximately 3 × 10^{6 }km^{2} (mean radius of approximately 10^{3 }km)^{71}, which covers approximately one third of the ice-covered area of the Arctic Ocean. We therefore assumed that one-third of all cyclones hits the ice pack at a certain location in the Arctic, that is, the ice pack experiences approximately 40 dynamic events per year. This corresponds to approximately 80–240 dynamic deformation events for the typical residence time of sea ice in the Arctic (2–6 years; Fig. 3a).

In addition, examples shown in Fig. 5 contain a simple thermodynamic term to mimic the effect of modal peak shift of thickness distribution due to thermodynamic ice growth:

$${X}_{i}={a}_{i-1}{X}_{i-1}+c/{X}_{i-1}$$

(9)

where *c* is the thermodynamic ice growth coefficient. This term comes from a simplified thermodynamic process without thermal inertia of sea ice and heat flux from the ocean^{72}:

$$\frac{{\rm{d}}H}{{\rm{d}}t}=\frac{{\kappa }_{{\rm{ice}}}({T}_{{\rm{f}}}-{T}_{{\rm{s}}})}{{\rho }_{{\rm{ice}}}{L}_{{\rm{f}}}H}$$

(10)

where *H* is the ice thickness, *κ*_{ice} is the heat conductivity of ice, *ρ*_{ice} is the ice density, *L*_{f} is the latent heat of freezing, *T*_{f} is the freezing temperature of sea ice and *T*_{s} is the temperature at the ice surface. We applied this formula with a simplification, Δ*H* *=* *c/H*, where Δ*H* is the ice thickness change due to a thermodynamic process, *c* is a thermodynamic ice growth coefficient corresponding to Δ*t**κ*_{ice} (*T*_{f} − *T*_{s})/(*ρ*_{ice}*L*_{f}). As the model does not include a process that forms new thin ice by lead opening, inclusion of the thermal forcing term without a compensating term makes the modal peak very steep after few years, that is, no ice exists in thickness ranges thinner than the thermal equilibrium thickness. To alleviate such an excessive modal peak generation and to take into account the insulating effect of the snow pack that substantially delays thermodynamic ice growth, we applied a moderate value, *c* = 0.015, which is about one-third of the value estimated from *c* *=* Δ*t**κ*_{ice} (*T*_{f} − *T*_{s})/(*ρ*_{ice}*L*_{f}), where Δ*t* ≅ 9 d (corresponding to 40 dynamic events per year) and annual mean surface air temperature of *T*_{s} = 263 K.

The initial condition of the ice thickness distribution in Fig. 5 is given by a thermodynamically grown sea ice without dynamic deformation, *X*_{0}. This is also a stochastic variable, having a normal distribution for simplicity:

$$g({X}_{0})=\frac{1}{\sqrt{2{\rm{\pi }}}{\sigma }_{0}}\exp \left(\frac{-{(x-{\mu }_{0})}^{2}}{2{\sigma }_{0}^{2}}\right)$$

(11)

where *μ*_{0} = 1.0 and *σ*_{0} = 0.25 are applied in the examples (that is, 1 m mean ice thickness with 0.25 s.d., shown by *m* = 0 in Fig. 5), which roughly corresponds to the thickness of new ice three months after its formation (based on Anderson’s freezing degree days law^{25}, with an assumption of *T*_{s} = 253 K). Figure 5 shows examples of ice thickness distribution after 60, 120 and 180 dynamic events, roughly corresponding to 1.5, 3 and 4.5 years of residence time of sea ice.

The current formulation contains three parameters to describe the dynamic ice thickening process: *b* (the proportionate thickening constant); *α* (the areal probability of dynamic thickening); and *m* (the number of dynamic events). In this study, we briefly describe the sensitivity of the ice thickness distributions to these parameters. In general, a smaller (larger) thickening constant *b* decelerates (accelerates) the dynamic thickening process, that is, a smaller *b* gives a smaller variance and steeper modal peak of thickness distribution if *α* and *m* are fixed. However, a large value of *b* (for example, *b* = 0.8, indicating that ridged ice can be 1.8 times thicker than the ice before an event at maximum) makes the distribution bimodal because the possible thickness gain at each dynamic event is far from the modal thickness and the ridged/rafted ice tends to generate another peak apart from the mode. Therefore, the possible and realistic range of *b* should be examined further together with the probability density function of the thickening increment *r* by high-resolution observations in the future. The areal probability of dynamic thickening, *α*, also affects the evolution of the dynamic thickening process. A larger *α* promotes dynamic thickening because a larger fraction of pack ice can be deformed at one event. The thickness dependency of the probability, equation (8), decelerates further thickening of thick ice. Although values of *b* and *α* affect the progress of dynamic thickening in the model, we obtained similar ice thickness distributions with a log-normal form sooner or later, that is, smaller *b* and *α* can be compensated by a large *m*, the number of deformation events, indicating a robustness of the formulation. The resulting shape of the distribution, its temporal evolution (Fig. 5) and its comparison with the observed change in distribution (Fig. 1b) suggest that the proposed stochastic ice thickening model captures the essence of the dynamic thickening process that resulted in the observed changes in ice thickness distribution.