Discrete element modeling of planetary ice analogs: mechanical behavior upon sintering

Potentially habitable icy Ocean Worlds, such as Enceladus and Europa, are scientifically compelling worlds in the solar system and high-priority exploration targets. Future robotic exploration of Enceladus and Europa by in-situ missions would require a detailed understanding of the surface material and of the complex lander-surface interactions during locomotion or sampling. To date, numerical modeling approaches that provide insights into the icy terrain’s mechanical behavior have been lacking. In this work, we present a Discrete Element Model of porous planetary ice analogs that explicitly describes the microstructure and its evolution upon sintering. The model dimension is tuned following a Pareto-optimality analysis, the model parameters’ influence on the sample strength is investigated using a sensitivity analysis, and the model parameters are calibrated to experiments using a probabilistic method. The results indicate that the friction coefficient and the cohesion energy density at the particle-scale govern the macroscopic properties of the porous ice. Our model reveals a good correspondence between the macroscopic and bond strength evolutions, suggesting that the strengthening of porous ice results from the development of a large-scale network due to inter-particle bonding. This work sheds light on the multi-scale nature of the mechanics of planetary ice analogs and points to the importance of understanding surface strength evolution upon sintering to design robust robotic systems.


Introduction
One of the main challenges in developing robots for in-situ space exploration is the uncertainty associated with the surface properties of the planetary bodies. Properties like ultimate bearing capacity and penetration resistance affect the ability of space robots to land, sample, and explore their environment. Hence, estimating these properties is critical for the design and optimization of the robots' hardware and control system and the development of suitable sampling systems.
Several robotic mission concepts are under development to explore Ocean Worlds [10,34]. Ocean Worlds are bodies in our solar system that are considered potentially habitable or inhabited. Saturn's moon Enceladus and Jupiter's moon Europa recently gained significant interest among the scientific community, and are now regarded as two of the most likely Ocean Worlds to harbor extraterrestrial life. Indeed, they appear to contain all essential ingredients of life, namely liquid water, energy, and nutrients.
Enceladus and Europa are likely to harbor an internal ocean of liquid water underneath their thick ice shells [66,82]. Enceladus' internal ocean spews material from large surface fractures around its southern pole known as the Tiger Stripes. The plume formed by the confluence of these geysers is primarily composed of micron-sized particles dominated by water ice [83] and are partly deposited back on the surface to form the outermost crust [44,50,81]. This unique geological context provides direct access to 12 Page 2 of 20 subsurface material, making Enceladus a high-priority target for planetary missions in the search for biosignatures [54]. Europa is also thought to be geologically active with plumes, though definitive data remain elusive [41]. The upcoming Europa Clipper mission [65] may confirm the presence of such features.
After the plume material is deposited on the surface, the size and microstructure of the ice particles are expected to evolve via sintering, transforming initially unconsolidated deposits into consolidated porous ice [50]. Sintering of ice is a metamorphism that describes the diffusion of water molecules within and between the ice particles, as well as the evolution of inter-particle bonds. Sintering leads to an increase in the inter-particle neck size and the densification of the aggregates, resulting in the material as a whole becoming stronger over time [5]. Sintering is a temperature and particle-size dependent process in which the evolution rate increases monotonically with increasing temperature and decreasing particle size. While earlier work provided an improved understanding of sintering timescales in planetary environments [60], planetary ice sintering remains a poorly understood process that is now the subject of active research.
In a recent study, laboratory icy plume deposit analogs were produced and left to sinter over extended periods of time and at several temperatures [16]. Cone penetration tests were performed at frequent intervals to investigate the mechanical strength of the samples. [16] showed that the observed temperature dependence of the strength evolution is commensurate with a mechanism dominated by diffusion of water molecules on the surface of ice particles. Their study revealed a link between thermodynamic processes and sample strength, but was not suited to understand and quantify the mechanical properties of planetary ice analogs at length scales ranging from particle-to macro-scale.
In this work, we develop a numerical model that aims to unravel the roles of particle-scale properties of the porous ice in the overall mechanical behavior. We use the Discrete Element Method (DEM) to represent the sintered ice particles and simulate cone penetration tests on ice plume deposit analogs at different consolidation levels. DEM allows studying complex behavior via simple particle-particle interactions. Alternative approaches to model this ice are presented in the Supplementary Material (Text S1). Our model, composed of homogeneous particles that interact via contacts, mimics the material's microstructure and reflects the ice properties both at the particle and at the contact level. As a result, the macroscopic mechanical behavior emerges from the microstructure and micromechanical interactions. Our model also describes the contacts and the sintering process based on the physics of the ice.
DEM has seen rapid development both in academic research [64,86] and in engineering to solve problems involving granular material, such as in the mining, food, and pharmaceuticals industries. DEM has also been recognized as an accurate and computationally effective tool to simulate terrain-tool interactions [55,79]. The method has been used to simulate ice, particularly for applications involving interactions with offshore structures and ship hulls [45,68,85]. [46] developed a DEM-based model but for sintered snow. However, to the best of our knowledge, no previous study has attempted to model sintered porous granular ice.
In this work, the ice particles are treated as rigid bodies that interact via contacts. The contact parameters are calibrated with experimental data of cone penetration tests. The inter-particle solid bonds are modeled with a cohesive force following a Simplified Johnson-Kendall-Roberts (SJKR) model [48]. Unlike other bond models [2,7,23,26,33,52,55,67,74], the SJKR model allows the bonds to fracture and reform dynamically throughout the simulation, allowing for indefinitely long simulations with large deformations and arbitrary interactions.
The paper is organized as follows. In Sect. 2, we briefly describe the ice preparation and testing methods. We also describe the fundamentals of DEM and detail the contact models used in our study. In Sect. 3, we introduce the simulation setup, the modeling procedure, and the parameters of our model. In Sect. 4, we investigate how to reduce the model complexity in an optimal and rigorous manner. In Sect. 5, we calibrate the model to fit experimental data following a proposed probabilistic method. Finally, we present and interpret the results in Sect. 6.

Laboratory ice analog preparation and testing
To create a laboratory analog of the ice plume deposits, deionized liquid water was atomized into liquid nitrogen. The water droplets instantly crystallize and form ice particles. This formed fine-grained ice was stored in sealed containers and left to sinter at four different temperatures (193 K, 223 K, 233 K, and 243 K) for time periods up to 14 months. The starting ice particles had a log-normal size distribution with a median diameter of 12 μm. The samples exhibited a porosity of 51.5% ± 1.6% , which remained sensibly constant throughout the aging process. The cone penetration resistance has been routinely measured using an in-house-developed experimental setup illustrated in Fig. 2a. The cone penetration tests consisted of driving a rod with a conical tip of 10mm diameter into the sample at a constant speed of 10 mm.s −1 . The vertical force acting on the cone tip is measured in relation to the penetration depth and used to derive the ice strength. A complete description of the material preparation, the measurement procedure, and a summary of the experimental results are presented in [16].

DEM modeling
In the DEM framework, as first described by [21], discrete rigid-body particles (or elements) interact with each other, transferring forces and torques via contacts. One appealing feature of DEM is that no assumptions about the material behavior at the continuum scale are required; instead, the material continuum behavior emerges from the large number of interactions at the particle scale [64]. In that sense, it is a faithful representation of the underlying physics. Nonetheless, DEM has a prohibitively high computational cost, limiting the simulation time and the number of particles in a given simulation.
In DEM, the translational and rotational motion of each individual particle is calculated by solving Newton's second law of motion at each time step, as summarized in the following equations: where m i , I i , v i , i are respectively the mass, the moment of inertia, the translational velocity vector, and the angular velocity vector of a particle i. g is the gravity vector, F n ij , F t ij , and r ij are respectively the normal force, the tangential force, and the rolling resistance applied on particle i from the interaction with particle j. R ij is the vector between the center of particle i and the contact point with the particle j. In this work, Eqs. (1) and (2) are solved using the Velocity Verlet scheme. More details on integration schemes are presented in the Supplementary Material (Text S2).
Prior to solving Eqs. (1) and (2), the forces and torques exerted on each particle by neighboring particles or boundaries are calculated using the three contact models depicted in Fig. 1. The normal and tangential forces at the contact point, respectively F n and F t , are modeled using the Hertz-Mindlin no-slip contact model with a linear spring-dashpot [59], as shown in Fig. 1a. In this framework, the normal and tangential forces have an elastic and a viscous component, denoted by the superscripts e and v, and are decomposed in the following equations: Additionally, the tangential force is governed by the Coulomb law of friction |F t | ≤ |F n | , where is the friction coefficient. Beyond this limit, slippage between particles occurs. The elastic force in the normal direction F e n is based on the classical Hertz's theory of contact between two spheres [36,47]. It assumes small elastic strains, small contact surfaces, and an ellipsoidal distribution of contact stresses. The elastic force in the tangential direction F e t is based on Mindlin's theory [59]. It is determined by both the normal and tangential overlap, respectively n and t . The tangential overlap t stems for the tangential velocity mismatch of two particles at their contact point. These models describe non-linear normal elastic force-displacement relationships which are expressed as follows: where k n and k t are the contact's normal and tangential elastic constants, R * = R i R j R i +R j the equivalent radius of the two bodies in contact, R i the radius of the i-th particle, ) the equivalent shear mod- ) the equivalent Young's modulus, i the Poisson's ratio of the i-th particle, and E i its Young's modulus. The viscous components in Eqs. (3) and (4) allow the system to dissipate energy to reach a steady state packing in reasonable time. In their original formulation, [21] proposed an expression of the critical damping ratio based on the critical damping time of a single degree-of-freedom spring mass dashpot system. [84] proposed that the critical damping ratio derives from the coefficient of restitution e -a physical parameter of the particles that characterizes the energy lost during collision and plastic straining -such that = −ln(e)∕ √ ln(e) 2 + 2 . When = 1 , the system is considered critically damped and reaches steady-state packing in the shortest achievable time. The normal and tangential viscous forces are expressed as follows: where n and t are the normal and tangential viscoelastic damping constants, v n and v t are the normal and tangential components of the relative velocity at the contact point, m * = m i m j m i +m j is the equivalent mass. S n and S t are, respectively, the normal and tangential stiffness and are defined by S n = 2E * √ R * � n � and S t = 8G * √ R * � n �. The rolling resistance represents the torques transmitted at the particles' contact region. The contact region formed when two particles are pushed together under a normal load can transfer torque owing to frictional forces distributed over the contact region. Thereby, the magnitude of the rolling resistance is proportional to the normal force-the source of the contact surface-and to a friction coefficient. It is calculated using a rolling friction model, as shown in Fig. 1b. While several rolling models exist [1,43,93], the Constant Directional Torque (CDT) is the most commonly used in DEM research owing to its relatively accurate and efficient calculations. In this model, the resistive torque r is proportional to the normal force F n and is oriented in the direction of the relative rolling motion ij | ij | . It is expressed as follows: where r is the coefficient of rolling friction and ij is the relative angular velocity of particle i with respect to particle j. Although the rolling friction has been somewhat successfully used to model particle shape [88], it remains unclear whether it is a physical property [88,93].
To represent inter-particle bonding due to sintering, we consider a supplemental cohesion model shown in Fig. 1c. Several bond models have been proposed like the dual spring [23], the Euler-Bernoulli [7,26], the cohesive beam [2,74], the bonded particle [33,55], the vector-based [52], and the parallel bond [67] models. These models typically describe a rigid bond between particles that forms initially and remains until breakage occurs. As a result, they can only be used to run short simulations with restricted interaction types.
On the other hand, the Johnson-Kendall-Roberts (JKR) model [48] acts by adding an additional normal force, also known as pull-out force, which tends to maintain the contact. The pull-out force corresponds to the force needed to separate two particles with zero overlaps. It can vanish and reappear throughout the simulation depending on the particle contacts, allowing the bond to fracture and reform dynamically. The JKR model was first developed from the observation of necks forming around the contact area of adhesive solids. It accounts for the electrostatic interactions such as the Van der Waals forces and other physical and chemical surface interaction effects, and, thus, its formulation depends on the material's surface energy [48].
In our work, we employ the simplified version of the JKR model (SJKR), in which the area increment caused by the formation of the neck is ignored, and the effective contact area is simply calculated as the intersection of the two particles. The cohesive force in the SJKR model is calculated as follows: where Ω is the specific cohesion energy density per unit volume (J/m 3 or Pa) and a 2 is the circular contact area. This cohesive force also implicitly represents the resistance of the inter-particle bond to shear and bending. Because the maximal tangential force F t and resistive torque r are proportional to the normal force, their magnitudes increase with the additional normal force due to the cohesion force |F coh | . The contacts configuration used in this study is summarized in the Supplementary Material (Table S1).

Simulation setup
The DEM simulations are performed using the open-source code LIGGGHTS © [51]. The cone geometry is replicated at true scale. The cone has a diameter of D cone = 10 mm . The ice samples are cylinders having a diameter of D cyl and a height of H = 130 mm , which is approximately (10) |F coh | = Ω a 2 equal to the height of the experimental samples. The cone was meshed using a non-uniform grid, where the average mesh size is 0.2 mm around the cone tip and 2 mm around the cone base and at the shaft. The mesh size was chosen to be approximately an order of magnitude smaller than the particle size to ensure an accurate transfer of forces between the particles and the solid. The particle diameter in our study ranged from 1.5 mm to 8 mm . The visualization of the simulations was performed using ParaView (version 5.2.0) and the postprocessing using a custom code implemented on MATLAB © (version r2018b). The simulations were run on Lonestar5 -a high-performance supercomputer at the TACC facilities (Texas Advanced Computing Center)-and took advantage of LIGGGHTS's parallel processing capabilities. To this end, an MPI (Message Passing Interface) framework was used to compile a parallel version of LIGGGHTS. Up to four nodes, each with two 12-core (Xeon E5-2690 v3 -Haswell 2.60 GHz) processing cores, were used for the various simulations. On this setup, a typical simulation takes about two hours.
To simulate the cone penetration test, the ice particles are poured into a cylinder from a height approximately equal to the final height of the sample. The particles are allowed sufficient time to settle under terrestrial gravitational conditions. The cone is then driven at a constant velocity of 100 mm.s −1 towards the bottom of the sample. It should be noted that, due to computational limitations, the penetration velocity was increased by a factor of ten relative to the experimental penetration speed. The simulation is terminated when the cone reaches 20 mm before the bottom of the sample. Fig. 2 shows a snapshot of the simulation at t = 0 s and t = 0.5 s cut along the middle of the sample.

Model input parameters
In our study, we varied six parameters, namely the particle's Young's modulus E p , the friction coefficient , the rolling friction coefficient r , the cohesion energy density Ω , the particle radius R, and the sample size D cyl . The parameters' range and all other model input parameters are summarized in Table 1. , r , Ω , E p are the degrees of freedom of the model. The range of particle properties is taken to be 1000 to 10 times lower than the bulk ice modulus owing to the peculiar microstructure of our porous ice [25,63,76]. The model's sensitivity to these parameters is investigated extensively in Sect. 4.3. The particle radius R and the domain size D cyl are the hyper-parameters of the model and are investigated extensively in Sect. 4.1.
To reduce the complexity of the model, we equate the particle-cone friction coefficients to the particle-particle friction coefficients, i.e. pp = pc and r,pp = r,pc . The particle's density was derived from the density of homogeneous polycrystalline isotropic ice at phase Ih . Even though the density of ice decreases slightly with temperature by ∼ 0.1 kg.m −3 ∕K [39], for simplicity, we assume a constant density at the reference temperature (i.e., T = 233 K). The Poisson's ratio p was taken as the average over our experimental temperature range [24,25,61,63,76,80]. The particle-particle coefficient of restitution e pp and particle-cone coefficient of restitution e pc were obtained from previous studies on collision between, respectively, two ice particles [38], and an ice particle and an ice block [37].  [16] to measure the cone penetration resistance of sintered ice samples. Dimensions are in mm b (left) DEM simulation of a cone penetration test of sintered ice at t = 0s (right) Cut-off view at the sample median at t = 0.5 s showing the particles in contact with the cone along with the cone dimensions. The cone geometry is replicated at scale. In the presented case, the sample diameter is D cyl = 100 mm and the particle radius is R = 2 mm Table 1 also reports the range of time step values used in the Velocity Verlet integration scheme. As the simulation time step is affected by both the particle properties and the particle size, the time step is adjusted for each simulation to match the critical time step to ensure algorithmic stability. The calculation of the critical time step is discussed in detail in the Supplementary Material (Text S3).

Model response variables
To study the response of the ice sample, we measure the force reaction on the cone upon penetration, which, according to Newton's third law, is equal to the force required to drive the cone at the given velocity. The force is extracted along three directions, with the Z-axis aligned with the cone axis. The X-and Y-axes are orthogonal to each other but arbitrarily oriented relative to the sample due to the symmetry of the problem. The stress components x , y and z are derived by dividing the components of the force vector by the largest cross-sectional area of the cone ( A cone = 78.5 mm 2 ). Figure 3a shows a typical stress response of a cone penetration simulation. We note that the stress components along the X-and Y-axes are not precisely equal, meaning that the problem is not perfectly symmetrical. This is primarily due to the granular nature of the material and the non-symmetry in the particle configuration.
The stress results from the top three centimeters starting of the samples are excluded from the analysis. It is suspected that the derived strength in this region is not representative of the sample strength since the tip influence zone has not yet fully developed [73]. Similarly, the stress results from the bottom three centimeters of the sample are excluded from the analysis. In this region, the tip influence zone interacts with the bottom of the sample, which distorts the strength measurements. The stress results in the remaining proof region are the most representative of the true strength of the specimen. The tip influence zone in this region is considered to be nearly undistorted and uniform. Figure 3a summarizes the extent of these regions.
To reduce the dimensionality of the model response, we extract a number of summary statistics from the stress profiles, namely the cone penetration resistance i , the dispersion factor s i , the penetration energy density E i along the three axes, and the strength-depth correlation factor . These metrics can be calculated as follows: where i,k is the stress component along the i-axis measured at depth z k . Figure 3b provides a visual depiction of the extracted summary statistics. i is indicative of the average sample strength, s i of the dispersion of the stress, E i is of the total energy required for the cone to penetrate through the sample, and the coefficient is indicative of any depth dependency in the stress profile.
= Corr( z , z) The coefficient takes values between − 1 and 1. A zero value corresponds to a plateau in the cone penetration strength profile. A positive value reflects a depth-strengthening behavior, while a negative reflects a depth-weakening behavior. The higher the absolute value of , the more pronounced the depth dependence.

Sample dimensions
Due to computational limitations, DEM samples are typically downscaled compared to their experimental counterpart, while particles are typically upscaled relative to the true particle size [31]. However, the selection of the downscaling and upscaling factors is not trivial and is usually performed heuristically. Here, we propose a method for rigorously selecting the optimal factors such that the numerical sample represents the true boundary conditions while entailing a relatively low computational cost. We simulate different sample configurations with a range of relative sample diameters D cyl ∕D cone and relative particle diameters D cone ∕d , as shown in Fig. 4. The sample diameter D cyl is expressed relative to the cone diameter D cone as it is indicative of the extent of the boundary influence zone. An infinitely large bed effectively mimics the true sample boundary conditions and, thus, large sample diameters are desired. The cone diameter D cone is expressed relative to the particle diameter d as it is indicative of the number of active particle-cone contacts. The total number of particles N needed for each configuration can be estimated using the formula N = 3 2 (1 − )D 2 cyl H∕d where is the sample porosity and H is the sample height. In this work, the sample porosity and height are taken to be equal to their experimental values, i.e., = 0.42 and H = 130 mm . With the available computation resources, several particles ranging from ∼ 300 to ∼ 65, 000 can be simulated in 3 minutes to 100 hours.
The computational time per node is defined as the time taken to complete a simulation multiplied by the number of nodes used. provides a measure that can be compared across simulations under some assumptions of equivalence, i.e., constant simulation performance, linear speed-up, little overhead time, and 100% parallelizable instructions. tends to increase with the number of particles, which itself increases with both the particle diameter and the sample diameter. Hence, there is a trade-off between the particle diameter and sample size, as illustrated by Fig. 4. We note that the physical configuration corresponds to an estimated particle count of N = 170.10 9 and an extrapolated computation time of = 3.10 11 h -a prohibitively expensive cost that reaffirms the need for scaling factors. To select the optimal configuration, we qualitatively analyze the sample displacement fields. Figure 5 shows the displacement map of the particles in the vertical direction z (left) and in the radial direction r (right), averaged azimuthally around the main axis. We observe that lower relative sample diameter leads to higher vertical displacement, which the law of mass conservation can explain. No significant change in the vertical displacement map can be seen for values higher than D cyl ∕D cone = 10 , indicating that this value might be an appropriate threshold. We also observe that finer particles lead to finer radial displacement maps, although no significant change can be seen beyond D cone ∕d = 2.5 , indicating that this ratio might be an appropriate choice.
To formally study the trade-off between sample size and computational time, we perform a Pareto analysis relating the displacement at the boundary to the computational cost, as shown in Fig. 6. The Pareto front connects the efficient configurations selected in such a way that no one objective can be improved without sacrificing the other. The objectives in our study are minimal displacement at the boundary and minimal computational cost. The displacement at the boundary is calculated as the average total displacement in the outer 10mm layer of the sample. The computational cost, measured in service-units (SUs), represents the computational resources needed in terms of time and number of cores. The Pareto optimal configuration correspond to the point closest to the origin of the graph shown in Fig. 6. By inspection and considering the available computational resources, we select D cyl ∕D cone = 10 and D cone ∕d = 2.5 . Henceforth, a sample diameter D cyl = 100 mm and a particle radius R = d∕2 = 2 mm are used in all simulations.
Finally, while varying the particle size in the above analysis, we observed a decrease in penetration energy density with decreasing particle size. Nevertheless, this effect does not impact our results and an in-depth study of its causes goes beyond the scope of this work.

Sample bed randomness
The sample bed is produced by pluviation, which has an inherent random component. The particles can assume a virtually infinite number of packing configurations while statistically retaining consistent macroscopic properties. This randomness implies a non-unique stress profile despite holding the same input parameters. To evaluate the effect of packing randomness on the model response, we compare the results of ten sample beds having arbitrarily different particle arrangements but the same input parameters. The results are summarized in Fig. 7.
We observe that the spread in the model outputs, measured by the interquartile range IQ, is comparable. In particular, we note that the interquartile ranges are smaller than the experimental confidence bounds derived from the laboratory tests [16], allowing for the use of a single representative numerical sample with given parameters for each experiment. Fig. 4 Matrix of relative sample dimensions D cyl ∕D cone as a function of the relative particle dimensions D cone ∕d . represents the computational time per node. Starred values are the estimated computation time derived from an exponential fit of the executed simulations. N refers to the number of particles in the sample, and starred values are estimated numbers. The case corresponding to the physical condition is reported in the top-right case (particle diameter ∼ 25 μm, sample diameter ∼ 160 mm). The color scale represents the computational time, where red cases correspond to simulations needing considerable computational resources that exceed the available resources for this study. Green to orange cases represents the computationally feasible cases with a varying need for resources To obtain a translatable measure of the dispersion, we calculate the coefficient of variation, which is expressed as C v (x) = s x ∕x , where s x is the standard deviation and x is the mean of the metric x. We obtain C v ( z ) = 4% , C v (s z ) = 14% , and C v (E z ) = 9% , indicating that the effect of sample bed randomness is limited. The coefficient of variation is not applicable to as it is a bounded measure (i.e., −1 < < 1 ). We take IQ as an estimate of the uncertainty associated with the sample bed randomness, which also indicates a limited spread relative to the experimental confidence bounds. We can conclude that the stochasticity introduced by the sample bed configuration will not significantly affect the subsequent model calibration and the use of a single representative sample bed is sufficient.
This limited effect of sample bed randomness has the added benefit of considerably reducing the simulation computational cost. Only one simulation is needed for each parameter set, and the associated confidence bounds accounting for sample bed randomness can be easily extrapolated. This reduces the number of simulations by a factor of approximately ten. Furthermore, creating the sample bed can take as long as running the virtual cone penetration test. Fig. 5 Matrix of the particle displacement field as a function of relative sample dimensions D cyl ∕D cone and relative particle dimensions D cone ∕d . In each panel, the left map represents the vertical displacement z while the right map represents the radial displacement r averaged azimuthally around the main sample axis Fig. 6 Computational cost versus average total displacement in the outer one centimeter layer of the sample. The displacement is calculated as = The dashed line is the estimated Pareto front, where optimal solutions correspond to the joint minimization of the objectives; displacement vs. computational cost Hence, using the same sample bed for all simulations further reduces the total computation time by a factor of two, allowing for a total reduction by a factor of about 20.

Sensitivity analysis
Before proceeding with the sensitivity analysis, we need to derive a lower-dimensional representation of the model output, i.e., a limited set of metrics that comprehensively capture the output characteristics. As seen in Sect. 3.3, ten variables are output from the DEM simulations. The degree of interdependence between the different model output parameters is measured by studying the correlation matrix, as shown in the Supplementary Material (Fig. S1). Output parameters with a high degree of correlation may be considered redundant, while output parameters with a degree of correlation close to zero may be considered orthogonal. We find that the penetration energy density E z is poorly correlated with the strength-depth correlation factor ( r = −0.17 ), indicating that and E z express very different characteristics of the ice. Henceforth, the output of the numerical test of any sample is described by the pair ( E z , ).
As described in Sect. 3.2, the DEM model has six input parameters (i.e., six degrees of freedom). In this section, we assess the impact of each input parameter on the model response in an attempt to reduce the model complexity without significantly affecting its flexibility and predictive capability. To this end, we conduct a sensitivity analysis in which the output response is determined by sequentially varying one input parameter while fixing all other parameters to nominal values. Fig. 8 shows the variation of the penetration energy density E z for a range of particle input parameters E p , , Ω , and r . Fig. 7 Box plots representing the distribution of the summary metrics of ten beds having arbitrarily different particle configurations but the same particle and interaction properties. The spread of the data points illustrates the variance in the model due to variations in particle arrangements.  Fig. 8a, we observe that the penetration energy density E z decreases with increasing Young's modulus E p , which can be explained by the fact that stiffer particles are harder to deform and thus have a lower contact area. Furthermore, since the cohesive forces are proportional to the contact area (see Eq. (10), they are also reduced, making it easier for the cone to penetrate the sample. Figure 8b shows that the penetration energy density E z increases monotonically with the friction coefficient , which can be explained by considering that more friction implies more energy for the particles to slide past each other. The strength appears to reach a plateau at low and high friction values. At low friction values ( ≈ 0.2), the strength is only a few kPa. In this regime, the particles can easily slide past each other and the strength is primarily determined by the resistance to rolling, which simulates the interlocking effect [88]. The relatively low contribution of the rolling resistance (few kPa) to the total sample strength (few MPa for a typical sample) supports the relative unimportance of the rolling friction parameter in our case. At high friction values ( > 2 ), the Coulomb threshold is high and particles cannot easily slide past each other. Instead, the particles' bonds rupture and the contact points are continuously reforming as the simulation progresses. In this case, the strength is limited by the bond strength and thus the cohesion value, which points to the relative importance of this parameter. Figue 8c shows a large influence of the cohesion energy density Ω on the model response, although the dependence is not monotonic. A higher inter-particle cohesion implies higher energy for the cone to penetrate through the sample. However, a peak strength is observed at Ω ∼ 80 MPa . Visualization of the simulations beyond this threshold indicates the formation of cracks which could explain the apparent decrease in strength beyond this limit, as shown in the Supplementary Material (Fig. S2). Finally, Fig. 8d shows the influence of the rolling friction coefficient r on the penetration energy density E z . The increase in strength for r < 1 can be explained by the fact that the particles require an increasing amount of torque to roll over each other. The maximum strength is reached when the rolling resistance becomes equal to the Coulomb friction limit i.e., r R * |F n |∕R = |F n | . In the case of homogeneous particle sizes, the equation leads to r = 2 . Beyond this limit, the particles will slip past each other before exceeding the torque limit, and thus the strength will decrease.
To quantify the influence of each parameter, we compute the parameter sensitivity matrix S k , whose formula, adapted from [89], reads as follows: where y = E z is the response variable and  Fig. 9.
Yan et al. [89] proposed to use the norm of the parameter sensitivity matrix norm(S k ) in DEM sensitivity studies to determine the most influential parameter. However, this metric is biased towards higher values due to the quadratic nature of the norm. Moreover, it does not provide information on the directional dependence of the response. As an alternative metric, we could take the mean value of the distribution < S k > which accounts for the directional dependence and gives all simulations the same weight, effectively solving two of the pitfalls of the norm metric mentioned above. The mean value metric, however, is sensitive to extreme values. As shown in Fig. 9, the distributions of E p , , and Ω are highly skewed, which biases the mean value. To mitigate these points, we propose an analysis based on the median value of the parameter sensitivity matrix median(S k ) . The median values, as shown in Fig. 9, are consistent both in direction and magnitude with the findings showcased earlier in Fig. 8. The friction coefficient entails a positive relationship where a 100% increase in leads to about 153% increase in penetration energy density. Similarly, a 100% Fig. 9 Quantitative sensitivity analysis of the model to its input parameters. The distribution of the sensitivity to each parameter experiment-wise is shown in the histograms. The X-axis represents the parameter sensitivity values S k , and the Y-axis represents the bin count 12 Page 12 of 20 increase in the cohesion energy density Ω leads to an average increase of about 36% in the penetration energy density. For Young's modulus E p , a 100% increase leads to about a 12% decrease in penetration energy density. The coefficient of rolling friction r seems to be the least influential parameter and is almost symmetrically distributed around − 8%, consistent with the non-monotonic relationship observed in Fig. 8.
We conclude from our sensitivity study that the coefficient of friction and the cohesion energy density parameter Ω are the most influential parameters and can, henceforth, be considered as the model parameters. Both Young's modulus E p and the coefficient of rolling friction r have a relatively low influence and can thus be held constant. For the subsequent simulations, we choose a value of Young's Modulus E p = 1 GPa, which is close to Young's modulus of snow [27] but 5-10 times lower than Young's modulus of isotropic polycrystalline ice [25,63,90]. While we expect Young's modulus of this type of granular ice to differ from snow and isotropic ice due to differences in microstructure, the low sensitivity values observed suggest that an exact estimate is not required. The coefficient of rolling friction is fixed to the value that maximizes the strength results, i.e. r = 1 , to emulate the interlocking effect of particle aggregates.

Model response maps
The model response is estimated by simulating a wide range of model input parameters. This step corresponds to the evaluation of the implicit function F that represents our DEM model: where i is the vector of model input parameters for the i-th simulation ( i , Ω i ) . We evaluate the function F at 110 parameter combinations distributed on a uniform grid. Then, the results are linearly interpolated to obtain a continuous map over the whole parameter space as shown in Fig. 10. We note that this interpolation is considered a "model of a model" and falls under the general framework of surrogate models [4,6,12,33,62,69,71,92]. We choose a linear interpolation for its simplicity and the relatively good approximation it provides for fine grids. Figure 10a shows the response map of the cone penetration energy density E z while Fig. 10b shows the response map of the strength-depth correlation factor . We see that the model can represent a broad range of ice conditions with varying strength values and depth-strengthening behavior. The modeled strength map, as seen in Fig. 10a, is regular and smooth with values ranging from ∼ 1 MPa (unconsolidated ice) to ∼ 13 MPa (heavily consolidated ice). The cone penetration energy density increases with increasing values of cohesion energy density Ω and friction coefficient , consistent with the expectation that stronger particle-particle interactions yield a higher macroscopic strength. We also note a decrease in penetration energy density beyond a critical value of Ω ∼80 MPa, which, as touched upon in Sect. 4.3, might be due to the formation of large cracks in the sample. As shown in Fig. 10b, the strength-depth correlation map is less smooth. In particular, at low Ω and high values, the model predicts a high depth-strengthening behavior. This observation is similar to the behavior exhibited by viscous liquids and is consistent with a loose frictional granular material. To better understand the depth-strengthening and depth-weakening behavior observed in our study (i.e. the value of the correlation factor ), future work will investigate the scaling laws of the penetration force with depth, similar to other work on impact in granular material [8]. These model response maps will subsequently be used in the calibration process to solve the inverse problem of finding the pair of model input parameters ( Ω, ) for a given output value ( E z , ) obtained experimentally.

Model calibration
DEM parameters are typically difficult to measure in the laboratory and cannot be easily related to measurable physical material parameters [35]. As a result, calibration is required to select the appropriate parameters for use in simulations.
Calibration is considered one of the most challenging steps in DEM modeling, and is now an active area of research [20]. It corresponds to solving the inverse problem of finding the model parameters ̃i that match the experimental bulk properties G exp i for the i-th experiment, as shown in the following equation: In previous DEM studies, G exp i have been typically obtained from biaxial and triaxial tests [42], angle of repose tests [72], and Brazilian tests [33,62]. Cone penetration tests are less commonly used for DEM calibration. They have only been used to calibrate non-cohesive materials in chambers with well-controlled boundary conditions [9,18,58]. Hence, calibrating this cohesive ice with cone penetration tests is an unprecedented task.
Multiple methods have been developed to solve the inverse problem in DEM [20,71]. However, they either require numerous experiments [30,35], a good a priori knowledge of the parameter dependencies [15,19,49,56,91], tend to find local minima rather than global optima [69], or tend to be computationally intensive [4,22]. Probabilitybased calibration methods like the sequential quasi-Monte Carlo [14], iterative Bayesian filtering [15], Transitional Markov Chain Monte Carlo [32], and Sequential Monte Carlo [13,14] are state-of-the-art and have been successfully used. However, they require expert knowledge and remain difficult to implement for many DEM users.
We propose a novel Bayesian parameter calibration method that optimally incorporates experimental uncertainties and uses a single calibration experiment. This method leads to a set of candidate DEM input parameters, together with their likelihood to yield the given experimental outcome. Considering the complete parameter space ensures that our method satisfies global optimality. To obtain a unique parameter set, we use a weighted mean approach where each parameter set is weighted by its likelihood. We run our method on each experiment independently, and estimate the corresponding posterior distribution on the model parameter space. The method makes use of the Bayes' theorem such that: where p(̃i|E exp z,i , exp i ) is the joint posterior probability distribution, p(̃i|E exp z,i ) is the posterior probability distribution on the parameter space given the experimental strength value E z , and p(̃i| exp i ) is the posterior probability distribution on the parameter space given the experimental strengthdepth correlation factor . This method assumes that the strength value and the strength-depth correlation factor are conditionally independent of the model input parameters. In other words, if the model input parameters were known, knowledge of the sample strength would not provide any information on its depth-strengthening behavior.
The uniqueness of our method lies in the estimation approach. We estimate the posterior distributions by solving the inverse problem described in Eq. (17) for a selected set of values within the experimental confidence bounds {G k |G k = G exp + k G exp , k = 1..N} , where G exp and G exp are, respectively, the mean and standard deviation of the experimental measurement of the property G, N the number of selected values in the set, and is a variable. We then search for each value G k on the response maps shown in Fig. 10, and record the associated parameter set {̃i ,k } . To each parameter ̃i ,k in the parameter set {̃i ,k } , we associate a probability p(̃i ,k ) equal to the corresponding value of the probability density function of a standard normal distribution ( k ) . For example, if we consider values at both ends of the experimental confidence interval, their associated parameters will be attributed a low probability value (i.e. a low weight), whereas for values at the center of the experimental confidence interval their associated parameters will be attributed high probability values (i.e. high weights). The algorithmic implementation for this procedure is presented in details in Supplementary Material (Algorithm S1).
The experimental strength measurements E exp z are represented by a Gaussian distribution with a mean and standard deviation equal to the experimental values, i.e., = E exp z and = s E exp z . For the correlation factor , we cannot assume a Gaussian distribution since is a bounded measure. We can transform the bounded correlation values to unbounded values using the Fischer transformation function F(r) = 1 2 ln( 1+r 1−r ) = arctanh(r), commonly used for correlation factors. The mean and standard deviation of the Gaussian distribution is equal to the mean and standard To assess the normality assumption, we compare the quantiles of the experimental and theoretical distributions in a Q-Q plot shown in the Supplementary Material (Fig.  S2). The results indicate a good match between the experimental and theoretical distributions overall, which suggests that the data reasonably satisfies the normality assumption and justifies the choice of the Gaussian distribution.

Comparison between experimental and simulation results
We first examine the cone penetration test results of ice samples that sintered at four different temperatures (193 K, 223 K, 233 K, and 243 K) for time periods up to 14 months [16]. Figure 11 shows the experimental strength-depth correlation factor as a function of the penetration energy density E z . We observe that the ice samples exhibit strength values up to 14 MPa and correlation factor values spanning almost the whole [−1, 1] interval. The data points are mostly scattered throughout the range of possible values, indicating that the ice can take a large number of physical configurations and have a rather large, unrestricted state space. A closer inspection of the data shows a cluster of very high positive strength-depth correlation factor values at penetration energy density below 2 MPa. Such high values of correspond to a linear increase in resistance with depth. For penetration energy density values higher than 2 MPa, the strength-depth correlation factor is uniformly distributed, which indicates an absence of apparent relationship between cone penetration resistance and depth. In addition, for E z > 2 MPa, the cone penetration profiles (e.g., Fig. 3) exhibit a jerky behavior. These findings are indicative of a brittle mechanical behavior for samples with resistance above 2 MPa, a value comparable to the brittle compressive strength of ice at strain rates greater than approximately 10 −2 s −1 [17]. The purple region in Fig. 11 corresponds to the ice states that are numerically reproducible with our current DEM model. It is derived from a smoothed representation of the combination of the model response maps shown in Fig. 10 and indicates that our DEM model can represent approximately 70% of the observed experimental states.
For each experiment presented as a data point in Fig. 10, we follow the calibration procedure described in Sect. 5.2. Figure 12 shows an example of the posterior probability distribution for an ice sample aged for 4 days at T = 243 K. The sample response is characterized by a cone penetration energy density E z = 3.06 ± 0.24 MPa (95% confidence interval) and a strength-depth correlation factor = −0.29 ( CI = [−0.62; 0.14] 95% confidence interval). Figure  In this particular example where the sample sintered for 4 days at T = 243 K, the most probable input parameters are in the lower-end of friction coefficient ≈ 0.33 and cohesion Ω ≈ 42 MPa, which is consistent with the low level of consolidation of the sample and its relatively young age. In the case of more consolidated ice, for example, a sample that aged for 28 days at T = 243 K developed a cone penetration energy density E z = 9.33 ± 2.28 MPa (95% confidence interval) and a strength-depth corrrelation factor = −0.69 ( CI = [−0.93; − 0.03] 95% confidence interval). Our analysis yields a most probable friction coefficient of ≈ 0.69 and cohesion parameter of Ω ≈ 87 MPa.  Figure 13 shows the evolution of the calibrated set of parameters as a function of sintering time and temperature for all tests performed. Fig. 13a and b show, respectively, the evolution of the calibrated cohesion energy density parameter Ω and the calibrated friction parameter , which both exhibit a clear temperature dependence. We constrain each parameter evolution at each temperature to the first order with a linear regression while weighting each parameter set by its probability. The dashed lines represent the best-fit trend lines.

Time evolution of particle-scale parameters
We observe from Fig. 13a and b that the evolution of the particle-scale parameters matches very well the evolution of the continuum-scale sample cone penetration resistance presented in [16]. The model parameters increase monotonically with time, and the evolution rates increase with temperature. This correspondence points to a link between the micro-and macro-mechanical properties and suggests that the evolution of the particle-scale interactions upon sintering is an important strengthening mechanism of the ice. Similarly to [16], we notice that the linear fits do not have the same intercept values with the Y-axis, which represents the parameters associated with a fresh ice sample. This may be an artifact of the simplistic linear evolution model assumed or due to a different strengthening rate at earlier sintering times.
To further investigate the temperature dependence, we extract the evolution rates of the model input parameters. Fig. 13c and d show the natural logarithm of the rates at which the cohesion energy density dΩ dt and the friction coefficient d dt evolve as a function of inverse temperature, also known as the Arrhenius plot. Such representation is widely used to study thermally activated processes. The evolution rates for each parameter seem to follow a linear trend, which supports that the evolution follows an Arrhenius law. The slope of the linear regression is given by −Q∕R , where Q is the activation energy of the underlying process and R = 8.314 J∕mol.K is the ideal gas constant. The linear regression yields an activation energy of For comparison, [16] derived an activation energy of Q = 24.3 ± 3.3 kJ/mol for the rate of strengthening of the experimental samples and found that it is consistent with diffusion of water molecules on the surface of ice particles. Within their confidence intervals, the activation energies representing the temperature evolution of both the cohesion energy density and the friction coefficient are broadly consistent with the experimentally-derived value for the strengthening rate. This suggests that the particlescale interaction parameters in the DEM model, despite larger uncertainties, reflect the underlying thermodynamics and represent well the effect of temperature that had been derived experimentally at the macroscopic scale.

Implications for mechanical properties from particle-to macro-scale
The correspondence between the calibrated model parameters and the physical properties of ice supports that the model parameters and Ω are representative of the physical micro-mechanical properties of the ice. The calibrated friction coefficient of a sample that sintered for 4 days at T = 243 K (i.e., ≈ 0.33 ) is in good agreement with published experimental measurements. [77] measured the kinetic friction coefficient of fresh ice to be = 0.29 ± 0.03 for a sliding velocity of V = 10 −3 ms −1 at a temperature of T = 223 K. Futhermore, the increase in the calibrated friction coefficient, as seen in the example above (from ≈ 0.33 to ≈ 0.69 in 24 days at T = 243 K), is consistent with Fig. 13 Evolution of the calibrated model input parameters a Cohesion energy density Ω and b Friction coefficient as a function of sintering time and temperature. The mean and standard deviation of the probability distribution of the calibrated parameters are shown. The dashed lines represent the weighted linear regression and are colorcoded with respect to temperature. Arrhenius plot of the calibrated model input parameters c Cohesion energy density Ω and d Friction coefficient as a function of inverse temperature increased inter-particle interactions due to sintering. The evolution of the model friction coefficient, as showcased in Fig. 13b, exhibits a clear positive trend and is reminiscent of static strengthening -a phenomenon observed in compact ice held under constant normal stress [77,78]. The cohesion energy density Ω is tightly linked to the material's surface energy , which has been established as the primary driver of ice sintering [5]. In fact, the JKR model, which represents various electrostatic, physical, and chemical surface interactions, expresses the cohesion force as a function of the material's surface energy (in unit surface). The SJKR model, which makes a geometrical simplification to the JKR model, expresses the cohesion force as a function of Ω (in unit volume). Furthermore, since increasingly sintered material require an increasingly high pull-out force to separate the particle, and since F pull−out = F coh = Ω a 2 , enhancement of sintering implies higher Ω values -a feature that is showcased in our results on the evolution of the model cohesion energy densities (Fig. 13a).
This modeling study along with prior works suggest interrelationships between the particle-scale interaction parameters, the evolution of a mesoscopic network of ice, and the macroscopic mechanical properties that govern the brittle compressive failure of ice. As discussed in Sect. 6.1, the evolution in cone penetration resistance of macroscopic samples of ice microspheres suggested a transition in mechanical behavior from loose unconsolidated ice particles at low cone penetration resistance values ( < 1 MPa ) to brittle compressive failure for values around 2 MPa or greater. This consolidation is not accompanied by any discernible change in bulk porosity, implying that sintering may be responsible for the development of an inter-particle network that binds the starting ice particles together. Once this network is established, the ice particles form consolidated aggregates which behaves as a coherent medium and fails in the brittle regime when compressed at strain rates greater than approximately 10 −2 s −1 [17].
Interestingly, the combination of friction and cohesive energy and its apparent dominant role in the brittle failure of porous ice has also emerged from experiments on the brittle compressive failure of fully-dense ice. There, the frictional-sliding wing crack mechanism [3,29,40,57] accounts quantitatively for the behavior [11,28,70,75,87]. Accordingly, sliding across the opposing faces of closed, parent/primary cracks that are inclined to the principal loading direction induces tensile stress at the crack tips. When sufficiently high, tension is relieved through the initiation of out-of-plane secondary cracks termed wings (or extensile cracks). As sliding continues, the wing-crack mouths open, thereby increasing the mode-I stress intensity factor at their tips until a critical level is reached, at which point the wings begin to grow in a stable, albeit jerky, manner. As sliding continues, the secondary/wing cracks lengthen, interact with other secondaries and eventually form a fault at which point the material collapses. In this model, resistance to sliding is set by the coefficient of friction and resistance to wing crack growth is set by fracture toughness to which surface energy, and hence cohesive energy, is a major contributor in ice [76, pp. 207-208].
That two different approaches-DEM of sintered ice and experiment-cum-physical modeling of pore-free ice-point to the underlying role of the same physical processes suggests that an intermediate level of porosity, while increasing microstructural complexity, may not change the underlying physics of brittle compressive failure.

Conclusion
In this study, we present a physics-based numerical model of planetary ice plume deposit analogs that explicitly represents the microstructure and its evolution upon sintering. We calibrated our DEM model using 100 published experiments of cone penetration tests on ice plume deposit analog samples that sintered at different temperatures for time periods up to 14 months. We tuned the sample dimension and particle size ( D cyl , d) following a Pareto-optimality approach and investigated the effect of the particle and bond parameters ( E p , , r , and Ω ) using a sensitivity analysis. We found that the friction coefficient and the cohesion energy density Ω are the primary model parameters. Furthermore, we proposed a novel easy-to-implement Bayesian probabilistic calibration method for numerically replicating experimental conditions while optimally incorporating experimental uncertainties.
Our DEM model has shown capable of reproducing experimental porous ice strengthening results and representing the physical micromechanics of ice and their evolution upon sintering. The evolution of the bond strength matches that of the macroscopic sample, suggesting that ice sintering is responsible for an increased interaction at the particlescale which manifests itself by the formation of a mesoscale network structure that mechanically behaves in a manner akin to fully-dense ice.
Our findings show that this methodology can provide a critical link between theoretical and experimental studies, allowing us to better understand the effect of sintering on the mechanical properties of plume deposits on icy worlds. In the future, this model could be used to study ice features that are impossible or extremely difficult to observe experimentally, such as the force network, ice fabric, micro-dynamics, and fracture behavior. In this regard, it would enable a more detailed exploration of ice mechanics and, potentially, discovering new fundamental insights about this unique icy material.
This numerical model approach can also be expanded to robot-terrain applications. Such models are beneficial for engineers optimizing the design of planetary exploration robots and sampling systems. For this purpose, it would need to be interfaced with multi-body dynamics software to simulate complex robotic interactions such as traversing or sampling. As our model is tuned to cone penetration tests, its performance should be first evaluated on other loading conditions, then generalized to accommodate realistic robotic systems, such as landing pads. Finally, the model can be extrapolated to Ocean Worlds surface conditions, accounting for the reduced gravity, reduced pressure, plume deposition rate, and differential sintering with depth. This extrapolation would allow the prediction of surface conditions on icy worlds and the identification of suitable landing and sampling sites.