Xue Chen

Xun Wang

Paul G Chen

Qiusheng Liu

Jens Hauslage

Ruth Hemmersbach

Valentina Shevtsova

Determination of Diffusion Coefficient in Droplet Evaporation Experiment Using Response Surface Method

Keywords: Droplet evaporation, Diffusion coefficient, Response surface, Marangoni flow, Microgravity

Evaporation of a liquid droplet resting on a heated substrate is a complex free-surface advection-diffusion problem, in which the main driving force of the evaporation is the vapor concentration gradient across the droplet surface. Given the uncertainty associated with the diffusion coefficient of the vapor in the atmosphere during space evaporation experiments due to the environmental conditions, a simple and accurate determination of its value is of paramount importance for a better understanding of the evaporation process. Here we present a novel approach combining numerical simulations and experimental results to address this issue. Specifically, we construct a continuous function of output using a Kriging-based response surface method, which allows to use the numerical results as a black-box with a limited number of inputs and outputs. Relevant values of the diffusion coefficient can then be determined by solving an inverse problem which is based on accessible experimental data and the proposed response surface. In addition, on the basis of our numerical simulation results, we revisit a widely used formula for the prediction of the evaporation rate in the literature and propose a refined expression for the droplets evaporating on a heated substrate.

Introduction

Evaporation of a liquid droplet deposited on a solid substrate has been studied intensively over the past decades due to its technological importance in a wide variety of applications, such as micro-cooling devices, ink-jet printing and biochemical assays [START_REF] Sartre | Effect of interfacial phenomena on evaporative heat transfer in micro heat pipes[END_REF][START_REF] Park | Flow and heat transfer characteristics of the evaporating extended meniscus in a micro-capillary channel[END_REF][START_REF] Bonn | Wetting and spreading[END_REF][START_REF] Erbil | Evaporation of pure liquid sessile and spherical suspended drops: a review[END_REF][START_REF] Chen | Spreading of annular droplets on a horizontal fiber[END_REF]. The evaporation of a sessile droplet is a complex free-surface advection-diffusion problem where heat and mass transfer processes occur simultaneously. Though the flow field within an evaporating droplet on a heated substrate is largely determined by surface-tension-induced (Marangoni) flow, the main driving force of the evaporation is the vapor concentration gradient taken between the droplet surface and its ambient gas phase, and by definition, the rate at which the liquid droplet evaporates into the ambient gas scales linearly with the (molecular) diffusion coefficient of the vapor. Therefore, an accurate determination of the diffusion coefficient under various environmental conditions is of great importance to better understand the evaporation process. [START_REF] Maxwell | The Scientific Papers of James Clerk Maxwell[END_REF] was the first to propose a diffusion coefficient model to describe the evaporation rate in a uniform and infinite gaseous medium. Since then, many researchers have followed the same approach to study the evaporation of a spherical drop. [START_REF] Morse | Diffusion and Supersaturation in Gelatine[END_REF] and [START_REF] Langmuir | The evaporation of small spheres[END_REF] performed the experiments of suspended iodine drops and obtained a linear scaling of the evaporation rate with the droplet radius. [START_REF] Picknett | The evaporation of sessile or pendant drops in still air[END_REF] analyzed the problem of a droplet resting on a solid substrate by using an analogy between diffusive flux and electrostatic potential. The linear scaling expression included a diffusive term (diffusion coefficient, concentration distribution) and a term representing a coupling between the liquid and the solid phases (initial contact radius and a function of contact angle). However, the question of how the function is related to the contact angle and other parameters remains controversially debated [START_REF] Birdi | A study of the evaporation rates of small water drops placed on a solid surface[END_REF][START_REF] Rowan | Evaporation of microdroplets and the wetting of solid surfaces[END_REF][START_REF] Hu | Evaporation of a sessile droplet on a substrate[END_REF][START_REF] Popov | Evaporative deposition patterns: spatial dimensions of the deposit[END_REF]; [START_REF] Hu | Evaporation of a sessile droplet on a substrate[END_REF] put forward an intuitively simple expression through a contact-angle function, which is valid for angles between 0 and π/2.

In all of the above-mentioned studies, the thermal properties of the substrate, as well as the effect of evaporative cooling, have not been taken into account. Recent studies [START_REF] Dunn | The strong influence of substrate conductivity on droplet evaporation[END_REF][START_REF] Brutin | Evaporation of ethanol drops on a heated substrate under microgravity conditions[END_REF][START_REF] Bouchenna | Generalized formulation for evaporation rate and flow pattern prediction inside an evaporating pinned sessile drop[END_REF]Chen et al. 2017a, b) highlighted those effects on the structure of bulk fluid flow in evaporating drops. [START_REF] Ristenpart | Influence of substrate conductivity on circulation reversal in evaporating drops[END_REF] showed that the Marangoni flow is determined by both relative thermal conductivities and relative substrate thicknesses; however, the study was limited to the case of a non-heated substrate. Very recently, Chen et al. (2017a) obtained simple scaling laws for the evaporation rate that scales linearly with the drop radius but follows a power-law with the substrate temperature. Combining numerical simulations with response surface method, Chen et al. (2017b) further investigated the thermal effects of the substrate on the Marangoni flow and clearly identified three characteristic bulk flow structures. These recent studies called for an improved empirical relation of [START_REF] Hu | Evaporation of a sessile droplet on a substrate[END_REF] for the prediction of the evaporation rate.

Besides, it is not feasible in space evaporation experiments to insert the instruments allowing a precise measurement of the diffusion coefficient. Given the uncertainty associated with the diffusion coefficient in a specific experiment in space, the question then arises as to what values should be used. The failure of some values used in the numerical and theoretical predictions, when compared with the experimental data, indicates that the diffusion coefficient is largely affected by the environmental conditions [START_REF] Erbil | Determination of diffusion coefficient-Vapor pressure product of some liquids from hanging drop evaporation[END_REF][START_REF] Carle | Contribution of convective transport to evaporation of sessile droplets: empirical model[END_REF]. Thus, predicting reliably the diffusion coefficient is primarily of importance for a thorough understanding of the evaporation process. The present work is a first attempt to address this issue. To this end, we construct a continuous function of output using response surface method based on a Kriging-based response surface [START_REF] Krige | A statistical approach to some basic mine valuations problems on the Witwatersrand[END_REF][START_REF] Matheron | Principles of geostatistics[END_REF], which allows to use numerical results as a black-box with a small amount of input and output data. This approach has been effectively used to deal with model parameter uncertainties in fluid dynamics and acoustic problems [START_REF] Echard | AK-MCS: an active learning reliability method combining Kriging and Monte Carlo simulation[END_REF][START_REF] Gayton | CQ2RS: a new statistical approach to the response surface method for reliability analysis[END_REF][START_REF] Wang | Shallow water sound source localization using the iterative beamforming method in an image framework[END_REF]. With the Kriging-based response surface, output result for any input model data can be easily estimated without further numerical computation. Relevant values of the diffusion coefficient can then be estimated by solving an inverse problem which is based on the experimental data and the proposed response surface. This response surface approach helps to interpret the experimental results and has also been applied in recent drop evaporation experiment under microgravity conditions onboard the Chinese scientific satellite SJ-10 (Hu et al. 2014).

The rest of the paper is organized as follows. The problem formulation and the methodology are described in "Problem Formulation". Numerical results are presented and discussed in "Results and Discussion", followed by conclusions drawn from this work.

Problem Formulation

We consider a millimeter-sized drop of an initial contact angle θ 0 and a constant contact radius R 0 resting on a heated substrate in still air, as shown in Fig. 1. The liquid is an Newtonian, incompressible fluid of constant density ρ, dynamical viscosity μ and thermal conductivity k. The surface tension σ is given by σ = σ 0γ (T -T 0 ), where σ 0 and γ are the surface tension evaluated at the reference temperature T 0 and its negative surface-tension coefficient, respectively. The bottom of the substrate is maintained at a temperature T s , and the gas far from drop is kept at a constant ambient temperature T ∞ . It is assumed that the axisymmetric sessile droplet maintains a sphericalcap shape and forms a constant contact line. It should be noted that, in practice, any chemical interactions, roughness and defects of the substrate, and the extend of wetting can cause contact angle hysteresis. We consider a pure liquid droplet on an ideal (smooth, perfectly flat and completely homogeneous) solid, a unique value of equilibrium contact angle is then determined by Young's equation. In other works, the contact angle hysteresis is neglected in our simulations. The assumption of constant-radius model is consistent with the evaporation process that we observed in experiments. In our ground and space experiments, we used a solid substrate and injected sufficient liquid volume to make sure that the contact line is well anchored to the substrate with its radii being at the edge of the substrate. Under such conditions, the evaporation process observed follows essentially the constant-radius mode during most of the drop's lifetime (90%). Therefore, it can be reasonably assumed that the contact radius R remains uncharged in our numerical simulations. Notice also that the lifetime of an evaporating drop is obtained by extrapolation to vanishing

Governing Equations

The physical problem is governed by the incompressible Navier-Stokes equations and the energy equation for the velocity vector u, hydrodynamic pressure p and temperature T , in the liquid phase :

ρ(∂ τ u + u • ∇u) = ∇ • -pI + μ ∇u + (∇u) T , (1) ∇ • u = 0, (2)
∂ τ T + u • ∇T = α∇ 2 T , ( 3
)
where I is identity tensor. Note that α = k/ρc p denotes the thermal diffusivity, c p is the specific heat capacity.

In the gas phase, vapor and heat transports in the atmosphere are solely by diffusion:

∂ τ c = D 0 ∇ 2 c, ( 4
)
∂ τ T = α∇ 2 T , (5)
where c stands for the molar concentration of the liquid vapor and D 0 is the vapor-in-air diffusion coefficient. The convection in the gas phase is neglected in this work although we do account for Marangoni flow along the droplet surface. The reason for that is twofold from the point view of fluid dynamics: a very small ratio of mass density between the gas and the liquid (i.e., about 1/800) and negligible buoyancy effects under microgravity conditions.

Boundary and Initial Conditions

The dynamic boundary condition at the droplet surface is the balance of the normal and tangential stresses:

n • (T • n) = -(2σ K + p ∞ ), (6) τ • (T • n) = γ ∇ T • τ (7)
where n and τ are, respectively, the unit normal and tangential vector to the liquid-gas interface (denoted by ), T = -pI + μ(∇u + (∇u) T ) stands for the full stress tensor of the fluid, K is the mean curvature of the liquidgas interface, and ∇ denotes the gradient tangent to the interface.

The heat flux across the interface experiences discontinuity due to the latent heat of vaporization H vap :

k g (∇T ) g • n -k l (∇T ) l • n = j m H vap , ( 8
)
where H vap denotes the latent heat of evaporation. The subscripts l and g represent the liquid phase and the gas phase, respectively, j m is the local evaporation flux, which can be calculated by Fick's law:

j m = -D 0 ∇c • n.
For the concentration of vapor c sat at the drop surface, it assumes to be saturated, defined by

c sat (T ) = p sat (T ) RT , ln p sat p ref = - H vap R 1 T - 1 T ref , ( 9
)
where p sat and p ref are the saturated vapor pressure at temperature T and T ref , respectively, R = 8.31451 J/(mol K) is the universal gas constant. The condition of mass conservation across the interface leads to:

ρ l (u • n -u ) = j m , ( 10
)
where u is the normal velocity of the interface. Given the total mass evaporation flux J m = j m d , the normal velocity can be geometrically determined:

u = - J m πρ l R 2 0 ⎛ ⎝ 1 - r sin θ R 0 2 -cos θ ⎞ ⎠ 1 + cos θ 1 -cos θ . ( 11
)
By assuming a spherical-cap shape of the droplet with the contact angle θ 0 and the contact radius R 0 , the initial height h 0 can be deduced by h 0 = R 0 tan(θ 0 /2). At initial time t = 0, the bulk flow inside the drop is at rest (i.e. u = 0), whereas the pressure within the drop is set equal to the ambient pressure augmented by the Laplace pressure

( p = 2σ 0 /R s with R s = (h 2 0 + R 2 0 )/(2h 0 )).
The initial temperature in the whole domain is set at the ambient temperature T ∞ (= 20 • C). The droplet evaporates into a non-saturated surrounding air with a relative humidity H . Then, the initial concentration in the vapor equals to c ∞ = Hp ∞ / RT ∞ .

We performed numerical simulations of transient fluid flow together with heat and mass transfer using COMSOL Multiphysics. Numerical procedure and model validation can be found in Chen et al. (2017a, b). Simulation results are subsequently used to construct a response surface as described below.

Response Surface

Response surface methodology is to model the behavior of a physical phenomenon versus the influential parameter variation using a certain number of numerical simulation results. The response surface, which is a function of concerned input parameters, is constructed using Kriging method in this paper. The reliability of the Krigingbased response surface is evaluated via the leave-one-out (LOO) cross-validation strategy [START_REF] Geisser | The predictive sample reuse method with applications[END_REF][START_REF] Saltelli | Variance based sensitivity analysis of model output: design and estimator for the total sensitivity index[END_REF]. Using uncertainty quantification analysis can help us numerically predict the change of concerned variables and estimate the unknown input. Here, we simply outline the methodology of the response surface method; more details about its principle and numerical procedure can be found in Chen et al. (2017b). In short, this analysis includes the following three steps: -First, we select some discrete points of inputs (diffusion coefficient and substrate temperature) and calculate the corresponding evaporation rates via numerical simulation. -Second, based on the numerical simulation results, a function of evaporation rate with respect to the diffusion coefficient and the substrate temperature is constructed using the Kriging-based response surface. -Finally, based on this response surface, with a measured output (e.g., the evaporation rate) and the substrate temperature, the diffusion coefficient can be estimated by solving the inverse of the response surface function.

Results and Discussion

Numerical results are presented for an ethanol droplet evaporating on a heated substrate covered with a thin film Teflon-PTFE. Simulation conditions correspond to our space experiments, in which the substrate temperature, as well as the environmental temperature is measured using a thermocouple. In order to maintain the substrate at a constant temperature T s , we inserted a thin film heater and a thermocouple fitted with a proportional-integralderivative (PID) controller to make sure that the temperature is uniformly distributed on the substrate. It is difficult, however, to measure precisely the environmental humidity; for an ethanol droplet evaporating in still air, the initial humidity is assumed to be 0, meaning that there is no ethanol vapor at the beginning of the experiment. We first present the determination of diffusion coefficient based on response surface method, and then propose an improved formula for the rate of evaporation.

Determining the Diffusion Coefficient Using Response Surface Method

The diffusion coefficient can be predicted by an indirect method with the known numerical results and the Krigingbased response surface. We first explain the numerical procedure to obtain the unknown parameter D 0 , and then compare the best estimated values with those obtained by a commonly used formulation.

It is realized that the diffusion coefficient is very difficult to measure in experiments, so we need to find a better method to determine its value. In the previous literature [START_REF] Gatapova | Evaporation of a sessile water drop on a heated surface with controlled wettability[END_REF][START_REF] Carle | Contribution of convective transport to evaporation of sessile droplets: empirical model[END_REF], the diffusion coefficient D 0 mainly depends on temperature T and pressure p ∞ :

D 0 (T , p ∞ ) = (D 0 • p) ref p ∞ T T ref 3/2 , ( 12
)
where (D 0 • p) ref = 1.337 Pa m 2 /s and T ref = 298.13 K for ethanol vapor in air. Equation 12 reveals that the diffusion coefficient is proportional to T 3/2 and inversely proportional to the atmosphere pressure p ∞ . Strictly speaking, the interfacial temperature T int should be used in Eq. 12. However, the interfacial temperature of an evaporating droplet exhibits a non-monotonic spatial distribution along the droplet surface due to several thermal effects at play, namely Marangoni effect and evaporative cooling on the droplet surface, heat conduction across the substrate and the liquid. In particular, the temperature drop at the liquid interface due to substantial evaporative cooling of a rapid evaporating droplet was found to increase with increasing evaporation rates. Instead of using the interfacial temperature which can only be obtained through a numerical simulation, we use a simplifying approach in two different ways. First, the substrate temperature T s is used in Eq. 12, which represents a limiting case for the diffusion coefficient (i.e., its upper bound), and secondly, we use an effective temperature [START_REF] Semenov | Boundary conditions for a one-side numerical model of evaporative instabilities in sessile drops of ethanol on heated substrates[END_REF]) to compute the diffusion coefficient. This effective temperature which lies between the ambient temperature T ∞ and the substrate temperature T s may be regarded to some extent as a "mean" interfacial temperature. We note that [START_REF] Gatapova | Evaporation of a sessile water drop on a heated surface with controlled wettability[END_REF] presented an iterative procedure by which the temperature dependence of the diffusion coefficient is taken into account to compute the evaporation rate of a sessile water drop on a heated substrate. While this is an interesting numerical procedure, our aim here is to infer the diffusion coefficient involved in a droplet evaporation experiment using response surface method.

T eff = [ √ T s T ∞ ( √ T s + √ T ∞ )/2] 2/3
There is also the added difficulty that the ambient pressure might not be measured or the measurements are not reliable, such that the estimation of the diffusion coefficient is imprecise. Consequently, the numerical simulations are inconsistent with the experimental results due to the improper input value. For this reason, the response surface method should be applied, in order to give a new way to obtain the unknown value of diffusion coefficient.

To construct response surface, a grid with 25 samples in (D 0 , T s ) is used, the values of the samples are shown by the circles in Fig. 2. The corresponding value of p ∞ , which is used in the numerical simulation, follows the relation of the parameters in Eq. 12. The required number of grid samples depends on the dimension of input parameters, the physical range of each input parameter and the complexity of response surface. In the considered numerical example, the grid with 25 samples is proven to be reasonable such that the Kriging response surface converges; we have tested cases with additional samples but they do not significantly improve the reliability of response surface. The lifetime of the evaporation t vap can be achieved, giving an average rate of evaporation | ṁ| = ρV 0 /t vap . Here, V 0 represents the initial volume of the droplet. Then, a 3D response surface of | ṁ|, denoted as

| ṁ| = f (D 0 , T s ), ( 13
)
is constructed based on the Kriging method with the numerical simulations. The constructed response surface which represents the evaporation rate variations with different substrate temperatures T s and diffusion coefficient D 0 , is shown in Fig. 3. The Kriging-based response surface is validated by the LOO method: the maximum crossvalidation is max s CV s (| ṁ|) = 3%, which means the response surface is reliable. Note that for a given substrate temperature T s , where Eq. 13 can be denoted by

| ṁ| = f T s (D 0 ), (14)
if f T s is a monotone function, the diffusion coefficient can be uniquely determined by

D0 = f -1 T s (| ṁ|). (15)
Therefore, by solving an inverse problem, i.e., Eq. 15 using the experimental results of substrate temperature and evaporation rate, we can infer the diffusion coefficient. Six samples of D 0 are assumed to be unknown and estimated by the measurements of T s and | ṁ|. The estimates of D 0 obtained from the response surface and the corresponding values computed using Eq. 12 (given T s and p ∞ ) are listed in Table 1. It should be noted that once we get the experimental results of | ṁ| and the imposed substrate temperature T s , the response surface method dos not require ambient pressure p ∞ to determine D 0 . Comparing the results in Table 1 indicates that Eq. 12 systemically underestimates diffusion coefficient. For example, if the substrate temperature is imposed at T s = 40 • C and the evaporation rate is | ṁ| = 2.028 × 10 -4 g/s, we obtain a value of D 0 = 30.07 mm 2 /s from the response surface, while Eq. 12 gives a value of D 0 = 28.43 mm 2 /s with T s and 27.06 mm 2 /s with T eff .

Revisiting the Evaporation Rate

Over the past decades, many efforts have been made to achieve a simplified, empirical formula for the (average) rate of evaporation of a sessile droplet. Indeed, such an empirical formula has been found to be beneficial to drop evaporation experiments as it provides a quick guide to designing and interpreting the experiments. The simplest and thus most widely used formula was proposed by [START_REF] Hu | Evaporation of a sessile droplet on a substrate[END_REF], based on a theoretical model of quasisteady diffusion-driven evaporation; it consists primarily of a contact-angle function representing the wettability of substrate. The formula also includes the concentration difference between the droplet surface and the surrounding environment, defined as follows:

| ṁ| = πR 0 D 0 cf (θ), ( 16
)
where the concentration difference c = c sat (T s ) -H c sat (T ∞ ) and the (dimensionless) contact-angle function f (θ) = 1.3 + 0.27θ 2 (θ in rad). It turns out that using Eq. 16 significantly overestimates the average rate of evaporation when compared with the space experimental results. The explanation for this follows by considering a liquid droplet where the temperature at the droplet interface is equal to the heating temperature of substrate. That means the temperature profile assumed by the model is instantly established in the whole droplet volume at the very beginning of evaporation process, which is obviously not the case, as an ambient, non-heated droplet is placed on the substrate during the experiments, and thus it would take some time for the flow of heat to develop. In other words, [START_REF] Hu | Evaporation of a sessile droplet on a substrate[END_REF] model underestimates the lifetime of the droplet evaporation, leading to a higher average rate of evaporation. A quantitative analysis is provided below. On the basis of our numerical results, we attempt to refine this formula taking into account the thermal effects such as substrate heating, evaporative cooling and Marangoni flow, which could allow us to quickly predict the evaporation rate of the space experiments in SJ-10. The key idea was to obtain a similar function which not only varies with the contact angle θ but also with the substrate temperature T s . To this end, we performed a series of numerical simulations for different initial contact angles θ , diffusion coefficient D 0 and substrate temperatures T s , resulting in the average rates of evaporation | ṁ| as a function of θ , D 0 and T s .

We first checked the relationship between the evaporation rate | ṁ| and the diffusion coefficient D 0 under different heating conditions. Numerical results are shown in Fig. 4. As expected, a linear scaling, i.e., | ṁ| ∝ D 0 , is recovered, in accordance with Eq. 16. It is further noted that the evaporation rate increases remarkably with the substrate temperature. This result is consistent with our previous findings (Chen et al. 2017a) that the evaporation rate of droplets could be well fitted by a power law with the substrate temperature, i.e., | ṁ| ∝ T b s with the power exponent b in the range of about 2.5. Now defining the ratio | ṁ|/(πR 0 D 0 c) (≡ f (θ, T s )), we plot its value as a function of the contact angle for the substrate temperatures ranging from 25 to 40 • C (Fig. 5). Indeed, the value is substantially lower than that given by the contact-angle function proposed by [START_REF] Hu | Evaporation of a sessile droplet on a substrate[END_REF], and the ratio decreases as the Finally, by curve fitting to the numerical results, we arrive at a function f (θ, T s ) = 1.305 + 0.12θ 2 -0.00175T s (θ in rad and T s in • C). At first sight, the proposed fitting function f (θ, T s ) = a + b θ 2 + c T s is physically meaningless owing to the two absolutely independent parameters (the contact angle θ and the substrate temperature T s ), which, of course, are measured in absolutely different units. However, it should be understood that the function is dimensionless, meaning that the coefficients b and c have, respectively, the units of rad -2 and • C -1 . With this new function, Eq. 16 could now be extended its application to a droplet evaporating on a heated substrate. Nevertheless, we caution that additional comparisons with other numerical and/or experimental results are needed to assess the validity of this formula under a wide range of conditions.

Conclusions

In this paper, we have presented a numerical model to investigate the evaporation process of a pinned, millimetersized drop deposited on a heated solid substrate. A limited number of numerical simulations have been performed to construct a continuous function of the evaporation rate with respect to the vapor diffusion coefficient in the atmosphere and the heating temperature using the Krigingbased response surface, from which relevant values of the diffusion coefficient involved in recent drop evaporation space experiment are obtained by solving an inverse problem. Our numerical results show that the conventional method systematically leads to a lower value of the diffusion coefficient than that predicted by the present approach. This methodology relies on a combination of numerical simulations and readily accessible experimental data (namely, heating temperature and overall evaporation rate) to determine precisely the values of some key physical parameters related to the experiments, and the proposed numerical procedure is general, so that the methodology could be a useful tool for situations requiring repeated simulations. We have also revisited a widely used empirical formula for the prediction of the evaporation rates of a sessile droplet, and proposed a new, improved empirical relation for the droplets evaporating on a heated substrate. The proposed formula, inspired from [START_REF] Hu | Evaporation of a sessile droplet on a substrate[END_REF] work, incorporates a combined thermal effect, namely substrate heating, evaporative cooling and Marangoni flow. Overall, this study provides a new way to simply and accurately determine the diffusion coefficients involved in the drop evaporation experiments and could help better devise future experiments.

  Fig.1A sessile droplet evaporating on a heated substrate in still air in a cylindrical coordinate system with radial coordinate r and axial coordinate z

Fig. 2

 2 Fig. 2 Sample points in (D 0 , T s ) for constructing response surface

Fig. 3

 3 Fig. 3 Response surface of | ṁ| as a function of D 0 and T s

Fig. 4 Fig. 5

 45 Fig. 4 Evaporation rate | ṁ| as a function of diffusion coefficient D 0 for various substrate temperatures T s . Initial conditions: R 0 = 6 mm, θ 0 = 27 • , T ∞ = 20 • C, and H = 0

Table 1

 1

	Diffusion coefficient D 0 (mm 2 /s) obtained from the	(T s , | ṁ|)	Determination of D 0
	Kriging response surfaces, for
	given experimental T s ( • C) and		Using response surface	Using Eq. 12 with T s	Using Eq. 12 with T eff
	| ṁ| (g/s), compared with its
	values from Eq. 12

Acknowledgments This work was financially supported by the National Natural Science Foundation of China (Grant No. 11532015), the Strategic Pioneer Program on Space Science, Chinese Academy of Sciences (Grant No. XDA 04020202-02) and by the CNES (Centre National d'Etudes Spatiales). X. Chen has benefited from financial support from the Guangxi's Key Laboratory Foundation of Manufacturing Systems and Advanced Manufacturing Technology (Grant No. 17-259-05-002Z) and China Postdoctoral Science Foundation (2018M633113).


Normah Mohd-Ghazali

Patrice Estellé
email: patrice.estelle@univ-rennes1.fr

Salma Halelfadl

Thierry Maré

Choon Tng

Siong

Ummikalsom Abidin

Thermal and hydrodynamic performance of a microchannel heat sink with carbon nanotubes nanofluids: Effect of concentration and channel section

Keywords: Microchannel heat sink, Optimized, Nanofluids, Surfactant, Carbon nanotube, Genetic algorithm

Comparison of optimized thermal and hydrodynamic performances of a circular and square microchannel heat sink cooled with carbon nanotube nanofluid, using experimentally obtained data.

• Analysis of surfactant and viscosity effects on CNT nanofluid performance.

• Utilization of multi-objective genetic algorithm (MOGA) to simultaneously minimize the thermal resistance and pressure drop.

Introduction

As microelectronics devices design continue towards miniaturization and becoming more complex, the power density is increasing with effective cooling of the active components needed to ensure their optimum operation. Sensitive components such as integrated circuits, which are the main generators of heat, are vulnerable to temporary or permanent failure if overheated. Thus, the landmark study by Tuckerman and Pease in 1981 [START_REF] Tuckerman | High-performance heat sinking for VLSI[END_REF] of using a microchannel heat sink (MCHS) in very large-scale integration (VLSI) circuits have attracted great interests from many researchers due to its proven effective and compact heat removal from a small area.

Efforts into the enhancement of the cooling capacity of a MCHS have much focussed on the structural materials [START_REF] Kleiner | High performance forced air cooling scheme employing microchannel heat exchangers[END_REF][START_REF] Gamrat | Conduction and entrance effects on laminar liquid flow and heat transfer in rectangular microchannels[END_REF][START_REF] Peng | Heat transfer characteristics of water flowing through microchannels[END_REF][START_REF] Choquette | Optimum design of microchannel heat sinks[END_REF][START_REF] Zhimin | Optimum thermal design of microchannel heat sinks[END_REF][START_REF] Kosar | Effect of substrate thickness and material on heat transfer in microchannel heat sinks[END_REF], different channel shapes [START_REF] Perret | Microchannel integrated heat sinks in silicon technology[END_REF][START_REF] Hestroni | Heat transfer in microchannels: comparison of experiments with theory and numerical results[END_REF][START_REF] Mchale | Heat transfer in trapezoidal microchannel of various aspect ratio[END_REF][START_REF] Ghazali-Mohd | Comparison of the Optimized Thermal Performance of Square and Circular Ammonia-cooled Microchannel Heat Sink with Genetic Algorithm[END_REF], effect of entrance channel shapes [START_REF] Moradikazerouni | Comparison of the effect of five different entrance channel shapes of a micro-channel heat sink in forced convection with application to cooling a supercomputer circuit board[END_REF] and various types of coolants [START_REF] Tuckerman | High-performance heat sinking for VLSI[END_REF][START_REF] Kleiner | High performance forced air cooling scheme employing microchannel heat exchangers[END_REF][START_REF] Adham | Optimization of an ammonia-cooled rectangular microchannel heat sink using multi-objective non dominated sorting genetic algorithm (NSGA2)[END_REF][START_REF] Hung | Heat transfer enhancement in microchannel heat sinks using nanofluids[END_REF][START_REF] Halelfadl | Optimization of thermal performances and pressure drop of rectangular microchannel heat sink using aqueous carbon nanotubes based nanofluid[END_REF]. To date, it seems that research into the structural material and common channel geometries have almost reached its saturation stage [START_REF] Adham | Thermal and hydrodynamic analysis of microchannel heat sinks: A review[END_REF]. However, the introduction of fins within the channels can lead to heat transfer improvement as shown in [START_REF] Li | Flow and heat transfer characteristics and optimization study on the water-cooled microchannel heat sinks with dimple and pinfin[END_REF]. This numerical analysis based on optimization procedure demonstrated that heat transfer of water in laminar regime within microchannel heat sink with pin-fin and dimple is enhanced with both the increase in pin diameter and the reduction in stream-wise spacing. The effect of triangular cavities and rectangular ribs on heat transfer performance of microchannel heat sink under laminar regime was also investigated numerically in [START_REF] Li | Characteristics of laminar flow and heat transfer in microchannel heat sink with triangular cavities and rectangular ribs[END_REF]. These authors evidenced uniform and lower temperature at the substrate of the heat sink and reported better heat transfer performance with cavities and ribs than without. They also proposed the better combination of cavities and ribs to achieve higher performance. Moreover, thermal performance of air-cooled heat sink was obtained from simulation by Moradikazerouni et al. [START_REF] Moradikazerouni | Investigation of a computer CPU heat sink under laminar forced convection using a structural stability method[END_REF] who also studied spatial variation of plate fins in forced convection flow. Vinoth and Senthil Kumar [START_REF] Vinoth | Channel cross section effect on heat transfer performance of oblique finned microchannel heat sink[END_REF] demonstrated from experiments that trapezoidal cross-section oblique finned micro-channel appear to be more suitable for heat transfer enhancement than square and semicircle cross section. In addition, Al2O3/water nanofluid of volume fraction 0.25% also increase heat transfer in comparison to water. Convective heat transfer and pressure drop of microchannel heat sink with multiple zigzag flow channels were experimentally studied by Duanthongsuk and Wongwises considering SiO2 -water nanofluids [START_REF] Duangthongsuk | An experimental investigation on the heat transfer and pressure drop characteristics of nanofluid flowing in microchannel with multiple zigzag flow channel structures[END_REF]. They compared the performance of a continuous zigzag flow channel and a zigzag flow channel with a single crosscutting at the middle position. An higher Nusselt number was obtained with nanofluids compared to water while better heat transfer performance was obtained with the presence of single cross-cutting. Few time ago, the fluid flow and heat transfer of γ-AlOOH-50/50 waterethylene glycol nanofluid flowing through a sinusoidal wavy channel was investigated by Vo et al. [START_REF] Vo | Numerical investigation of γ-AlOOH nano-fluid convection performance in a wavy channel considering various shapes of nanoadditives[END_REF] showing also the influence of nanoparticles shape in heat transfer performance.

Moreover, improvement of heat transfer capabilities with variations in the thermophysical properties of fluids have been continuously of interest. The thermal conductivity and heat capacity -two properties responsible towards a coolant ability to remove heat -of air, water, and refrigerants could be manipulated through the addition of moisture in the case of gas, and particles for the liquids. The latter has become the trending coolant with continuous studies into different particles with various concentrations dispersed in various base fluids [START_REF] Byrne | Experimental thermal-hydraulic evaluation of CuO nanofluids in microchannels at various concentrations with and without suspension enhancers[END_REF][START_REF] Adham | Optimization of nanofluid-cooled microchannel heat sink[END_REF][START_REF] Mahian | Recent advances in modeling and simulation of nanofluid flows-Part I: Fundamentals and theory[END_REF]. The enhanced heat transfer attainable with nanofluids have made them as a popular choice of coolants [START_REF] Maré | Consideration of Carbon Nanotubebased Nanofluid in Thermal Transfer[END_REF], as also evidenced in some of the previous references.

Unlike the standard nanofluids with spherical oxide or metallic particles immersed in based fluid, fewer studies have been completed on the use of carbon nanotubes (CNT) based nanofluids in microchannel heat-sink. Recently, Arani et al. [START_REF] Arani | Heat transfer improvement of water/single-wall carbon nanotubes (SWCNT) nanofluid in a novel design of a truncated double-layered microchannel heat sink[END_REF] numerically investigated heat transfer rate of CNT nanofluids flowing through a truncated double-layered microchannel heat sink. From theoretical models for nanofluid properties, they showed the increase in heat transfer with nanoparticle content and the reduction in thermal resistance with the increase in Re or nanoparticle loading. An experimental investigation on thermal performance and fouling of CNT/water nanofluid inside a heat sink with rectangular microchannel was performed by Sarafraz et al. [START_REF] Sarafraz | Fouling formation and thermal performance of aqueous carbon nanotubes nanofluid in a heat sink with rectangular microchannel[END_REF]. Temperature profile along the channel length was reported to decrease with nanoparticle content that also lead to the increase in heat transfer compared to water. Lastly, the laminar heat transfer of kerosene nanofluid/multi-walled carbon nanotubes in the microchannel heat sink was numerically studied by Arabpour et al. [START_REF] Abedin | The study of heat transfer and laminar flow of kerosene/multi-walled carbon nanotubes (MWCNTs) nanofluid in the microchannel heat sink with slip boundary condition[END_REF]. Increasing Reynolds number and volume fraction of nanoparticles produced the increase in convection heat transfer that reduced the maximum temperature. The reduction in thermal resistance that causes the enhancement of heat transfer amount was also reported to be dependant of slip velocity at walls. Different results have been reported in the literature [START_REF] Su | The effect of carbon nanotubes on the physical properties of a binary nanofluid[END_REF][START_REF] Murshed | Superior thermal features of carbon nanotubesbased nanofluids -A review[END_REF][START_REF] Murshed | A state of the art review on viscosity of nanofluids[END_REF] on the analytical models of the thermophysical properties on CNT. Thus, only experimentally obtained properties of the CNT nanofluid, specifically for thermal conductivity and viscosity, have to be preferably used in modelling of its performance in any cooling system as reported by Halelfadl et al. [START_REF] Halelfadl | Optimization of thermal performances and pressure drop of rectangular microchannel heat sink using aqueous carbon nanotubes based nanofluid[END_REF] and Nik Mazlam et al. [233]. It has been established that dispersion and stability of the nanoparticles can be prevented with the addition of surfactant though their specific effects have hardly been studied. As reported earlier in laminar pipe flow [START_REF] Minea | Numerical study on CNT nanofluids behavior in laminar pipe flow[END_REF], experimentally obtained thermophysical properties of CNT nanofluids are here used to determine the optimized thermal and hydrodynamic performance of a MCHS. Two types of CNT nanofluids are investigated and mainly differs in types of surfactants used, lignin and sodium polycarboxylate, respectively.

Although circular geometry is better in terms of fluid flow due to the absence of corners, fabrication of a circular MCHS has been an issue till late. Advances in the manufacturing technology could make it feasible and advantageous in heat transfer considering the physics of the flow. Thus, the objective of this study is also to compare the performance of circular and square geometry through simultaneous minimization of the thermal resistance and pumping power of a MCHS cooled with CNT nanofluids. The thermal resistance model is utilized with multi-objective optimization algorithm (MOGA), a search tool the authors found to be capable of a global search in the presence of multiple conflicting objectives [START_REF] Ghazali-Mohd | Comparison of the Optimized Thermal Performance of Square and Circular Ammonia-cooled Microchannel Heat Sink with Genetic Algorithm[END_REF][START_REF] Moradikazerouni | Comparison of the effect of five different entrance channel shapes of a micro-channel heat sink in forced convection with application to cooling a supercomputer circuit board[END_REF][START_REF] Halelfadl | Optimization of thermal performances and pressure drop of rectangular microchannel heat sink using aqueous carbon nanotubes based nanofluid[END_REF][START_REF] Adham | Optimization of nanofluid-cooled microchannel heat sink[END_REF][START_REF] Nik-Mazlam | Thermal and hydrodynamic performance of a microchannel heat sink cooled with carbon nanotubes nanofluid[END_REF]. It has been known that the pressure drop tends to increase with the increase in the thermal resistance of a MCHS, MOGA has been used by the authors in their previous studies with a MCHS. The model developed by Ghazali-Mohd et al. [START_REF] Ghazali-Mohd | Comparison of the Optimized Thermal Performance of Square and Circular Ammonia-cooled Microchannel Heat Sink with Genetic Algorithm[END_REF] is used here to look into how different surfactant affecting the thermophysical properties of CNT nanofluids contributes on the thermal and hydrodynamic performance of a CNT nanofluid-cooled MCHS.

Mathematical modeling

A schematic structure of the aluminium circular and square MCHS used is shown in Fig.

1 [START_REF] Ghazali-Mohd | Comparison of the Optimized Thermal Performance of Square and Circular Ammonia-cooled Microchannel Heat Sink with Genetic Algorithm[END_REF]. The MCHS has the dimension L×W of 1×1 cm with a height of 0.213 mm. The size is based on previous studies [START_REF] Tuckerman | High-performance heat sinking for VLSI[END_REF][START_REF] Halelfadl | Optimization of thermal performances and pressure drop of rectangular microchannel heat sink using aqueous carbon nanotubes based nanofluid[END_REF]. A series of parallel channels run the length of the MCHS. Each microchannel has a hydraulic diameter, Hc. A uniform constant heat flux is assumed to be applied from the bottom of the MCHS which is connected to the electronics chip. Working fluid at the flow rate of G = 4.7×10 -6 m 3 /s enters the series of parallel channels along the length of the MCHS to remove the heat flux generated. On top of the MCHS, it is covered by an adiabatic cover plate. The entire chip area can be cooled by this system. The dimensions of the MCHS are given in Table 1.

Table 1. Geometrical properties of the model [START_REF] Ghazali-Mohd | Comparison of the Optimized Thermal Performance of Square and Circular Ammonia-cooled Microchannel Heat Sink with Genetic Algorithm[END_REF] Length, L(m) 0.01

Width, W (m) 0.01 Substrate thickness, t (m) 213 10
Thermal conductivity, k (W/m•K) 238

The total thermal resistance is contributed by the substrate below the channels, Rcond, the convective heat transfer, Rconv, and the coolant heat capacity, Rcap, where the resistances are assumed to be in series. In the optimization process of the MCHS, Hc and the ratio of the wall width, Ww, to Hc are the variables to be optimized to achieve simultaneous minimization of the thermal resistance and pumping power. The total thermal resistance according to the thermal resistance model can be determined from,

!" $% (1)
where $% takes the value of 1/3 and 1/ &, for the square and circular MCHS respectively, while takes the value of unity and 4/ & for the square and circular MCHS respectively. The hydraulic diameter, Hc, and channel wall width to channel width ratio, β, are regarded as variables,

' ( (2)
Similar to the thermal resistance, the hydraulic resistance is represented by the total pressure drop used in [START_REF] Ghazali-Mohd | Comparison of the Optimized Thermal Performance of Square and Circular Ammonia-cooled Microchannel Heat Sink with Genetic Algorithm[END_REF]:

∆* + , -1 ' (3)
In this investigation of the effect of surfactants on thermphysical properties and heat transfer performance of nanofluids flowing in MCHS, m for the circular and square channels used are , respectively. Equation (1) and equation ( 3) have conflicting effects. Both are desired to be as low as possible but the lowering of Equation ( 1) results in an increase in equation ( 3), and vice versa. Thus, simultaneous minimization of both functions is sought. This can be achieved with MOGA [START_REF] Ghazali-Mohd | Comparison of the Optimized Thermal Performance of Square and Circular Ammonia-cooled Microchannel Heat Sink with Genetic Algorithm[END_REF].

Optimization with MOGA

Generally, optimization problems consist in discovering the best solution from all feasible solutions. The optimization procedure applied here is based on past studies [START_REF] Ghazali-Mohd | Comparison of the Optimized Thermal Performance of Square and Circular Ammonia-cooled Microchannel Heat Sink with Genetic Algorithm[END_REF], [START_REF] Moradikazerouni | Comparison of the effect of five different entrance channel shapes of a micro-channel heat sink in forced convection with application to cooling a supercomputer circuit board[END_REF], [START_REF] Hung | Heat transfer enhancement in microchannel heat sinks using nanofluids[END_REF], mathematically represented as,

+345 ̅ 7̅ 85 9̅ , 5 . 9̅ , … , 5 < 9̅ = (4)
for m = 1, 2, …, M for the number of objective functions, in this study M = 3. Meanwhile, 7̅ is a vector containing the design variables which in this case are Hc and β.

The optimization solver of multi-objective optimization using Genetic Algorithm is chosen as the current problem exhibited a multi-objective nature. The two objective functions to be minimized are selected to describe the overall performance of the microchannel heat sink.

The first objective function is associated with the thermal performance that is the total thermal resistance, Equation (1), while the second is related to the hydrodynamic performance that is the pressure drop, Equation (3). Next, the bounds are created for the design variables. Two design variables used in the optimization procedure include the hydraulic diameter, Hc, and channel wall width to channel width ratio, β, with their lower and upper bounds limit specified in Table 2, based on past studies. In the optimization process of a MCHS, Hc, and β were regarded as variables while the other parameters were fixed. After that, the optimization options such as population, selection, reproduction, mutation, crossover, migration and stopping criteria were also set in order to generate a Pareto optimal front for both objective functions. Several solutions were also generated to determine the most appropriate combinations of the optimization options. Those options are specified in Table 3.

Carbon Nanotube Nanofluids

The water based-CNT nanofluids presently considered were previously investigated [START_REF] Estellé | Thermal conductivity of CNT water based nanofluids: Experimental trends and models overview[END_REF]. In short, they consisted in carbon nanotubes (density 1800 kg/m C ; purity 90%) of 9.2 nm and 1.5 µm in average diameter and length respectively, dispersed in a mixture of water and surfactant. It is noted that the weight ratio between surfactant and nanoparticle content was kept constant for each studied concentration in nanotubes. Two types of surfactant were considered namely lignin and sodium polycarboxylate. They were used to improve the dispersion and stability of CNT within water and reduce clogging and sedimentation with time. Their preparation and the evaluation of their thermophysical properties at 20°C were reported earlier for volume fraction ranging from 0.0055% to 0.55%. In brief, density and heat capacity were obtained from mixing rules using equations developed by Pak and Cho [START_REF] Pak | Hydrodynamic and heat transfer study of dispersed fluids with sub -micron metallic particles[END_REF] and O'Hanley et al. [START_REF] Hanley | Measurement and model validation of nanofluid specific heat capacity with Differential Scanning Calorimetry[END_REF] respectively as described in [START_REF] Estellé | Thermophysical properties and heat transfer performance of carbon nanotubes water-based nanofluids[END_REF]. Thermal conductivity and viscosity were experimentally determined from well-designed procedures in [START_REF] Estellé | Thermal conductivity of CNT water based nanofluids: Experimental trends and models overview[END_REF][START_REF] Estellé | Lignin as dispersant for water-based carbon nanotube nanofluids: Impact on viscosity and thermal conductivity[END_REF][START_REF] Halelfadl | Caractérisation des propriétés thermo-physiques et d'échanges de chaleur des nanofluides à base de nanotubes de carbone[END_REF]. All these results are shown in Figure 3. In the following, nanofluids produced with lignin are labelled N2. N3 is used for nanofluids obtained with sodium polycarboxylate as surfactant.

So, Figure 3 shows that density and heat capacity of nanofluids increases, respectively decreases, with volume fraction as generally reported in literature. In addition, both properties are similar for both nanofluids. Regarding the thermal conductivity of nanofluids, it obviously increases with nanoparticle content with similar trend and values for both nanofluids. While the type of surfactant used in the CNT nanofluid preparation has no effect on the previous thermophysical properties, a difference is observed for viscosity. The effect is graphically represented in Figure 3 where, as the volume fraction increases just beyond 0.03%, the viscosity for N2 and N3 diverged in value, N3 being highly affected with a sharp increase in viscosity up to 0.55%. These graphs already indicate the necessity to appraise the types of surfactant used and its effect on the desired properties looked for, the thermal and hydrodynamic performances.

They also show the interest of using experimental data rather than models that could not be relevant for nanofluids.

All thermophysical properties at 20°C are modeled through the experimental data reported in Figure 3 in order to make easier their use for the optimization procedure. These equations are presented in Table 5. Based on previous observations, similar equations were used with both nanofluids for density, heat capacity and thermal conductivity.

Results and discussion

The effects of the particle volume fraction on the thermal resistance and pressure drop is presented first for N2 in a circular MCHS, for the volume fractions of 0.1%, 0.2%, 0.3%, and 0.5% at 20°C. The pattern and trend exhibited are expected for the thermal and hydrodynamic performance of a MCHS where a decrease in the thermal resistance is followed by an increase in the pressure drop as shown in Fig. 4. Up to 0.3%, the difference between the volume fractions does not seem to affect much on the performance of the MCHS. A previous study [START_REF] Perret | Microchannel integrated heat sinks in silicon technology[END_REF] showed that the increased nanoparticles volume fraction increased the thermal conductivity of the coolant as well as increased the convective heat transfer coefficient which in turn reduces the convective thermal resistance.

Since both N2 and N3 nanofluids have the same nanotubes and thermophysical properties except for the viscosity, a similar trend is expected with their thermal and hydrodynamic performances.

Thus, a comparison of the performance in a circular and square is presented next in Fig. 5, Fig. 6 and Fig. 7 for the volume fraction of 0.1%, 0.3% and 0.5% respectively. For each case of a circular or square MCHS, at any particular thermal resistance, Rth, the pressure drop is higher for N3 than for N2, showing the effect of the surfactant used and its influence on dispersion state of nanotubes and viscosity increase. For a circular MCHS at 0.63 K/W, the pressure drop for N3 is 21% higher than that for N2, while for the square MCHS at 0.63 K/W, the pressure drop for N3 is 28% higher than that for N2. Hence, N2 nanofluid is more preferable if a lower pumping power is desirable for the same thermal performance. The difference in the performance becomes more significant as the volume fraction increases as seen in Fig. 6 and Fig. 7. The pressure drop also increases as the volume fraction is increased at any thermal resistance. As the hydraulic diameter increases, the expected decrease in pressure drop is observed. N3, however, showed a significant increase in the pressure drop which subsequently will result in a higher pumping power for the heat sink to operate with the coolant. The type of surfactant used has no effect on the thermal performance as listed in Table 6 for 0.1% volume fraction, selected at the hydraulic diameter of 201 mm. The wall width, Ww, and number of microchannels, n, in each MCHS are also listed in Table 6. An increase of 29% and 28% in pressure drop is expected when N3 is used instead of N2, for a circular and square MCHS respectively. For the same channel height and number of channels per MCHS, the thermal performance is the same for N2 and N3 but hydrodynamically, N2 outperformed N3 because of the change in viscosity. This applies for both the circular and square geometries. The higher pumping power required for the circular MCHS is due to the smaller cross-sectional area involved, 20% smaller than that of the square geometry. The same effect is observed at other volume fractions. With a lower thermal resistance at a higher pumping power, flow in circular channels eliminates the potential secondary flow that could occur in the square MCHS. The phenomenon could reduce the effectiveness of the channel wall and consequently affect the thermal performance of the MCHS.

Finally, the outcomes of this study showed that serious consideration is needed with the selection of surfactant when energy use is an issue. Actually, while stability and dispersion of nanoparticle is improved, the thermophysical properties can be modified and leads to drastic changes in hydraulic and thermal performance. Utilization of a MCHS does not only involve 1 unit of 1 cm by 1 cm, but many which translate into a high pumping power requirement.

Conclusion

We report in this study the optimization of a microchannel heat sink (MCHS) cooled with carbon nanotubes nanofluids flowing under laminar regime. This was done using multiobjective genetic algorithm (MOGA) to minimize simultaneously the thermal resistance and pressure drop of the MCHS. This optimization analysis was performed from viscosity and thermal conductivity values of carbon nanotubes nanofluids at 20°C that were experimentally obtained. Two types of CNT nanofluids were studied and mainly differs in types of surfactants used, lignin and sodium polycarboxylate, respectively. The effect of nanotube volume content and shape of channel section on the thermal and hydrodynamic performances was also investigated.

It was reported that nanofluids with lignin as surfactant produced a better thermal and hydrodynamic performance for the MCHS compared to the ones containing sodium polycarboxylate. This is attributed to higher viscosity of these later nanofluids as all other thermophysical properties are similar. As an example, a 29% and 28% increase in pressure drop is expected for sodium polycarboxylate based nanofluid at a volume fraction of 0.1% for a circular and square MCHS respectively. The pattern is observed at other volume fraction as well, even higher pressure drop with increasing volume fraction. Also, circular channel geometry presents a lower thermal resistance than square channel geometry. Thus, it is imperative to properly determine the properties experimentally or use relevant models to predict the expected performance of nanofluid in flowing energy systems like microchannels before deciding on the nanofluids use in real applications.

Figure 1 .

 1 Figure 1. Microchannels considered in the study.

Figure 2 .

 2 Figure 2. Optimization procedure (adapted from [14])

Figure 3 .

 3 Figure 3. Thermophysical properties of CNT nanofluids at 20°C -effect of CNT concentration and surfactant

Figure 4 .

 4 Figure 4. Pressure drop and thermal resistance of a circular MCHS at various volume fractions.

Figure 5 .Figure 6 .

 56 Figure 5. Comparison of the pressure drop and thermal resistance of the circular and square MCHS at 0.1%.

Figure 7 .Fig. 8 andFigure 8 .

 788 Figure 7. Comparison of the pressure drop and thermal resistance of the circular and square MCHS at 0.5%.

Figure 9 .

 9 Figure 9. Effect of hydraulic diameter on the pressure drop of the circular and square MCHS at 0.5%.

Table 2 .

 2 Design variables limitsDesign variables Lower Limits Upper limits

	Hc	0.0005	0.003
	β	0.1	0.2

Table 3 .

 3 Optimization options set for the current study

	Population size	120
	Selection function	Tournament
	Crossover fraction	0.8
	Mutation function	Constraint Dependent
	Crossover function	Intermediate
	Migration direction	Forward
	Stopping criteria (function tolerance)	1 10 >
	Lastly, Fig. 2 describes the complete optimization procedure employed in the current
	research.

Table 5 .

 5 Models for the thermophysical properties of CNT nanofluids.

		0.8				12
	Thermal conductivity (W/m.K)	0.55 0.6 0.65 0.7 0.75	N2	N3	Viscosity (mPa.s)	2 4 6 8 10	N2	N3
		0.5				0
		0	0.2	0.4	0.6	0		0.2	0.4	0.6
			volume fraction (%)				volume fraction (%)
		D E/+. G		y = 0.0072ln(x) + 0.6405
		H DI/+ C		y = 0.0019x + 1.0011
		J K/DI. G		y = -36.794x + 4178.7
		µ (mPa.s)		y = 2.7933x 2 -0.0924x + 1.0539 (N2)
					y = 28.98x 2 -0.8218x + 1.1761 (N3)

Table 6 .

 6 Pressure drop (DelP) and pumping power (Pp) at the lowest thermal resistance, 0.1%.

		Hc (µm) Rth (K/W) DelP (kPa) Pp (W) β	Ww (β×Hc) (µm) n
	N2C (0.1%)		0.27	25.8	0.121	0.0154	3.1	49
	N3C (0.1%)	201	0.27	33.4	0.157	0.0101	2.03	49
	N2S (0.1%)		0.34	18.0	0.085	0.0106	2.13	49
	N3S (0.1%)		0.34	23.1	0.109	0.0115	2.31	49

Acknowledgements

The authors would like to thank the Universiti Teknologi Malaysia for the research grant Vot19H60 for the funding to do the research.


Rafael Ferreira Da Silva

Anne-Cécile Orgerie
email: anne-cecile.orgerie@irisa.fr

Henri Casanova
email: henric@hawaii.edu

Ryan Tanaka
email: ryanyt@hawaii.edu

Ewa Deelman
email: deelman@isi.edu

Frédéric Suter
email: frederic.suter@cc.in2p3.fr

Accurately Simulating Energy Consumption of I/O-intensive Scientific Workflows

Keywords: Scientific workflows, Energy-aware computing, Workflow profiling, Workflow scheduling

While distributed computing infrastructures can provide infrastructure-level techniques for managing energy consumption, application-level energy consumption models have also been developed to support energy-efficient scheduling and resource provisioning algorithms.

In this work, we analyze the accuracy of a widely-used application-level model that have been developed and used in the context of scientific workflow executions. To this end, we profile two production scientific workflows on a distributed platform instrumented with power meters. We then conduct an analysis of power and energy consumption measurements. This analysis shows that power consumption is not linearly related to CPU utilization and that I/O operations significantly impact power, and thus energy, consumption. We then propose a power consumption model that accounts for I/O operations, including the impact of waiting for these operations to complete, and for concurrent task executions on multi-socket, multi-core compute nodes. We implement our proposed model as part of a simulator that allows us to draw direct comparisons between real-world and modeled power and energy consumption. We find that our model has high accuracy when compared to real-world executions. Furthermore, our model improves accuracy by about two orders of magnitude when compared to the traditional models used in the energyefficient workflow scheduling literature.

Introduction

Computational workloads that require from a few hours to a few months of execution are commonplace in scientific simulations. These simulations often comprise individual computational (but often I/O-intensive) tasks with some dependency structure, which is why many scientists today formulate their computational problems as scientific workflows [START_REF] Taylor | Workflows for e-Science[END_REF]. To obtain simulation results within acceptable time-frames, large scientific workloads are executed on distributed computing infrastructures such as grids and clouds [START_REF] Romanus | The anatomy of successful ECSS projects: Lessons of supporting high-throughput high-performance ensembles on XSEDE[END_REF]. The need to manage energy consumption across the entire suite of information and communication technology has received significant attention in the last few years [START_REF] Baliga | Green cloud computing: Balancing energy in processing, storage, and transport[END_REF][START_REF] Beloglazov | A taxonomy and survey of energy-efficient data centers and cloud computing systems[END_REF]. As a result, large data-centers have developed techniques for managing cooling and energy usage at the infrastructure level. Concurrently, researchers have investigated application-level techniques and algorithms to enable energy-efficient executions [START_REF] Orgerie | A survey on techniques for improving the energy efficiency of large-scale distributed systems[END_REF]. In the context of scientific workflows, researchers have proposed a range of energy-aware workflow task scheduling or resource provisioning algorithms [START_REF] Ghose | Energy efficient scheduling of scientific workflows in cloud environment[END_REF][START_REF] Li | Cost and energy aware scheduling algorithm for scientific workflows with deadline constraint in clouds[END_REF][START_REF] Pietri | Energy-aware workflow scheduling using frequency scaling[END_REF][START_REF] Shepherd | Workflow scheduling on power constrained vms[END_REF][START_REF] Wu | Soft error-aware energy-efficient task scheduling for workflow applications in DVFS-enabled cloud[END_REF]. Results therein are obtained based on a model of power consumption that is easy to instantiate but that makes strong assumptions: power consumption is considered to be linearly correlated with CPU utilization, and equally divided among virtual machines or CPU cores within a computational node. An interesting question is whether this model is accurate in practice, and whether it can be applied to I/O-intensive workflow executions.

Our broad objective in this work is to characterize the energy consumption behavior of complex workflow applications that execute on distributed platforms. We profile real scientific workflow applications on a distributed platform that comprises multi-socket, multi-core compute nodes equipped with power meters. We select two widely scientific workflows, each of which has some I/Ointensive tasks, for which we conduct a comprehensive analysis of the power and energy consumption of their execution. Via this analysis, we quantify the accuracy, or lack thereof, of power consumption models commonly used in the energy-efficient workflow scheduling literature. We then propose a more accurate power consumption model. More specifically, this work makes the following contributions:

1. The power and energy consumption profiles of two real I/O-intensive scientific workflow applications; 2. A comprehensive analysis of these profiles with respect to resource utilization and I/O operations; 3. An evaluation of the accuracy of the power model that is widely used in workflow scheduling research; 4. A power consumption model for I/O-intensive workflows that accounts for the allocation of cores to sockets, CPU utilization, and I/O operations; 5. An experimental evaluation of the proposed model that shows that it can produce nearly accurate energy consumption estimates, with improvements over traditional models by almost two orders of magnitude.

Workflow Characterization

The analysis presented in this work is based on the execution of two production scientific workflow applications on the Grid'5000 [START_REF] Balouek | Adding virtualization capabilities to the grid5000 testbed[END_REF] platform. Grid'5000 is a testbed for experiment-driven research, which provides resource isolation and advanced monitoring and measurement features for the collection of power consumption traces. We consider these two I/O-intensive workflows:

-Epigenomics [START_REF] Juve | Characterizing and profiling scientific workflows[END_REF]: A bioinformatics workflow that maps the epigenetic state of human cells on a genome-wide scale by processing multiple sets of genome sequences in parallel. We consider an Epigenomics instance with 577 tasks. -SoyKB [START_REF] Joshi | Next Generation Resequencing of Soybean Germplasm for Trait Discovery on XSEDE using Pegasus Workflows and iPlant Infrastructure[END_REF]: A bioinformatics workflow that re-sequences soybean germplasm lines selected for desirable traits such as oil, protein, soybean cyst nematode resistance, stress resistance, and root system architecture. We consider a SoyKB instance with 676 tasks. We profiled these workflows when executed with Pegasus [START_REF] Deelman | Pegasus, a workflow management system for science automation[END_REF], a state-of-the-art workflow management system. Pegasus monitors and logs fine-grained profiling data such as I/O operations, runtime, memory usage, and CPU utilization [START_REF] Juve | Practical resource monitoring for robust high throughput computing[END_REF].

The workflows were executed on the taurus cluster at the Grid'5000 Lyon site, which is instrumented at the node level with power meters. We used a single node to run the workflow tasks and collect power measurements (although not efficient, it allowed us to collect non-biased measurements). Each node is equipped with two 2.3GHz hexacore Intel Xeon E5-2630 CPUs, 32GB of RAM, and standard magnetic hard drives. Power measurements are collected every second from power meters (with an accuracy of 0.125 Watts5 ) that are connected to a data collector via a serial link. We are interested in identifying relationships between power consumption, task duration, CPU utilization, and volume of I/O. Detailed execution profiles (but without power/energy data) and performance analysis for both workflows can be found in [START_REF] Ferreira Da Silva | Online task resource consumption prediction for scientific workflows[END_REF].

Table 1 shows the execution profiles of Epigenomics and SoyKB tasks, with one row per task type. Since most Epigenomics tasks require 1 CPU core, power measurements were collected from a resource where only a single core was enabled (i.e., only 1 CPU slot is advertised by the resource manager). Only the pileup task requires 2 cores, but there is only one such task in the workflow. For SoyKB, many tasks require 2 CPU cores. Therefore, we collected power measurements from a resource configured with two cores. The last two columns in Table 1 show the average power consumption per task and the energy consumption to compute all tasks of that type in sequence. As power measurement were collected every second, tasks with very short runtimes (e.g., sol2sanger in the Epigenomics workflow) may not allow accurate power measurements, and are not emphasized in our upcoming analyses.

Workflow Energy Consumption Analysis

Energy-aware workflow scheduling studies [START_REF] Ghose | Energy efficient scheduling of scientific workflows in cloud environment[END_REF][START_REF] Li | Cost and energy aware scheduling algorithm for scientific workflows with deadline constraint in clouds[END_REF][START_REF] Pietri | Energy-aware workflow scheduling using frequency scaling[END_REF][START_REF] Shepherd | Workflow scheduling on power constrained vms[END_REF][START_REF] Wu | Soft error-aware energy-efficient task scheduling for workflow applications in DVFS-enabled cloud[END_REF] typically assume that the power consumed by the execution of a task at time t, P (t), is linearly related to the task's CPU utilization, u(t), as:

P (t) = (Pmax -Pmin) • u(t) • 1 n , (1)
where P max is the power consumption when the compute node is at its maximum utilization, P min is the idle power consumption (i.e., when there is no or only background activity), and n is the number of cores on the compute node. Therefore, the energy consumption of the task, E, is defined as follows:

E = r • Pmin + r 0 P (t)dt, (2)
where r denotes the task's runtime. To determine the idle power consumption P min , we collected power measurements on one node of our cluster at every second whenever no activity was performed on that node over a 2-month period (for a total of 216,000 measurements). The average idle power consumption from these measurements is 98.08W (standard deviation 1.77W). The power model in Eq. 1 does not consider the energy consumption of I/O operations, and hereafter we quantify the extent to which this omission makes the model inaccurate. Fig. 1 shows scatter plots of the power consumption versus CPU utilization for all task types of both workflows. We observe very low Pearson's correlation coefficient values between power consumption and CPU utilization (0.38 for Epigenomics, -0.02 for SoyKB). This means that no linear increase is observed in the power consumption as CPU utilization increases. For example, the align reference SoyKB task has an average CPU utilization at about 108% and consumes about 135W, while the sort sam task from that same workflow has a CPU utilization at about 124% but consumes only 102W. This difference in power consumption is mostly explained by volumes of I/O (reads and writes). analysis, Pearson's correlation coefficient values are 0.86 for Epigenomics, and 0.64 for SoyKB. These results show that power consumption is not strictly dependent, or even mainly influenced, by CPU utilization u(t) (Eq. 1), but that it depends significantly on I/O volumes. Hence, we conduct a principal component analysis (PCA) to evaluate the variance of each parameter (CPU utilization, I/O reads, and I/O writes) and their impact on the power consumption. In this analysis, we aim to understand how CPU utilization and I/O operations are influencing (positively or negatively) power consumption, and consequently quantify the weight of each parameter. From the principal components, set of values of linearly uncorrelated variables, we obverse the loadings (the weight by which each standardized original variable should be multiplied to get the component score), which contain the data variance.

Table 2 shows the principal component (PC) loadings (rotations) for each parameter. For Epigenomics, the first two PCs explain most of the variability (85.3%). All parameters present similar variance for PC1, with the I/O read and I/O write parameters dominating, while CPU utilization has greater impact on PC2. Since PC1 explains most of the variance (64.3%), the power consumption of the Epigenomics workflow is also significantly impacted by the number of I/O operations (in particular I/O reads) as shown in Fig. 2-left. Similarly, the first two PCs for SoyKB explain most of the variability (85.4%). I/O read has greater impact on PC1, while PC2 is mostly impacted by CPU utilization and I/O write. Although I/O read has significant impact on PC1, this component only explains 49% of the variance, thus I/O read has less influence on the power consumption for SoyKB (Fig. 2-right). Note that the impact of I/O read on PC2 is minimal. Overall, these results provide motivation and quantitative bases for developing a more accurate power model that captures the implications of I/O operations on power consumption in addition to that of CPU utilization.

(1,0) (1,1) (1,2) (1,3) (1,4) (1,5) (0,0) (0,1) (0,2) (0,3) (0,4) (0,5) (1,0) (1,1) (1,2) (1,3) (1,4) (1,5)

Analysis of Power and Energy Consumption for Concurrent Task Execution

The power consumption model in Eq. 1 assumes that the consumed power is simply the CPU utilization divided by the number of cores. To evaluate the validity of this assumption we collected and analyzed power measurements for solitary and concurrent workflow task executions. Since our cluster nodes are all equipped with dual, hexacore CPUs, we performed task executions with two schemes for core allocation (see Fig. 3): (1) unpaired -cores are enabled in sequence on a single socket until all cores on that socket are enabled, and then cores on the next socket are enabled in sequence; and (2) pairwise-cores are enabled in round-robin fashion across sockets (i.e., each core is enabled on a different socket than the previously enabled core). We report on results for only a subset of workflow tasks because (1) some tasks are unique; (2) some task runtimes are very short and overheads in Pegasus, such as releasing the next task, make the benefit of running these tasks in parallel negligible; or (3) energy measurements may not be accurate for tasks with very short runtimes due to the measurements interval of 1s. Finally, all our results report average runtime, power and energy measurements for concurrent executions of instances of the same task type. Epigenomics: Fig. 4 shows the average task runtime, average task power consumption, and total energy consumption (i.e., to run all 140 tasks) when running map tasks concurrently for different numbers of CPU cores. Task performance is significantly impacted when multiple cores are used within a single socket. For example, when 2 cores are enabled in different sockets (pairwise), no performance decrease is observed. However, a performance degradation of about 25% occurs when both cores are within a single socket (unpaired ). The above is due to the fact that each socket has a single L3 cache shared between its cores. q q q q q q q q q q q q q q q q q q q q q q q q Energy Consumption (kWh) Average Task Power Consumption (W) Average Task Runtime (s) While the use of multiple cores within a single socket limits performance, it consumes less power per unit of time: on the order of 10% (Fig. 4-center). According to Eq. 1, power consumption should grow linearly. Instead, we observe that power consumption is not equally divided among the number of cores per CPU. Eq. 1 thus underestimates the energy usage per unit of time-root mean squared error (RMSE) for pairwise is 10.64, and 4.92 for unpaired.

The energy profile shown in Fig. 4-right accounts for the execution of all 140 map tasks. Although power consumption is lower when using a single socket, the total energy consumption is higher due to higher task runtimes. Workflow task executions may benefit from single socket CPU usage if task runtimes are very short. In this case, the performance loss is negligible and the difference of power consumption may save energy (e.g., the filterContams task in Epigenomics). The energy consumption for the set of map tasks presents a logarithmic decrease as a function of the number of cores. This logarithmic behavior is due to the increase in power consumption. The estimation errors propagated by Eq. 1 into Eq. 2 leads to energy consumption estimation errors up to 23% (RMSEs are 0.02 for pairwise and 0.03 for unpaired ). SoyKB: Fig. 5 shows the average task runtime, the average task power consumption, and total energy consumption (i.e., to run all 500 tasks) when running haplotype caller tasks concurrently using 2 up to 8 CPU cores. Due to disk space quota on Grid'5000, we were unable to run workflow instances that used more than 8 cores concurrently. We only report on results for more than 2 cores because the workflow cannot be executed on a single core. Task runtime differences between unpaired and pairwise is minimal regardless the number of cores used. A small degradation in runtime is observed when the number of cores increase from 2 to 4. However, there is a significant performance decrease when the number of cores exceeds 4. This is because haplotype caller performs substantial I/O operations (it only has 67% of CPU utilization on average). The performance degradation is due to simultaneous I/O operations, which cause tasks to idle due to I/O resources being unavailable and/or saturated. This idle time (IOWait) is reported in the logs generated by Pegasus. q q q q q q q q q q q q q q Energy Consumption (kWh) Average Task Power Consumption (W) Average Task Runtime (s) Similar to Epigenomics, the unpaired scheme consumes slightly less power (about 5%, as see in Fig. 5-center. The power consumption estimated by Eq. 1 lies between the real-world consumption with the two schemes, with prediction errors up to 10% (RMSE up to 4.85 for pairwise). In Fig. 5-right, we see that the actual energy values are well above the estimated values (up to 22% higher). The main factor for this discrepancy is I/O operations, including the time spent waiting for I/O to complete (as indicated by IOWait values in the Pegasus logs).

Modeling and Simulating Energy Consumption of I/O-intensive Workflows

In this section, we present an augmented model for power consumption that accounts for I/O in addition to CPU utilization. This model also accounts for the number of cores and the way in which they are activated (unpaired or pairwise schemes), as well as for the time spent waiting for I/O operations to complete (IOWait).

Model

We model P (t), the power consumption of a compute node at time t, as:

P (t)=PCPU(t)+P I/O (t), (3)
where P CPU (t), resp. P I/O (t), is the power consumption due to CPU utilization, resp. I/O operations. In what follows, we detail the model for both these terms.

CPU -Let s denote the number of sockets on the compute node, and n the number of cores per socket, so that the total number of cores on the compute node is s • n. Let K the set of tasks that use at least one core on the compute node. We have:

PCPU(t)= k,i,j PCPU(k,i,j,t), (4)
q q q q q q q q q q 4 6 8 q pairwise unpaired socket 1 unpaired socket 2 q q q q q q q q q q q q q q q q 25 50 75 100 40 60 80 100 I/O-intensiveness (MB/s) Dynamic Power (W) # cores q q q q 2 3 4 5

q pairwise unpaired Fig. 6. Linear regression models. Left: power consumption increase in function of number of cores enabled for the Epigenomics map tasks. Right: dynamic power consumption vs. I/O-intensiveness for the SoyKB realign creator, indel realign, haplotype caller, and genotype gvcfs tasks.

where P CPU (k, i, j, t) is the power consumption of CPU utilization at time t due to the execution of task k (k ∈ K) on socket i (0 ≤ i < s) at core j (0 ≤ j < n).

In previous sections, we examined the impact of cores/socket allocation on power consumption in addition to CPU utilization. We have seen that the power consumption does not have constant increase as cores are enabled on sockets, and the behavior depends on the scheme used to enable further cores (pairwise or unpaired ). Fig. 6-left shows a scatter plot of power consumption increase for each additional enabled core for the map task of the Epigenomics workflow. The increase for the unpaired scheme can be approximated by linear regression with negative slope. For the pairwise scheme, an approximation by linear regression leads to nearly constant increase (noting that the RMSE is relatively high). Although this figure is for a particular task of the Epigenomics workflow, very similar results are obtained for all tasks for both production workflows considered in this work. Therefore, we derive a model that is only dependent on the task's CPU utilization and the hardware platform.

Based on the above, we now define our model for P CPU (k, i, j, t) as:

PCPU(k,i,j,t)=    (Pmax-Pmin)• u(t) s•n if j = 0 (first core on a socket)
0.881•PCPU(k,i,j-1,t) if j > 0 and pairwise 0.900•PCPU(k,i,j-1,t) if j > 0 and unpaired [START_REF] Deelman | Pegasus, a workflow management system for science automation[END_REF] where u(t) is the task's CPU utilization at time t (which can be computed by benchmarking the task on a dedicated compute node). The model is written recursively as the power consumption due to enabling a core on a socket depends on the power consumption due to previously enabled cores on that socket. The 0.881 and 0.900 constants above are obtained from the aforementioned linear regressions. Finally, note that P CPU (k, i, j, t) does not depend on i since only the rank (j) of a core in a socket matters.

I/O -Similarly to the above model for power consumption due to CPU utilization, we have:

P I/O (t)= k,i,j P I/O (k,i,j,t), (6)
where P I/O (k, i, j, t) is the power consumption of I/O operations at time t due to the execution of task k (k ∈ K) on socket i (0 ≤ i < s) at core j (0 ≤ j < n) on the compute node. Fig. 6-right shows dynamic power consumption (i.e., power consumption beyond P min ) vs. I/O-intensiveness for 4 tasks of the SoyKB workflow (realign creator, indel realign, haplotype caller, and genotype gvcfs). We define the I/Ointensiveness as the I/O volume (for reads and writes) in MB divided by the time the task spends performing solely computation (i.e., the runtime minus the time for performing and waiting for I/O operations). A higher value indicates a more I/O-intensive task, as it represents I/O overhead per second of CPU usage. We are able to compute the I/O-intensiveness of each task based on profiling data in Pegasus logs. The four task types in Fig. 6 exhibit a range of CPU utilizations, with relatively high volumes of data read/written. As for the results in Fig. 6, similar results are obtained for all tasks in the workflows we consider. We use a linear regression, shown in the figure, which has positive slope regardless of the core allocation scheme (with a steeper slope for the pairwise scheme). Based on these results, we model P I/O (k, i, j, t) as follows:

P I/O (k,i,j,t)= 0.486•(1+0.317•ω(t))•PCPU(k,i,j,t) if pairwise 0.213•(1+0.317•ω(t))•PCPU(k,i, j,t) otherwise (7)
where the 0.486 and 0.213 values above come from the linear regressions, and ω(t) is 0 if I/O resources are not saturated at time t, or 1 if they are (i.e., idle time due to IOWait). ω(t) is equal to 1 whenever the volume of I/O requests placed by concurrently running tasks exceeds some platform-dependent maximum I/O throughput. When using this model, e.g., to drive simulations of workflow task executions so as to evaluate energy-efficient workflow scheduling algorithms, it is then necessary to keep track of simulated I/O requests so as to set the ω(t) value accordingly. It turns out that, in our result, the impact of IOWait does not show any strong correlation with the features of different task types. This is why ω(t) in Eq. 7 is weighted by a single factor (0.317). We computed this factor as the average of the most accurate such factor values we computed individually for each task type. Our evaluation of the model (see Section 5.2) shows that it achieves high accuracy across task types. It is thus tempting to claim that the impact of the IOWait effect on power consumption can be captured reasonably well using a single, application-independent value for the above factor. Providing a definitive answer as to whether this claim is general would require a comprehensive set of experiments with more workflow applications running under this condition.

Experimental Evaluation

To evaluate the accuracy of our model, we extended a simulator [START_REF]WRENCH Pegasus Simulator[END_REF] of the state-of-the-art Pegasus [START_REF] Deelman | Pegasus, a workflow management system for science automation[END_REF] workflow management system (WMS), which is the WMS we used to perform the experiments described in Section 2. This simulator is built using the WRENCH simulator framework [START_REF]The WRENCH Project[END_REF], which can be used to q q q q q q q q q q q q q q q q q q q q q q q q q q map haplotype_caller indel_realing q estimation real-pairwise real-unpaired wrench-pairwise wrench-unpaired q q q q q q q q q q q q q q q q q q q q q q q q q q map haplotype_caller indel_realing build simulators of WMSs that are accurate, can run scalably on a single computer, and can be implemented with minimal software development effort [START_REF] Casanova | WRENCH: A Framework for Simulating Workflow Management Systems[END_REF]. We extended the simulator by replacing its simulation model for power consumption (the traditional model in Eq. 1) by the model proposed in Section 5.1. We provide the simulator with a description of the hardware specifications of the taurus Grid'5000 cluster and with traces from individual Epigenomics and SoyKB workflow task executions. As a result, our simulator can simulate the exact procedure used for obtaining all real-world experimental results described in previous sections, making it possible to draw direct comparisons between real-world and simulated results. The simulator code, details on the simulation calibration procedure, and experimental scenarios used in the rest of this section are all publicly available online [START_REF]WRENCH Pegasus Simulator[END_REF]. Fig. 7 shows the simulated power and energy consumption measurements as well as with the traditional model based on Eqs. 1 and 2 (shown as estimation) and with our proposed model (shown as wrench-*). Due to space constraints, we only show results for the map Epigenomics task, and the haplotype caller and indel realign SoyKB tasks. For the map tasks, the RMSE for pairwise is 4.24, and 3.49 for unpaired, which improves over the traditional model by about two orders of magnitude for the former, and a half for the later. Similarly, RMSEs for the haplotype caller tasks are 2.86 and 2.07 for pairwise and unpaired respectively, or improvements of about two orders of magnitude for both schemes. Last, RMSEs for the indel realign tasks are 0.59 for pairwise and 0.47 for unpaired, or improvements by about an order of magnitude. Predicted energy consumption based on our proposed model nearly match the actual measurements for both schemes for all task types (RMSEs 0.01).

Related Work

In the past few years, green computing has become a major topic of discussion in the scientific computing community. Many recent studies have addressed green solutions, in particular on distributed computing platforms. Research efforts in this field commonly include powering off or putting idle machines into low power states based on predictions of future workloads. On the application side, efforts are mainly focused on the optimization of resource provisioning and workload scheduling constrained by budgets and application deadlines. A recent survey [START_REF] Orgerie | A survey on techniques for improving the energy efficiency of large-scale distributed systems[END_REF] of techniques for improving energy-efficiency describes methods to evaluate and model the energy consumed by resources on distributed systems. The survey presents taxonomies of compute node and network energyaware techniques classified according to the technology employed. These techniques include adjustment of the processor's frequency and power consumption through DVFS [START_REF] Guérout | Energy-aware simulation with dvfs[END_REF], workload consolidation by running multiple tasks on the same physical machine in order to reduce the number of nodes that are powered on [START_REF] Lefevre | Towards energy aware reservation infrastructure for large-scale experimental distributed systems[END_REF], energy-aware task scheduling [START_REF] Kliazovich | Dens: data center energy-efficient networkaware scheduling[END_REF][START_REF] Wang | Energy-aware parallel task scheduling in a cluster[END_REF], virtual machine migration [START_REF] Beloglazov | A taxonomy and survey of energy-efficient data centers and cloud computing systems[END_REF][START_REF] Liu | Performance and energy modeling for live migration of virtual machines[END_REF], the coordination of network protocols [START_REF] He | Energy-efficient capture of stochastic events under periodic network coverage and coordinated sleep[END_REF], etc. These strategies often model energy consumption as a function of runtime, or do not consider the performance loss of running multiple tasks within a socket.

Several models have been developed to predict the power consumption of distributed system workloads. Most of them focus on measuring the resource utilization of distributed systems [START_REF] Enokido | An extended power consumption model for distributed applications[END_REF][START_REF] Lee | Energy efficient utilization of resources in cloud computing systems[END_REF][START_REF] Samak | Energy consumption models and predictions for large-scale systems[END_REF]. In [START_REF] Enokido | An integrated power consumption model for distributed systems[END_REF], an integrated power consumption model, which incorporates previous approaches into a single model, describes a distributed system where several clients issue requests to a central storage server. Most of these models are limited to single-core and energy consumption is related to CPU usage. A few models consider data transfers, but as a separate operation (I/O operations during task execution are not considered).

In the context of scientific workflows, several works [START_REF] Ghose | Energy efficient scheduling of scientific workflows in cloud environment[END_REF][START_REF] Li | Cost and energy aware scheduling algorithm for scientific workflows with deadline constraint in clouds[END_REF][START_REF] Pietri | Energy-aware workflow scheduling using frequency scaling[END_REF][START_REF] Shepherd | Workflow scheduling on power constrained vms[END_REF][START_REF] Wu | Soft error-aware energy-efficient task scheduling for workflow applications in DVFS-enabled cloud[END_REF] have proposed energy-aware algorithms for task scheduling or resource provisioning. These algorithms are often designed to meet energy budget or deadline constraints. Their model assumes that the total energy usage is equal to the integral of the consumed power, which is linearly related to the resource utilization. In this work, we have shown that I/O operations also have significant impact on the power consumption, and thereby the energy. To the best of our knowledge, this is the first work that profiles and analyzes power and energy consumption of real scientific workflow applications at a fine-grained level, and proposes a model that also accounts for cores/sockets allocation and I/O usage.

Conclusion and Future Work

In this work, we have profiled and analyzed the power consumption of two production scientific workflow applications executed on a distributed platform. We have investigated the impact of resource utilization and I/O operations on the energy usage, as well as the impact of executing multiple tasks concurrently on multi-socket, multi-core compute nodes. In contrast to traditional power consumption model used in the energy-efficient workflow scheduling literature, we find that power consumption is impacted non-linearly by the way in which cores in sockets are allocated to workflow tasks. Furthermore, our experimental results show that I/O operations have significant impact on power consumption. Based on these results, we proposed a power model for I/O intensive workflows that accounts for the above behaviors. Experimental evaluation of this model shows that it accurately captures real-world behavior, with order of magnitude improvement over the traditional model.

In future work, we plan to instantiate and validate our proposed model for other workflows and platform configurations. In particular, we hope to use powermetered platforms in which compute nodes have SSDs instead of HDDs. With SSDs, the impact of I/O on power consumption may exhibit different behaviors that could mandate augmenting our model. The power consumption of I/O could also be smaller relative to that of computation, but note that platforms that target extreme-scale computing also often employ low-power compute nodes (i.e., equipped with ARM processors). Another future work goal is to extend the synthetic workflow generator in [START_REF] Ferreira Da Silva | Community resources for enabling and evaluating research on scientific workflows[END_REF], which produces realistic synthetic workflow configurations based on profiles extracted from workflow execution traces. The objective is to extend the generated workflow descriptions to include data obtained from real-world power profiles that is sufficient to instantiate the power consumption model proposed in this work.

Fig. 2 Fig. 1 .Fig. 2 .

 212 Fig. 1. Task power consumption vs. CPU utilization for the Epigenomics (left) and SoyKB (right) workflows.

  0) (0,1) (0,2) (0,3) (0,4) (0,[START_REF] Deelman | Pegasus, a workflow management system for science automation[END_REF]

Fig. 3 .

 3 Fig. 3. Example of CPU core usage for the unpaired (left) and parwise (right) schemes when 6 cores are enabled.

Fig. 4 .

 4 Fig. 4. Average task runtime (left), average task power consumption (center), and energy consumption to run all 140 map tasks from Epigenomics (right). Power and energy consumption computed using Eqs. 1 and 2 are shown as estimation.

Fig. 5 .

 5 Fig.5. Average task runtime (left), average power consumption (center), and energy consumption to run all 500 haplotype caller tasks from SoyKB (right). Power and energy consumption computed using Eqs. 1 and 2 are shown as estimation.

Fig. 7 .

 7 Fig.7. Per-task power (top) and total energy (bottom) consumption measurements for the Epigenomics map task and the SoyKB haplotype caller and indel realign, as well as estimated with Eq. 1 and 2 (estimation) and our proposed model (wrench-*).

Table 1 .

 1 Execution and energy profiles of the Epigenomics (top) and SoyKB (bottom) workflow tasks. Energy measurements are for running all tasks of that type in sequence. Runtimes are shown in seconds, I/O operations in MB, and power in W. (µ is the mean, and σ the standard deviation.)

	Task	Count #cores	Runtime CPU util. I/O Read I/O Write µ σ µ σ µ σ µ σ	Power µ σ	Energy (Wh)
	fastqSplit		7	1 5.8 1.9 99.8% 0.0 508.1 173.2 254.1 86.6 126.9 5.48	1.4
	filterContams		140	1 1.2 0.2 99.1% 0.0 25.4	3.7 12.7	1.8 100.9 5.6	4.6
	sol2sanger		140	1 0.4 0.1 95.7% 0.2 66.9	9.8 29.0	4.3 98.5 3.8	1.4
	fast2bfq		140	1 0.8 0.1 97.8% 0.1 35.5	5.2	6.4	0.9 98.3 3.6	2.9
	map		140	1 57.9 5.0 99.9% 0.0 437.9	2.4	2.6	0.6 126.8 0.9	285.7
	mapMerge		8	1 5.9 6.9 99.5% 0.0 171.2 205.6 84.0 103.4 113.5 7.7	1.5
	maqIndex		1	1 33.5	-99.9% -511.7	-338.3	-125.1	-	1.2
	pileup		1	2 38.4	-80.8% -559.3	-264.1	-135.5	-	1.4
	Task		Count #cores	Runtime µ σ	CPU util. I/O Read I/O Write Power Energy µ σ µ σ µ σ µ σ (Wh)
	align reference		25	2		1.8 0.0 53.9% 0.0 2609.7 0.0 186.6 0.01 134.8 4.7	1.6
	sort sam		25	2		1.3 0.1 61.9% 0.0 901.5 0.0 187.2 1.6 101.7 1.9	0.9
	dedup		25	2		2.0 0.0 60.7% 0.0 901.9 0.0 186.9 0.2 106.2 4.3	1.5
	add replace		25	2		1.3 0.0 62.0% 0.0 901.5 0.0 186.9 0.0 102.6 1.7	0.9
	realign creator		25	2	133.1 2.6 75.9% 0.0 3230.8 8.7 189.6 2.8 135.3 0.3 125.1
	indel realign		25	1		34.3 0.0 18.9% 0.0 953.8 5.8 187.0 0.0 123.2 0.6	25.9
	haplotype caller	500	1		79.3 6.9 66.7% 0.0 1149.8 24.2 186.9 0.0 130.8 1.0 1329.5
	genotype gvcfs		20	1	263.8 29.6 95.9% 0.0 1058.0 16.2 187.6 0.1 126.6 0.3 185.5
	comb variants		1	1		35.5	-26.5% -958.0	-186.9	-108.9 -	1.1
	variants indel		1	2		48.6	-23.7% -1699.5	-454.4	-114.0 -	1.5
	filtering indel		1	1		34.7	-20.3% -955.2	-186.9	-109.1 -	1.0
	variants snp		1	2		48.6	-23.2% -1699.5	-454.4	-115.4 -	1.5
	filtering snp		1	2		34.7	-10.2% -955.3	-186.9	-109.6 -	1.0
	merge gcvf		1	1 46804.5	-99.9% -3061.2	-238.8	-128.9 -1675.3

Table 2 .

 2 Principal component (PC) loadings (rotations) for the Epigenomics and SoyKB workflows.

	Parameter	PC1	Epigenomics PC2	PC3	PC1	SoyKB PC2	PC3
	CPU Utilization	0.53	0.84	-0.03	-0.55	-0.62	0.56
	I/O Read	0.59	-0.35	0.71	-0.73	0.04	-0.67
	I/O Write	0.59	-0.40	-0.69	-0.39	0.78	0.47

Manufactured by OMEGAWATT: http://www.omegawatt.fr/gb/index.php

Acknowledgments. This work is funded by NSF contracts #1642369 and #1642335, "SI2-SSE: WRENCH: A Simulation Workbench for Scientific Workflow Users, Developers, and Researchers", and CNRS under grant #PICS07239; and partly funded by NSF contract #1664162, and DOE contract #DE-SC0012636.