Uncertainty analysis of single-and multiple-size-class frazil ice models

. The formation of frazil ice in supercooled waters has been extensively studied, both experimentally and numerically, in recent years. Numerical models, with varying degrees of complexity, have been proposed; these are often based on many parameters, the values of which are uncertain and difﬁcult to estimate. In this paper, an uncertainty analysis of two mathematical models that simulate supercooling and frazil ice formation is carried out within a probabilistic framework. The two main goals are (i) to provide quantitative insight into the relative importance of contributing uncertain parameters, to help identify 5 parameters for optimal calibration, and (ii) to compare the output scatter of frazil ice models with single and multiple crystal size classes. The derivation of single-and multi-class models is presented in light of recent work, their numerical resolution is discussed, and a list of the main uncertain parameters is proposed. An uncertainty analysis is then carried out in three steps. Parameter uncertainty is ﬁrst quantiﬁed, based on recent ﬁeld, laboratory and numerical studies. Uncertainties are then propagated through the models using Monte Carlo simulations. Finally, the relative inﬂuence of uncertain parameters on the 10 output time series - i.e. the total frazil volume fraction and water temperature - is assessed by means of Sobol indices. The inﬂuence of input parameters on the long-term asymptotic as well as short-term transient evolution

The main drivers for the formation of frazil are the water cooling rate resulting from heat exchanges with the atmosphere, the initial seeding of frazil nuclei and the turbulent mixing.Then, the thermal growth process deriving from the heat exchange between water and primitive crystals, allows a fine description of the balance between growth frazil ice and water supercooling.
Previous mathematical models describing the evolution of frazil ice and water temperature, were based on ideas pioneered by Daly (1984Daly ( , 1994)).Omstedt (1985) developed a model based on a turbulent channel-flow boundary-layer theory, in which a mean particle is considered of constant geometric properties and constant Nusselt number is considered.Mean particle-size models have been incorporated in numerical hydraulic tools such as CRISSP2D (Shen and Wasantha Lal, 1993;Shen et al., 1995;Shen, 2010).More complex models are based on the ice-number continuity equation introduced by Daly (1984), taking into account the complexity of crystal distribution.In contrast with single-size-class models, the crystal number continuity equations allow introduction of secondary nucleation and flocculation processes.When combined with thermal growth, these complement the modeling of frazil crystal size evolution.Svensson and Omstedt (1994) proposed a numerical model for solving the equations introduced by Daly (1984) by considering discrete radius intervals.Conceptual models for secondary nucleation and flocculation were introduced and calibrated to fit the experimental data of Michel (1963) and Carstens (1966).Hammar and Shen (1991) derived a variable particle-size model in two dimensions, and later improved it with secondary nucleation and flocculation (Hammar and Shen, 1995).Variable particle-size models has also been integrated in TELEMAC-MASCARET (Souillé et al., 2020).Numerous simplified implementations have been made that consider a well-mixed water body.For example, Wang and Doering (2005) worked with the same model as Svensson and Omstedt (1994) and further discussed calibration of influential parameters in initial seeding, secondary nucleation, and flocculation by comparing them with the experimental data of Clark and Doering (2004).Implementation of multiple size class frazil dynamics in the context of sea ice and ice shelf water plumes has also been proposed by Smedsrud (2002), Smedsrud and Jenkins (2004) and Holland and Feltham (2005).In the same context, Rees Jones and Wells (2018) have recently shed new light on multiple-size-class models and discussed crystal growth rate, and the occurrence of frazil explosion.They also identified and characterized steady state crystal size distribution.Henceforth, single-and multiple-size-class models will be referred to as (SSC) 1 and (MSC) 2 models respectively.
Common to all numerical model of frazil dynamics is their use of a large number of parameters, making calibration difficult.
The modeling studies mentioned above show that frazil ice models can be fitted to reproduce evolution of temperature and frazil volume fraction.However, given the uncertainty of fitting parameters, it is questionable whether these models are predictive.
As Rees Jones and Wells (2018) point out, the consistency between experiments and models does not necessarily means parameterization is correct.This is particularly true when different processes compete and have similar impact on the crystal size distribution.Besides, introducing new processes in the models increases the number of parameters that need calibration, which : .:::: This comes at the cost of introducing new uncertain parameters , and raises the question of models trustworthiness.
In this paper, we intend to bridge that gap by proposing an uncertainty and sensitivity analysis of single and multipleclass frazil ice models, using a of full description of parameters uncertainty (by probability density functions) and variance decomposition of model outputs.First, we introduce the models, and main hypotheses and discuss their implementation in section 2. The uncertainty assessment methodology is presented in section 3. We then rely on numerous field and experimental measurements to quantify the uncertainties of input parameters of the two frazil numerical models in section 4. Using Monte Carlo experiments, we next study propagation of these uncertainties on model outputs in section 5: in other words, the evolution of the frazil ice volume fraction and the water temperature.The evolution of output scatter is discussed and compared to the asymptotes of the dynamic systems.Finally, sensitivity analysis based on Sobol indices (Sobol, 2001) and aggregated Sobol indices enables us to propose a selection of the most influential parameters in both models.We conclude about this study's perspectives in section 6.

Frazil ice dynamics
This section introduces the continuum equations used to describe evolution of water temperature and frazil volume fraction, as well as their discrete counterparts and their numerical resolution under the assumption of well-mixed water bodies.The parameters of the models, that can be considered uncertain, are also introduced.

Mathematical description
Frazil crystals are supposed to be dics ::::::: assumed ::: to :: be ::::: discs : of radius r and thickness e : λ.The ratio between diameter and thickness, denoted R = 2r/e :::::::: R = 2r/λ : is supposed to be constant as crystals grow.Let us define n as the crystal number density, which corresponds to the number of crystals per unit volume per unit length in radial space.The total number of particles per unit volume is then defined as N = ∞ 0 ndr ::::::::::::: N = ∞ 0 n(r)dr.We also introduce the frazil volume fraction density as c = nV in which V = πr 2 e :::::::: V = πr 2 λ : is the crystal volume.The volume proportion of the water-ice mixture occupied by frazil ice, i.e., the total volume fraction of frazil, is then defined as C = ∞ 0 cdr ::::::::::::: C = ∞ 0 c(r)dr.An incompressible water-ice mixture of velocity u and depth h is considered.Following Daly (1984) the number density balance equation can be defined as (1) in which ::::: where ν c is the frazil particles diffusivity, w r is the buoyancy velocity, and remaining right hand side terms are source terms described hereafter: -The source term (1) represents the crystal size evolution due to thermal growth, resulting from the heat exchange between the supercooled water and ice particles.The crystals can be considered to grow mainly on their peripheral area, defined by a = 2πre ::::::: a = 2πrλ : (Daly, 1994).Frazil evolution is then driven by the radial growth rate G, defined as :::: equal :: to : in which ρ i is the frazil ice density, L i is the latent heat of ice fusion and the right hand side is the heat exchange between the crystal of temperature T i and the surrounding water of temperature T .The latter is modeled as a function of the temperature delta ∆T = T i − T , the thermal conductivity of water k w , the Nusselt number N u :: N u , which represents the ratio between turbulent heat transfer and conduction heat transfer, and a thermal boundary layer length scale, denoted δ T ::::: (which :: is ::::: either :::::: chosen :: as ::: the :::::: radius :: or :::::::: thickness :: in ::::::::: literature).The crystal temperature T i is assumed to be equal to the freezing temperature denoted T f .The Nusselt number can be described as a function of ratio m * = r/η, where η is the Kolmogorov length scale, which can be defined as a function of the turbulent dissipation rate ε as η = (ν 3 /ε) 1/4 , in which ν is the molecular viscosity of water (Daly, 1984).The Nusselt number formulation described by Holland et al. (2007) is used.For small particles, :: i.e. ::::::: m * ≤ 1, : heat transfer is governed by diffusion and convection, and the Nusselt number can therefore be written as : where P r :: P r is the Prandlt number, defined as the ratio between molecular and thermal diffusivity.For larger particles (m * > 1), heat transfer is governed by turbulent mixing of the boundary layer around the crystal, and the Nusselt number is defined by where α T is the turbulent intensity.
-The source term (3) of Equation ( 1) is a flocculation source term supposed to represent the net effect of both flocculation and breakup, as introduced by Svensson and Omstedt (1994), who chose F = F 0 r 2 , where F 0 is a constant.mechanisms.::::::::: Moreover, :: it : also depends on frazil concentration and hydrodynamic conditions, of which the effects on flocculation are still poorly characterized in literature.To be consistent with models in the literature, the formulation proposed by Svensson and Omstedt (1994) is chosen for this paper.
-The source term (4) of Equation (1) represents the buoyancy of frazil ice crystals, where w r is the rise velocity.In the present paper, it is set to w r = 32.8(2r) 1.2 , following Svensson and Omstedt (1994), who simplified the model proposed by Daly (1984).For consistency with previous modeling studies, we retain this simple formulation, even if many other formulations have been proposed as summarized by Morse and Richard (2009) and McFarlane et al. (2014).
Finally, the thermal balance of the water-ice mixture complements Equation (1).Supposing :::::::: Assuming : C 0, the water fraction temperature equation is obtained : where T is the temperature of the water fraction of the water-ice mixture, ν t : is : the thermal diffusivity, ρ the density of water, c p the specific heat of water and φ the net heat source resulting from heat exchanges with the atmosphere in W.m −3 .
In the following sections, a well-mixed water body is considered.Equations ( 1) and ( 6) are written in terms of spaceaveraged quantities.This assumption allows us to neglect the convection and diffusion operators and to focus on solving the average temperature and average frazil volume fraction.Therefore, the partial differential equations of frazil and temperature are simplified to a coupled set of Ordinary Differential Equations (ODEs) including only the source terms.

Radial space discretization
In this section we introduce a multiple-size-class (MSC) frazil ice model, which is a discrete version of Equations ( 1) and ( 6) in radial space.It consists of considering m classes of constant radius chosen between a minimum and a maximum radius (Svensson and Omstedt, 1994).For each class i, the radius and the thickness are supposed to be equal to r i and e i = 2r i /R ::::::::: The peripheral area as well as the surface and volume of frazil crystals are defined as a i = a(r i ) = 2πr i e i , s i = s(r i ) = 2π(r i e i + r 2 i ) and V i = V (r i ) = πr 2 i e i ::::::::::::::::: a i = a(r i ) = 2πr i λ i , ::::::::::::::::::::::: ::: and ::::::::::::::::: V i = V (r i ) = πr 2 i λ i , respectively.A log-Uniform discretization of the radial space is chosen as in previous studies (Svensson and Omstedt, 1994;Wang and Doering, 2005;Rees Jones and Wells, 2018).The number of crystals of class i is noted n i , and similarly, the volume fraction of crystals of class i is noted c i .The total volume fraction of frazil, and total number of particles, are then defined as C = m i=1 c i and N = m i=1 n i with c i = n i V i , respectively.The discrete version of Equation ( 1) written in terms of volume fraction balance reads: with the boundary conditions The discrete versions of the source terms introduced in Equation 1 are defined following previous work (Svensson and Omstedt, 1994;Wang and Doering, 2005;Holland and Feltham, 2005): -The thermal growth source term (1) is composed of thermal growth and melt functions Γ i and Λ i , which are defined as : where G i = G(r i ), ∆V i = V i+1 −V i and H = He(T f − T ), with He :::::::::::::: H = H e (T f − T ), :::: with ::: H e : the Heaviside function allowing to a switch between melting or freezing.We suppose :::::: assume that frazil crystals grow from their peripheral area a i but melt from their surface s i (Holland and Feltham, 2005).
-The flocculation source term (3) is defined as -The buoyancy of crystals is simplified into a gravitational removal sink term (4) defined as γ i = −w r (r i )a d /h, in which h is the water depth.We introduce a coefficient a d , to account for the uncertainty of the rise velocity and gravitational removal process.
Finally, writing the thermal growth source term as , one can write the discrete version of Equation ( 6) as A general discrete frazil ice model is implemented in the present work (Equations 7 and 12) so that one can easily retrieve the formulations developed in previous papers.For example, supposing H = 1, a d = 1, τ s = 0 and writing c i = n i V i , Equation ( 7) is equivalent to Equation (4) in Wang and Doering (2005) and Equation (1) in Svensson and Omstedt (1994) with ζ = 1.
By also neglecting flocculation (a f = 0), it is equivalent to Equation (10) in Rees Jones and Wells (2018).

Numerical methods to solve governing equations
Let us denote t k = t 0 + k∆t the time at iteration k, t 0 the initial time and c k = c(t k ).A semi-implicit theta scheme is implemented for the MSC model, with a constant time step ∆t .It consists of an Euler explicit time discretization of the left-hand side derivative, while taking c i = θc k+1 i + (1 − θ)c k i in the right-hand side.Choosing θ = 0 leads to a fully explicit time scheme while choosing θ = 1 is equivalent to an implicit scheme on c.The non linear terms are treated semi-implicitly i.e.
An important aspect of the numerical frazil ice models presented in this paper is the need to provide a non-zero initial condition for the frazil volume fraction (in absence of seeding: i.e. τ s = 0).In the case of MSC models, Svensson and Omstedt (1994) assumed a uniform initial number of particles n 0 in each radius class i.e. n i (t 0 ) = n 0 (1 ≤ i ≤ m).We followed the same principle to initialize the system, but fixed the number of initial particles at zero for classes with a radius exceeding a threshold r 0 : i.e. n ri≤r0 (t 0 ) = n 0 ::: like ::::::: proposed ::: in :::: other :::::: works ::::::::::::::::::::::::::::::::::::::::::::::: (Holland and Feltham, 2005;Rees Jones and Wells, 2018).To be able to compare SSC and MSC models, we worked in terms of initial volume fraction of frazil C(t 0 ) = C 0 .The system was then initialized with As the initial setup has a big :::::::: significant influence on the results (Holland and Feltham, 2005), C 0 and r 0 are considered in the following sections as uncertain parameters. :: (a) To test convergence in terms of the number of radius classes, we make it vary from m = 10 to m = 1000 ::::::: m = 200.The results, presented in Figure 1, are consistent with observations by Rees Jones and Wells (2018) that the system requires m 100 to converge.Note that Svensson and Omstedt (1994) and Wang and Doering (2005) took m = 20 and m = 40, respectively.
It should also be noted that convergence in the number of classes depends on the initialization method of the system, as highlighted by Holland and Feltham (2005).To perform all the numerical simulations required for a complete sensitivity analysis, m = 100 was chosen as a trade-off between numerical convergence and computational cost.

Study cases and model parameters
In the present study, we focus on the evolution of water temperature and total frazil volume fraction in a supercooled, well mixed water body of depth h = 1 m.To focus on the frazil modeling rather than heat budget, uncertainties deriving from exchanges with the atmosphere are neglected and the cooling rate φ is considered deterministic and constant over time in all experiments.Nonetheless, different values of φ were tested, ranging from −50 to −1000 W.m −3 , to test the variability of the results in different cooling rate situations.As described in previous sections, frazil ice models are driven by many parameters, some of which are subject to a significant degree of uncertainty.The list of parameters considered in the present study for conducting the uncertainty analysis is shown in Table (1).Some physical properties are considered constant and are taken at Taking the SSC model as a reference, one can simply visualize the main feature of frazil dynamic systems.When time is close to zero, the initial temperature decrease rate is defined by lim t→0 Ṫ = φ/ρc p .In the absence of seeding and gravitational removal (i.e.τ s = 0 and a d = 0), the temperature decreases to the maximum supercooling point, characterized by a maximal temperature depression, denoted θ = max(T f − T ), and a critical time, denoted t c .After the maximum supercooling, the frazil production rate releases more heat compared to that released in the atmosphere, leading to an increase in water temperature which tends toward freezing point.Finally, when time tends to infinity, there is a balance between the frazil growth source term and the cooling rate.This leads to a linear increase in frazil volume fraction at a constant rate i.e.
In the present study, we analysed the uncertainty of input parameters in the main transient and steady-state features of the frazil ice models with and without gravitational removal.We started the simulation with a water temperature equal to the freezing temperature (which is supposed to be equal to zero in the absence of salinity : i.e.T f = 0 • C).We checked that all the simulations ran beyond the critical time.For the range of parameters considered there, it was found that choosing a final time of t f 1h was enough to capture the transient features of the ODE systems.
Note that both the SSC and MSC models are able to reproduce the experiments of Michel (1963), Carstens (1966) and (Clark and Doering, 2004) as shown by Svensson and Omstedt (1994) and Wang and Doering (2005).In Figure (2) we give an example of the results of both SSC and MSC models compared to Carstens (1966) experimental results (Case I).

Probabilistic framework for uncertainty analysis
The objective of this section is to describe the mathematical tools necessary to study the uncertainties of SSC and MSC frazil ice models.An uncertainty study of a numerical model can be performed in three steps (Sudret, 2007).The first step consists in identifying the uncertain parameters, and characterizing the probabilities of occurrence for their values through Probability Density Functions (PDFs), which is referred to as quantification of uncertainty sources.The second step concerns the propagation of uncertainties through the interest models, generally using sampling techniques (e.g., Monte Carlo) to obtain possible values of the target outputs (here frazil ice concentration, water temperature) and compute their statistical estimates (mean, standard deviation, percentiles, etc.).The third step, called sensitivity analysis, focuses on ranking the uncertain parameters in terms of influence on the target output.In the following, formulations for the statistical estimates and sensitivity analysis indices are given, with a summary of elements essential to understand this paper.Curious readers can find theoretical details ::::::::: Theoretical :::::: details ::: can ::: be ::::: found : in Soize (2017) and Sudret (2007).All uncertainty analysis computations were performed with the OpenTURNS library (version 1.18) developed by a partnership between Airbus Group, EDF R&D, IMACS, ONERA and Phimeca Engineering (Baudin et al., 2016b).

Uncertainty quantification
Let X = X 1 , . . ., X n X be the vector of uncertain parameters of the frazil ice model, where n X is the number of uncertain parameters.Let f : g be the interest frazil model, either SSC or MSC models presented in Equations ( 14) and ( 12) for temperature and ( 13) and ( 7) for frazil volume fraction.The output of the interest models are the frazil volume fraction and water temperature discrete time series i.e.T (t k ) and C(t k ) (1 ≤ k ≤ n t ), which we will generally denote Y = Y 1 , . . ., Y nt below for the sake of simplicity.Variables ::: The ::::::: random :::::: vectors X and Y are linked through frazil models f such that Y = f (X, d) : g :::: such ::: that ::::::::::: Y = g(X, d), where d is a deterministic vector, i.e. fixed parameters in contrast with the uncertain set of inputs X.
To undertake the uncertainty studies, both the inputs X and the interest outputs Y are considered to be random variables.
This means, taking the example of :::::: vectors.::: We ::::::: assume :::: that X , that they are characterized with probabilities of occurrence for given values x denoted P(X = x).These probabilities link to the PDF denoted p X through ::: has : a ::::::: density ::: p X , ::: so :::: that , where E is a subset from the space of all possible values D X .Each element X i and Y k of the input and output vectors is hence characterized by a PDF.
The results of the uncertainty analysis are directly linked to the UQ ::::::::::: (Uncertainty :::::::::::: Quantification) : study specification and consequently to the description of the uncertain input parameters X.Thus, special attention is paid to propose a relevant quantification of uncertainty sources.This particular point is addressed in section (4), in which frazil literature is explored to provide adequate bounds and PDFs for each parameter identified in Table (1).

Uncertainty propagation of random variables
In the following, we suppose that PDFs of inputs X are known, in contrast to those of Y.The objective is therefore to characterize the latter by estimating statistical indicators as the PDF's mean and standard deviation.To achieve this, a Monte Carlo sampling method can be used.
Statistical estimators of the response can then be computed.The Monte-Carlo estimation of the mean, denoted µ Y k and standard deviation, denoted σ Y k at time t k reads: These estimations converge to the true values following the law of large numbers (conditioned by the existence of corresponding PDF's first and second moment, i.e. expectation and variance (Soize, 2017)).The convergence order of the Monte-Carlo method is given by the central limit theorem leading to decrease of the confidence intervals' sizes proportional to √ N :::::: 1/ √ N .In the results of this paper, confidence intervals of the estimated mean and standard deviation are computed using a bootstrap method.
In that case, we only estimate first-order and total indices.Total indices quantify the part of variance of the output explained by an input and its interactions with all the other inputs parameters .Total Sobol indices are defined as follows, for k ∈ {1, . . ., n t } and i ∈ {1, . . ., n X }: where To compute Sobol indices, the method proposed by Saltelli (2002) is used, in which two independent samples of size N , denoted A and B are generated.Both can be written as matrices (cf.Equation B4) in which each line is a realization of the random vector X.A third matrix, denoted C i is then created by replacing only the column i of the matrix A by :::: other :::::::: uncertain ::::::::: parameters :: is ::::::::: determined ::: by ::::::::: subtracting ::: the :::: total ::::: Sobol ::::::::: sensitivity :::: index :::::: (ST k i ) ::: and : the column i of the matrix B (cf. Equation B4). .
First order and total Sobol indices are computed using estimations of V i (Y k ) and V −i (Y k ) computed using samples A, B and C i and noted Vi (Y k ) and V−i (Y k ) respectively.These estimations are defined as follows: where f is the centered model defined by f = f − f in which f is the mean of the combined output samples f (A) and f (B).To compute the second-order Sobol indices, an additional matrix is used, noted C j which is created by replacing only the column i of the matrix B by the column i of the matrix A. Then the estimation Vij (Y k ) is computed as: For a sample size N , estimation of the first-order and total Sobol indices requires (n X + 2) × N simulations.

Uncertainty source quantification
Many laboratory experiments have been carried out to better understand frazil ice dynamics as summed up by Barrette (2020Barrette ( , 2021)).These experiments, as well as field measurements, help us define parameter variability.In this section, we consider the uncertain parameters listed in Table (1), and deduce adequate variability and PDFs for uncertainty propagation and sensitivity analysis using the aforementioned data.The main geometrical properties of crystals and radial space discretization are first discussed, as well as initial concentration.We then review all the parameters involved in frazil source terms in the same chronology as presented in section (2).

Crystals' geometrical properties
Description of the frazil crystals' shape has been the subject of several field and laboratory measurements (Arakawa, 1954;Daly, 1984Daly, , 1994)).The crystals have been described as thin discs that grow mainly from their peripheral area.More recently, photos have brought valuable confirmation of the disc shape of frazil crystals, but also highlighted the complexity in larger flocs' geometry (Clark and Doering, 2006;Ghobrial et al., 2012;McFarlane et al., 2014McFarlane et al., , 2015;;Schneck et al., 2019).Let us recall that, in both SSC and MSC models, discretization is done in radial space.Therefore, either the ratio R or the thickness e : λ must be considered constant to fully describe disc-shaped particles.In the present study, frazil crystals are assumed to have the same aspect ratios, which is consistent with previous studies hypothesis (Svensson and Omstedt, 1994;Holland et al., 2007).Note that Rees Jones and Wells (2018) considered constant thickness instead, but highlighted a weak dependency of the thermal growth to aspect ratio.Let us discuss the uncertainties associated with the choice of the radius discretization as well as the diameter to thickness ratio.
Either a mean radius or a radius space discretization must be specified for SSC and MSC models, respectively.In both cases, observed radius distribution and particle size, reported in a many studies, are taken into consideration.For the sake of synthesis, we analyzed the corresponding publications from 1950 to 2019, and report the available observations in Figure 3.
This figure is complementary to Table 7 of McFarlane et al. (2017), who summarize particle sizes from both laboratory and field measurements.Daly and Colbeck (1986) and Clark and Doering (2006) reported log-normal distributions of particles in laboratory experiments, with a mean radius ranging from 0.12 to 0.25 mm and 0.49 to 1.4 mm, respectively.This was later confirmed by Ghobrial et al. (2012) and McFarlane et al. (2015) in the University of Alberta Cold Room Facility, with mean size ranging from 0.66 to 0.94 mm.Schneck et al. (2019) recently published similar results under different salinities, with mean radius ranging from 0.45 to 0.52 mm.Log-normal distributions were also observed in field studies, with the mean ranging from 0.59 to 0.94 mm in a report by McFarlane et al. (2017).Gathering together all reported mean radius ranges, an uncertainty interval of 1.2 × 10 −4 to 2.1 × 10 −3 m was chosen for the mean radius in the SSC model.When their full variation range is considered, frazil crystal sizes are known to follow a log-normal distribution as previously argued, but there is no evidence or reason for the mean radius being log-normal as well.The number of field and laboratory observations is still too small to fit an empirical PDF on the mean radius with reasonable confidence.Hence, a log-uniform PDF is chosen with the previously described bounds.Using Log-uniform approximate distributions allows us to explore by means of evenly distributed values, parameters that vary over several orders of magnitude.
For the MSC model, both a minimum and a maximum radius are to be determined for the bounds of the radial discretization.
The aspect ratio ranged from 6 to 13 in the study by Daly and Colbeck (1986), while there were greater values, from 12.2 to 16.33, in Clark andDoering (2004, 2006).Arakawa (1954) reported an even wider aspect ratio range, from 5 to 100.More recently, McFarlane et al. (2014) reported aspect ratios ranging from 11 to 71 with a mean of 37 and a standard deviation of 11.Considering all these observations, an uncertainty range of 5 to 100 was chosen for the diameter-to-thickness ratio.As not enough data to properly fit a PDF on aspect ratios was provided, a uniform PDF is chosen for the sensibility analysis (Maximum Entropy Principle Soize (2017)).
Weak dependency of the thickness to the radius was reported by McFarlane et al. (2014), who suggest assuming a increasing aspect ratio as discs grow instead of a constant aspect ratio.Not :::: Note also that, as frazil ice forms larger flocs, the disc shape hypothesis commonly accepted in models does not hold anymore.This could lead to erroneous estimations of the thermal growth rate of larger flocs.However, neither the variability of the aspect ratio in the crystals' distribution, nor that of the shape, are considered in the present study ::: and ::: we ::::: make ::: the ::::::::: assumption :::: that : r :::: and :: R ::: and :::::::::: independent.

Initial concentration
As previously mentioned, either a non-zero initial volume fraction or a non-zero seeding rate is required to trigger thermal growth in the model.In experimental facilities, frazil ice nuclei are sometimes artificially introduced in water to initiate frazil ice growth (Muller, 1978;Tsang and Hanley, 1985) but the initial concentration is rarely addressed and depends on the method of initial seeding.As such, a large uncertainty domain is considered in this study for the initial volume fraction of frazil, ranging from 10 −8 to 10 −4 , to account for the lack of data on this parameter.Additionally, a log-uniform PDF is retained since the parameter range variation vary over several orders of magnitude.
For the MSC model, it is necessary to choose how the initial volume fraction is distributed over size classes.One possibility is to initialize the system with log-uniform distributions, similar to the one observed in laboratory experiments.However, observations of nuclei close to r c in size are limited.Consequently, authors have preferred simpler initialization methods, whereby a constant number of crystals is distributed equally over all classes (Svensson and Omstedt, 1994;Wang and Doering, 2005; Rees Jones and Wells, 2018).Holland et al. (2007) argued that distributing the initial volume fraction over a range of radii, thus changing the initial number of particles per class, significantly impacts results.They found that filling only the first size class, as in Hammar and Shen (1995) has less impact on results and should therefore be the preferred initialization method.Given poor evidence of initial predominance of each class in nature, we decided to test both initialization methods i.e.
distributing the initial volume fraction over a range of small radii as presented in section (2.4) then feeding only the first class as suggested by Holland et al. (2007).For the first method, ice nuclei are supposed to be initially spread between the minimum radius r min and a threshold radius r 0 , which we suppose can only be smaller than the mean radius (cf. Figure 3).Therefore, the threshold r 0 is considered uncertain within the bounds [1.2 × 10 −4 , 2.1 × 10 −3 ].For the second method, we set r 0 = r min , so that only the first receives the initial concentration.

Thermal growth and turbulent parameters
As shown in Equation ( 2), thermal growth (1) is mainly affected by two uncertain parameters, which are the Nusselt number and the thermal boundary layer thickness.In this section, we discuss the uncertainty of both parameters.Since the Nusselt number is being modeled via turbulent parameters, namely the turbulent dissipation rate and the turbulent intensity, their uncertainty is also addressed.
As discussed by Rees Jones and Wells (2018), the thermal boundary layer thickness δ T is not constant around a crystal and in recent studies there has been an inconsistent scaling of thermal growth.Svensson and Omstedt (1994) chose the length scale as δ T = e :::::: δ T = λ : while others have chose δ T = r (Smedsrud and Jenkins, 2004;Holland et al., 2007), leading to a serious underestimation of the growth rate, ::::: since ::::: λ < r ::: for ::::: frazil :::: discs.Rees Jones and Wells (2015) have shown that there is a logarithmic dependency of the thermal growth on the aspect ratio that favors the δ T = e ::::: δ T = λ : scaling.To account for the variability of the thermal boundary layer thickness, δ T should be taken as an uncertain parameter.However, it should be mentioned that, with the variation range being δ T ∈ [e, r] ::::::::: δ T ∈ [λ, r], the ANOVA methodology (section 3.3) cannot be applied, since the set of uncertain inputs, which contains δ T and e/r ::: λ/r, cannot be considered independent anymore.Sophisticated methods such as the Rosenblatt transformation (Soize and Ghanem, 2004) can be applied to transform the dependent parameter spaces into independent Gaussian inputs.This is however out of the scope of this paper.To overcome this difficulty, the sensitivity analysis is conducted with the most appropriate scaling: i.e., δ T = e ::::: δ T = λ, for both the SSC and MSC models.
A similar method was described by McFarlane et al. (2015) to estimate the dissipation rate in rivers leading to values ranging from 4.2 × 10 −4 to 1.4968 m 2 .s−3 .Schneck et al. (2019) summarized dissipation rates observed in oceans, ranging from 10 −9 to 10 −2 m 2 .s−3 and from 10 −10 to 10 −3 m 2 .s−3 in polar regions.In the present study, we consider dissipation rates varying between 10 −9 and 1.5 m 2 .s−3 , which cover most flows encountered in rivers and oceans.Similarly, we consider a wide range of flows from low turbulence to high turbulence intensity which include most of the work presented in Figure (3), leading to a turbulent intensity ranging from 1% to 20%.Log-uniform PDFs are considered for both the dissipation rate and the turbulent intensity.Note that turbulence influences not only the rate of growth of frazil crystals but also the rate of secondary nucleation, thus impacting the frazil size distribution.

Other source terms
Following the same order as in Equation ( 1), let us discuss uncertain parameters involved in frazil source terms other than thermal growth: secondary nucleation (2), flocculation (3) and gravitational removal (4).It should be mentioned that almost no direct observations of these processes have been reported in the literature.Therefore, the definition of the uncertainty intervals is mainly based on expert knowledge and past numerical experiments in which parameters were determined by comparison to observed radius distributions.
-The uncertainty interval of the seeding rate is set after Daly (1994) to [3 × 10 −1 , 10 −4 ] m −2 .s−1 , and to simplify the uncertainty analysis, is considered constant over time.For the secondary nucleation, Svensson and Omstedt (1994) proposed setting a common value for n max that would allow focus on calibration of the initial seeding.They found that a value of n max = 4 × 10 6 , along with initial seeding of the order of magnitude of n 0 = 10 4 , gave satisfactory results compared to the experimental results of Michel (1963) and Carstens (1966).Wang and Doering (2005) found that fitting a single n max value for all experiments was unsatisfactory.Instead, they proposed fitting both the initial seeding and n max , leading to values ranging from 2 × 10 4 to 2 × 10 5 for Clark and Doering ( 2004) experiments, and from 2 × 10 4 to 2 × 10 5 for Carstens (1966) experiments.Smedsrud (2002) proposed a different value of n max = 10 3 , which was also used by Smedsrud and Jenkins (2004) and Holland and Feltham (2005).
-Similarly, the flocculation parameter a f was calibrated to a value of 10 −4 s −1 by Svensson and Omstedt (1994), who compared the results of their simulation to the expected size distribution spectra (other parameters in their model, :::: like :: the ::::::: number ::: of ::::: class, ::: the :::: heat :::: sink, ::: the :::::::: turbulent ::::::::: dissipation :::: rate ::: and :::::: initial :::::::: condition : were set to m = 20, φ = −1000 . From all previous calibration efforts :::::: studies, one could choose n max from 10 3 to 10 7 and a f = 10 −4 .But given the fact that secondary nucleation and flocculation are still very poorly understood and modeled, we choose large uncertainty intervals for these parameters.We added an order of magnitude for the bounds leading to n max :::: n max : ranging from 10 2 to 10 8 and we suppose that flocculation can, depending on the conditions, be negligible, so that a f was varied from 10 −8 to 10 −3 .n max ::: The :::::::::: parameters :::: n max : and a f are considered independent from hydrodynamic parameters, and Log-uniform PDFs were chosen for both. -The gravitational removal is affected by both buoyant rise velocity of frazil particles, and deposition process once particles reach the surface of the water column.Several attempts to measure the frazil rise velocity were done (Osterkamp and Gosink, 1983a;Wueben, 1984;McFarlane et al., 2014)  (2014), with velocities ranging from approximately 0.7 to 16 mm/s for a radius of 1 mm.Models exhibit significant differences as well (cf.figure 4).To take into account uncertainties inherent to the choice of the rise velocity model, a buoyancy parameter a d is introduced.A rise velocity envelope, combining the results of all models, can be defined by upper and lower bounds denoted w max (r, R) and w min (r, R) (cf.figure 4).The interval of a d is defined from the mean of the gap between the simplified formulation used the present study for w r (Svensson and Omstedt, 1994) and upper and lower bounds of the rise velocity envelope.To avoid the modeling of dependencies, a constant value is taken for a d , even if the envelope depends on the radius and on the diameter to thickness ratio.Finally, we propose in which a + = mean(w max (r, R)/w r (r)) and a − = mean(w min (r, R)/w r (r)).By considering r ∈ [10 −5 , 10 −2 ]m and R ∈ [5, 100], the following interval is obtained: Finally, it should be noted that this uncertainty quantification is rather imprecise since the rise velocity dispersion depends on the shape of the particles and flow conditions.
Therefore, our analysis could be refined by taking into account dependencies.Also note that low values of a d could also be justified from the uncertainty of the deposition process, which is not modelled in the present study.

Uncertainty propagation and sensitivity analysis
In this section, we present the different Monte Carlo simulation cases considered for the uncertainty analysis, as synthesized in Table (3).Then we discuss the results of the uncertainty propagation and the sensitivity analysis obtained for both SSC and MSC models.

Propagation cases
Monte Carlo simulations are first carried out without seeding and gravitational removal (τ s = 0 and a d = 0) for both (SSC) and (MSC) models, which correspond to cases ( 1) and ( 2), respectively.Seeding and gravitational removal processes are subsequently considered in cases ( 3) and ( 4).For each case, the set of uncertain parameters is described in Table (2).When not in X, parameters are considered deterministic with the following default values: Statistical estimators are evaluated every 10s of physical time, leading to n t = 400.To cope with the computational cost of multi-class experiments we used clusters to run all simulations in parallel.The computation time of 4.5 × 10 5 multi-class simulations with m = 100 is ∼ 24 hours with 960 processes.A total of 4 million simulations were carried out for uncertainty propagation.

Results without gravitational removal
By neglecting the seeding rate and gravitational removal source terms, the steady state corresponds to the constant frazil production rate that only depends on the heat sink φ, as shown in section (2.5).As expected, the statistical estimator time series for temperature and frazil volume fraction, presented in Figure (5) for cases (1) and ( 2), show a very narrow scatter of the output PDF at the start of simulation and at steady state.The results converge towards the two asymptotes (T 0 and C ∞ respectively) of the mono-class ODE system: i.e., the constant cooling rate when t → 0 (initial supercooling phase) and the constant frazil production rate when t → ∞ (recovery phase).However, for both SSC and MSC models, the transition between the two asymptotes of the ODE system is spread out, and there is a significant difference in maximum supercooling between the 5 th , 25 th , 75 th and 95 th percentiles.For the median, a maximum supercooling of T (t c ) −0.018 • C is reached at t c = 180s for the SSC model.The gap between the 5 th and 95 th percentiles maximum supercooling is ∆t c = 14.8 min and ∆θ = 0.1 • C :::::::: ∆t c = 890 :: s ::: and ::::::::::::: ∆θ = 0.097 • C. Note that the envelope obtained with the standard deviation gives a poor description of the output since its PDF is not normal, :: as :: it ::: can ::: be :::: seen :::: from ::: the ::::: μ ± σ::::: lines :: in ::: the ::::: figure : 5.
Similar results were obtained with the MSC model, which recovers the same asymptote at steady state, however, there is a slight residual scatter at recovery.For case (2), we observed slightly less scatter at the maximum supercooling than with the SSC model, and the gap between the 5 th and 95 th percentiles maximum supercooling was ∆t c = 13.7 min ::::::::: ∆t c = 820 : s : and ∆θ = 0.08 • C. :: (a)  Time series of the first-order Sobol indices are presented in Figure ( 6) (see Appendix E and F for details).In Appendix C, we also presente the first-and total-order Sobol indices at times t c , 2 × t c and t f with 95% confidence intervals, as well as the aggregated Sobol indices, computed via Equation (B7).The time evolution of Sobol indices of temperature and frazil were similar for both the SSC and MSC with the exception of the initial concentration, which, as expected had more impact on frazil volume fraction at the start of simulation.
At the recovery, the parameters of secondary nucleation and flocculation processes ::::: (n max :::: and ::: a f ), both impacting steadystate crystal distribution, become more influential.However, the hierarchy of parameter is less obvious than for the (SSC) model.In addition, we observed strong interactions between parameters when the model reaches steady state ::: (cf.::::: blank ::::: space :: on :::::: While the approach adopted by Svensson and Omstedt (1994): i.e., tweaking the values of n max and a f , was adequate to calibrate the frazil distribution in (MSC) models, the present sensitivity analysis shows that it might not be the best option to calibrate water temperature and frazil total volume fraction.Results suggests that more focus should be on the initial condition to calibrate supercooling, by modifying initial seeding like it was done by Wang and Doering (2005) or by modifying the initial distribution itself.We therefore suggest calibrating the supercooling curve with the help of the initial distribution, along with secondary nucleation and flocculation parameters to calibrate the evolution of size distribution over time.Hopefully, recent observations of the transient evolution of frazil size distribution (McFarlane et al., 2015;Schneck et al., 2019), will provide the necessary data to carry out an optimal calibration of the identified parameters. :: (a)  2), (3) and (4).

Influence of gravitational removal and seeding
Long-term evolution that does not take account of gravitational removal leads to infinite increase in frazil concentration as long as the cooling rate remains constant (cf.Equation 15).This asymptotic behavior of the models has never been observed in experiment nor in nature.Clark and Doering (2006) observed a peak in the number of particles per image they recorded, located shortly after maximum supercooling, after which there was a small decrease in the number of particles and a stagnation at residual supercooling.Similar observations were also reported by McFarlane et al. (2015) and Schneck et al. (2019).The models in which only thermal growth is considered do not incorporate the required physics to properly reproduce what is observed.However by introducing gravitational removal, as shown in section (2.5), the models converge towards a constant frazil volume fraction (cf.Equation 16).In this section we analyze the results of the uncertainty propagation and sensitivity analysis for cases (3) and (4) which include the seeding rate and gravitational removal source terms.
At steady state, the introduction of gravitational removal leads to wide scatter of both temperature and frazil as shown in :: by :: (a)  15).For the MSC model, the most influential parameters at recovery are n max and a d (S i -C = 0.37 and 0.26 ::::::::::: S nmax = 0.37 ::: and :::::::::: S a d = 0.26, respectively at t = t f ), which is coherent :::::::: consistent with Equation ( 27) and ( 29) of Rees Jones and Wells (2018).The hierarchy of the most influential parameters is similar to cases (1) and (2) prior to maximum suprcooling.However accurate modeling of the buoyancy velocity of frazil crystals is essential, as it has a very important influence on the long-term evolution of the system, and therefore merits particular attention.

Maximum supercooling scatter
The results discussed above were obtained with a cooling rate of −500 W.m −3 .Several cooling rates, ranging from −50 to −1000 W.m −3 , were tested with case (1) to assess variations in maximum supercooling predictions.We found that the higher the cooling rate, the greater the scatter of the predected :::::::: predicted maximum supercooling temperature, as presented in Figures (10) and (11).This is the opposite for the time until maximum supercooling peak, where the higher the cooling rate, the lower the scatter in supercooling time.δ T = r.The results highlight the significant influence of the choice of scaling for the boundary layer.Choice of scaling also 620 explains inconsistencies in calibrated parameters in literature.
The median maximum supercooling is only −0.001 • C and is reached in only 20 seconds.This confirms the sensitivity test carried out by Holland and Feltham (2005), who suggested distributing initial concentration on one class.However, the results show that this type of initialization might not be the best option, as it almost totally does away with transient evolution of the model.
In this paper, we considered a well-mixed water body and a simplified gravitational removal sink term.Spatial variations for temperature and frazil may result in different conclusions.A poorly mixed water body, where the cooling rate at higher layers is more severe than at the bottom, would cause an heterogeneity in the initial formation of frazil.Additionally, the meteorological seeding of frazil nuclei occurs mainly at the free surface, increasing the heterogeneity.Furthermore, as the rise velocity is may be removed from choosing a mean radius.We have shown, however, that scatter is similar somehow in both models, and derives from new uncertain parameters inherent in radial space discretization.Note that, in many numerical tools, modeling frazil distribution requires the resolution of multiple advection-diffusion equations.Given the number of classes required for a model convergence, one can easily grasp the high numerical cost of using the MSC model for large scale, multi-dimensional applications.This makes the SSC model a relevant ::::::: suitable candidate for multi-dimensional frazil ice modeling, and the present study shows that it is still a very good compromise between uncertainty and model complexity.However, it should be noted that such uncertainty in the MSC model could be overcome in future :::::::: laboratory :::::::::: experiments : by a better estimation of the initial crystal size distribution.
The sensitivity analysis, allowed us to address with confidence the choice of calibration parameters.Relying on first-and total-order Sobol indices, we quantified the relative influence of each uncertain parameter on the output distribution, for both SSC and MSC models, and proposed a selection of parameters to be used for calibration.For the SSC model, the most influential factor for both temperature and frazil is the mean radius.Initial concentration played a secondary role although it was initially identified as a predominant factor.We therefore suggest using the average radius as the main calibration parameter.
The turbulent dissipation rate also plays a major role and as such should be specified with care.As it is often appropriately quantified ::: can :: be :::::::: estimated ::: via ::::::::: turbulence :::::: models, we suggest using initial concentration and diameter-to-thickness ratio as secondary calibration parameters.With the MSC model, we showed that the dispersion is somehow similar to what we observed for the SSC model but originates from new uncertain parameters.Thus, the most influential parameters on the transient phase are the parameters specific to the initial condition.However, once the steady state was reached, we observed an increasing influence of the secondary nucleation and flocculation parameters.The long-term evolution of the system also showed increasing interactions between parameters, which can : .:::: This ::::: could : be explained by the balance in the physical processes involved in class interactions : , ::: and :::::: further ::::::::::: investigation ::::: using :::: class ::::::: volume ::::::: fraction :: as :::::: model :::::: outputs :::::: would :::: help.When gravitational removal was introduced in the models, the stationary state was modified and the concentration converged towards a finite limit instead of diverging.In the case of the SSC model, the asymptotic limit is a function of the ratio between the gravitational removal term and the heat flux while in the case of the MSC model, the stationary state is also a function of the steady state radius distribution (which depends on the balance between secondary nucleation, flocculation and gravitational removal).Our study, confirming previous asymptotic analyses, showed that both secondary nucleation and gravitational removal parameters are the most influential on total frazil volume fraction.The buoyant rise velocity, the uncertainty of which was rarely taken into account in previous modeling studies, should therefore be one of the main foci of future efforts to calibrate frazil ice models.
Contrary to frazil, water temperature was mostly influenced by the initial condition, even at steady state.This should impel us to use both water temperature and frazil volume fraction measurements to calibrate the models.Fortunately, recent laboratory and field studies (Schneck et al., 2019;McFarlane et al., 2015McFarlane et al., , 2017)), particularly on the evolution of frazil distributions over time, offer precious data that can assist in developing appropriate calibration of the models.In this regard, using optimal calibration techniques, that allow consideration of both modeling and data uncertainties, would be a natural and complementary extension of the present study.

Figure 3 .
Figure 3. Particle size ranges (thick grey line), mean range (thin dark line) and log-normal distributions' mean radius (red vertical ticks) and many formulations were proposed as summarized by McFarlane et al. (2014).Significant scatter can be observed in the data as shown on Figure (10) of McFarlane et al.
Figure :: 6 ::: and ::::: Total ::::: Sobol :::::: indices ::: in Appendix F : ), which might lead to difficulties in the calibration process.Aggregated Sobol indices summarized in Figure (7), confirm the relative influence of each parameter over the whole duration of the simulation, taking into account both the transient phase and steady state.

Table 1 .
Description of uncertain parameters of the frazil ice models.
represents a possible configuration of the frazil model.The output realizations i.e. y i = y 1 i , . . ., y nt i , are generated by N successive deterministic simulations with corresponding inputs.For a sampling :::::: sample : of size N , the output of the model is a N × n t matrix: