Size controls on the crossover from normal to self-inhibited sintering of ice spheres

Anomalies in the sintering of ice spheres were observed in previous studies as a growth of protrusions that evolve into porous bonds which eventually retreat. These anomalies were discussed in terms of Mullins-Sekerka theory which predicts the self-amplified growth of small perturbations. This interpretation would imply that the anomalous (mechanically weakening) sintering should crossover to normal (mechanically strengthening) sintering if the sphere size is reduced below a critical size. Here we show that this is indeed the case with sintering experiments in random sphere packs. Six samples, varying in particle diameter from 1 to 2.3 mm were sintered over ten days at −10 °C and then imaged by 3D microcomputed tomography. For two additional samples with small (1 mm) and large (2.1 mm) particles, we monitored the sintering evolution. Surface and bond properties were quantified by image analyses, which revealed solidly sintered bonds in small, and porous bonds in large spheres. Microstructure-based finite element simulations of the temperature field revealed higher temperature gradients between the particles in the bond region, as a driver for the protrusion growth. The temperature gradients, as well as the bond porosity, decreased after an initial increase up to five days and signal the retreat of the bonds. Our work details microstructural controls on anomalous sintering in ice as a self-inhibited bond growth process if only slight deviations from isothermal conditions occur. Besides the insight into fundamental sintering mechanisms, the results may have practical consequences for piste preparation, food processing, or research on icy planetary bodies. © 2021 The Author(s). Published by Elsevier Ltd on behalf of Acta Materialia Inc. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
Sintering is an important process for high-temperature materials that alters and dictates their properties through the growth of intergranular bonds. A "hot" material, in terms of the homologous temperature, is close to a phase transition where thermal energy is available to drive atomic or molecular mass fluxes, a precondition for sintering. Water-ice is a typical hot temperature material in ambient conditions and ubiquitous in our daily lives, for example as snow or in food, where sintering plays a central role.
During sintering a solid connection is formed between contacting ice crystals, often called neck or bond, by transporting mass to the contact region. Various mass transport mechanisms on molecular scales contribute to the sintering process, including vapor diffusion, surface and grain boundary diffusion, and plastic flow [1,2] . Mathematically, these processes can be described by a generalized * Corresponding author.
E-mail address: loewe@slf.ch (H. Löwe). power-law [1] , which characterizes the sintering mechanism by two integer exponents, n and m , as well as by the temperature and material dependent rate constant B . This powerlaw describes the bond growth in time t in the ratio of bond radius r b and particle radius r p , assuming a circular neck geometry. In various experimental sintering studies with spherical ice particles, necks that conform this growth law were observed and used to examine the transport mechanisms [3][4][5] . This law was also used to investigate sintering of snow, e.g. [6,7] .
However, deviations from this sintering theory have been reported during sintering of large spherical ice particles [8] as well as for non-spherical snow particles [9] . Both studies reported the evolution of needle like structures on the particle surfaces, referred to as protrusions, which formed a porous bond in the contact region. Moreover, they report from visual observations that some of the protrusions retreated after several ( ∼seven) days, if they did not reach the other particle's surface. In a more detailed investigation of the protrusions on 4 mm spheres [10] , the protru-sion length scale was related to surface perturbations for which Mullins-Sekerka theory predicts an instable interface. The theory of interface instability during diffusion or heat controlled crystal growth by Mullins and Sekerka describes the fate of a small perturbation on a crystal surface [11,12] . A perturbation that protrudes into a diffusion field is exposed to a higher concentration gradient than its surroundings and is therefore prone to growth. At the same time, capillary forces act on the perturbation, which tend to minimize the surface and reduce the surface energy. These competing mechanisms lead either to a stable interface when a perturbation decays, or to an unstable interface, on which a protrusion may grow. Accordingly, the authors [10] concluded that an unstable ice-air interface facilitates the growth of small perturbations on the ice crystal surface, which led to the evolution of the observed protrusions. A hitherto unmentioned but important corollary of this interpretation is that the transition to unstable neck growth should actually depend on the size of the spheres.
A central result of the Mullins-Sekerka theory for spherical particles is the link between the length scale of the perturbation and the underlying particle radius for unstable growth. For a fixed supersaturation, the instability only occurs for particle sizes larger than a critical radius [11] , R c = 14 k (c ∞ −c 0 ) /c 0 . 1 The critical particle radius R c depends on the capillary constant k , as well as the degree of saturation, which is described by the concentrations in the far field c ∞ and at the interface c 0 . However, the value Mullins and Sekerka obtain for R c , seven times the critical radius of nucleation theory ( = 2 k / (c ∞ − c 0 ) /c 0 ), may only be seen as lower limit, as for prominent instabilities orders of magnitude larger values for R c are required [11] . The radius quantifies the surface curvature H ( = 1 /R for spheres), which strongly affects the local diffusion field and as such crystal growth. Theoretical instability studies [13,14] highlight the important roles of temperature, saturation, and thermal conductance of the system, which has been detailed for a planar vapor-ice interface by Style and Worster [15] . Overall, these theoretical investigations are crucial for understanding the evolution of surface instabilities, but a prediction of their occurrence is complicated by their sensitivity to complex growth conditions in sphere assemblies [10] .
This dependence on particle size may explain the different outcomes of sphere sintering in the experiments described above, where solid necks were observed for sphere diameters of 0.02 to 0.7 mm [4,5] while porous bonds were found for 4 mm spheres [8] . Moreover, the weakly sintered ice spheres with a diameter of 2.1 mm and their peculiar mechanical behavior in the compression experiments [7] can be explained with porously sintered bonds. Overall, it can be expected that the critical particle size D c ( = 2 R c ) lies in the range between 0.7 and 4 mm. However so far, no systematic experiments are available that neither confirm an explicit size dependence of this anomalous sintering mechanism nor constrain the transition to normal sintering. A potential transition in that range would however be relevant for ice crystals in natural snow (for example in a natural snow pack with advanced metamorphism, on ski pistes and roads in spring), or in the food industry [16,17] . Since the quality of the sintered bonds highly affects physical properties (thermal, optical, and mechanical), better knowledge about sintering in this particle size range is important for applications.
To contribute to this problem and investigate the size dependence on the nature of sintering and growth conditions between ice spheres, we have conducted ice sphere sintering experiments geometry can be described by superposition of spherical harmonics, viz. with particle sizes between 1 and 2.3 mm in diameter. Sintering proceeded over ten days at −10 • C and the results were captured via 3D micro-computed tomography ( μCT) imaging. For six samples with different particle sizes, we only captured the final state (referred to as size series). For two additional samples with 1 and 2.1 mm particles, we observed the sintering process in time-lapse μCT scans, with images after 0, 3, 5, 7, and 10 days (time series).
Using various image analysis methods, we determined the specific surface area and surface curvature, an estimate of the porosity of the bonds, and the bond size. The bond size increased in time for all samples, indicating an apparently monotonic sintering process. In contrast, other bond metrics revealed contrary evolutions between small and large spheres, indicating the growth of solid bonds in the sample with 1 mm particles and porous bonds for 2.1 mm particles. The sintering conditions in the 2.1 mm sample of the time series were further investigated using Finite Element Method (FEM) simulations of the temperature field in the sample, as a key driver for vapor transport. The combination of microstructural and thermal information illustrates that even small fluctuations in the temperature of the sample environment, which can hardly be avoided experimentally, induce significant temperature gradients within the sample that are concentrated in the contact region between two particles. Besides the localization of the temperature gradient, the FEM simulations revealed a decrease in the local temperature gradient after a certain time. The observed nonmonotonic evolution of the simulated local temperature gradient in the bond suggests an interpretation of the anomalous sintering as a self-inhibited bond growth process. The paper is organized as follows. The next section describes the sample preparation, sintering conditions, and the image analysis methods. In the subsequent sections, the results are presented and discussed in view of particle size dependence, retreating protrusions, and practical implications.

Sample preparation
The ice spheres were produced with the tumbling method described in [7] . The latter publication provides a comprehensive characterization of the ice spheres, including an analysis of their sphericity and crystallographic structure. By varying the tumbling time, we obtained particles at different sizes. The sizes were selected by sieving: either as a size range that is defined by the sieve mesh size of two consecutive sieves used in series, or only by one sieve, where we only took the particles that remained in the mesh but passed through the sieve when shaking it. The particles were sieved directly after tumbling before they could sinter to make sure they passed the sieves one by one. Using such a simple method of segregation, different size classes could be obtained as listed in Table 1 . The smallest and the largest sizes were used in time-lapse experiments, which focus on the continuous time evolution. The respective experiments are referred to as time series ( tS) in the following. The size classes in between the smallest and largest spheres are used in an experiment to continuously assess the diameter dependence of the final sintering state, this experiment is referred to as size series ( dS). In total, eight samples were investigated. A sample is specified by the particle size, e.g. tS2.1 Fig. 1. Touching particles were extracted from the sample volume to examine the sintered bonds in a smaller subvolume (A). The porosity was quantified in 2D slices, in which air between protrusions appears as enclosed pores. By iterating through all slices along the three main axes (B, here z-axis), the air bubbles within the particle volume were removed from each particle (C), to only count the pores in the bond (D). and dS2.1 for both samples with spheres of 2.1 mm diameter. Here, and in the calculations, we used the middle value of the size range as particle size d p .
For the size series ( dS), six samples were prepared with the particle sizes d p listed in Table 1 . The particles were filled into cylindrical μCT sample holders with a diameter of 10 mm to a height of approximately 10 mm. The sample holders were sealed with a lid to avoid evaporation. For isothermal sintering, all six samples were placed in a Styrofoam box, which was stored for ten days in a small climate chamber. The chamber was set to −10 • C and remained closed for the ten days. After ten days, the six samples were consecutively taken out for scanning in the μCT.
For the time series ( tS), two samples were prepared in the same way, with particles at the upper and lower end of the particle size range (2.1 and 1 mm). For both samples, the sintering process was captured by μCT time-lapse scans. The samples remained in the μCT scanner for ten days and were scanned initially (0), after 3, 5, 7, and 10 days. The temperature in the scanner was set to −10 • C .

Temperature logging
For all samples, the temperature was logged every 10 min during the sintering process. We used iButtons from the company Dallas Semiconductor (now Maxim Integrated) with a resolution of 0.5 • C . In order to derive a potential temperature gradient, the temperature was measured at the top and bottom of the sample. One iButton was fixed on top of each sample holder and a second was placed beside the sample, either in the μCT scanner for tS, or in the Styrofoam box for dS.

3D μCT image analysis
The samples were scanned in a Scanco Medical μCT 40 scanner. With a sampling time of 1.2 s, a current and peak energy of 145 μA and 55 kVp for the X-ray tube, the scanning duration was 2.8 h and the voxel size was 6 μm. From the resulting grey scale images, a volume of interest of (1100 ×1100 ×550) voxels was then segmented, by fitting the sum of three Gaussian distributions to the grey scale histogram [18] . The obtained binary image of the two phases, ice and air, was used for the subsequent evaluation. For dS1.2, the section had to be reduced to (1100 ×110 0 ×40 0) voxels, to eliminate a scanning error at the top of the sample.

Specific surface area, surface curvature, and pore space
The specific surface area SSA, the ratio of ice surface area to ice mass, was used for the investigation of the sample surface and the growth of small structures. It was calculated using a triangulationbased method that is provided by the vendor software [19] .
The pore space of the tS samples was evaluated regarding the mean pore sizes SP. To this end, the pore space is fitted with spheres of different radii and their average is SP. The method [21] is provided by the vendor software.
The specific surface area, curvature, and pore space are characteristics of the entire sample and will be used to characterize the anomalous features of the ice-air surface during sintering. To detect also anomalous, volumetric features (like porous bonds) a non-standard analysis of the bonds is required.

Bond porosity and bond size
The quantification of the bond properties relies on automated bond-specific image analysis that was performed on the 3D binary images using the programming language Python. The analysis of the bonds is based on the identification of particles and their neighbors, which is achieved by watershed segmentation, similar to the method described in [7] . After the particles were identified and labeled, pairs of touching particles were extracted from the volume to examine their bond in a smaller subvolume.
For the quantification of the porosity, the subvolumes were analyzed in 2D slices, by iterating through 2D slices along the main axes (sketched in Fig. 1 ). This was required, since the bond porosity that formed during sintering, is not necessarily due to enclosed pores in the ice phase, but rather formed by several column-like structures that connect the two particles. Thus, in 3D the air in the bond remains connected to the main air phase that surrounds the particles. In 2D, the air phase between two columns is separated from the main air phase and appears as enclosed porewhich makes it countable. To only investigate the pores in the bond, air bubbles in the ice volume of the particles should be excluded from the evaluation and had to be removed. This was done for each particle in a separate subvolume ( Fig. 1 (C)) by identifying the largest connected air region and assigning all other voxels to the ice phase. Afterwards, the two particles were rejoined and enclosed air pores remained only in the bond between the particles ( Fig. 1 (D)). A 2D slice of a subvolume containing two particles after elimination of the air bubbles is depicted in Fig. 2 . Then the closed pores in the ice volume of each 2D slice were labeled and counted as n pores , as well as the number of slices n slices in which the two 2D-particles were in contact, i.e. the slices that contain the actual bond. From this, we defined the bond pores χ by normalizing the number of pores in a bond by the bond extension: The identified pores are presented in different colors in the inset of Fig. 2 . The entire evaluation was performed in all three main directions of the randomly oriented spheres, i.e. the volume was sliced along the x -, the y -, and the z-axis, to show the reliability of the results. The results of the three evaluations may deviate from each other for a single bond, but are highly correlated in the average over all bonds in the volume of interest. The presented results show the mean of all three evaluations and their standard deviation.
The estimate of the bond size r b relies also on the watershed segmentation. The algorithm separates two identified particles by a one voxel thick area, the watershed, which we take as the bond between the particles. Due to the porosity of the bond, this contact area is not necessarily connected and may consist of several patches. By again separating two particles from the volume, all contact patches that belong to the bond between these two particles were identified. The total number of voxels within the contact area n vox , was converted with the voxel size s vox to a radius r b = (n vox /π ) 1 / 2 · s vox of an area-equivalent circular contact area.

FEM simulations of the temperature field
In addition to the structural information from the binary images, we obtained insight into the thermal conditions from voxelbased FEM simulations of the temperature field in the digitized sample structures. Local temperature gradients play a key role in crystal growth during metamorphism (and sintering) in snow [15,22,23] and provide us information about the local growth conditions. Based on the assumptions that 1) the vapor density in the pore space is saturated for a given temperature [24] (equilibrium condition) and that 2) the ice-air interface evolves slow compared to vapor and heat diffusion [25,26] (stationary diffusion approximation), it has been shown that the velocity of the ice-air interface during diffusion-limited crystal growth can be approximated as a function of only the local temperature gradient [23] . The leading order role of temperature gradients (in solid and vapor phases) has also been highlighted in an analysis of interface stability [15] . Accordingly, we used the temperature field within the sample to obtain insight into the evolution of the growth conditions and investigate the anticipated non-monotonic growth of the bonds and protrusions [7][8][9] .
The FEM simulations were performed on all five time-lapse μCT images of the tS2.1 sample, similarly as described by Pinzer et al. [22] . All voxels (ice and air) were converted to eight-node brick finite elements, on which the steady-state temperature distribution equation was solved. The boundary condition (Dirichlet) was defined with fixed temperatures at the top and bottom nodes, according to the temperatures measured in the experiments with the iButtons at the sample top and bottom. The thermal conductivities of ice and air were assumed to be κ ice = 2 . 34 WK −1 m −1 and κ air = 0 . 024 WK −1 m −1 . Heat transfer by advection was neglected.
From the resulting temperature field within the sample, the temperature gradients in all three main directions ( x , y , z) were determined. Thus, the simulation results consist of 3D images of the temperature gradient and the heat flux q . The 3D images of the vertical temperature gradient ( z-direction) were used for further evaluation. We determined the average ∇T and the maximum ∇T max of the temperature gradient in the entire sample and also in subvolumes of the bond regions, denoted as ∇T b and ∇T b max . The respective subvolumes for calculating these values in the bond regions were identified in the watershed segmented images that were used for the evaluation of the bond properties. By evaluating the contours of adjacent particles using the contours function in the Python Skimage library, the contact between the particles and its extension in the x -, y -, and z-directions could be determined. This defines a volume around the bond, which was augmented by 10 voxels in each direction for the evaluation of the temperature gradient.

Results
A visual overview of all investigated samples is presented in Fig. 3 . It shows a 2D slice of each 3D μCT scan: the five images of the two time series (upper and middle row) show the same slice after 0, 3, 5, 7, and 10 days; the lower row shows one slice of each sample of the size series. The slices were taken from an x − z plane and are oriented such that the z-axis points upward, matching the vertical orientation of the samples in the scanner.
While the sphere surfaces in tS1 remain apparently smooth (top panel), the images of the tS2.1 sample (middle panel) reveal protrusions that can be seen on the top surface of the upper particle, as well as at the contact with the particle below where a porous bond is formed. These features, protrusions and porous bonds are most distinct in this sample. In the other samples, some protrusions (e.g. in dS2.3) and porous bonds (green circle in tS1) can be observed. In the following, we first present the results of the overall sample, then of the bond specific metrics.

Specific surface area
The measured specific surface area SSA CT is evaluated by calculating the ratio with the SSA ideal = 3 /r p of a size-equivalent sphere, which was calculated with the known particle radius r p ( = d p / 2 ) for all investigated samples. The ratio emphasizes the evolution of features on the particle surfaces. It is shown in Fig. 4 (A) for all samples against the sintering time and reveals contrarily evolving curves of the time series: from zero to ten days, the ratio increases for tS1 and decreases for tS2.1. For the dS samples, the ratio is shown at ten days, which decreases with increasing particle size d p and lies consistently in between the two time series. This size dependence becomes more clear in Fig. 4 (B), which shows the same data as a function of particle diameter d p . The errorbars in the dS samples denote the size ranges of the sieves and show that the entire range from 1 to 2.3 mm is covered, with the largest uncertainty in the center.
A decreasing SSA-ratio, for tS2.1 in time and dS with size, implies an increase of the surface area. This increase can be caused by the growth of small structures on the particle surface, such as  the protrusions observed on some particle surfaces, as well as the porous bonds. In the increasing SSA-ratio of tS1, we observe a decrease of the sample surface area in time. This indicates the formation of solid particle contacts, which reduces the surface area as it is known from normal sintering.

Surface curvature H
Further evidence for the emergence of protrusions can be obtained from interfacial curvatures. An illustrative example of the curvature evaluation is given in the 3D presentation of the surface in tS2.1 after seven days in Fig. 5 . It shows the neck region between two particles colored by the mean curvature. Dominant protrusions are on the surface of the lower particle, which protrude towards the particle above (the z-axis indicates the sample top similarly as visualized in the 2D representation in Fig. 3 ). This vertical direcionality is due to a temperature gradient, as detailed further below. The protrusions seem to emerge from general (low amplitude) surface patterns which are visible on the remaining parts of the surface. This variety of features gives rise to a dis-tribution of mean curvatures which is different from the one of a sphere. The distribution of the main surface curvature H is shown in Fig. 6 , in (A) for both time series and in (B) for all samples after ten days. Negative values correspond to concave curvatures, and positive values correspond to convex curvatures. In all samples, the surface roughness gives rise to a maximum around 0 mm −1 and the characteristic curvature scale that corresponds to the particle size ( H = 1 /r p ) cannot be derived. Aside from surface roughness, negative curvatures stem mainly from concave air bubbles within the ice and from the contact region between two particles. Positive curvatures are attributed to protrusions, besides surface roughness.
In time, the curvature of tS1 changes little. This can be seen more clearly in the average curvatures listed in Table 2 . The most obvious changes are observed in the tS2.1 sample for small positive curvatures (inset in Fig. 6 (A)) that also lead to an increase of the average curvature in time ( Table 2 ). Initially, the small positive curvatures are equally high in number as for tS1, but are distinctly higher at later times. Similarly increased positive curvatures were measured in the dS2.3 sample. This increase in the restricted range of small positive H between 15 and 40 mm −1 is attributed to the  protrusions. These values correspond to the curvature of spheres with a radius of 0.025 to 0.07 mm, which gives us an estimate for the length scale of the protrusions that matches well the qualitative observations in the images.

Pore space SP
The time evolution of the average pore size SP, listed in Table 2 , reveal contrary evolutions in the two tS samples as well. While SP in tS1 increased only slightly, in tS2.1 it decreases about 10% during ten days. This is another evidence of the small structures that evolve on the particle surfaces and grow into the pore.

Bond pores χ and bond size r b
The bond pores χ and the bond size r b were derived as bondspecific metrics for each bond. The average of χ for all bonds and the three directions is shown against sintering time for all samples in Fig. 7 (A). The errorbars (very narrow, only visible in tS2.1 at seven days) show the standard deviation of the three independent evaluations along the main directions. The value of χ decreases monotonically for the tS1 sample with an initial maximum of 0.38, which means that about 2/3 of the bonds sintered solidly, on average. In the tS2.1 sample, the average number of pores per bond increases from 1.3 to 2. The size series reveals increasing χ with increasing particle size. Almost no pores were measured in the smallest sample, while in the largest sample on average 3/4 of the bonds contain a pore.
The two time series tS show again opposite behavior for the evolution of χ. The constant decrease of χ in tS1 can be explained by filling of pores that might form initially due to surface roughness and non point-like contacts. An example of this is shown in the bond with a green circle in Fig. 3 . The development of χ in the tS2.1 sample is more interesting. An increase of χ corresponds to the formation of porous bonds, the growth of protrusions that form the contact. This is observed in tS2.1 up to seven days, as already indicated in previous metrics (SSA, H). Interestingly, χ decreases from seven to ten days. This could either be explained by a closing of the pores, as described for tS1, or it could imply that the protrusions actually retreat, invert the formation of a contact, and then open the pores again. The contact size is shown in Fig. 7 (B) against sintering time. The contact size is expressed as the radius r b of an area-equivalent circular contact area. It increases for both time series and shows no clear tendency with particle size.

Temperature measurements
The logged temperatures are shown in Fig. 8  For both samples, the top temperature is higher than the bottom temperature. A larger temperature difference within the tS sample is obvious. From this differences, we calculated the temperature gradients within the samples and obtained an average of 0.015 K/m for the dS samples and 0.87 K/m for the tS sample. The   latter was used as input for the FEM simulations of the temperature field within the tS2.1 sample.

FEM simulations
The spatial pattern of the z-component of the temperature gradient around the spheres is visualized with colors in a 2D slice of the ten days image for sample tS2.1 in Fig. 9 (A).
It shows largest temperature gradients between the particles with local maxima at the protrusion tips. In the ice phase, the temperature gradient is minimal, indicating a rather homogeneous temperature distribution within each particle. The time evolution is shown in Fig. 9 (B) with the average ∇T and maximum ∇T max of the vertical gradient for all five images of the tS2.1 sample (normalized by the initial value). The average is constant over time as it has to be due to the imposed simulation boundary conditions, but the maximum decreases towards ten days. The same evaluation of average and maximum values was performed on the subvolumes of all bonds and is shown in the same figure. The maximum ∇T b decreases in time, similarly as ∇T max in the entire sample. In contrast, a non-monotonic behavior was found for the average temperature gradient around the bonds ∇T b max , which initially increases and then decreases after five days.
The temporary increase of ∇T b max in the bond volumes indicates a localization of the gradient in the contact regions. Its decrease after five days, as well as the decrease of the maximum, can be attributed to improved thermal contact during ongoing sintering, which is evidenced by the increasing heat flux q in the sample ( Fig. 9 (B)).
For one bond, the time evolution of the temperature gradient and the related structural changes are shown explicitly: 2D slices of the five FEM images are shown together with the equivalent 2D slices of the binary images in Fig. 10 , and the corresponding ther-mal ( ∇T b max and ∇T b ) and structural parameters ( χ and r b ) from this 3D bond subvolume are plotted against the sintering time in Fig. 11 (normalized to the initial values). This analysis reveals the concurrent evolution of the maximum temperature gradient, the extension of the protrusion at which it is located, and the bond pores χ. The 2D FEM images ( Fig. 10 upper row) show the localization of the maximum temperature gradient at the protrusion tip that formed between zero and three days. The maximum temperature gradient in the entire subvolume ∇T b slightly increases during the first three days and decreases to the end of the experiment ( Fig. 11 ). The respective protrusion grows in time but retreats again, without contacting the other particle ( Fig. 10 lower row). Its maximum length seems to be reached after seven days. The evolution of χ in this bond, with its maximum at seven days, is in accordance with the average χ in this sample. The bond size reduces after three days again. In this bond, we observe the interrelation between growth and local temperature gradient. It seems in this particular case that the decrease in χ is caused by dissolving protrusions, as seen visually and in the measured r b . In all investigated parameters, except for r b , this bond is representative for the bonds in the tS2.1 sample.

Sintering by unstable crystal growth
The two series of size and time dependent sintering experiments with ice spheres investigate sintering anomalies, their growth conditions and dependence on particle size. The anomalies, which have also been found in previous experiments [8][9][10] , consist of the evolution of protrusions on the particle surface and the formation of porous bonds between particles, instead of solidly sintered bonds. These features are visible in the 3D (and 2D) μCT images and were quantified via surface metrics (SSA, H) and bond metrics ( χ, r b ) that characterize the evolution during sintering.
An estimate of the protrusion length scale (typical width of the protrusions) was obtained from the surface curvature distribution. If the range of curvatures that are affected by the protrusion evolution ( H between 15 and 40 mm −1 , Fig. 6 (A)) is converted to an equivalent radius R H , we estimate the protrusion width with 2 R H . For comparison with previous works [10] , where the protrusion length scale was expressed as wavelength λ of a sinusoidal perturbation on a growth interface, we consider λ = 4 R H (since the protrusion shape is covered by λ/ 2 ). The values we obtain for λ are in the range of 0.1 to 0.27 mm. For these length scales, instability was predicted [10] , using the dispersion relation of Mullins-Sekerka theory [12] . This theory of surface instability during diffusion limited crystal growth predicts the evolution of small surface perturbations ( λ) for particles larger than a critical particle size D c . The size of the particle dictates the local curvature. The local cur-   vature, in turn, dictates the concentration field above the surface, which governs the local crystal growth.
Crystal growth or sintering in snow is mediated by vapor diffusion, which is highly sensitive to temperature gradients. Experimentally, a temperature gradient is hardly avoidable, as the top and bottom temperature measurements in the samples show ( Fig. 8 ). In fact, temperature gradients below 10 K/m are commonly considered "isothermal", e.g. [27,28] . Even if the average temperature gradient across the sample is very small (as for the dS2.1 sample), considerably higher gradients occur in the contact region due to the geometry (initial spheres in point contact) and large difference in thermal conductivities between air and ice [29] . The main heat flux occurs through the ice phase and the temperature in a particle is quasi homogeneous, except for the parts where it contacts other particles that have slightly different temperatures (cf. Fig. 9 and [29] ). These contact points are bottlenecks of the heat flow and therefore regions of high temperature gradients. They are the regions where sintering takes place. Therefore, although we aimed for isothermal sintering at −10 • C , the unavoidable small experimental temperature fluctuations lead to significant temperature gradients at the particle contacts. The influence of the downward directed temperature gradient can be directly seen in the images in Fig. 3 : the gradient determines the growth direction, upwards, and amplifies the protrusion growth (cf. tS2.1 and dS2.3). This fact, however, provides interesting insights into the relation of temperature gradients and critical particle size, since we measured different temperature gradients in the two sets of experiments, discussed in the next section.

Particle size dependence
The two samples of the time series bound the investigated particle size range and provide the upper and lower limits in most of the evaluated metrics. For both samples, the bond size r b increases in time, as expected during sintering. In contrast, the other investigated parameters reveal a contrary evolution of the two samples. This is consistent with the growth of solid and porous bonds for t S1 and t S2.1, respectively. The specific surface area decreases for tS1 as the increasing contact area between touching particles vanishes from the surface, and the mean of the bond pores χ decreases due to the filling of initial pores (green circle in Fig. 3 ).
The mean curvature ( Table 2 ) of tS1 slightly decreases in time as negative curvatures in the concave contact area gain weight in the H distribution ( Fig. 6 ). In tS2.1, the opposite is observed. The SSA increases about 20% during the ten days, which can be ascribed to the growth of protrusions and porous bonds. Chen et al. [10] measured the SSA explicitly in the bond region between two particles and observed an increase of the SSA with the evolution of a porous bond as well. In their case, the SSA increased even faster (from ∼2.5 to ∼4 mm −1 within two days; while we observed an increase from 4.4 to 5.7 mm −1 during ten days), which can be explained by the fact that they evaluated only the SSA around the contacts, whereas we considered the entire sample volume. We observed in χ an increasing bond porosity in time and the mean of χ is about ten times higher than in tS1 ( Fig. 7 (A)). The mean curvature increases in time due to the small positive curvatures of the convex protrusions. In the H distribution, the positive curvatures of the initial images were similar for both tS samples, but clearly increased in tS2.1 for all other times. This confirms the visual observation of the protrusion evolution in the 2D images ( Fig. 3 ).
For the dS samples, the values of all parameters lie consistently in between those of both tS samples after ten days, mainly in order of particle size. For most parameters ( r b , SSA, χ), the values of the smallest sample are very close to those of tS1, whereas a larger difference was measured for the large particles. This can be explained by the higher temperature gradient in the tS sample ( Fig. 8 ), which amplifies the growth of the protrusions. The results of both series show the roles temperature gradient and particle size play in the evolution of protrusions: a temperature gradient accelerates the process (comparison of tS and dS samples) and dictates the growth direction ( Fig. 3 ), but as seen in both series, despite similar temperature conditions only the large particles evolve protrusions while the small ones do not. This highlights the role of the particle size and indicates that the critical particle size D c must lie between 1.2 and 2.1 mm for sphere packs which are subject to quasi-isothermal conditions with temperature fluctuations up to 1 K/m. An exact value is not resolved by the size series. The plot of SSA against particle size in Fig. 4 (B) shows the uncertainty of the particle sizes (sieve sizes) and the largest uncertainty occurring in the center of the investigated range. Nevertheless, the results allow for a significant reduction of the range of D c from previously (0.7, 4) mm [5,8] to (1.2, 2.1) mm.
It is important to point out here that D c theoretically also depends on the concentration fields [11] , in our case the excess vapor, which is influenced by the exact value of the temperature gradient. This cannot further be quantified by our experiments. Even in the FEM simulations, the limited accuracy of the applied voxel-based method raises a substantial error (up to a factor of 12.5 [23] ) in the numerical solution at the interface, due to the two-phase nature of the system. However, we measured different temperature gradients in the two sets of experiments, but the results (SSA, H, χ, r b ) agree with the observations of solid bonds in small and porous for large particles. Thus, the range found for D c does not seem to be heavily influenced by the precise value of the temperature gradient below 1 K/m. If we look at previous studies in the literature, we can assume: 1) similar conditions in the 4 mm spheres sintering experiments under "quasi-isothermal" conditions at -10 • C , where porous bonds were observed [8] ; 2) a saturated environment for the d p ≤0.7 mm spheres that sintered in a sealed container with ice-coated walls, where solid bonds developed [5] . Both are in accordance with our D c -range. In contrast, the protrusions that were observed on non-spherical snow particles [9,30] most likely grew on surfaces with higher curvatures compared to the curvature of 1 mm spheres (which would correspond to spheres with even smaller diameter than 1 mm). However, the particles were most likely exposed to higher supersaturation in both cases, which reduces the D c [11] . In one case [9] small, fractured ice par-ticles with sharp edges that are prone to evaporation were mixed into the sample by sieving of the snow. In the other case [30] , the samples were exposed to extreme temperature changes during short transportation phases such that even frost was observed in the samples, which must have been induced by a high temperature gradient.

Retreating protrusions
The decrease in χ after seven days in the tS2.1 sample could be explained by retreating protrusions and by filling of pores. Both retreating protrusions [8,9] and filling of pores by vapor condensation [8] was observed on similar time scales before. The predominance of the two processes cannot be deduced from the structural metrics ( χ, r b ), but the observed thermal conditions provide complementary information. In sintering experiments with 3 mm spheres and an externally applied sign-alternating temperature gradient, retreat was observed when the gradient was inverted [31] , and explained by this obvious change in thermal conditions, which is crucial for protrusion growth. The existence of a temperature gradient in our experiments and its amplifying influence on protrusion growth has been discussed above. The time series of the 2.1 mm sample with the FEM simulations of the temperature field allow us to extend this discussion and point out the possibility of retreating protrusions with unidirectional temperature gradients.
The visualization of the FEM simulations reveals the localization of the temperature gradient in the bond region and especially at the protrusion tips ( Fig. 9 (A)), which drives the protrusion growth, as pointed out before. Its magnitude decreases over the course of the ten days, as we observed for ∇T max in the entire sample, as well as for ∇T b in the bond subvolumes ( Fig. 9 (B)). In a closeup analysis of one bond, we observed the interrelation of temperature gradient and protrusion growth ( Fig. 10 ), where ∇T b first increases as the protrusion forms, but decreases with time as the protrusion also decreases. With the decreasing temperature gradient, the driving force for protrusion growth weakens and may be overcome by the disposition of the system to reduce its surface area and energy by dissolving the protrusion (or filling the pore space). Although this visual observation is restricted to a 2D representation, the supporting metrics ( χ, r b , ∇ T b , ∇ T b max ) were taken in the 3D bond subvolume. The analysis of the entire sample and all other bond subvolumes, in terms of in χ, ∇T b , and ∇T b max , is consistent with this particular observation.
The question remains why the magnitude of the temperature gradient decreases after five days. The key to the explanation lies in the dynamics of the system. While the bonds grow during sintering (porous or not porous), the thermal contact improves, as evidenced by the increasing heat flux q ( Fig. 9 (B)). From an analytical investigation of this problem [29] , we estimate a q that is five times higher for a contact area of r b = 0 . 4 mm ( Fig. 7 (B)) compared to a point contact (with κ ice = 0 . 024 WK −1 m −1 and κ air = 2 . 34 WK −1 m −1 ). The latter study concludes that the system's conductivity dominates over the material's conductivity and emphasizes the importance of the contact geometry. As a consequence of increased r b and q , the local temperature gradient in the bond region decreases, which reduces the driving force for protrusion growth. Yet, a conclusion of whether the reduction in χ is induced by dissolving protrusions or filling of pores cannot be drawn. A combination of both is possible, as both would satisfy the system's desire to decrease surface area and energy.
Overall the results suggest the following mechanistic picture of anomalous sintering in ice spheres packs: Due to the initial singular geometry of large spheres with point contacts, small temperature gradients across the sample significantly increase the gra-dients in the vicinity of the bonds which are initially the bottlenecks of heat flow. The high values of local temperature gradients locally induce high supersaturations that fulfill the conditions for unstable growth. This situation is more favorable for large spheres. The onset of unstable protrusion growth leads to improved particle contacts which are porous in nature due to the underlying growth patterns. The growing particle contacts now improve the heat flow and thereby reduce the temperature gradients in the vicinity of the necks. This feedback acts as a self-inhibitor of the growth process as it reduces the driving force for unstable growth. The established, high curvature protrusions that form the neck relax thermodynamically and thereby retreat by minimizing their curvature and surface energy.

Implications
For various applications, it is of practical relevance to know that sintered contacts might differ from the theoretically predicted solid necks and monotonic evolution. Better knowledge about the size dependence, growth conditions, and retreat of bonds is important for their usage, design, and interpretation. The nature of the bonds governs the physical properties of the system: the mechanical [7] and thermal [29,32] properties are affected by the bond porosity, and the optical properties are directly related to the SSA [33] , which is considerably affected by the sintering anomalies. Theoretical investigations of surface instabilities are crucial to understand the growth mechanisms and estimate the occurrence, but exact predictions are hindered by unknown parameters [10,11] . Previous experimental reports are very limited and leave a wide range for the critical particle size.
Our experiments provide help for selecting an appropriate particle size for future experiments, since ice spheres are a popular remedy for experimental investigations of simplified and well defined structures [4,5,7,8] . The experiments presented here explain, with the measured bond pores, the weak sintering of the 2.1 mm ice bead samples in the compression experiments and the decreasing compressive strength between seven and ten days [7] . These ice bead samples sintered similarly to the dS samples at -10 • C in a Styrofoam box. Since these compression experiments were designed for model calibration, and control of the bonds is crucial, we now know that particles with diameters less than 1.2 mm would be favorable. Yet, the compression experiments, together with the sintering experiments presented here, show the mechanical relevance of porous and even retreating bonds. This is significant for other systems of similarly sized particles, for example in the food industry while designing or preserving food [16,17] , or for the preparation of ski pistes or roads. The desired strength of a ski piste can often not be achieved in spring, when snow mainly consists of large melt form particles that do not sinter sufficiently over night [34] . Sintering of ice is also relevant for extra-terrestrial applications. The surface of planetary bodies often consists of sintered ice particles and the necessity of knowledge about sintering mechanisms and anomalies for investigating these remote objects is obvious. Sintered contacts and their physical properties are crucial for investigating planetary surfaces from telescopic or thermal observations, or for developing landing technologies of space craft [32,35] .
Our experiments further constrained the critical particle size for which surface instabilities may occur and provide evidence for a change in thermal conditions that allow for the retreat of bonds. The discussed image analysis methods are capable of resolving the protrusions as structures on the surfaces (SSA, H), and the pores of the bonds ( χ). These methods and metrics are applicable in further experiments for investigation of sintering anomalies. The discussion of protrusion retreat is based on FEM simulations of the temperature field within the tS2.1 sample. An automated method that explicitly resolves the protrusions would be necessary to further investigate this interesting process.

Conclusions
Two series of ice sphere sintering experiments monitored by μCT imaging have detailed the size controls of sintering anomalies (protrusions and bond porosity) in ice sphere sintering. A critical sphere size D c is required to observe the anomalies. Only slight deviations from isothermal conditions, which can be hardly avoided in any experiment, is the driving mechanism for protrusion growth and formation of porous bonds. The growth is self-inhibited when the local temperature gradient decreases due to the improved thermal contact in the course of sintering. The achieved reduction of the possible range of the critical size D c from (0.7, 4) mm [5,8] to (1.2, 2.1) mm is of practical relevance for applications of sintering. It is worthwhile to further investigate D c and retreating protrusions to explore the dependence on temperature on longer time scales. Locally, the observed retreat of a bond and decrease of the thermal contact may eventually bring the temperature gradient back to conditions of unstable growth. This non-linear behavior thus raises the fundamental question if even oscillatory modes of bond growth in sintering under slight temperature gradients exist in nature.

Declaration of Competing Interest
None.