2.1 Flow field solution: Hess-Smith panel method

Two-dimensional Hess-Smith panel method is used to determine the flow velocity components required for droplet trajectory calculations and the velocity distribution around the aerofoil for boundary-layer calculations. The aerofoil is divided into quadrilateral panels, each having a constant strength source singularity element plus a vortex singularity that is constant for all panels. Using the flow tangency boundary condition at the collocation points of the surface panels, the singularity strengths are obtained. This provides the velocity potential at any point in the flow field. The derivatives of the velocity potential are the velocity components at a given point in the flow field. Moreover, integral boundary-layer method is used in TAICE for calculating heat transfer coefficients when they are not supplied with another tool. Since equations which are used for calculating heat transfer coefficients depend on the flow being laminar or turbulent, transition prediction should be considered. The method which is used for the calculation of these coefficients are given in well-known references^{(Reference Gent, Dart and Cansdale13,Reference Schlicting27)} .

2.2 Particle trajectory calculations

The equations governing droplet and crystal trajectories are identical. However, the models used for drag force and heat transfer predictions differ depending on the phase of cloud particles. The differences are pointed out below, where appropriate. The following assumptions are made for the particles:

• • The presence of cloud particles has no effect on the flowfield,

• • Gravity and aerodynamic drag are the only forces acting on the cloud particles.

Thus, the governing equations for particle motion in Cartesian coordinates are^{(Reference Myers22)}:

(1) $$\begin{equation} m_p\ddot{x} = -D \cos {\gamma }, \\ \end{equation}$$

(2) $$\begin{equation} m_p\ddot{y} = -D \cos {\gamma }+mg, \\ \end{equation}$$

(3) $$\begin{equation} \gamma = \arctan \frac{\dot{y}_p-V_y}{\dot{x}_p-V_x}, \\ \end{equation}$$

(4) $$\begin{equation} D = \frac{1}{2} \rho V^2 C_D A_p, \\ \end{equation}$$

(5) $$\begin{equation} V = \sqrt{\left(\dot{x}_p-V_x\right)^2+\left(\dot{y}_p-V_y\right)^2}. \end{equation}$$

In the above expressions, C _D is drag force coefficient, and V _x and V _y are velocity components of the particle. Calculation of the drag coefficient depends on the particle type. Ice crystals are formed generally within the tropical clouds and may have various shapes. They differ from liquid droplets in terms of their shape and usually are not exactly spherical. Thus, for accurate prediction of the trajectories, drag force coefficient, heat transfer and phase change models should be modified for non-spherical particles in order to be used for ice crystals. Required models are depicted along the trajectory shown in Fig. 2. Brief explanations about the trajectory models are given in section A of the Appendix. More detailed information can be found in the studies of Trontin et al, Villedieu et al^{(Reference Trontin, Blanchard and Villedieu28,Reference Villedieu, Trontin and Chauvin30)} and Ayan et al^{(Reference Ayan, Özgen and Canıbek2)}.

Figure 2.

Models related to mixed phase icing starting from the release point till accretion on the surface.

2.3 Impingement model for ice crystals

Impinging ice crystals on a surface may lead to three different scenarios which are sticking, bouncing and shattering. Impact model which is discussed in this section is mainly based on HAIC experiments and literature data^{(Reference Trontin, Blanchard and Villedieu28,Reference Villedieu, Trontin and Chauvin30)} . A brief information will be given here about these regimes.

In bouncing and sticking regimes, particle size and shape do not change after impact. However, it may split into smaller particles in the shattering regime. In the bouncing regime, particles just bounce off from the wall and are re-emitted to the flow. On the other hand, in sticking regime, particles adhere to the wall completely. The threshold between these regimes depends on the kinetic energy of the particle before impact. Below, a dimensionless number L, which is the ratio of the kinetic energy to the surface energy of the particle, is introduced^{(Reference Trontin, Blanchard and Villedieu28,Reference Villedieu, Trontin and Chauvin30)} .

(6) $$\begin{eqnarray} L=\frac{1}{12}\frac{\pi \rho _p d_p^{3}v_{pn}^{2}}{\pi e_\sigma T_p d_p^{2}}=\frac{1}{12}\frac{\rho _p d_p v_{pn}^{2}}{e_\sigma T_p} , \end{eqnarray}$$

In the above equation, v _pn is the normal velocity of the particle and e _σ is the surface energy per unit area and can be calculated with the following relation:

(7) $$\begin{eqnarray} e_\sigma =e_{\sigma 0}\exp \left(\frac{Q_s}{RT}-\frac{Q_s}{RT_0}\right),\vspace*{-5pt} \end{eqnarray}$$

where Q _s is the activation energy. As it is given in Higa et al^{(Reference Higa, Arakawa and Maeno18)}, Q _s = 4.82 × 10 ⁴ J K⁻¹.mol⁻¹ and according to the Vidaurre and Hallet^{(Reference Vidaurre and Hallett29)} e _σ0 = 0.12 J.m⁻² at T ₀ = 253 K.

Thus, when L ⩾ 90, high loss of kinetic energy and high fragmentation occurs upon impact; when L ⩽ 0.5, no fracturing is expected. On the other hand, when L is in between these values, inelastic bouncing with some particle fragmentation is expected. However, less kinetic energy loss will occur in this case. These values are determined based on both literature data^{(Reference Hauk, Roisman and Tropea17)} and Hauk’s experimental results^{(Reference Hauk, Roisman and Tropea17)}.

Predicting sticking ratio is the challenging part of this model and a simple relation is suggested by Trontin et al^{(Reference Trontin, Blanchard and Villedieu28)}. In their study, sticking ratio depends only on the mean melt ratio of the impinging ice crystals and total water content of the flow. The model is calibrated using the National Aeronautics and Space Administration and National Research Council of Canada (NASA - NRC) experiments^{(Reference Currie, Fuleki, Knezevici and MacLeod6-Reference Currie, Struk, Tsao, Fuleki and Knezevici8)}. However, in this study, a more simplified form of their relation is used:

(8) $$\begin{eqnarray} \epsilon _s&=&K_l\psi _l\vspace*{-5pt} \end{eqnarray}$$

In the above equation, sticking ratio is related to the dimensionless parameter, K _l , and the liquid ratio of the impinging particles, ψ _l . According to Trontin et al^{(Reference Trontin, Blanchard and Villedieu28)}, in rime ice conditio,n ice crystals do not stick on the surface and K _l should be equal to zero for consistency. Moreover, after calibration study, K _l is fixed to 2.25 for glaze ice conditions.

Impinging liquid and solid ratios in total amount of impinging particles are given below.

(9) $$\begin{equation} \psi _l = \frac{\dot{m}_{imp,l}}{\dot{m}_{imp,c}+\dot{m}_{imp,l}} \\ \end{equation}$$

(10) $$\begin{equation} \psi _c = \frac{\dot{m}_{imp,c}}{\dot{m}_{imp,c}+\dot{m}_{imp,l}} \end{equation}$$

Moreover, it is observed from the experiments that the sticking ratio of the ice crystals increases with increasing ψ _l up to a threshold value. However, further increase does not cause additional increase in the sticking ratio of the ice crystals. Thus, a limitation is introduced here for ψ _l for this purpose. Its value could not exceed 0.2, which is determined during the calibration study.

In addition, the model predicts behaviour of secondary particles which are produced in case of bouncing or shattering regimes. The model is mainly empirical and the adjustable constant has been fitted thanks to literature and HAIC experimental data.

Since bouncing ice crystal particles may re-impinge and accumulate in downstream, their trajectory calculations should be continued after impact. Relations which are used for modelling these secondary particle behaviours can be found in Refs 2 and 28.

2.4 Erosion Model

It is observed in experiments that in mixed phase icing conditions, upcoming solid particles cause reduction in the amount of accreted ice on the surface^{(Reference Al-Khalil, Irabi and Miller1)}. Generally, this reduction is much more noticeable in glaze ice conditions. Wright et al^{(Reference Wright, Jorgenson and Veres32)} proposed a model for prediction of erosion effect. In this study, this model is modified by using additional available experimental data. A modified version of the model is given below which is used to predict erosion effect on accreted ice.

(11) $$\begin{eqnarray} \epsilon _{er}=\epsilon _{ct}\left(\frac{d_p}{d_{p,ref}}\right)^3\left(\frac{v_{tang}}{v_{ref}}\right)^{2.25}\left(\frac{IWC}{IWC_{ref}}\right)^{0.2} \left(\frac{exp\left(14-\frac{Q}{RT_{surf}}\right)}{exp\left(14-\frac{Q}{RT_{ref}}\right)}\right)^{1.25} \end{eqnarray}$$

In the above relation, ε _er by definition refers to erosion rate and its value should be between 0 and 1, ε _ct is erosion constant, d _{p, ref} is particle reference diameter, V _ref is particle reference velocity, IWC _ref is reference ice water content and T _ref is reference equilibrium temperature of the panel. Q is activation energy (Q = 69000 J.mol⁻¹) and R is perfect gas constant (R = 8.314 J.mol⁻¹.K⁻¹). According to Wright et al and after calibration study, ε _ct = 0.17, d _{p, ref} = 150 µm, V _ref = 51.44 m/s, IWC _ref = 0.7 g/m³ and T _ref = 265.14 K are taken. Moreover, erosion efficiency in ice accretion calculation relations are modified according to the suggestion of Trontin et al^{(Reference Trontin, Blanchard and Villedieu28)} for easy implementation and is further calibrated in this study.

(12) $$\begin{eqnarray} \dot{m}_{er}&=&\epsilon _{er}\dot{m}_{imp,c} \end{eqnarray}$$

In the above relation, $\dot{m}_{er}$ is the amount of reduced mass rate due to erosion and can be calculated by multiplying erosion efficiency with impinging mass rate. This reduced mass includes both water and ice parts as it is shown below:

(13) $$\begin{eqnarray} \dot{m}_{er}&=&\dot{m}_{er,i}+\dot{m}_{er,w}, \end{eqnarray}$$

where $\dot{m}_{er,i}$ and $\dot{m}_{er,w}$ are ice and water erosion mass rates, respectively, and they are taken into consideration while total accretion rates are calculated^{(Reference Trontin, Blanchard and Villedieu28)}.

2.5 Ice Accretion: Further Extended Messinger Model

The ice accretion problem is based on the phase change or Stefan problem. According to the Extended Messinger Model, the problem is governed by four equations: energy equations in the ice and water layers, the mass conservation equation and the phase change condition at the ice/water interface^{(Reference Myers22)}:

(14) $$\begin{equation} \frac{\partial T}{\partial t} = \frac{k_i}{\rho _i C_{pi}}\frac{\partial ^2 T}{\partial y^2},\\ \end{equation}$$

(15) $$\begin{equation} \frac{\partial \theta }{\partial t} = \frac{k_w}{\rho _w C_{pw}}\frac{\partial ^2 \theta }{\partial y^2}, \\ \end{equation}$$

(16) $$\begin{equation} \rho _i \frac{\partial B}{\partial t}+\rho _w \frac{\partial h}{\partial t} = r_l \rho _l \beta _l V_\infty + r_c \rho _c \beta _c V_\infty + \dot{m}_{in} - \dot{m}_{e,s}, \\ \end{equation}$$

(17) $$\begin{equation} \rho _i L_f \frac{\partial B}{\partial t} = k_i \frac{\partial T}{\partial y}-k_w\frac{\partial \theta }{\partial y} \end{equation}$$

1. (i) It is assumed that ice is in perfect contact with the substrate (y = 0), which has high conductivity and a thermal mass much greater than ice. Thus, in this model, temperature variation of the substrate is assumed as negligible and it has a fixed and known value. Moreover, for ice accretion on the surface, the substrate temperature should not be greater than T _f . Otherwise, only water can accumulate on the surface.

   (18) $$\begin{eqnarray} T(0,t) = T_{subs} = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}T_{rec} ,& \text{if }\ T_{rec} \le T_f \\ T_f ,& \text{if }\ T_{rec} > T_f \end{array}\right. \end{eqnarray}$$

2. (ii) The temperature is continuous at the phase change boundary (y = B) and equal to the freezing temperature T _f . This condition is not imposed in rime ice conditions:

   (19) $$\begin{equation} T(B,t)=\theta (B,t)=T_f. \end{equation}$$

3. (iii) At the air/ice interface (y = B) in rime ice conditions, heat flux is determined by convection Q _con , sublimation Q _sub , cooling by incoming droplets Q _cl , cooling by incoming crystals Q _cc (notice that ice crystals may be partially melted, and may contain both liquid and solid phase, taken into account by terms Q _ccl and Q _ccs , respectively, and Q _cc = Q _ccl + Q _ccs in Equation (20)), kinetic energy of incoming droplets Q _kl , kinetic energy of incoming crystals Q _kc , aerodynamic heating Q _aero , latent heat release Q _lat and energy brought in by runback water Q _in :

   (20) $$\begin{eqnarray} -k_i \frac{\partial T}{\partial y}&=&\left(Q_{con}+Q_{sub}+Q_{cl}+Q_{cc}\right)-\left(Q_{aero}+Q_{kl}+Q_{kc}+Q_{lat}+Q_{in}\right)\nonumber\\ && \hspace{11.38092pt} {\rm at} \hspace{11.38092pt} y=B \end{eqnarray}$$
   In rime ice conditions, the heat terms are expressed as:
   (21) $$\begin{equation} Q_{con} = h_c\left[T\left(B,t\right)-T_a\right],\\ \end{equation}$$

   (22) $$\begin{equation} Q_{sub} = \chi _{sub} \left[p_vT(B,t)\frac{T_o}{T(B,t)}\left(\frac{p_i}{p_o}\right)^{-1/\gamma }-(RH)p_v(T_a)\frac{p_o}{p_a}\right],\\ \end{equation}$$

   (23) $$\begin{equation} Q_{cl} = r_l \rho _l \beta _l V_\infty C_{pw} \left[T(B,t)-T_l \right], \\ \end{equation}$$

   (24) $$\begin{equation} Q_{ccs} = (1-\eta _{melt}) r_c \rho _c \beta _c V_\infty C_{pi} \left[T(B,t)-T_c \right],\\ \end{equation}$$

   (25) $$\begin{equation} Q_{ccl} = \eta _{melt} r_c \rho _c \beta _c V_\infty C_{pw} \left[T(B,t)-T_c \right],\\ \end{equation}$$

   (26) $$\begin{equation} Q_{kl} = r_l \beta _l V_\infty \frac{V_{imp,l}^2}{2}, \\ \end{equation}$$

   (27) $$\begin{equation} Q_{kc} = r_c \beta _c V_\infty \frac{V_{imp,c}^2}{2},\\ \end{equation}$$

   (28) $$\begin{equation} Q_{aero} = r h_c \frac{V_\infty ^2}{2C_p}, \\ \end{equation}$$

   (29) $$\begin{equation} Q_{lat} = \rho _r L_f \frac{\partial B_l}{\partial t}, \\ \end{equation}$$

   (30) $$\begin{equation} Q_{in} = \dot{m}_{in}C_{pw}\left[T_f-T(B,t)\right] \end{equation}$$

   At the air/water interface (y = B + h) in glaze ice conditions, heat flux is determined by convection Q _con , evaporation Q _ev , cooling by incoming droplets Q _cl , cooling by incoming crystals Q _cc (notice that ice crystals may be partially melted, and may contain both liquid and solid phase, taken into account by terms Q _ccl and Q _ccs , respectively and Q _cc = Q _ccl + Q _ccs in Equation (31)), kinetic energy of incoming droplets Q _kl , kinetic energy of incoming crystals Q _kc , aerodynamic heating Q _aero and energy brought in by runback water Q _in :

   (31) $$\begin{equation} -k_w \frac{\partial \theta }{\partial y}=(Q_{con}+Q_{ev}+Q_{cl}+Q_{cc})-(Q_{aero}+Q_{kl}+Q_{kc}+Q_{in}) \hspace{11.38092pt} {\rm at} \hspace{11.38092pt} y=B+h. \end{equation}$$
   In glaze ice conditions, the heat terms are expressed as:
   (32) $$\begin{equation} Q_{con} = h_c[\theta (B+h,t)-T_a], \\ \end{equation}$$

   (33) $$\begin{equation} Q_{ev} = \chi _{ev} \left[p_v\left(\theta (B+h,t)\right)\frac{T_o}{\theta (B+h,t)}\left(\frac{p_i}{p_o}\right)^{-1/\gamma }-RHp_v(T_a)\frac{p_o}{p_a}\right],\\ \end{equation}$$

   (34) $$\begin{equation} Q_{cl} = r_l \rho _l \beta _l V_\infty C_{pw} \left[\theta (B+h,t)-T_l \right], \\ \end{equation}$$

   (35) $$\begin{equation} Q_{ccs} = (1-\eta _{melt}) r_c \rho _c \beta _c V_\infty C_{pi} \left[\theta (B+h,t)-T_c \right], \\ \end{equation}$$

   (36) $$\begin{equation} Q_{ccl} = \eta _{melt} r_c \rho _c \beta _c V_\infty C_{pw} \left[\theta (B+h,t)-T_c \right], \\ \end{equation}$$

   (37) $$\begin{equation} Q_{kl} = r_l \beta _l V_\infty \frac{V_{imp,l}^2}{2}, \\ \end{equation}$$

   (38) $$\begin{equation} Q_{kc} = r_c \beta _c V_\infty \frac{V_{imp,c}^2}{2}, \\ \end{equation}$$

   (39) $$\begin{equation} Q_{aero} = r h_c \frac{V_\infty ^2}{2C_p}, \\ \end{equation}$$

   (40) $$\begin{equation} Q_{in} = \dot{m}_{in}C_{pw}\left[T_f-\theta (B+h,t)\right] \end{equation}$$

4. (iv) Surface (y = 0) is initially (t = 0) devoid of ice and water:

   (41) $$\begin{equation} B=h=0 \hspace{11.38092pt} {\rm at} \hspace{11.38092pt} t=0 \end{equation}$$

2.5.1 Rime ice growth and temperature profile

Setting h = 0 in Equation (16) yields:

(42) $$\begin{equation} B(t)=\frac{r_l \rho _l \beta _l V_\infty + r_c \rho _c \beta _c V_\infty +\dot{m}_{in}-m_{sub}}{\rho _r}t \end{equation}$$

The temperature profile in the ice layer is assumed to be linear:

(43) $$\begin{eqnarray} \frac{\partial ^2 T}{\partial y^2}=0 \end{eqnarray}$$

Integrating twice:

(44) $$\begin{eqnarray} T(y,t)=a_1y+a_2 \end{eqnarray}$$

Applying Equation (18) yields:

(45) $$\begin{eqnarray} a_2=T_{subs} \end{eqnarray}$$

In order to find a ₁, Equation (20) is employed:

(46) $$\begin{eqnarray} -k_{i}a_1=(Q_{con}+Q_{sub}+Q_{cl}+Q_{cc})-(Q_{aero}+Q_{kl}+Q_{kc}+Q_{lat}+Q_{in}) \hspace{11.38092pt} {\rm at} \hspace{11.38092pt} y=B \nonumber\\ \end{eqnarray}$$

Then, using Equations (21) to (30):

(47) {\fontsize{9}{11}\selectfont{$$\begin{eqnarray*} a_1&=&\frac{-h_c\left[T(B,t)-T_a\right]-r_l \rho _l \beta _l V_\infty C_{pw}\left[T\left(B,t\right)-T_l\right]+\dot{m}_{in}C_{pw}\left[T_f-T(B,t)\right]}{k_i}\nonumber\\ &&+\,\frac{-r_c \rho _c \beta _c V_\infty \left(\left(1-\eta _{melt}\right)C_{pi}+\eta _{melt}C_{pw}\right)\left[T(B,t)-T_c\right]+\left(Q_{aero}+Q_{kl}+Q_{kc}+Q_{lat}-Q_{sub}\right)}{k_i}\nonumber\\ \end{eqnarray*}$$}}\vskip-45pt$$\begin{equation} \end{equation}

Evaluate Equation (44) at y = B and by substituting a ₁ and a ₂ from Equations (45) and (47)

(48) $$\begin{eqnarray} T(B,t)=\frac{k_iT_{subs}+\left(aT_a+bT_l+cT_c+eT_f+d\right)B}{k_i+\left(a+b+c+e\right)B}, \end{eqnarray}$$

where

(49) $$\begin{equation} a = h_c, \\ \end{equation}$$

(50) $$\begin{equation} b = r_l \rho _l \beta _l V_\infty C_{pw}, \\ \end{equation}$$

(51) $$\begin{equation} c = r_c \rho _c \beta _c V_\infty \left(\left(1-\eta _{melt}\right)C_{pi}+\eta _{melt}C_{pw}\right),\\ \end{equation}$$

(52) $$\begin{equation} d = Q_{aero}+Q_{kl}+Q_{kc}+Q_{lat}, \\ \end{equation}$$

(53) $$\begin{equation} e = \dot{m}_{in}C_{pw} \end{equation}$$

For the calculation of Q _sub in Equation (22), vapour pressure needs to be calculated. This is a function of T(B, t). Therefore, T(B, t) needs to be calculated through an iterative process until convergence in Q _sub is reached. Substituting T(B,t) from Equation (48) back into Equation (47) yields:

(54) $$\begin{eqnarray} a_1=\frac{-a(T_{subs}-T_a)+[d+b(T_l-T_{subs})+c(T_c-T_{subs})+e(T_f-T_{subs})]}{k_i+(a+b+c+e)B} \end{eqnarray}$$

Substituting a ₁ from Equation (54) and a ₂ from Equation (45) yields:

(55) $$\begin{eqnarray} T(y,t)=a_1 y+T_{subs} \end{eqnarray}$$

2.5.2 Glaze ice growth and temperature profile

While calculating ice growth in the rime ice condition, since all upcoming water droplets and liquid parts of the partially melted ice crystals freeze, a water layer does not appear above the ice surface. However, when the surface temperature of the ice layer reaches the melting temperature, water starts accumulating on the surface. In this case, the physical phase of ice crystals needs more attention. At this point, three different scenarios may occur:

• • Depending on the energy transfer between ice crystals and water layers on the surface, upcoming ice crystals may completely melt.

• • Depending on the energy transfer between ice crystals and water layers on the surface, upcoming ice crystals may partially melt.

• • If surface temperature of the water layer is below melting temperature, they may remain solid.

Thus, the melting ratio is an additional unknown and can be predicted by solving the energy balance equation. Energy and mass terms which are included in the control volume are depicted in Fig. 3.

Figure 3.

Mass and energy terms in the balance equations.

The energy terms that contribute to the energy balance are as follows:

(56) $$\begin{equation} \lambda _{h,w} = r_l \rho _l \beta _l V_\infty C_{pw}\left[T_f-T_l\right]+\dot{m}_{in}C_{pw}\left[T_f-T_{in}\right],\\ \end{equation}$$

(57) $$\begin{equation} \lambda _{h,i} = r_c \rho _c \beta _c V_\infty (1-\eta _{melt})C_{pi}\left[T_f-T_c\right]+r_c \rho _c \beta _c V_\infty \eta _{melt}C_{pw}\left[T_f-T_c\right], \\ \end{equation}$$

(58) $$\begin{equation} \lambda _{in} = r_l \beta _l V_\infty \frac{V_{imp,l}^2}{2}+r_c \beta _c V_\infty \frac{V_{imp,c}^2}{2}, \\ \end{equation}$$

(59) $$\begin{equation} \lambda _{out} = h_c\left[T_{rec}-T_f\right]+\dot{m}_{ev}L_{ev} \end{equation}$$

The first term in previous equations, λ _{h, w} , refers to the heat transfer due to impinging droplets and runback water in. Heat transfer due to impinging solid and liquid parts of ice crystals are calculated by the second relation, λ _{h, i} . On the other hand, λ _in refers to the total kinetic energy of the upcoming particles and λ _out shows heat transfer due to convection and evaporation. Inclusion of the heat loss by melting part of the incoming ice crystals causes an imbalance in the energy equation and this must be compensated by additional freezing. Thus, at melting temperature, the net energy before computing the melting ratio can be calculated as follows:

(60) $$\begin{equation} \lambda _{net} = \lambda _{in}-\lambda _{out}-\lambda _{h,w}-\lambda _{h,i} \\ \end{equation}$$

(61) $$\begin{equation} \alpha _{melt} = \frac{\lambda _{net}}{r_c \rho _c \beta _c V_\infty L_f} \end{equation}$$

Three different scenarios are possible after calculating the net energy terms:

• • When λ _net > 0 and α _melt < 1, some energy is available in the control volume and it will melt some part of the ice crystals. Surface temperature of the water layer will be equal to the melting temperature.

• • When λ _net > 0 and α _melt > 1, some energy is available in the control volume and all ice crystals will melt. Remaining energy will increase the temperature of the water layer and the surface temperature of the water layer will be greater than the melting temperature.

• • When λ _net < 0, all upcoming ice crystals will stay as solid, even after impact. The surface temperature of the water layer will be less than the melting temperature.

Next, it is assumed that the temperature profile in the ice layer is linear:

(62) $$\begin{eqnarray} \frac{\partial ^2 T}{\partial y^2}=0 \end{eqnarray}$$

Integrating twice T(y, t) can be obtained:

(63) $$\begin{eqnarray} T(y,t)=a_3 y+a_4 \end{eqnarray}$$

Employing the condition defined in Equation (18) yields:

(64) $$\begin{eqnarray} a_4=T_{subs} \end{eqnarray}$$

Thus, at y = B:

(65) $$\begin{eqnarray} T_f=a_3B+T_{subs} \end{eqnarray}$$

Assume that the temperature profile in the water layer is linear as well:

(66) $$\begin{eqnarray} \frac{\partial ^2 \theta }{\partial y^2}=0 \end{eqnarray}$$

Integrating twice:

(67) $$\begin{eqnarray} \theta (y,t)=a_5 y+a_6 \end{eqnarray}$$

Employing Condition 14 yields:

(68) $$\begin{eqnarray} T_f=a_5B+a_6 \end{eqnarray}$$

In Equation (68), the calculation of a ₅ is not necessary when the temperature of the water layer surface is equal to the melting temperature, which corresponds to the 0 < α _melt < 1 case. However, for other cases, the condition defined in Equation (18) needs to be imposed:

(69) $$\begin{eqnarray} -k_wa_5=(Q_{con}+Q_{ev}+Q_{cl}+Q_{ccs}+Q_{ccl})-(Q_{aero}+Q_{kl}+Q_{kc}+Q_{in}) \end{eqnarray}$$

Substituting Q _con , Q _ev , Q _cl , Q _ccs , Q _ccl and Q _in from Equations (33) to (40):

(70) $$\begin{eqnarray} a_5=\frac{-(a+b+c+e)\theta (B+h,t)+aT_a+bT_l+cT_c+eT_f+d}{k_w}, \end{eqnarray}$$

where

(71) $$\begin{equation} a = h_c, \\ \end{equation}$$

(72) $$\begin{equation} b = r_l \rho _l \beta _l V_\infty C_{pw}, \\ \end{equation}$$

(73) $$\begin{equation} c = r_c \rho _c \beta _c V_\infty \left(\left(1-\eta _{melt}\right)C_{pi}+\eta _{melt}C_{pw}\right), \\ \end{equation}$$

(74) $$\begin{equation} d = Q_{aero}+Q_{kl}+Q_{kc}+Q_{in}, \\ \end{equation}$$

(75) $$\begin{equation} e = \dot{m}_{in}C_{pw} \end{equation}$$

Evaluating Equation (67) at y = B + h, substituting a ₅ from Equation (70):

(76) $$\begin{eqnarray} \theta (B+h,t)=\frac{k_wT_f+(aT_a+bT_l+cT_c+eT_f+d)h}{k_w+(a+b+c+e)h} \end{eqnarray}$$

For the calculation of Q _ev in Equation (33), vapour pressure needs to be calculated. This is a function of θ(B + h, t). Therefore, θ(B + h, t) needs to be calculated through an iterative process until convergence in Q _ev is reached.

Substituting θ(B + h, t) from Equation (76) back into Equation (67):

(77) $$\begin{equation} a_5=\frac{a(T_a-T_f)+b(T_l-T_f)+c(T_c-T_f)+d+e(T_in-T_f)+f(T_a-T_f)}{k_w+(a+b+c+e)h} \end{equation}$$

Therefore, Equation (67) becomes:

(78) $$\begin{eqnarray} \theta (y,t)=a_5 y+T_f \end{eqnarray}$$

Integrating Equation (16) yields:

(79) $$\begin{equation} h=\frac{r_l \rho _l \beta _l V_\infty + \eta _{melt} r_c \rho _c \beta _c V_\infty + \alpha _{melt} r_c \rho _c \beta _c V_\infty + \dot{m}_{in}- \dot{m}_{ev}}{\rho _w}(t-t_g)-\frac{\rho _g}{\rho _w}(B_l-B_g) \end{equation}$$

This integration involves B _g and t _g constants which are the ice thickness and time at which water first appears in glaze ice. Substituting Equation (55) and (78) into Equation (17) yields:

(80) $$\begin{eqnarray} \rho _gL_f \frac{\partial B_l}{\partial t}&=&k_i\frac{T_f-T_{subs}}{B_l}-k_wa_5 \end{eqnarray}$$

Here, a ₅ can be calculated from Equation (77) and h from Equation (79). In order to find B _g (hence t _g ), substitute ∂B _l /∂t from Equation (42) into Equation (80). However, only the ice thickness due to liquid water coming from droplets or the liquid part of partially melted ice crystals is taken into account:

(81) $$\begin{eqnarray} L_f\left(r_l \rho _l \beta _l V_\infty + \eta _{melt} r_c \rho _c \beta _c V_\infty + \dot{m}_{in}- \dot{m}_{sub}\right)=k_i\frac{T_f-T_{subs}}{B_g}-k_wa_5 \end{eqnarray}$$

Rearranging:

(82) $$\begin{eqnarray} B_g=\frac{K_i(T_f-T_{subs})}{(r_l \rho _l \beta _l V_\infty + \eta _{melt} r_c \rho _c \beta _c V_\infty + \dot{m}_{in}- \dot{m}_{sub})L_f+k_wa_5} \end{eqnarray}$$

Note, a ₅ has to be calculated for h = 0 here because at B = B _g , h = 0:

(83) $$\begin{equation} a_5=\frac{a(T_a-T_f)+b(T_l-T_f)+c(T_c-T_f)+d}{k_w} \end{equation}$$

The time at which glaze ice first appears is found from Equation (42):

(84) $$\begin{eqnarray} t_g&=&\frac{\rho _r}{r_l \rho _l \beta _l V_\infty + r_c \rho _c \beta _c V_\infty + \dot{m}_{in}- \dot{m}_{sub}}B_g \end{eqnarray}$$

When we collect all unfrozen water in one term, by rewriting $\dot{m}_{i}$ :

(85) $$\begin{eqnarray} \dot{m}_{w}=\frac{\rho _wh}{t-t_g} \end{eqnarray}$$

(86) $$\begin{eqnarray} \dot{m}_{i}=\rho _g \frac{\partial B_l}{\partial t} \end{eqnarray}$$

As discussed in the erosion model part, water and ice erosion rates can be related with the total erosion rates as follows:

(87) $$\begin{eqnarray} \dot{m}_{er,i}&=&\dot{m}_{er}\frac{\dot{m}_{i}}{\dot{m}_{i}+\dot{m}_{w}} \end{eqnarray}$$

(88) $$\begin{eqnarray} \dot{m}_{er,w}&=&\dot{m}_{er}\frac{\dot{m}_{w}}{\dot{m}_{i}+\dot{m}_{w}} \end{eqnarray}$$

Water and ice mass accumulation rates can be updated after calculating the erosion rates:

(89) $$\begin{eqnarray} \dot{m}_{i,upd}&=&\dot{m}_{i}-\dot{m}_{er,i} \end{eqnarray}$$

(90) $$\begin{eqnarray} \dot{m}_{w,upd}&=&\dot{m}_{w}-\dot{m}_{er,w} \end{eqnarray}$$

Moreover, Trontin suggested that some of the accreted water does not run downstream and is trapped inside the accumulated ice layer^{(Reference Trontin, Blanchard and Villedieu28)}.

(91) $$\begin{eqnarray} \dot{m}_{out}&=&2\zeta \dot{m}_{w,upd} \end{eqnarray}$$

(92) $$\begin{eqnarray} \zeta &=&\frac{\psi _l+\eta _{melt}\epsilon _s\psi _c}{\psi _l+\epsilon _s\psi _c} \end{eqnarray}$$

The above relation is applicable only when ζ ⩽ 0.5. Otherwise, all the accreted water is assumed to run downstream.