Eulerian–Lagrangian CFD-microphysics modeling of a near-field contrail from a realistic turbofan

Aircraft contrails contribute to climate change through global radiative forcing. As part of the general effort aimed at developing reliable decision-making tools, this paper demonstrates the feasibility of implementing a Lagrangian ice microphysical module in a commercial CFD code to characterize the early development of near-field contrails. While engine jets are highly parameterized in most existing models in a way that neglects the nozzle exit-related aspects, our model accounts for the geometric complexity of modern turbofan exhausts. The modeling strategy is based on three-dimensional URANS simulations of an aircraft nozzle exit involving a bypass and a core jet (Eulerian gas phase). Solid soot and ice particles (dispersed phase) are individually tracked using a Lagrangian approach. The implemented microphysical module accounts for the main process of water-vapor condensation on pre-activated soot particles known as heterogeneous condensation. The predictive capabilities of the proposed model are demonstrated through a comprehensive validation set based on the jet-flow dynamics and turbulence statistics in the case of compressible, turbulent coaxial jets. Simulations of contrail formation from a realistic nozzle-exit geometry of the CFM56-3 engine (short-cowl nozzle delivering a dual stream jet with a bypass rate of 5.3) were also carried out in typical cruise flight conditions ensuring contrail formation. The model provides reliable predictions in terms of the plume dilution and ice-particle properties as compared to available in-flight and numerical data. Such a model can then be used to characterize the impact of nozzle-exit parameters on the optical and microphysical properties of near-field contrails.


Introduction
The continuous growth in air passengers by about 5% per year has increased scientific concerns about the climate impact of commercial aviation. 1 Cruise flights contribute to global warming through soot and watervapor emissions from aircraft engines that are responsible for visible white condensation trails (contrails) behind aircraft. 2,3 Persistent contrails composed mainly of supersaturated ice crystals survive in the atmosphere for several hours and even evolve to form large cirrus formations referred to as aircraft-induced clouds (AIC). 4 The spreading and persistence of this artificial cloudiness in the atmosphere determine their impact on the earth's radiative balance. 5 Global assessments 6 show that the impact of linear contrails is about 3-9Á10 23 WÁm 22 , while the impact of AIC is 10 times higher, 6 which corresponds to three times the impact of aviation CO 2 alone. 7 For this reason, the contrail formation process has been a topic of increasing interest in the literature. Owing to both complexity and onerous costs of in-flight measurements, the development of innovative mitigation strategies has to be based on reliable decision-making tools. As such, several modeling approaches with different levels of microphysics/ dynamics representations have been used in literature to investigate contrails. The following is a short review of main approaches and findings.
Modeling contrail life cycle involves several scales ranging from local scales in the engine jet phase (t~0-10 s after the exhaust emission) to the global scale in the dispersion phase (few minutes to hours after the exhaust emission; see Figure 1). From the thermodynamic point of view, the revised Schmidt-Appleman criterion 8 helps to determine atmospheric regions of contrail occurrence based on water liquid saturation conditions. This criterion was used to initiate contrail formation in global models, such as in the Contrail Cirrus Prediction Tool (CoCiP), 9 and made it possible to assess the global environmental impact of contrails. To allow coverage of the full evolution cycle of a large population of contrails, the physics modeling was relatively simplified, such as for initialization of the species concentration with a Gaussian profile. Even though the results of regional air-traffic simulations 9,10 were quite consistent with in-flight measurements and satellite observations, the microphysical processes of ice production remain misunderstood. A more precise characterization of ice particles was achieved using 0-D microphysical trajectory box models. 11,12 This class of models allows detailed chemical kinetics and microphysics in the engine jet phase. For instance, simulations of young contrails 11,13 highlighted the domination of soot-induced heterogeneous freezing due to the large amount of emitted soot condensation nuclei from kerosene-fueled engines. In this respect, the reduction of soot emissions with the use of alternative jet fuels can help alter the optical properties of contrails. 14 The box model results for plume dilution, however, showed some discrepancies compared to in-flight measurements due to (1) the absence of a detailed turbulence modeling 11,14 and (2) the simplified level of dynamics using averaged bulk exhaust (provided by semiempirical models). 11,14 Some authors remedied the shortcomings through offline-coupling strategies with an external fluid-flow solver. 15,16 In the last two decades, computational costeffectiveness has increased interest in 3D CFD-based models to refine the early plume dynamics and associated mixing process. Early plume properties are known to influence later contrail stages (optical depth, radiative forcing) as shown by the studies on the engine jet and vortex phase. 7,17 As such, they are usually used as inputs in mesoscale models. 18,19 Solid particles motion can be tracked in CFD using either an Eulerian 20 or a Lagrangian. [21][22][23] The former can capture the motion of higher particle concentrations 24 with either bulk parametrization 25,26 or bins 20 for microphysics modeling. The simplest representation of particle-size distribution considers a single bin (monodisperse distribution) as used by Khou et al. 27,28 who investigated the near-field contrails (engine jet phase) behind realistic commercial aircraft geometries using Reynolds-averaged Navier-Stokes (RANS) equations for turbulence modeling. Results revealed that wingtip vortices drained and cooled jets and highlighted the role of ambient relative humidity on ice-crystal properties. Using a single bin could, however, prevent accurate estimates of icecrystal concentration. 28 The most physically accurate treatment of ice particles in Eulerian terms considers size-resolved (binned) microphysics, 20 which avoids the problems of numerical diffusion inherent in bulk models 29,30 that could lead to slight discrepancies in the dispersion of ice and water-vapor fields. The computational cost, however, is excessive (see the detailed review 31 ).
In contrast, Lagrangian tracking allows for an explicit process-oriented modeling of ice microphysics, 30 which accounts for its popularity in contrail research for a while. 15,22,30 Most if not all Lagrangian studies are highly parameterized in a way that neglects nozzle exit-related aspects. For example, Paoli et al. 32 modeled free turbulent jets with classical expansion laws to investigate contrail formation for two-and four-engine configurations. Unterstrasser and Go¨rsch 33 also used Gaussian-like plume profiles to study the vortex phase and to characterize the impact of aircraft type on contrail evolution. These simplified representations of the nozzle exhaust jet do not allow for accurately taking into consideration the features of modern turbofan nozzles. As the flow structure and mixing in the early plume depend on engine exhaust parameters, the impact of high bypass rates and complex geometries (e.g. short and long cowls and chevrons) on the nearfield microphysics remain unclear. This particularly motivated the development of a CFD-microphysical model that integrates the nozzle-exit geometry and avoids using parameterizations to describe the exhaust plume.
Therefore, the main objective of this study was to implement a Lagrangian microphysical module in a CFD code to simulate ice-particle growth in near-field contrails starting from a realistic turbofan nozzle exit. To the authors' knowledge, no contrail study has integrated the actual nozzle-exit geometry in a Lagrangian approach. A high-resolution spatial unsteady RANS method was therefore used to demonstrate the feasibility of the present approach. The paper is organized as follows. Section 2 begins with a description of the modeling approach, followed by the mathematical equations of fluid dynamics, particle motion, and ice microphysics. Section 3 describes the model and presents how the microphysical module was coupled to Simcenter Star-CCM + CFD software. Section 4 discusses the validation test results for jet dynamics, plume dilution, and ice-particle properties.

Model description
Approach for modeling the formation of near-field contrails Figure 2 shows that the nozzle exit flow from a turbofan engine at cruise flight consists of two jets: a hot core flow and a surrounding cold bypass flow. Both jets are mixed with the freestream moist air. The bypass flow simply contains ambient moist air that has been slightly accelerated by the engine fan, while the core jet drives the fuel-combustion products composed of gaseous emissions laden with soot particles (dispersed or solid phase). Water vapor in the dilution plume condenses onto soot particles, which serve as condensation nuclei when activated with sulfuric acid (heterogeneous nucleation 11 ). Under saturation conditions, ice particles form and grow by taking up more condensable matter -that is, water vapor -from the plume and ambient air, leading to contrail formation.
In this study, the Eulerian-Lagrangian approach was used to model the formation of contrail ice particles in the near field of an aircraft engine. The continuous gas phase was modeled with the Eulerian approach, while the dispersed-phase of solid particles (soot and ice particles) were tracked with the Lagrangian approach. Lastly, a one-way microphysical coupling was considered to account for ice-particle growth due to watervapor condensation.

Gas-phase equations
The compressible flow of the gas phase (''g'' index) is considered to be a mixture of two miscible species: air (''a'' index) and water vapor (''w'' index). Assuming an ideal gas, the mixture density r g , temperatureT g , and pressure p g are linked with the equation of state: where r g is the specific gas constant. The overbar signdenotes the Reynolds (ensemble) averaging, whiler efers to the Favre (density-weighted) averaging defined as rF = rF and F =F + F 00 . The latter is considered to account for compressibility effects. As the growth of ice particles is a time-dependent process, the flow is assumed to be unsteady. Therefore, the unsteady Favre-averaged Navier-Stokes (FANS) equations expressing the conservation of the gas-phase global mass, momentum, and energy were solved. In Cartesian coordinates (i, j, k), these equations read as follow with some modeling assumptions. 34 Solving equations (2)-(4) provides the primary variables of the mean flow r g , the velocity field u g =ũ g ,ṽ g ,w g À Á , and the total energyẽ tot g , such as c p, g denotes the specific heat m g the dynamic viscosity, m t, g the eddy viscosity, Pr g the Prandtl number, and Pr t, g the turbulent Reynolds number of the gas phase. The latter physical properties are all functions ofT g , which was derived from the total energy using equation (5). The averaging of nonlinear terms on the right-hand side of equations (3) and (4) yields the Reynolds stress e t turb ij and the turbulent heat fluxq turb j, g , in addition to their laminar (or viscous) counterparts e t visc ij .andq visc j, g , given by equations (6)-(9), respectively, for a Newtonian fluid.
q visc j, g = c p, g m g Pr g ∂T g ∂x j ð8Þ Figure 2. 2D Scheme of exhaust flow from a turbofan engine with ice-particle formation.
To compute the mixture composition in the gas phase (F g =Ỹ a F a +Ỹ w F w ), the mass fractions of air and water vapor (Ỹ a andỸ w , respectively) were derived with equations (10) and (11) defined as follow: Equation (10) translates the mass conservation of water vapor where the source term v ice on the right-hand side is the rate of condensed matter representing the mass transfer due to ice-particle growth. This coupling term is presented in greater detail in the next section. The variables s g and s t, g denote the Schmidt number and the turbulent Schmidt number, respectively taken as being equal to 1 35 and 0.9. 27,28 Lastly, the eddy viscosity m t, g term introduced by the Boussinesq approximation (equation (7)) was modeled using the Realizable k À e model 36 to close the set of equations. This model is known to provide improved predictions in terms of the dissipation rate for round jets. 37 It was chosen after comparing three eddyviscosity models to the available experimental data in the case of compressible coaxial jets. Validation tests are presented below in section 5.1. The k À e Realizable formulation uses two transport equations: one for the turbulent kinetic energy k and one for the dissipation rate e, such as: where S ij = 1 2 ∂ũ i, g =∂x j + ∂ũ j, g =∂x i À Á represent the deformation rate of the mean flow. C m , s k , s e , C 1e , and C 2e are the model constants (see Refs. 34,38 ).

Dispersed-phase equations and microphysics model
Particles motion. At cruise conditions, the number density of physical particles is so high that a modeling simplification is required to save simulation time. For example, a typical value for emitted soot particles from an engine nozzle is about 10 11 #Á m 23 . 39 Therefore, computational particles are considered as spheres of radius r p , with each representing N physical particles. 40 Numerically, soot and ice particles are treated as passive tracers since their size and relaxation time are relatively small compared to the characteristic scales of the gas-phase flow (for ice particles t p1 0 À5 s with the largest r p of a few microns). 41 Thus, they have the same velocity as the carrier phase, but the drag and gravity forces are neglected since they have neither mass nor volume. The equation governing particle movement reads as follows: where x p is the particle position vector.
Ice-particle growth. In contrails, ice-particle growth is due to condensation of water vapor onto suitable nucleation sites on pre-activate particles (mostly soot particles). 8 The ratio between the local water vapor partial pressure (P w ) and the saturated vapor pressure with respect to liquid water P 0 w, liq defines the saturation ratio 5 noted S: such as the local pressure P w is calculated according to Dalton's law: where X w is the mole fraction of water vapor and P tot is the absolute total pressure. When the saturation ratio is above 1, heterogeneous ice nucleation takes place and spherical-shaped ice particles form and grow in the dilution plume. In this study, the microphysics of icecrystal growth was modeled, while the processes of soot activation 42 and homogeneous ice nucleation 43 were not taken into account. Therefore, the transition from continuum to molecular regime is supposed as being instantaneous. 44 The mass change of a single-particle (computational) growth due to condensation/evaporation effects can be expressed by a diffusion law as described by the Fukuta and Walter 44 model, which describes the change of both particle mass and size as follows: where m p denotes the mass of a particle p, r p the particle density, and A the Kelvin effect (defined in Appendix A). The factors C T and C r control the watervapor temperature and air density just above the particle surface. The detailed formulations of both factors can be found in the Appendix A. Therefore, we obtain the mass variation of each ice crystal and a fortiori the vapor mass variation. This leads to the coupling term v ice in equation (10), which can be derived as follows where n is the number density of particles contained in a volume cell V, such that n = N p =V. The total number of computational particles is N com , each of them representing N p physical particles. More details about the model are available in Fukuta and Walter. 44 The latter model was implemented in a microphysical module and coupled to the CFD commercial code in Simcenter STAR-CCM + v13.04. A detailed description of the numerical configuration is provided below.

Numerical methodology
The CFD solver is based on a cell-centered finite-volume approach for unstructured grids. 20 The SIMPLE algorithm 45 was used to solve the flow governing equations. A second-order upwind scheme was used for space discretization, while time integration was achieved with a second-order Euler scheme. For the model of ice-particle growth, a fourth-order Runge-Kutta scheme was used to discretize the equation (20) of particle size change. The discretization of the equation (20) is presented in Appendix B. The microphysical module was implemented as a user-defined function in the commercial CFD code in Simcenter STAR-CCM + v13.04. The script of the microphysical module (coded in C) can be downloaded from Ref. 46 The methodology for coupling with the CFD solver is detailed in Figure 3. At each time step, the primary flow variables of both gas-and dispersed-phase are derived from the CFD solver. The pressure and temperature at each grid point help to compute the thermodynamic variables (P w , P 0 w, liq , C T , and C r ), allowing derivation of the local saturation ratio. The residence time (t p ) of each particle is computed at each time step. Radius of new particles (t p = 0) is initialized with a monodisperse distribution using an initial radius of r p, 0 = 20nm, as per numerical studies. 22,27 The number density of soot particles was set to n = 10 11 #Ám 23 , which corresponds to an emission index EI soot , of 4 3 10 14 particles per kg of fuel according to in situ measurements. 39 The radius of old particles (t p . 0) was recovered from a particleradius table. For the following time steps, the evolution of particle radius depends on the local saturation-ratio values. Nonsaturated particles keep the same radius, while saturated particles grow. The rates of change for their mass ( _ m p ) and size ( _ r p ) were computed with equations (19) and (20). As such, the particle-radius table was updated and the mass-transfer term v ice derived with equation (21). The coupling term was injected into the CFD solver to correct the water-mass fraction in the following time step.

Validation tests and results
Two validation tests were performed to assess the reliability of the CFD-microphysics model presented above for studying ice-particle growth in the exhaust jet of an aircraft engine. In the first test, the model's gas-phase resolution (i.e. jet dynamics and turbulence) was assessed with experimental results from the canonical case of compressible coaxial jets. This validation test made it possible to choose the most suitable turbulence model. In the second test, the model's microphysical resolution in terms of dilution and ice-particle properties (size, distribution, optical depth) were assessed in the exhaust jet of a CFM56-3 engine at realistic cruise flight conditions.

Compressible coaxial turbulent jets
Computational domain and grid. Guitton et al. 47 experimentally characterized the flow dynamics and turbulence in the case of isothermal compressible coaxial jets of air. This experiment was used as a validation test case for our computational model. As shown in Figure  4, the nozzle had a length L noz = 0:158m and was composed of two jets: an inner primary jet, whose diameter was D p = 0:055m, and an outer secondary jet with a diameter of D s = 0:1m. The mean flow was assumed to be axisymmetric. Indeed, it has been observed experimentally that the vortex structures forming from the rollup of the shear-layers remain axisymmetric over a distance of a few diameters. 48 Hence, the computational domain was the two-dimensional rectangular box measuring 40D s 3 15D s based on the half-nozzle model, and the Reynolds-averaged Navier-Stokes (RANS) equations were solved in cylindrical coordinates O, r ! ,ỹ À Á . Two velocity-inlet boundaries were used for the engine primary-and secondary-jet flows, with corresponding values of U p = 170m:s À1 and U s = 120m:s À1 , respectively (Figure 4). A symmetry condition was applied to the centerline (bottom boundary), while the remaining boundaries were treated as ambient-pressure outlet. A no-slip boundary condition was used on nozzle walls (5 mm thick).
A grid convergence study was carried out using three grids: coarse (290,110 cells), medium (618,651 cells), and fine (1,227,567 cells), respectively labeled grid 1, 2, and 3. The GCI method 49 was used to compute the discretization errors. Results of grid errors based on the axial velocity at a point probe located in the mixing layer between the two exhaust jets (r = 0:5D p , y = 10D s ) were sufficiently low: GCI 12 coarse = 1:45% and GCI 23 fine = 0:94%. As such, the medium grid was selected for computations. This mesh configuration had a total of 35 prismatic layers next to all wall boundaries, which ensured a y + 41.
Jet dynamics and turbulence statistics. The mixing process and associated turbulence effects have a significant impact on the dilution and ice-particle microphysics in an exhaust jet. 32 As such, the model's gas-phase resolution was first validated by studying the performance of three eddy viscosity turbulence models: k À e Standard, 36 k À e Realizable, 50 and k À v SST. 51 This comparative study covered first-and second-order statistics, that is, mean axial velocity, turbulent kinetic energy (TKE), and the Reynolds stress. Figure 5 presents the cross-sectional profiles of the exhaust flow against experimental data. 31 Results of axial velocity profiles ( Figure 5(a)) confirm that the predictions of the three tested turbulence models are satisfactory overall. The Standard and SST versions slightly overestimated the axial velocity on the external part of the jet (i.e. r ø 0:5D s ). Small discrepancies can also be observed near the primary-jet centerline (i.e. y ø 7D s ) for the k À e Standard, while the Realizable version seems to offer the best fit for experimental measurements. 32 The Reynolds-stress profiles ( Figure 5(b)) also confirm the latter conclusion in which the Standard and SST models overestimated the experimental profiles for y ø 5D s .
To better assess the discrepancies observed near the jet axis in Figure 5(a), the performance study on turbulence models was carried out along the primary-jet centerline (see the results in Figure 6). The three models tested yielded similar performance in the potential core region (Figure 6(a)). In contrast, both velocity decay and TKE increase (Figure 6(b)) were best captured by the k À e Realizable model. Quantitatively, the k À e Realizable provided the lowest mean relative error with respect to experimental data on the mean axial velocity: 2.2% compared to 3.8% and 4.9% for k À v SST and k À e Standard, respectively. Based on this comparative study, the k À e Realizable was selected.

Contrail formation in the CFM56-3 exhaust jet
Computational domain and grid. Ice-particle formation was modeled in the exhaust plume of the CFM56-3 engine at realistic cruise flight conditions that ensure contrail formation. 27 The engine speed was set at U cruise = 252 mÁs 21 and ambient conditions were taken at 35, 000ft. The static temperature and pressure values were T amb = 219K and P amb = 23, 800 Pa, respectively. The mass fraction of water vapor in ambient air was set at Y amb = 6:08 Á 10 À5 to account for 60% relative humidity for comparison with numerical results. 27 As the mean flow was assumed to be axisymmetric, only the 90°sector based on the engine nozzle was considered for URANS calculations, as shown in Figure 7. The geometry of the CFM56-3 short-cowl nozzle is composed of a core and a bypass flow whose diameters are D c = 0:646m and D b = 1:586m. The nozzle configuration is representative of modern commercial turbofan engines, 52 in which a central plug is included to encourage early mixing between the core and bypass flows (hereafter the inner shear layer). The exhaust conditions used by Garnier et al. 53 were used in this study (see Table 1). According to Paoli et al. 32 physical soot particles in the range 10 6 -10 8 preserves the statistical reliability of the simulation results. As such, 10 7 physical particles were used herein. This ensures 40,000 computational particles, and a total number of emitted soot particles at n = 10 11 #Ám 23 .
As shown in Figure 7, the computational domain had a radius R = 5D b and a length L = 76D b to ensure ambient pressure at far-field boundaries. Three inlet sections were used upstream: a freestream pressure inlet (P 0 amb ,T 0 amb ), a bypass pressure inlet (P b , T b ), and a core pressure inlet (P c , T c ). A no-slip adiabatic condition was imposed on engine walls, while a symmetry condition was used on both lateral planes A turbulence intensity of 10% was used for the exhaust jets, as per Garnier et al., 53 while a value of 0.1% was applied on both freestream and far-field boundaries (P amb ,T amb ). Remaining boundaries were treated with a far-field condition.
The geometry and mesh of the computational domain were generated with ANSYS ICEM CFD software. 54 A grid convergence study was performed on the gaseous phase using three grids: coarse (4.35M cells), medium (9.8M cells), and fine (20.9M cells), respectively labeled grid 1, 2, and 3. The mean errors computed in the inner shear layer (i.e. between the core and the bypass jets) for the static temperature, watervapor mass fraction, and core axial velocity are below 1% (see Table 2). Furthermore, discretization errors estimated using the GCI method 49 at a probe point (r = 0:7D c , y = 10D b ) were sufficiently low with the medium mesh. As such, computations were only performed on the medium grid.
The medium grid is mainly composed of hexahedral elements, as shown in Figure 8. Three cell numbers n r , n u , and n z characterize the global mesh in cylindrical coordinates (r, u, z). Refinement blocks were used in the exhaust jet and shear layer, where gradients are expected to be relatively high. Table 3     Number of computational particles, N com (#) 4 10 4 Initial soot-particle diameter, r p, 0 (nm) Jet flow and plume dilution. Figure 9(a) and (b) present the stream-wise contours of the axial mean velocity and Mach number, respectively, for the CFM56-3 plume.
Since both dual-stream jets were at sonic conditions (M = 1), weak shock waves were observed near the exhaust sections and edges (Figure 9(b)); similar behavior was reported at nozzle-exit flows. 55,56 A flow separation on the central plug can be seen in the core flow due to the over-expansion shock. The velocity contours in Figure 9(a) show that the potential core length was about 5D b ; this distance is comparable to both experimental (5.4D b ) 52 and numerical (5D b ) 57 values. The dilution ratio defined by Schumann et al. 58 was used to characterize the evolution of flow variables (such as temperature and contrail diameter) along the Table 2. Grid study results: mean relative error e, standard deviation s, and GCI index 53 between grids 1 (coarse), 2 (medium), and 3 (fine) in the inner shear layer (r = 0:7D c ).
where N exit is the dilution factor at the exhaust section 58 andZ is the ratio of the amount of fluid from the jet region to that from the freestream. 53 The latter was calculated from whereZ y, plume is the passive scalar concentration in the exhaust plume,Z i, atm is the initial atmospheric concentration, andZ i, jet denotes the initial jet concentration. The concentration measurement was taken along the plume centerline. Figure 10 shows the evolution of the dilution ratio N Z as function of the plume age compared to the bulk mean dilution interpolated from in-flight measurements. 59 For instance, the interpolation was based on a large set of data from more than 70 plume encounters with the DLR Falcon research aircraft and the NASA ER-2 aircraft. 59 It should be noted that the bulk dilution ratio corresponded to the radial average of the local dilutions ratios within the plume at a given distance from the nozzle exit. This explains why the predicted dilution ratio along the core jet as seen in Figure  10 was relatively smaller than the bulk interpolation of the mixed primary and secondary jets (as also highlighted by Schumann et al. 59 ), especially near the engine exhaust. Even with that, the variations of the dilution ratio in the simulated plume are of the same order of magnitude as the bulk mean interpolation. The study of ice-crystal properties during contrail formation is presented in the next section.
Contrail formation and ice-particle properties. Figure 11(a) presents the spatial distribution of ice-particle radii in the plume exhaust. Up to a distance of 10D b downstream, the jet temperature was high enough to prevent ice coating of soot particles with an initial radius of r p, 0 = 20nm. The mixing with ambient air downstream cooled the plume, and the first ice crystals can be observed within a thin layer of the plume edge at y = 10D b . More ice particles formed and grew subsequent to the thickening of the external shear layer by transverse diffusion up to y50D b where water-vapor condensation reached the plume core.
The rate of condensation along the exhaust plume was fully determined by the local ice saturation ratio S, which, in turn, depends on the local temperatureT g and water mass fractionỸ w . 32 Figure 11(b) and (c) compare two transverse cross-sectional distributions at   y = 30D b and 72D b from the core exhaust section, respectively. Both engine core and bypass sections are represented by black circles (solid line). At y = 30D b , the particles were spatially distributed over a circular area of diameter 2D b . Ice condensation was only observed on the external shear layer (r ø 0:5D b ) while soot particles transported along the plume core remain under-saturated. The transverse dispersion of ice crystals enhanced by mixing with ambient air downstream resulted in a radial expansion of the plume up to a diameter 3D b at y = 72D b . At this section, the center of the plume was fully saturated and, due to the relatively high water content in the core jet, ice particles seemed to grow faster along the plume core (r4D b ), up to a radius about 50r p, 0 , as compared to the plume edge. At the outer edge of the plume; ice particle radii reached a constant value around 500 nm for y=D b . 50. Indeed, the growth of ice-particles (i.e. condensation rate) is mainly governed by the local value of the saturation ratio S, which is below 1 in the outer edge of the plume. As a consequence, the ice-particles do not grow in the outer edge and reach a constant size of 500 nm.
The evolution of the saturation ratio along the exhaust jet shows that the plume was fully saturated, that is, S = 1, at a distance about 55D b behind the nozzle exit (see Figure 12). This is consistent with the numerical results of Khou 60 obtained behind a full aircraft geometry with an Eulerian approach for particle tracking. Since the growth rate of ice crystals is directly related to the local saturation-ratio value, the mean particle radius follows the same trend as the bulk mean saturation ratio. It should be noted that that the mean particle radius \ r p y ð Þ . corresponds to the arithmetic average of particle radii over the subdomain of length h centered on the position y, defined as 32 where N is the particle number within the considered subdomain and r p, i is the radius of a given particle i.
Overall, the results show that our model captures well the global dynamic of ice-particle growth compared to the numerical prediction by Khou. 60 The mean radius was initially constant near the engine exhaust and started growing around y = 10D b . At 0:5s of the plume age (y = 72D b ), the mean particle radius was in the good range of 500-1000 nm, as confirmed by several numerical studies. 15,16,61 To better assess the particle-size distribution along the exhaust plume, the probability density functions (PDF) were investigated at four distances: 20D b , 30D b , 50D b , and 72D b behind the engine exhaust (see Figure  13). Results show that the PDF distributions shifted toward large radii following the growth of ice crystals along the exhaust plume. Moreover, Gaussian distributions were obtained from y = 50D b indicating polydispersed distributions of ice-particle radii, as observed in other numerical studies. 22,32 Lastly, the visibility in the predicted contrail was assessed with the optical depth t. This parameter measures the attenuation of radiation passing through clouds and is often reported at a wavelength in the visible part of the radiation spectrum. 7 The optical depth t was evaluated as formulated by Naiman et al. 62 where N com, g denotes the number of computational particles in the grid cell, N g the total number of hysical particles in the grid cell, and DV its volume. Equation (26) approximates the Mie extinction efficiency Q ext at a given wavelength l w , such as the refractive index for ice of m r = 1:31 and the wavelength of visible light is considered as l w = 0:55mm. Figure 12. Evolutions of the mean particle radius and the saturation ratio along the jet plume. Figure 13. Probability density functions (PDF) of particle radius at four axial distances: 20D b , 30D b , 50D b , and 72D b .
The visibility threshold of jet contrails is defined by the criterion t ø 0:5. 61 Figure 14 shows that the predicted contrail from the CFM56-3 engine was visible at a distance about 26m from the engine exhaust. The latter value is good agreement with the in situ observations reporting visible contrails at a distance in the range 10-30 m. 59,63 The comprehensive analysis performed shows that this model's predictions are consistent with available numerical data. Such model can be used in future investigations to characterize the impact of complex nozzle geometries -such as chevron nozzles -on the optical and microphysical properties of near-field contrails. Research perspectives to be explored next involve the interactions of soot and sulfur species in the early plume. This requires the implementation of a gas-phase chemistry with the microphysics of volatile particles.

Conclusion
This paper demonstrated the feasibility of Eulerian-Lagrangian CFD coupling with ice microphysics to predict the near-field contrails starting from the nozzleexit geometry of the CFM56-3 engine. The microphysical module implemented accounts for the main process of ice-particle growth based on heterogeneous nucleation. Two validation tests were performed to assess the reliability of the proposed CFD-microphysics model. In the first validation, the canonical case of compressible coaxial jets made it possible to choose the Realizable k À e turbulence model based on its good performance in terms of the jet dynamics and turbulence properties. In the second test case, the model's microphysical resolution was assessed in the exhaust jet from the CFM56-3 engine at realistic cruise flight conditions. The predicted plume dilution showed good agreement with available experimental data. Qualitatively, the spatial distribution of ice particles revealed that the first ice crystals were formed at the plume edge at a distance of about 10D b from the engine exhaust. Further downstream, ice particles formed and grew subsequent to the thickening of the external shear layer due to transverse diffusion. Results of the saturation ratio show that the plume core axis was saturated at a distance of about 55D b . Polydispersed distributions of particle radii were obtained in the simulated contrail with the largest ice particle reaching a radius about 1 m m in the plume core. As the jet evolved downstream, the particle-size distributions shifted toward larger particles. Lastly, the optical-depth analysis shows that the predicted contrail was visible around 26 m behind the engine exhaust, which is consistent with contrail observations.
The model presented herein is intended to serve as the foundation for more advanced numerical simulations to investigate the impact of complex nozzle-exit geometries, such as chevron nozzles. The implementation of gas-phase chemistry and volatile-particle microphysics are planned to extend the model capabilities.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The operation of this supercomputer is funded by the Canada Foundation for Innovation (CFI), the Ministe`re de l'É conomie, de la science et de l'innovation du Que´bec (MESI), 531, and the Fonds de recherche du Que´bec-Nature et technologies (FRQ-NT). The authors would also like to acknowledge the financial support of Safran Aircraft Engine for this work. This work is part of the Safran Industrial Research Chair on the Development of Sustainable Aero-Propulsion Systems at É TS.

ORCID iDs
The saturation pressure of ice (Pa) was calculated from: