Heat and water vapor transfer in the wake of a falling ice sphere and its implication for secondary ice formation in clouds

TL;DR

This study uses direct numerical simulations to analyze heat and water vapor transfer behind falling ice spheres, shedding light on supersaturation zones that may promote secondary ice formation in clouds.

Contribution

It introduces a detailed simulation framework for the coupled heat and vapor transfer around falling graupel, considering various settling regimes and their effects.

Findings

Wake supersaturation zones can form behind falling graupel.

Settling regime influences scalar field structures.

Modest Richardson numbers may affect wake dynamics.

Abstract

We perform direct numerical simulations of the settling of an ice sphere in an ambient fluid accounting for heat and mass transfer with the aim of studying in a meteorological context the case of falling graupel in humid air. The study is motivated by the fact that falling graupels in clouds are heated by the latent heat released during the accretion of liquid water droplets. They may therefore be considerably warmer than their surrounding and evaporate water vapor, which mixes with the surrounding air in the wake of the graupel, thereby creating transient zones of supersaturation there. The problem of a falling graupel is modeled as that of a heated sphere falling in a quiescent ambient fluid under the action of gravity. The coupling between the temperature and velocity fields is accounted for by the Boussinesq approximation. This problem can be parameterized by four parameters: the particle/fluid density ratio, the Galileo number, the Prandtl number, and the Richardson number. A separate scalar transport equation accounts for the vapor transport. Typical cloud conditions involve small temperature differences between the sphere and the surrounding, yielding relatively small Richardson numbers for both heat and mass transport. We give a special emphasis to the Galileo numbers 150, 170, 200 and 300 in order to analyze the specificities of each settling regime. The questions addressed in this study are mainly methodological and concern the influence of the settling regime and the mobility of the sphere on the structure of the scalar fields, the possible influence of modest Richardson numbers on the structure of the wake, and the possible application of this simulation framework to the investigation of the saturation in the wake of a falling graupel.

Tables

Table 1. Table 1: Maximal relative error between S w subscript 𝑆 𝑤 S_{w} and S w T superscript subscript 𝑆 𝑤 𝑇 S_{w}^{T} for different sets of boundary conditions.

$T_p$ ($°C$)	$T_{\infty }$ ($°C$)	Max relative error
$0$	$-15$	$7.02\times {10}^{-2}$
$0$	$-10$	$3.53\times {10}^{-2}$
$0$	$-5$	$1.38\times {10}^{-2}$
$-5$	$-15$	$3.41\times {10}^{-2}$
$-5$	$-10$	$1.23\times {10}^{-2}$
$-10$	$-15$	$1.08\times {10}^{-2}$

Equations

$$\tilde{\nabla }\cdot \tilde{\mathbf{u}}$$

$$\frac{\partial \tilde{\mathbf{u}}}{\partial \tilde{t}}+\left[(\tilde{\mathbf{u}}-{\tilde{\mathbf{u}}}_{\mathbf{p}})\cdot \tilde{\nabla }\tilde{\mathbf{u}}\right]$$

$$\frac{\partial \tilde{T}}{\partial \tilde{t}}+\left[(\tilde{\mathbf{u}}-{\tilde{\mathbf{u}}}_{\mathbf{p}})\cdot \tilde{\nabla }\tilde{T}\right]$$

$$\frac{\partial {\tilde{n}}_v}{\partial \tilde{t}}+\left[(\tilde{\mathbf{u}}-{\tilde{\mathbf{u}}}_{\mathbf{p}})\cdot \tilde{\nabla }{\tilde{n}}_v\right]$$

$$\tilde{T}$$

$${\tilde{n}}_v$$

$$R{i}_T$$

$$\frac{{\rho }_p}{{\rho }_{\infty }}\frac{\mathrm{d}{\tilde{\mathbf{u}}}_{\mathbf{p}}}{\mathrm{d}\tilde{t}}$$

$$\frac{{\rho }_p}{{\rho }_{\infty }}\frac{\mathrm{d}{\tilde{\mathbf{\omega }}}_{\mathbf{p}}}{\mathrm{d}\tilde{t}}$$

$$n_v$$

$${\mathbf{e}}_{\parallel }$$

$${\mathbf{u}}_h$$

$${\mathbf{e}}_h$$

$${\mathbf{e}}_{\perp }={\mathbf{e}}_{hz\perp }\times {\mathbf{e}}_{\parallel }.$$

$$N{u}_{loc}$$

$${\dot{Q}}_{cond}$$

$$Nu=\frac{4}{\pi D^2}{⟨{\int }_{\mathcal{S}}N{u}_{loc}\mathrm{d}S⟩}_t,$$

$$N{u}_{c1}$$

$$N{u}_{c2}$$

$${\tilde{u}}_d({\tilde{r}}_{\perp },{\tilde{z}}_{\parallel })$$

$$T_p=\{-10,-5,0\}°C,$$

$$T_{\infty }=\{-15,-10,-5\}°C,$$

$$S_w(\mathbf{x},t)$$

$$=$$

$$\frac{\mathrm{d}{S}_w}{S_w}$$

$$S_w^T=\frac{k_b\times \left(n_{v,\infty }+\tilde{T}\times (n_{v,p}-n_{v,\infty })\right)\times \left(T_{\infty }+\tilde{T}\times (T_p-T_{\infty })\right)}{e_{sat,w}(T_{\infty }+\tilde{T}\times (T_p-T_{\infty }))}.$$

$$S{t}_{mp}$$

$${\nabla }^{\ast }\cdot \mathbf{u}$$

$$\frac{\partial \mathbf{u}}{\partial t^{\ast }}+\mathbf{u}\cdot {\nabla }^{\ast }\mathbf{u}$$

$$\frac{\partial T}{\partial t^{\ast }}+\mathbf{u}\cdot {\nabla }^{\ast }T$$

$$\frac{\partial n_v}{\partial t^{\ast }}+\mathbf{u}\cdot {\nabla }^{\ast }{n}_v$$

$$\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot {\nabla }^{\ast }\mathbf{u}$$

$$A(\mathbf{x}({\mathbf{x}}^{\ast },t^{\ast }),t)$$

$${\nabla }^{\ast }A({\mathbf{x}}^{\ast },t^{\ast })$$

$$\frac{\partial A({\mathbf{x}}^{\ast },t^{\ast })}{\partial t^{\ast }}$$

$$\frac{\partial A({\mathbf{x}}^{\ast },t^{\ast })}{\partial t^{\ast }}+\mathbf{u}({\mathbf{x}}^{\ast },t^{\ast })\cdot {\nabla }^{\ast }A({\mathbf{x}}^{\ast },t^{\ast })$$

Full text

1   Introduction

Convection associated with mass transfer around a settling sphere is of relevance for many industrial and environmental systems such as clouds or combustion chambers. In the case of a graupel falling in a cloud an additional effect becomes important. The ice particle can be considerably warmer than the surrounding due to the latent heat released by the riming process, i.e. the freezing of liquid cloud droplets colliding with the falling ice particle. The diffusion of heat and vapor in the wake of the ice ball can lead to local vapor supersaturation that may even lead to secondary ice nucleation in its wake through heterogeneous ice nucleation mechanisms, if ice nuclei are present. Heterogeneous ice nucleation, in contrast to homogeneous ice nucleation, requires the presence of another substance usually called ice nucleating particle (INP) and is more likely to occur in the present context since homogeneous ice nucleation requires very low temperatures (below -38$\degree C$) and large relative humidity (above 145$\%$). The conditions up to which this heterogeneous nucleation can occur depend on both temperature and water content as well as the type of ice nuclei (Hoose and Möhler,, 2012). Heterogeneous ice nucleation can take place under various modes depending on the atmospheric condition and the type of INP (Vali et al.,, 2015): ice can directly form from water vapor (deposition) or from supercooled liquid (freezing) with different possible freezing mechanisms (immersion, contact, condensation freezing). A description of the advection-diffusion of both temperature and water content in the wake of a particle is therefore to be seen as a first step for a better understanding of the fundamental processes involved in secondary ice nucleation induced by a falling graupel.

The case of a single heavy sphere falling in a quiescent environment, without accounting for heat or mass transfer, could seem to be simple but is known to feature complex properties (Ern et al.,, 2012; Jenny et al.,, 2004). It can be parametrized by two dimensionless quantities, as described in Jenny et al., (2004), namely the ratio between particle and fluid density $\rho_{p}/\rho_{\infty}$ and the Galileo number $Ga=u_{g}D/\nu$, where $u_{g}=\sqrt{\lvert(\rho_{p}/\rho_{\infty}-1)g\rvert D}$ is the gravitational velocity, $D$ the diameter of the particle, $g$ the gravitational acceleration and $\nu$ the kinematic viscosity of the fluid. Figure 1 summarizes the different settling regimes to be expected for a density ratio $\rho_{p}/\rho_{\infty}=10$: for the lowest Galileo number the particle falls on a straight vertical path and the wake is steady axisymmetric. When $Ga$ increases, above a critical value $Ga_{c1}$, the particle follows an oblique path and the wake loses its axisymmetry to become planar reflectional symmetric. In this steady oblique regime the value of $Ga_{c1}$ is independent of the density ratio (Fabre et al.,, 2012). Steadiness is lost when $Ga$ becomes larger than $Ga_{c2}$ where a Hopf bifurcation occurs: the wake oscillates in time while the particle continues to fall in a plane. This is the oblique oscillating regime (Jenny et al.,, 2004). The flow becomes then chaotic when $Ga$ further increases above $Ga_{c3}$ and the particle motion is fully three-dimensional. The density ratio has been shown to have negligible effect on the transition scenario for values above 2.5 (Zhou and Dušek,, 2015). In the case of the flow around a fixed sphere, which can be thought of as a case for which $\rho_{p}/\rho_{\infty}$ tends to infinity, the transition scenario of the wake structure is parametrized by the Reynolds number $Re=u_{\infty}D/\nu$, where $u_{\infty}$ is the relative velocity of the unperturbed fluid. This scenario is similar to that of free spheres of large density ratio, and it has been investigated in many experimental and numerical works (Ormières and Provansal,, 1999; Johnson and Patel,, 1999; Ghidersa and Dušek,, 2000).

Turning now to the case of flows with heat or mass transfer, the literature mostly focused on configurations where the ratios between scalar diffusivity and fluid viscosity (namely the Prandtl number $Pr=\nu/\mathcal{D}_{T}$ for the temperature and the Schmidt number $Sc=\nu/\mathcal{D}_{m}$ for the mass) is equal to $0.7$, $1$ or $7$. The numerical work of Bagchi et al., (2001) for passive scalar transport featuring $Pr=0.72$ and a fixed sphere with Reynolds number ranging from 50 to 500, showed that the structure of the scalar field is highly dependent on the structure of the flow field with a clear influence of three-dimensionality and vortex shedding. Ignoring 3D effects indeed leads to an underprediction of drag and an overprediction of the Nusselt number. Other studies from the literature explored the influence of density variation due to temperature inhomogeneities, with more attention given to the configurations of assisting and opposing flow (i.e. where gravity points in the same or in the opposite direction of the unperturbed flow respectively). This influence of temperature variations on the buoyancy term is usually parametrized by the Grashof number $Gr=\beta g(T_{p}-T_{\infty})D^{3}/\nu^{2}$ (where $\beta$ is the thermal expansion coefficient, $T_{p}$ the temperature of the sphere and $T_{\infty}$ the temperature of the incoming flow) or Richardson number $Ri$, that are linked by the relation $Ri=Gr/Re^{2}$. It has been observed in the literature that this can induce modifications of the structure of the flow and of the heat transfer coefficient (Kotouč et al.,, 2008, 2009; Bhattacharyya and Singh,, 2008). Bhattacharyya and Singh, (2008) studied the influence of buoyancy for a Reynolds number ranging from $1$ to $200$ and a Richardson number ranging from [math] and $1.5$. They observed the development of a recirculating eddy in the downstream of the sphere which further collapses when the temperature of the sphere increases, consequently leading to the development of a buoyant plume above the sphere. A main effect is the delay of the flow separation with $Ri$. Kotouč et al., (2008) also observed that convection tends to stabilize the flow in the assisting flow configuration by preventing detachment of the boundary layer and the formation of a recirculation zone, while in the case of opposing flow a destabilizing effect is observed with a decrease of the critical Reynolds number for the onset of recirculation (Kotouč et al., (2009)). The transition scenario remains unchanged for weakly heated spheres, while new regimes appear for larger values of the Richardson number (of the order of 0.3 Kotouč et al.,, 2008, 2009). One can cite for example steady flow regimes with four or six vorticity threads and respectively two and three symmetry planes. Even a rather weak level of buoyancy with $Ri=0.1$ appeared to have an influence on the critical Reynolds number at which the transition between the different regimes appears (Kotouč et al.,, 2009). As expected, buoyancy can also lead to a modification of the drag, lift and transfer coefficients. An increase of the drag coefficient with $Ri$ in the assisting configuration has been observed by Bhattacharyya and Singh, (2008) , Kotouč et al., (2009), while the lift coefficient seems to decrease with $\left|Ri\right|$ (Kotouč et al.,, 2009). Concerning the heat transfer, an influence is only visible for high values of the Grashof number.

The references mentioned above focused on the influence of buoyancy on fixed spheres and very little attention has been given to the case of freely falling objects. Dan and Wachs, (2010) compared the terminal velocity of a sphere for Reynolds numbers ranging between $10$ and $130$ and Grashof numbers $Gr=\{-100,\;0,\;100\}$ and observed that a cold sphere settles faster than a warm one, with larger modifications induced at lower Reynolds numbers. Such modifications on the wake structure can further influence particle trajectories as well as the way that they might interact. Indeed, the numerical work of Gan et al., (2003) showed that the equilibrium position and the trajectory of a particle settling in a vertical channel might change with the Grashof number. They further observed that two particles settling in a channel tend to separate if they are colder than the fluid while they would aggregate if hotter.

Numerical simulations have been employed more recently, in a meteorological context, to investigate the flow around falling graupels of different shapes (Wang and Kubicek,, 2013) and the ventilation coefficient of a spherical hailstone (Cheng et al.,, 2014). Both of them neglected the influence of buoyancy, and set the graupel to be fixed.

Based on this several questions can be raised, which will be the object of the current study:

• •

As will be seen in the next section, the typical Galileo numbers for falling graupels are relatively large (from $\mathcal{O}(10^{2})$ to $\mathcal{O}(10^{3})$), corresponding to the chaotic regime. Therefore the question arises: is it necessary to consider the chaotic regime when studying secondary ice formation? For this purpose we will detail the characteristics of the thermal wake, and will base our analysis on the description of the wake in a tilted coordinate system similar to Uhlmann and Dušek, (2014). This is, to our knowledge, the first time that this has been proposed for a thermal wake.

• •

The density ratio representative for falling graupel is typically large (several hundreds), meaning that the density might be large enough for the system to be equivalent to a fixed sphere (i.e. $\rho_{p}/\rho_{\infty}\rightarrow\infty$). The second question which will be addressed concerns thus the influence of the density of the sphere on the structure of the thermal wake.

• •

The third point concerns the influence of buoyancy effects: The typical thermal Richardson number of the system is expected to be small ($\mathcal{O}(10^{-4})$), but the large Galileo number implies that the corresponding Grashof number is not necessarily small. Since we consider lower Galileo numbers than the largest values observed in clouds the Grashof numbers investigated here are smaller than for clouds. We suggest therefore to investigate configurations featuring larger Richardson numbers ($\mathcal{O}(10^{-3})$) and above to test if buoyancy can modify the wake. It is not clear whether buoyancy would potentially affect the system, since the literature showed that weak buoyancy ($Ri=0.1$) can modify the flow structure (Kotouč et al.,, 2009). This issue will also be addressed in the current study.

• •

We will be then concerned more specifically by secondary ice nucleation and will focus the discussion on the structure of the saturation field. The first question which we will address is purely methodological and concerns the necessity to separate the transport equation of heat from the transport equation of mass. As it will be described in the next sections, it is important to take into account the variations of both temperature and water vapor in the flow induced by the presence of the graupel. The transport equations of both scalars can be formulated in non-dimensional form by the same type of advection-diffusion equation, with the same set of boundary conditions, the only difference lying on the difference of diffusivities. As those diffusivities are relatively similar, it might be unnecessary to solve two separate transport equations for the two scalars, which will be discussed in the current study.

• •

The last point concerns a qualitative description of the saturation that can be reached in the wake of the sphere: does the graupel increase the level of saturation in its vicinity such that it could have an effect on the ice content (according to the values indicated by Hoose and Möhler,, 2012)?

Let us finally emphasize that the goal of the current study is not to investigate in details the question of secondary ice nucleation that might be induced, but to discuss the methodological framework that could be used to address this question and provide a new database which can be further used to study the possibility of condensation/deposition in the wake (Prabhakaran et al.,, 2017).

The paper is organized as follows: the methodology followed in the current paper is presented in section 2, then we will describe the structure of the velocity and temperature fields for different settling regimes first by neglecting the buoyancy effect (in section 3) and then for Richardson numbers, which we will define here as $Ri_{T}=Gr/Ga^{2}$, ranging from $0.001$ to $0.1$ (in section 4). Both sections 3 and 4 concern purely non-dimensional scalar and velocity fields, while section 5 will be devoted to the description of the saturation field and the necessity to separate both heat and mass transport equations for the computation of the saturation with respect to ice in the wake of a graupel, for given sets of temperature conditions. We will finally draw conclusions with respect to the questions addressed in the paper.

2   Methodology

We propose to consider a system representing a sphere of diameter $D$ falling in quiescent moist air as sketched in figure 2. We follow the Boussinesq hypothesis and model the fluid with constant viscosity $\nu$, and a density $\rho$ which depends linearly on the local variations of temperature $T$. The density variations due to variable vapor concentration are negligible in the present context since they are about two orders of magnitude smaller.

The sphere is supposed to have a uniform and constant temperature $T_{p}$, the incoming flow has a constant velocity $\mathbf{u_{\infty}}$, a constant temperature $T_{\infty}$ and, consequently, a constant density $\rho_{\infty}$. The system is represented by the Navier-Stokes equations (eq. 1 and 2) coupled to one advection-diffusion equation for the transport of heat (eq. 3), under the Boussinesq approximation for the representation of the buoyancy effects. A second (passive) scalar field representing the water vapor concentration is transported according to the advection-diffusion (eq. 4). The system is described in a dimensionless form with reference quantities chosen as in Kotouč et al., (2009) and Jenny and Dušek, (2004), where $\mathbf{u_{p}}$ represents the velocity of the sphere and $\mathbf{u}$ corresponds to the fluid velocity with respect to a fixed frame $(0,x^{*},y^{*},z^{*})$ and evaluated on a grid $(x,y,z)$ moving with the sphere center as represented in figure 2. Here lengths are made non-dimensional with respect to the particle diameter, velocity components with the gravitational velocity $u_{g}$, timescale with the gravitational timescale $D/u_{g}$, pressure with $\rho_{\infty}u_{g}^{2}$. As a convention, henceforth we denote non-dimensional quantities with a tilde. The dimensionless temperature and vapor concentration are defined as in equations 5-6. Under these assumptions the governing equations read as follows (for more details on the formulation the reader is referred to the appendix A):

[TABLE]

where $\mathbf{k}$ is the normalized unit vertical vector such that $\mathbf{g}=-\lvert\mathbf{g}\rvert\mathbf{k}$. We introduce the dimensionless Richardson number $Ri_{T}$ defined as

[TABLE]

The sphere is free to move and its motion and rotation are governed by the following equations (with $\boldsymbol{\omega_{p}}$ the angular velocity of the sphere, $\tilde{P}$ the hydrodynamic pressure without the hydrostatic part and ${\boldsymbol{\tilde{\tau}}}$ the dimensionless viscous stress tensor whose components are $\tilde{\tau}_{ij}=\left[{\partial\tilde{u}_{i}}/{\partial\tilde{x}_{j}}+{\partial\tilde{u}_{j}}/{\partial\tilde{x}_{i}}\right]/Ga$):

[TABLE]

Equations 1-9 are solved in a strongly coupled fashion with the method of Kotouč et al., (2008) and Jenny et al., (2004), namely a spectral/spectral-element method on a cylindrical domain with axis parallel to the gravity vector, and the sphere center placed at the domain axis. We use a spectral element discretization in the axial/radial $(z,r)$ plane and a Fourier decomposition in the azimuthal direction which was shown to be an efficient choice for the simulation of the different transitions to be observed (Ghidersa and Dušek,, 2000). We have used the same computational domain as Zhou and Dušek, (2015) with an extension of $37D$ and $8D$ in the axial and radial directions, respectively, and a discretization into 245 elements with 6 collocation points per element in each direction. Figure 3 shows the structure of this mesh in a $(z,r)$ plane. We used a discretization in the azimuthal direction with a truncation at the 7th Fourier mode for steady cases and 10th for unsteady ones. Tests have been performed with 8 collocation points in each direction and a truncation up to the 16th mode and no significant difference has been observed.

The Prandtl number is taken equal to $Pr=0.72$ throughout this work, meaning that the thermal and mass boundary layers should be resolved with the same quality as the background flow and no extra refinement is required here.

A no-slip boundary condition is used at the particle’s surface and the incoming flow is set constant. A zero-gradient boundary condition is applied at the outflow face for the velocity field as well as for the temperature. At the lateral boundaries of the domain a no-stress, adiabatic and hermetic Neumann condition is imposed with zero pressure.

In a meteorological context, the transport of water vapor has to be taken into account. It is indeed important to consider the local saturation $S$ which is defined as the ratio between the local vapor partial pressure $e$ and the saturation pressure $e_{sat}$ (which in our case we assume to be only function of temperature). The local partial pressure and vapor concentration are linked according to the relation

[TABLE]

with $k_{b}$ the Boltzmann constant. The saturation will therefore depend on both variation of temperature and vapor concentration. The transport of $n_{v}$ is governed by an advection-diffusion (equation 4) with a diffusion coefficient yielding $Sc=\nu/\mathcal{D}_{m}\sim 0.63$ (Young,, 1993) and a non-dimensional number concentration number defined according to equation 6. The boundary conditions $n_{v,\infty}$ and $n_{v,p}$ are also set constant and uniform both at the inflow and on the sphere surface. The corresponding Richardson number is of the order $\mathcal{O}(10^{-6})$ such that density variations due to vapor transport can be neglected. Microdroplets will often collide with the falling ice particle, and freeze when they come in contact with the ice particle. This leads to an increase of the sphere’s surface temperature due to the release of latent heat during freezing. For this reason we propose to focus here on configurations for which $T_{p}>T_{\infty}$, i.e. a sphere which is warmer than the surroundings. Concerning the boundary conditions on the partial pressure (at infinity and on the particle) we assume that the saturation is equal to unity with respect to the liquid phase for the incoming flow and of unity with respect to ice at the sphere’s surface, yielding $n_{v,p}=e_{sat,i}(T_{p})/k_{b}T_{p}$ and $n_{v,\infty}=e_{sat,w}(T_{\infty})/k_{b}T_{\infty}$ (with $e_{sat,w}$ the saturation vapor pressure with respect to liquid water and $e_{sat,i}$ the saturation vapor with respect to ice). The boundary condition at infinity is chosen to reflect the presence of supercooled droplets in mixed phase clouds. Please note that the boundary condition at the particle is capable of describing dry growth as well as wet growth, since the latter occurs at $T_{p}\approx 0\degree C$, for which $e_{sat,i}=e_{sat,w}$. We used the correlations of Murphy and Koop, (2005) to set the evolution of the saturation partial pressure $e_{sat,j}$ for the phase $j$ as a function of the temperature. For more details on this step the reader is referred to appendix (B).

The results in section 3 and 4 are presented in a dimensionless form and do not include the evolution of the dimensionless concentration number. The difference between the Schmidt and Prandtl numbers are indeed too small to induce significant difference between the non-dimensional temperature and concentration fields. Contrarily, the results of section 5 dealing with the description of the evolution of the saturation field include the solution for the transport equation for $\tilde{n}_{v}$.

A part of the analysis in the wake of the sphere is made with the help of a coordinate system that is not only centered at position of the particle, but also follows the orientation of its settling velocity, as previously used by Uhlmann and Dušek, (2014). For this purpose we first introduce the velocity of the particle with the Cartesian coordinates $\mathbf{u}_{p}=(u_{px},u_{py},u_{pz})^{T}$ and the corresponding unit vector

[TABLE]

We also introduce the projection of the relative velocity onto the $(x,y)$ plane $\mathbf{u}_{h}$ and its associated unit vector $\mathbf{e}_{h}$, namely

[TABLE]

then we define another horizontal vector $\mathbf{e}_{hz\bot}=\mathbf{e}_{z}\times\mathbf{e}_{h}$, and finally a unit vector which is now perpendicular to both $\mathbf{e}_{\parallel}$ and $\mathbf{e}_{hz\bot}$ and which is defined by

[TABLE]

An example of the structure of this new Cartesian coordinate system is shown in figure 4 for a configuration in the oscillating oblique regime. We will preferentially consider this coordinate system to describe the structure of the wake. We denote the position along the $\mathbf{e}_{\parallel}$ axis as $\tilde{z}_{\parallel}$ and along the $\mathbf{e}_{\bot}$ axis as $\tilde{x}_{\bot}$. A tilted cylindrical coordinate system with $\tilde{z}_{\parallel}$ as a main axis will also be used, and we denote the corresponding radial position by $\tilde{r}_{\bot}$, as will be used in section 3. Statistics which imply time averaging of the temperature or velocity field, presented in this coordinate system, are performed as follows: the instantaneous velocity and temperature fields are first projected onto this instantaneous inclined frame for each snapshot and then averaged over time.

We focus here on a density ratio $\rho_{p}/{\rho_{\infty}}=10$ and investigate Galileo numbers ranging from 10 to 300 such that all different regimes to be observed with such density ratio will be covered. A special emphasis is given on $Ga=150$, $170$, $200$ and $300$ in order to describe the specificities of each regime. In the case of a falling graupel the density ratio would typically range between 330 and 650 with a diameter between 0.5 and 6mm, leading then to Galileo numbers ranging from 60 and 2800 (Pruppacher and Klett,, 2010). In the present study we propose to focus on intermediate Galileo numbers in order to test the influence of the wake regime on the structure of the scalar field. It would be necessary in the future to get also more insight on the influence of turbulence or collective effects but this would go beyond the scope of the current paper and we neglect those effects in the current study. We set first $\rho_{p}/{\rho_{\infty}}=10$ and not $\rho_{p}/\rho_{\infty}=600$ as it has the advantage of reducing the characteristic timescale of the particle and with it the duration of the simulations, but is large enough to behave in a similar way as $\rho_{p}/{\rho_{\infty}}\simeq 600$. A real density ratio of 600 can be seen as an intermediate point between $\rho_{p}/{\rho_{\infty}}=10$ and a fixed sphere (i.e. $\rho_{p}/\rho_{\infty}=\infty$). In order to test the necessity to reproduce the influence of the mobility in a meteorological context we investigated the density ratios $\rho_{p}/\rho_{\infty}=\{1.5;\;10;\;\infty\}$.

We focus on ambient temperature $T_{\infty}$ ranging between -15$\degree C$ and -5$\degree C$ and sphere temperature ranging between -10$\degree C$ and 0$\degree C$, the resulting thermal Richardson number for this system is therefore estimated to vary between $0.5\times 10^{-4}$ and $1.75\times 10^{-4}$. The corresponding Grashof number $Gr=RiGa^{2}$ would range between $Gr=\mathcal{O}(1)$ and $Gr=\mathcal{O}(10^{3})$. The highest Galileo number considered in the current study (300) combined with the highest expected Richardson number value yield $Gr\sim 15$. Considering the real values of $Ri_{T}$, with this Galileo number range, would not yield any visible effect on the flow. For this reason we propose to explore larger ranges and consider $Ri_{T}=(0;\;0.001;\;0.05;\;0.1)$, in order to cover larger values of the Grashof number.

3   The case of passive scalar transport

We propose in this section to focus on the structure of the scalar field in the configuration where no buoyancy effects are accounted for ($Ri_{T}=0$), for a Prandtl number $Pr=0.72$, with the aim of providing a first description on the influence of the settling regime as well as the density ratio on the structure of velocity and the scalar field. Figure 1 gives a visual impression of the velocity and temperature field for different Galileo numbers and for $\rho_{p}/\rho_{\infty}=10$. The main characteristics of the flow regime are easily recognizable on both velocity and temperature field: $Ga=150$ features an axisymmetric wake while for $Ga=170$ and $Ga=200$ the wake is oblique, with oscillations in space only visible for $Ga=200$. The structure of the velocity and thermal wakes are chaotic for $Ga=300$. Figure 5 gives also a visual impression of the different zones of the flow.

3.1 Recirculation length

We start our description with the evolution of the extent of the recirculation region behind the sphere as well as the evolution of the mean transfer coefficients as they provide important information for the following description. Figure 6 shows the evolution of the vertical component of the instantaneous velocity as well as the instantaneous temperature field in planes containing the tilted axis $x_{\perp}$ and $x_{hz\perp}$ for different Galileo numbers. It highlights the influence of the recirculation region on the temperature distribution in the wake. This region is indeed characterized by a toroidal vortex whose projection onto two perpendicular planes represents counter-rotating vortices that affect the temperature as sketched on figure 5: cold fluid is transported from the outer shear region towards the rear stagnation point. Conversely warm fluid is transported from the back of the sphere to the shear region, leading to local maximum of temperature in the core of the vortices and local minimum between them. The influence of the symmetry breaking is also visible. If one considers for instance the steady oblique case, projected into the $(x_{\perp},z_{\parallel})$ plane ($Ga=170$ in figure 6b,f), the recirculation is mainly located on the side $x_{\perp}>0$ and therefore cold ambient fluid is primarily entrained on this side, making the temperature globally lower in this region of the wake than on the side $x_{\perp}<0$.

Figure 7 shows the recirculation length which is defined as the largest distance to the sphere surface for which the time averaged streamwise velocity is equal to zero, estimated on the three dimensional contours of the mean velocity field (with this definition the position of the corresponding point can thus be located away from the wake axis as can be clearly seen in figure 6). We make the distinction here between two time averaged flow fields: The velocity field is either averaged in the original coordinate system $(x,y,z)$ or after rotation to the oblique coordinate system $(\mathbf{e}_{\parallel},\mathbf{e}_{hz\bot},\mathbf{e}_{\bot})$. It shows that the recirculation length first increases with the Galileo number until unsteadiness is reached. It then decreases in the oblique oscillating regime for further increasing in the chaotic regime. The increase of $L_{r}$ with $Ga$ is consistent with the observation, for a fixed sphere, of Magnaudet et al., (1995) and Bouchet et al., (2006) for the steady axisymmetric regime, but they observe a first decay once the steady non-axisymmetric flow develops while here the decay starts with the onset of unsteadiness. This discrepancy is most probably due to estimation of the recirculation length not on the wake axis but as the largest distance of the isocontour $\tilde{u}_{z}=0$ to the sphere surface. We indeed computed the length on the axis and found similar trend as in Bouchet et al., (2006). The decrease of $L_{r}$ of the time-averaged flow for the oblique oscillating case is a signature of the vortex shedding and the oscillation of the counter-rotating vortices observed in the plane $(x_{\perp},z_{\parallel})$. For more clarity we omitted the configurations featuring $\rho_{p}/\rho_{\infty}=1.5$ or a fixed sphere on the figure, for which we observed either a negligible influence on $L_{r}$ (fixed sphere), or a small decrease of $L_{r}$ ($\rho_{p}/\rho_{\infty}=1.5$).

3.2 Transfer coefficient

We use the same definition of the local Nusselt number as Kotouč et al., (2008), namely

[TABLE]

where $\dot{Q}_{cond}$ is the theoretical purely conductive heat flux:

[TABLE]

with $\lambda$ the thermal conductivity of the fluid. The mean Nusselt number is then defined as the integral of $Nu_{loc}$ over the surface of the sphere, which is time-averaged in unsteady flows:

[TABLE]

where the operator $\left<\;.\;\right>_{t}$ refers to time averaging.

Figure 8 gives the evolution of the mean heat transfer coefficient as a function of the Galileo number and a comparison with the correlation of Ranz and Marshall (Ranz and Marshall,, 1952), namely:

[TABLE]

which can be used for Reynolds number up to $Re=5\times 10^{4}$, as well as another correlation which has been validated against direct numerical simulation for $1<Re<400$ and $0.25<Pr<100$ (Clift et al.,, 1978):

[TABLE]

Applied to this study, the correlations would refer to a time-averaged Reynolds number, defined in two ways. The first one accounts for the three components of the sphere’s velocity $Re(t)={\left|\mathbf{u}_{p}\right|D}/{\nu}$, while the other accounts only for its vertical component $Re_{T}(t)={\left|u_{p,z}\right|D}/{\nu}$. Using the balance of external and hydrodynamical forces one can show that the mean Reynolds number $\left<Re_{T}\right>_{t}$ is linked to the Galileo number according to the relation $\sqrt{3C_{D}/4}\left<Re_{T}\right>_{t}\sim Ga$ (Jenny et al.,, 2004), so that the mean Reynolds number can be estimated as a function of the Galileo number once assuming that the drag coefficient follows a known correlation of Schiller and Naumann, (1935) with a difference of at most one percent. The combination of the correlations provides a good estimate of the Nusselt number as a function of the Galileo number, with almost no difference observed between the values computed with the actual time averaged Reynolds number or the with the one obtained from the drag correlation (figure omitted). The standard deviation of those coefficients are larger for the chaotic regimes than the oblique oscillating one with a standard deviation of approximately one percent for the chaotic regimes as can be for example observed on the time variations of $Nu$ and $Re$ at $\rho_{p}/\rho_{\infty}=10$ on figure 8. Interestingly, the time averaged value of the Nusselt coefficient shows little dependence on the density ratio, which is in accordance with the formulation of the correlations and the relation between the drag coefficient and the Reynolds number.

3.3 Structure of the scalar field and influence of the density ratio

We now move on to the structure of the scalar field. We characterize the amplitude of the wake with the decay of the scalar field downstream of the sphere. This decay is represented on figure 9 both for the scalar field and for the velocity deficit $\tilde{u}_{d}$. The latter is defined as

[TABLE]

where $u_{z_{\parallel}}=\mathbf{u}\cdot\mathbf{e}_{\parallel}$ and $\left<\;.\;\right>_{\theta t}$ refers to the average over the azimuthal direction and in time. The evolution on the axis $\tilde{r}_{\perp}=0$ of the velocity deficit as well as the mean temperature are represented on figure 9.

We observe a negligible influence of the wake regime in the rear of the sphere as both fields seem to be dominated by the influence of the recirculation: the velocity defect first increases because of the reverse fluid motion on the wake axis in the recirculation zone and then it further decreases. The corresponding temperature evolution is shown on figure 9(b). It features a rapid decrease within a distance of $0.2D$ to attain a value approximately equal to $0.4$ as already observed by Bagchi and Kottam, (2008). After this zone, both velocity deficit and temperature decrease, with decay rates that differ from the expected trend ($z^{-1}$ for a laminar wake and $z^{-2/3}$ for a turbulent wake), most probably because this asymptotic state would be expected at larger distances from the sphere.

The influence of the recirculation zone is also visible downstream of the sphere with the cross-stream evolution of the temperature at different positions as shown in figure 10. Due to the recirculation, cold fluid is transported towards the downstream stagnation point on the sphere which will induce the large gradient observed on the temperature decay (fig. 9(b)), and from there hot fluid is transported to the shear region creating a local peak of temperature, in accordance with the temperature contours observed by Bagchi et al., (2001), Bhattacharyya and Singh, (2008) and Bagchi and Kottam, (2008). The amplitude of this peak decreases for for $Ga=150,170$ and $200$ and then slightly increases for the last Galileo number. This first decrease is consistent with the increase of the recirculation length for $Ga=\{150,\;170,\;200\}$ as visual inspections of the temperature field have shown that recirculation vortices with large streamwise extensions lead to colder fluid in the vortex core. The small increase of the peak observed for $Ga=300$ compared to $Ga=200$, which is in contradiction with the increase of $L_{r}$, can be attributed to the decrease of the temperature in the inter-vortex region observed from $Ga=200$ to $Ga=300$ for the streamwise evolution of the temperature (fig. 9): fluid is getting colder in the intra-vortex region due to the increase of $L_{r}$, but it also gets colder in the inter-vortex region, decreasing with this the normalized amplitude of the temperature peak. After the recirculation zone, the temperature becomes a monotonic function of $\tilde{r}_{\bot}$ and attains a self-similar Gaussian profile for $\tilde{z}_{\parallel}>5$ as shown in figure 10, which is consistent with the observations made by Bagchi and Kottam, (2008) who found Gaussian profiles in the thermal wake at streamwise distances of $6.5D$ and $13D$.

The influence of the regime is more visible away from the centerline as can be clearly seen on the streamwise evolution of the half width of the wake (fig. 11). We introduce this quantity for the velocity deficit and the temperature, respectively $L_{hw}^{u_{d}}(\tilde{z}_{\parallel})$ and $L_{hw}^{T}(\tilde{z}_{\parallel})$. We define it as the radial coordinate $\tilde{r}_{\bot}$, for a given location $\tilde{z}_{\parallel}$ on the axis where the corresponding quantity is equal to $e^{-1/2}$ times the value on the centerline (this definition has also been used by Bagchi and Balachandar,, 2004; Legendre et al.,, 2006; Bagchi and Kottam,, 2008). The evolution of $L_{hw}$ is represented on figure 11 and shows a similar trend for the scalar quantities and the velocity deficit: the width first decreases with $\tilde{z}_{\parallel}$ which is consistent with Legendre et al., (2006) and then increases almost linearly as observed with most free shear flows. The width of the hydrodynamic and thermal wakes follow different trends with $Ga$ as $L_{hw}^{T}$ shows similar size for the first three Galileo numbers while $L_{hw}^{u_{d}}$ is lower at $Ga=150$ than $Ga=170$ and $Ga=200$. We have also included a case featuring $Ga=100$ to emphasize the effects of both diffusion and wake regime, since diffusion should tend to increase the growth of the half width with the distance $\tilde{z}_{\parallel}$ and the structure of the wake should affect this width at the end of the recirculation region. The evolution of $L_{hw}^{u_{d}}$ confirms the intuition that the half width of the velocity deficit should increase much faster for $Ga=100$ than $Ga=150$ due to diffusion effects, which is in accordance with the observations of Bagchi and Balachandar, (2004) who observed a larger width at $Re=107$ than at $Re=241$ and $Re=261$. The half width is surprisingly larger at $Ga=170$ than $Ga=150$ although diffusive effects should be expected to be slightly larger at $Ga=150$. The difference in the viscosities are indeed not large enough to induce a significant difference in the slope of $L_{hw}^{u_{d}}$ while the obliqueness of the wake at $Ga=170$ tends to increase its width in the recirculation region once averaging in the azimuthal direction. In the steady oblique regime, the wake projected onto the oblique coordinate system is not axisymmetric (see for example the visualizations in Uhlmann and Dušek,, 2014). The unsteadiness of the wake will also act as diffusive effects and tend to increase the width of the average wake, as we observe here. This last point is in accordance with the simulations of Legendre et al., (2006) in the case of a bubble for which a thicker wake is observed at $Re=500$ compared to $Re=200$. Considering now the half width of the temperature, this trend is less pronounced (fig. 11(b)). The influence of the viscosity is clearly visible with a larger thermal wake obtained at $Ga=100$ compared to the other Galileo numbers. But the configurations with $Ga=150,\;170$ and $200$ feature similar wake width. This can be attributed to the low difference in diffusivities between those cases. The enlargement of the wake attributed to the unsteadiness of the regime is mostly visible on the chaotic regime.

The density ratio does not have an important influence on streamwise evolution of the velocity deficit and on the temperature (figure omitted). The main difference can be attributed to the slow decrease of the recirculation length to be observed for $\rho_{p}/\rho_{\infty}=1.5$ compared to $\rho_{p}/\rho_{\infty}=10$, but which is not observed when comparing the larger density ratio to the fixed sphere configuration. The half width shows a larger sensitivity to the sphere mobility as represented on figure 12 which displays an increase of the half width as the density ratio decreases. This can be attributed to the tendency for the wake to align differently than the particle trajectory when the density ratio decreases. With this the angle between the wake and the trajectory might increase which causes an enlargement of the wake once averaging. Another difference at the lowest density ratio is the convergence of the statistics that has been more difficult to reach for the half widths, for similar observation time, as can be seen in figure 12. Therefore here the comparison should remain mostly qualitative. The scaled evolution of the temperature as a function of the radial distance shows little influence of the density ratio (fig 10), apart from slightly larger peak of temperature in the recirculation region obtained for $\rho_{p}/\rho_{\infty}=1.5$, as a consequence of the small decrease of the extent of this zone. The evolution of $L_{hw}^{T}$ confirms that density has less influence on the thermal wake at short distances than on the velocity deficit.

It appears therefore that in a meteorological context it is important to investigate configurations featuring a chaotic regime as it affects most of the statistics of the scalar field, starting from the characteristics of the recirculation zone. The mobility of the sphere can be neglected, or density ratio has to be taken as large.

4   The case of active scalar transport

Figure 13 gives a visual impression of the three-dimensional structure of the velocity and temperature fields for the same Galileo numbers as described in the previous section, here for $Ri_{T}=0.1$. It first appears that buoyancy effects tend to elongate the thermal wake in the direction of gravity, due to the transport of warmer fluid in the upward direction. Another major point is the modification of the regime of the wake for $Ga=170$ and $Ga=200$. We indeed observe that axisymmetry is preserved at $Ga=170$ while this case is steady oblique for $Ri_{T}=0$. In the same way at $Ga=200$ the wake does not feature space oscillations anymore and lost its unsteadiness. Buoyancy, in the assisting flow configuration, therefore tends to stabilize the wake, in accordance with the previous observations of Kotouč et al., (2008) and Kotouč et al., (2009). Figure 14 gives the map of the different regimes observed for $\rho_{p}/\rho_{\infty}=10$ and for $Ri_{T}$ ranging from 0 to 0.1. The stabilizing influence of buoyancy is clearly visible as the different transition thresholds are pushed towards larger values of $Ga$.

It is known from the literature that a buoyant plume will be induced at larger Richardson numbers, and the competition between the upward plume and the original downward motion tends to decrease the extent of the recirculation region. The evolution of this recirculation as a function of the Galileo number for those new Richardson numbers is included in figure 7. The trend of the evolution of $L_{r}$ with $Ga$ is conserved: it will first increase with $Ga$ as long as the regime remains steady, then unsteadiness will be associated with a decay of the recirculation length until the wake becomes chaotic and then $L_{r}$ increases with $Ga$. The recirculation length decreases with $Ri_{T}$, in accordance with the observations of Bhattacharyya and Singh, (2008). This decrease can be seen as the signature of the superposition of two effects: the initial recirculation induced in the rear of the sphere, and the plume with positive streamwise velocity induced by the buoyancy effect. It influences the heat and mass transfer by a decrease of the local Nusselt number in the region of the rear stagnation point. This is consistent with the observations of Bhattacharyya and Singh, (2008) and Kotouč et al., (2008) (figure omitted). The buoyancy induced plume reduces the transport of cold fluid from the shear region towards the rear stagnation point, and with this the local temperature gradient and therefore the local Nusselt number decrease. But this modification is relatively small at the Richardson numbers investigated here and no significant difference is observed on the mean transfer coefficient (figure omitted). This is consistent with the simulations of Bhattacharyya and Singh, (2008) for which only small variations of the Nusselt number are observed in the range $0.0<Ri_{T}<0.1$ and $1<Re<200$.

We consider now the influence of buoyancy on the evolution of the velocity and temperature field. We first consider the evolution on the centerline (fig 15) and observe that buoyancy not only decreases the size of the recirculation zone but also the magnitude of the velocity in this zone as a consequence of the upward flow induced by density variations. The influence of buoyancy on the velocity deficit on the centerline seems to be confined in this region as little influence of $Ri_{T}$ is to be observed at $\tilde{r}_{\bot}=0$ for larger streamwise positions. The temperature is more affected on the centerline and it increases with $Ri_{T}$. This can be attributed to the modification of the recirculation: as this zone is more confined and less intense, the fluid transported to the back of the sphere is hotter than with $Ri_{T}=0$. This tends to decrease the cooling to be observed in this region of the flow, leading to larger temperatures. If we now consider the evolution of the half width, we see that the influence of buoyancy is more pronounced as represented on figure 15, where it can be seen that buoyancy tends to decrease the width of both thermal and velocity wakes, meaning that the wake is more confined. A possible interpretation of this is that buoyancy tends to favor the alignment of the wake in the vertical direction and through this the alignment with the sphere’s trajectory leading to a thinner wake.

5   Saturation profiles in the wake

We now examine the evolution of the saturation in the wake of the sphere for physical conditions that are the most relevant in a meteorological context. For this reason the case with $Ga=300$ is chosen, which features a chaotic behavior, with no contribution from buoyancy ($Ri_{T}=0$). Since the scalar fields seem to be little affected by the particle mobility at high density ratios, we present the case of a fixed sphere, which also ensures consistency with former work (Wang and Kubicek,, 2013; Cheng et al.,, 2014). We will investigate the following set of temperatures of the sphere:

[TABLE]

and for the temperatures of the incoming fluid:

[TABLE]

with the condition that the particle will always be warmer than the fluid. In the current section we will briefly describe the structure of the saturation field with respect to ice defined as the ratio between the local partial pressure $e$ and the saturation water pressure with respect to liquid water $e_{sat,w}$:

[TABLE]

The aim in the analysis is twofold. First we will test whether the local saturation can increase up to values that might trigger further ice nucleation (Hoose and Möhler,, 2012). Second, a methodological issue will be addressed with a discussion on the necessity to separate the transport equations of heat and water vapor. The difference between Schmidt and Prandtl number is in the present case relatively small, and the boundary conditions in non-dimensional form are the same for the temperature and the vapor concentration, which could suggest to use the approximation $\tilde{n}_{v}=\tilde{T}$. In the current study, our goal is to have access to the detail of the saturation field in the wake of the particle in order to address these two questions, namely what is the range of saturation reached? And does the approximation $\tilde{n}_{v}=\tilde{T}$ have a negligible influence?As previously explained in section 2 we set the saturation to be equal to unity with respect to the liquid phase at infinity and equal to unity with respect to ice at the particle surface such as to model the presence of micro-droplets at the inflow. It could also represent a wet growth regime at $T_{p}=0\degree C$ since we would have then $e_{sat,w}=e_{sat,i}$. The saturation field with respect to water will depend on the set of temperatures taken at the boundaries and figure 16 gives a visual impression of the saturation $S_{w}$ for the smallest and largest temperature differences between the particle and the updraft. It shows that the largest temperature difference (i.e. a very warm graupel) will induce the largest saturation, which can be larger than 1.1 in a large portion of the wake. For the lowest temperature difference (i.e. configurations with the coldest hailstone) the level of saturation reached is much lower. This last configuration does therefore not appear to be favorable for secondary ice nucleation.

From the definitions 10, 24 and 46 one can show that the incremental change of $S_{w}$ is linked to the incremental change of $\tilde{n}_{v}$ and $\tilde{T}$ according to

[TABLE]

This shows that very close to the particle, the variations of $\tilde{n}_{v}$ and $\tilde{T}$ will act in different ways, since the coefficient in front of $\mathrm{d}\tilde{n}_{v}$ is always positive and the coefficient in front of $\mathrm{d}\tilde{T}$ always negative, due to the evolution of $e_{sat,w}$ in the range of temperature considered here. It explains the small saturation observed in the vicinity of the particle for very small $\tilde{z}_{\parallel}$ on figure 16. Indeed, in this region we will have a decrease of both $\tilde{n}_{v}$ and $\tilde{T}$. It leads on the one hand to a decrease of the local partial pressure, but on the other hand to the decrease of temperature $T$ and inducing a decrease of the threshold $e_{sat,w}(T)$ at which we attain supersaturation.

If we neglect the differences in the diffusivity of mass and heat (i.e. $\tilde{n}_{v}\approx\tilde{T}$), then equation 24 can be written in the form $S_{w}=f(\tilde{T},T_{p},T_{\infty})$.

We can introduce with this the surrogate saturation $S_{w}^{T}$ that would be obtained by computing $n_{v}$ from the non-dimensional temperature, hence providing information on the expected saturation. It is computed as a function of $\tilde{T}$ only and the boundary conditions according to:

[TABLE]

The evolution of this "ideal" saturation $S_{w}^{T}$ is represented as a function of the temperature and for different sets of boundary conditions on figure 17.

Equation 25 shows not only that partial pressure and temperature have opposite contributions but also the necessity to account for both variations. We now investigate to which extent the approximation $\tilde{n}_{v}=\tilde{T}$ would affect the results obtained here. For this we introduce the relative error, $\varepsilon(\textbf{x},t)=\left|S_{w}(\mathbf{x},t)/S_{w}^{T}(\mathbf{x},t)-1\right|$ to quantify the differences reached. Table 1 summarizes the maximal relative error between $S_{w}^{T}$ and $S_{w}$ obtained for different sets of boundary conditions. It shows that the error can attain up to roughly 7 percent for the largest temperature difference. Figure 18 depicts the evolution in a $(\tilde{r}_{\bot},\tilde{z}_{\parallel})$ plane of the relative error obtained on an instantaneous field. It shows that the largest error is reached in the mixing region and in the recirculation region next to the sphere, indicating that the local temperature might have an influence on this error.

Let us take a closer look at the typical saturation $S_{w}$ that is encountered for a given temperature. For this we proceed as follows: we first divide the temperature domain into bins of temperature $T_{c}$ and we denote by $\mathscr{S}_{w}(T_{c},t)$ the ensemble of saturation values obtained on the subvolume $\Omega(T_{c},t)$, itself defined by the condition $T(\mathbf{x}\in\Omega(T_{c},t),t)=T_{c}$. We then average, for each temperature $T_{c}$, over the samples collected in time by the ensemble of saturation $\mathscr{S}_{w}(T_{c},t)$. We denote the corresponding mean by $\left<\mathscr{S}_{w}\right>(T_{c})$. The evolution of this conditional-averaged saturation as a function of the conditional temperature $T_{c}$ is represented on figure 17, where the standard deviation is also indicated. It confirms the trend given by the evolution of $S_{w}^{T}$, namely that for the smallest values of $\tilde{T}$ (i.e. far away from the sphere) $S_{w}$ increases with $\tilde{T}$ until a non-dimensional temperature $\tilde{T}_{c}$ ranging between 0.2 and 0.4, after which the saturation $S_{w}$ decreases with $\tilde{T}$. It also explains the increase of the saturation observed on figure 16 in the vicinity of the particle, as $\tilde{T}$ will decrease when going from the particle towards the fluid (negative components of $\mathrm{d}\tilde{T}$). The evolution of $S_{w}$ emphasizes also the influence of the boundary conditions since it appears that a large temperature difference tends to increase $S_{w}$ in the regions $\tilde{T}<\tilde{T}_{c}$, leading to larger increase of the saturation for a given temperature difference. This last point explains the difference observed for two sets of boundary condition represented on figure 16. The case $(T_{p}=0\degree C,T_{\infty}=-15\degree C)$ appears indeed to be more favorable to the development of supersaturated regions than the configuration $(T_{p}=-10\degree C,T_{\infty}=-15\degree C)$. The comparison between the evolutions of $S_{w}^{T}$ and $\mathscr{S}_{w}$ shows that the "ideal" saturation $S_{w}^{T}$ provides a good approximation of the effective saturation, with the tendency to underestimate it.

Let us now finish with a brief discussion on the evolution of the saturation as a function of the temperature represented on fig. 17(b). It indeed shows that transport of both heat and mass increases the saturation above the level of saturation equal to unity with respect to the liquid phase, which was already known to be a value at which a marked increase in ice-nucleating activity is achieved (Beard,, 1992). This indicates that, as proposed in the literature (Pruppacher and Klett,, 2010; Dye and Hobbs,, 1968; Rosinski and Morgan,, 1991), a falling graupel could be at the origin of a local increase of the saturation that might trigger further secondary ice nucleation. A deeper analysis would be required to quantify more precisely the potential of secondary ice nucleation and will be the object of future analysis. At this state of the work several conclusions can be formulated:

• •

The presence of the graupel can induce an increase of the saturation compared to the ambient. This may have several implications: it could "trigger" the formation of secondary ice nuclei by increasing $S_{w}$ up to the threshold of ice nucleating particle (INP) activation (which depends on the type of INP). If this threshold is already reached in the ambient then the increase of $S_{w}$ would also increase the activation rate. Furthermore, it could finally favor the growth of an ice crystal that was already formed.

• •

The difference between the saturations $S_{w}(\mathbf{x},t)$ and $S_{w}^{T}(\mathbf{x},t)$ is more important in the zones of high saturations, with error up to 7 percent and the tendency to underestimate the saturation if only one transport equation is considered. This means that both INP activation and activation rate might be underestimated, as well as the growth rate of ice crystals or micro-droplets already present in the flow.

• •

It is difficult to generally state whether this underestimation could lead or not to substantial error as it requires to know the typical regions that could be occupied by the micro-particles (ice nucleating particles, micro-droplets or ice crystals). To our knowledge, the concentration of micro-particles in the wake of a sphere has been only studied in the work of Homann and Bec, (2015). A key parameter is the Stokes number of the micro-particles $St_{mp}$ which is defined as the ratio between the relaxation time scale of the micro-particles $\tau_{mp}=\rho_{mp}D_{mp}^{2}/(18\nu\rho_{\infty})$ (with $\rho_{mp}$ and $D_{mp}$ respectively the density and diameter of the micro-particles) and the time-scale of the flow $\tau_{f}=D/\langle\lvert u_{p,z}\rvert\rangle_{t}$. Based on this the Stokes number will be a function of the properties of the micro-particle and the graupel according to:

[TABLE]

A large range of Stokes numbers can be reached depending on the type of micro-particles (featuring a small or very small diameter ratio $D_{mp}/D$) and the type of graupel (featuring a large or very large Galileo number). The largest Stokes will be reached by micro-droplets since they feature large diameter ratios (ranging from 0.001 to 0.5 depending on the type of graupel). The density ratio is expected to be of the order of $\mathcal{O}(10^{3})$ for micro-droplets, and for aerosols similar to the ones investigated by Kanji et al., (2011). The size ratio can vary between $\mathcal{O}(10^{-6})$ and $\mathcal{O}(10^{-3})$, leading either to very small Stokes number ($St_{mp}=\mathcal{O}(10^{-6})$) or intermediate Stokes numbers ($St_{mp}\approx 0.1$). For the smallest Stokes number, $\tau_{mp}$ is so small that the micro-particles react instantaneously to the fluid and can sample the entire wake. But for larger Stokes number and at the Reynolds numbers of interest Homann and Bec, (2015) have shown the formation of a cylinder of high concentration around the sphere, with a low-concentration shadow behind the sphere, followed by a zone of high concentration for larger streamwise distances. The zone with the largest saturation and largest error due to the approximation $\tilde{n}_{v}=\tilde{T}$ in the shear region coincides with the shadow region, such that the largest error would not have a significant influence at this Stokes number. For the largest temperature difference the error features blobs with errors of the order of one percent for larger streamwise distances in the high concentration cone (fig. 18). It is difficult to estimate whether or not this error could have an influence on the activation fraction as the corresponding saturation is close to the onset saturation (Hoose and Möhler,, 2012). Indications on the history of the saturation seen by micro-particles would be necessary to properly determine the importance of the error. Simulations accounting for the transport of aerosols should therefore be considered to clearly state whether this underestimation could lead to substantial error.

6   Conclusion

Direct numerical simulations have been conducted to study the case of a falling sphere made out of ice in moist air with the aim of exploring the saturation field in the wake of this sphere and test whether it might be large enough to be at the origin of secondary ice nucleation. For this we used a spectral/spectral-element solver to simulate the Navier-Stokes equations coupled to two transport equations to account for variation of heat and vapor content in the flow, under the Boussinesq assumption. The sphere is assumed to be at constant temperature, warmer than the surroundings, and at saturation of unity with respect to ice. The surrounding fluid is assumed to have a constant velocity with incoming constant temperature and constant saturation of unity with respect to liquid water.

The goal of the present study is to address several methodological questions to build a numerical framework that can be used in the future to investigate the ability of a graupel, for a given set of conditions, to trigger secondary ice nucleation.

The first question concerns the influence of the settling regime on the structure of the scalar fields, in the absence of buoyancy effects. The analysis focused first on the temperature field. Since heat and mass in this system have similar diffusivities the conclusions hold also for the vapor partial pressure. As expected it shows that the regime plays an important role on the structure of the temperature field, and that a configuration with a chaotic wake should be chosen for the application aimed here. The recirculation region appears to play a significant role on the temperature field as it drives the transport of cold fluid from the wake toward the back of the sphere and conversely. The influence of the mobility is considered here either by varying the density ratio between the particle and the upcoming fluid or by fixing the sphere, and it appears necessary to represent the graupel by an object featuring large density or by a fixed object.

The relatively small temperature differences and large density ratio together imply in the current context to have small Richardson numbers (less than $0.001$). Simulations accounting for buoyancy effects, for larger but still small Richardson number (i.e. up to $0.1$) have been performed in order to test whether such small values might affect the results. It appears that light buoyancy effects at $Ri_{T}=0.1$ are sufficient to stabilize the wake and trajectory and the sphere by pushing the thresholds of appearance of each regime towards larger Galileo numbers, with qualitative accordance with the observation from the literature for fixed spheres (Kotouč et al.,, 2008, 2009). Simulations involving $Ri_{T}=0.05$ and $Ri_{T}=0.1$ confirmed that the recirculation zone is also affected by the buoyancy with a decrease of the extent of this region, at much smaller Richardson numbers than that of Bhattacharyya and Singh, (2008) and in non-axisymmetric and unsteady chaotic regimes. The actual Richardson number of $Ri_{T}=10^{-4}$ appears to barely affect the flow, meaning that buoyancy effects do not need to be accounted for in the meteorological context.

Finally, the saturation in the wake of the sphere is explored for different sets of boundary conditions and it shows that substantial supersaturation can indeed be reached in the wake of the sphere, and that large temperature differences between the particle and the surrounding are more favorable to the development of this zone. Here as well the recirculation region has an important impact on the level of the saturation. This is, to our knowledge, the first time that this increase of saturation has been observed for chaotic regimes with the aid of direct numerical simulation. The necessity to separate both heat and mass transport equations is discussed, and small differences in the saturation field are obtained, with the tendency to underestimate the saturation if mass diffusivity is approximated with the heat diffusivity. No general recommendation can be formulated on the necessity to separate both transport equations since the error introduced would depend on the explored system. Configurations with small ice nucleating particles featuring small Stokes numbers appear to be the most sensitive ones since the micro-particles might sample the regions of the flow with the largest error (up to 7$\%$) at saturation close to the activation threshold. Additional simulations involving the transport of micro-particles would therefore be required to discuss the necessity to simulate separately the transport of heat and water vapor. The main prospects of this work concern the more systematic analysis of the saturation field, for different sets of boundary conditions in order to discuss in detail the expectation to observe secondary ice nucleation in the wake of a graupel.

Acknowledgements

The simulations were partially performed at SCC Karlsruhe. The computer resources, technical expertise and assistance provided by this center are thankfully acknowledged. We acknowledge support by Deutsche Forschungsgemeinschaft and open access publishing fund of Karlsruhe Institute of Technology.

Appendix A Governing equations

The dimensional form of the Navier-Stokes equations under the Boussinesq approximation and the advection-diffusion equations for heat and mass, both formulated in an inertial frame of reference with independent variables $(x^{*},y^{*},z^{*},t^{*})$, read

[TABLE]

where $\mathbf{u}$ is the fluid velocity, $P$ the pressure and $T$, $n_{v}$ the two scalars. The Boussinesq approximation states that the variation in density due to variations in temperature and vapor concentration affects only the buoyancy term $\rho/\rho_{\infty}\mathbf{g}$.

We subtract the contribution of the hydrostatic pressure from the pressure gradient term by introducing a modified pressure field $p$ such that $\mathbf{\nabla}^{*}P=\mathbf{\nabla}^{*}p+\rho_{\infty}\mathbf{g}$. With this the momentum equation becomes

[TABLE]

where the buoyancy term $\mathbf{B}$ is defined as $\mathbf{B}=\left((\rho-\rho_{\infty})/\rho_{\infty}\right)\mathbf{g}$.

In order to evaluate our governing equations on a mesh which is translating with the particle motion, we introduce the following independent variables of the moving mesh $(x,y,z,t)=(x^{*}-x^{*}_{p}(t^{*}),y^{*}-y^{*}_{p}(t^{*}),y^{*}-z^{*}_{p}(t^{*}),t^{*})$. The origin is located at the center of the sphere and the axes are aligned with the frame of reference. Each field $A(\mathbf{x},t)$ is then given in the mesh by

[TABLE]

Please note, that the components of vector fields remain expressed in the inertial frame, and hence merely the indexing has changed. The transformation is constant in space, and thus

[TABLE]

The temporal derivative of $A$ on the moving mesh can be obtained by applying the multivariable chain rule to (33), which results in

[TABLE]

where $\mathbf{u}_{p}\equiv\frac{\partial\mathbf{x}^{*}_{p}}{\partial t^{*}}$ is the particle velocity in the inertial frame. The substantial derivative of quantity $A$ can therefore be expressed on the translating mesh by

[TABLE]

Thus, the formulation of the Navier-Stokes equations and transport equations of the scalar quantities, on the particle-attached grid read

[TABLE]

The system of equations can then be reformulated in non-dimensional form, after normalizing the spatial coordinates according to $\tilde{x}_{i}=x_{i}/D$, the time $\tilde{t}=tD/u_{g}$, velocity as $\tilde{\mathbf{u}}=\mathbf{u}/u_{g}$, pressure as $\tilde{p}=p/(\rho_{\infty}u_{g}^{2})$, and temperature and concentration according to eq. 5 and 6, leading to the set of equations 1-4.

We decompose the buoyancy term into two contributions with help of a first-order Taylor series approximation around $\rho_{\infty}$. The first contribution, $\mathbf{B}_{T}$, accounts for the modification of density due to variations in temperature and the second one, $\mathbf{B}_{n_{v}}$, due to variations in vapor concentration. Their definitions are

[TABLE]

We assume that the fluid behaves as a perfect gas, implying that the partial derivatives of $\rho$ with respect to $T$ and $n_{v}$ yield

[TABLE]

where $N_{A}$ is the Avogadro constant, $M_{w}$ the molar mass of water and $M_{d}$ the mixture molar mass of dry air. Both buoyancy contributions can therefore be written in non-dimensional form as $\tilde{\mathbf{B}}_{T}=Ri_{T}\tilde{T}\mathbf{k}$ and $\tilde{\mathbf{B}}_{n_{v}}=Ri_{n_{v}}\tilde{n}_{v}\mathbf{k}$, with the thermal Richardson number defined as in equation 7 and the vapor concentration Richardson number according to

[TABLE]

where the term $(\rho_{p}/\rho_{\infty}-1)^{-1}=D\lvert\mathbf{g}\rvert/u_{g}^{2}$ comes from the non-dimensionalization.

Appendix B Computation of the saturation field

The saturation $S_{j}$ with respect to the phase $j$ is defined as the ratio between the vapor partial pressure $e$ and the saturation vapor pressure with respect to this phase $e_{sat,j}$. We used the correlations of Murphy and Koop, (2005) for the computation of the saturation vapor pressure with respect to ice and liquid water, respectively $e_{sat,i}$ and $e_{sat,w}$, namely

[TABLE]

with temperatures expressed in degrees Kelvin and pressures expressed in Pascal. Figure 19 depicts the evolution of the saturation vapor pressure as a function of the dimensional temperature for both liquid and solid phases.