Abdelhamid Maali

Bharat Bhushan

Measurement of Slip Length on Superhydrophobic Surfaces

Keywords: fluid slip, atomic force microscopy, surface force apparatus, rheometry, particle image velocimetry, pressure drop

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Liquid flow close to the interface has received much attention recently, and it has been studied theoretically and experimentally by several groups. Using the newly developed techniques, it is actually possible to drive and probe the flow at the micro/nanoscales. For a smooth hydrophobic surface, there is an agreement that the fluid flow is subject to a slip close to the boundary [START_REF] Churaev | Slippage of liquids over lyophobic solidsurfaces[END_REF][START_REF] Lauga | Microfluidics: the no-slip boundary condition[END_REF][START_REF] Neto | Boundary slip in Newtonian liquids: a review of experimental studies[END_REF][START_REF] Maali | Nanorheology and boundary slip in confined liquids using atomic force microscopy[END_REF][START_REF] Bocquet | Nanofluidics, from bulk to interfaces[END_REF]. The amount of fluid slip is commonly represented by the slip b, which is typically of the order of tens of nanometres. Very large slip length has been reported on superhydrophobic surfaces [START_REF] Bocquet | Flow boundary conditions from nano-to micro-scales[END_REF][START_REF] Wang | Boundary slip and nanobubble study in micro/nanofluidics using atomic force microscopy[END_REF][START_REF] Rothstein | Slip on superhydrophobic surfaces[END_REF]. Large slip length indicates that the liquid flow between confining surfaces experiences a lower friction, and a solid moving in a liquid experiences a lower drag force. Therefore, to create surfaces that increase slip length becomes a hot topic in micro/nanofluidics devices.

In this paper, various techniques used to measure drag and slip length on micro/nanostructured superhydrophobic surfaces with large slip length are reviewed. Some theoretical models and calculated expressions for the effective slip length on the superhydrophobic surface are presented first. Next, the different techniques to measure the slip length are presented. This review is divided into *Authors for correspondence (a.maali@loma.u-bordeaux1.fr; bhushan.2@osu.edu).

One contribution of 10 a Theme Issue 'Biosensors: surface structures and materials'. Adapted from [START_REF] Cottin-Bizonne | Dynamics of simple liquids at heterogeneous surfaces: molecular-dynamics simulations and hydrodynamic description[END_REF]. (Online version in colour.) three sections. In §2, models and theoretical calculations of the effective slip length on structured superhydrophobic surface are presented. In §3, various techniques to measure large slip length, which include rheometry, particle image velocimetry, pressure-drop, surface force apparatus (SFA) and atomic force microscope (AFM) measurements, and some representative results are presented. In §4, a summary and outlook are provided.

Calculation of the effective slip length

As mentioned earlier, several groups have theoretically investigated liquid flow on superhydrophobic surfaces. They have calculated the effective slip length in different situations. Here, the simplest models are recalled, and their results are presented.

A planar surface with one-dimensional features, such as long grooves, usually aligned or perpendicular parallel to the flow directions has been studied by several authors. Philip [START_REF] Philip | Flows satisfying mixed no-slip and no-shear conditions[END_REF][START_REF] Philip | Integral properties of flows satisfying mixed no-slip and no-shear conditions[END_REF] studied the flow in a pipe of radius R with stripes parallel to the flow. In such situations, the pipe surface is patterned with a periodic no-slip and slip region. When R → ∞, the effective slip length is given by

b Philip = L p ln cos pf g 2 , ( 2.1)
where L is the pitch of the pattern and f g is the gas fraction area (fraction of the surface where the slip is infinity). f g is equal to (La)/L, where a is the stripe width with zero slip length. For similar surfaces, but with a pattern perpendicular to the flow, Lauga & Stone [START_REF] Lauga | Effective slip in pressure-driven stokes flow[END_REF] calculated the slip length, and when the radius of the pipe R → ∞, their expression is reduced to Cottin-Bizonne et al. [START_REF] Cottin-Bizonne | Dynamics of simple liquids at heterogeneous surfaces: molecular-dynamics simulations and hydrodynamic description[END_REF] conducted semi-analytical calculations to calculate the effective slip length on microstructured superhydrophobic surfaces. A schematic of their situation is presented in figure 1. They solved the Navier-Stokes equation and calculated effective slip length defined as b = V s / ġ, where V S is the effective velocity slip on the surface and ġ is the shear rate that is imposed at an infinite distance from the surface. They considered first a surface composed of a succession of stripes of width a characterized by a slip length (b 1 = 0) and stripes of width La having an infinite slip length (b 0 → ∞). The results are presented in figure 2 for the case of flow parallel to the stripes (figure 2a) and perpendicular to the stripes (figure 2b). As shown in figure 2, their approach is in excellent agreement with the analytical results of Philip [START_REF] Philip | Flows satisfying mixed no-slip and no-shear conditions[END_REF] (flow parallel to stripes) and Lauga & Stone [START_REF] Lauga | Effective slip in pressure-driven stokes flow[END_REF] (flow perpendicular to stripes). They also studied two other particular cases: one was a surface composed of alternating stripes of no-slip (b 1 = 0) and partial-slip length (b 0 ) and the second was a surface composed of alternating stripes of infinite-slip (b 0 → ∞) and partial-slip length (b 1 ). Figure 3a A surface with two-dimensional patterning, usually composed of a series of posts distributed on a regular lattice as illustrated in figure 4a, is expected to display lower friction than a surface with one-dimensional patterning. Ybert et al. [START_REF] Ybert | Achieving large slip with superhydrophobic surfaces: scaling laws for generic geometries[END_REF] were the first to calculate an expression that predicted the value of the effective slip length for the given surface. Their results are based on semianalytical calculations, as mentioned above. They computed the effective slip length as a function of the solid fraction (surface area where the slip length vanishes). Then, a linear regression was performed on the numerical calculated data that allows extraction of a useful and very simple expression for the effective slip length b eff as function of the solid fraction f s ,

b Lauga-Stone = L 2p ln cos pf g 2 . ( 2
b eff = L 0.325 √ f s -0.44 . (2.3)
In figure 4b, the computed effective slip length divided by L, where f s = 1f g = a/L, is plotted for different values of the solid fractions f s . Davis & Lauga [START_REF] Davis | Hydrodynamic friction of fakir-like superhydrophobic surfaces[END_REF] conducted analytical calculations to derive an expression of the effective slip length on the superhydrophobic surface composed of posts suitably distributed on a surface. For posts on a square lattice and in the limit of small coverage by the post (small solid fraction area), the value of the slip length can be obtained analytically. Their expression predicts that the relation between the slip length and the solid fraction is given by

b eff = L 3 16 p f s - 3 2p ln(1 + √ 2) ≈ L 0.332 √ f s -0.421 . (2.4)
Notice here the good agreement between the expression derived by Davis & Lauga [START_REF] Davis | Hydrodynamic friction of fakir-like superhydrophobic surfaces[END_REF] using analytical calculations and the expression derived by Ybert et al. [START_REF] Ybert | Achieving large slip with superhydrophobic surfaces: scaling laws for generic geometries[END_REF] using numerical calculations.

Measurement methods and results

(a) Rheometry experiments

In such experiments, the slip length was obtained through torque measurements with a rheometer [START_REF] Lee | Structured surfaces for a giant liquid slip[END_REF][START_REF] Lee | Maximizing the giant liquid slip on superhydrophobic microstructures by nanostructuring their sidewalls[END_REF][START_REF] Lee | Influence of surface hierarchy of superhydrophobic surfaces on liquid slip[END_REF][START_REF] Ming | Fluid drag reduction on superhydrophobic surfaces coated with carbon nanotube forests (CNTs)[END_REF]. In figure 5, the schematic of a cone-plate rheometer is shown; this is the most popular geometry because it produces a uniform shear rate over the sample. When a cone of radius R with a very small angle q 0 rotates at angular velocity U, the slip length is related to the measured torque by

M = 2phUR 3 3q 0 1 - 3b eff 2Rq 0 + 3b 2 eff R 2 q 2 0 - 3b 3 eff R 3 q 3 0 ln Rq 0 + b eff b eff , (3.1)
where h is the viscosity of the liquid. Lee et al. [START_REF] Lee | Structured surfaces for a giant liquid slip[END_REF] and Lee & Kim [START_REF] Lee | Maximizing the giant liquid slip on superhydrophobic microstructures by nanostructuring their sidewalls[END_REF] used a commercial rheometer (AR 2000, TA instruments) whose operational torque range was between 0.1 and 200 mN m, with an angular velocity between 10 -8 and 300 rad s -1 . They studied surfaces presenting well-regulated post and grate structures. The samples were created by photolithography and deep reactive etching on a silicon wafer. The grates were designed to be concentric so that they were parallel to the liquid flow in the rheometer system, and the posts had a circular cross section and were distributed in a square lattice. They studied the effective slip length versus the gas fraction (solid fraction) and also the value of the effective slip length versus the pitch value. In figure 6, the measured effective slip length for a fixed value of the gas . The effect of the pitch on the slip length with a gas fraction fixed at 98% [START_REF] Lee | Structured surfaces for a giant liquid slip[END_REF]. The solid line represents theory for grates [START_REF] Lauga | Effective slip in pressure-driven stokes flow[END_REF]; dashed line, theory for posts [START_REF] Ybert | Achieving large slip with superhydrophobic surfaces: scaling laws for generic geometries[END_REF]; triangles, experimental data for grates [START_REF] Lee | Structured surfaces for a giant liquid slip[END_REF]; circles, experimental data for posts [START_REF] Lee | Structured surfaces for a giant liquid slip[END_REF]. Reproduced with permission. (Online version in colour.) fraction (98%) is shown, but with varying pitches of 20, 50 and 60 mm for a post and 50, 200 and 250 mm for grates. The depth for each sample was maintained to be the same as the varying pitch. The slip length increases linearly with the pitch for both posts and grates, which is in agreement with the theoretical predictions. This result illustrates that as the pitch is larger, the produced slip length is larger with a limitation fixed by the condition for the wetting transition.

Bocquet et al. [START_REF] Bocquet | Comment on 'Large slip of aqueous liquid flow over a nanoengineered superhydrophobic surface[END_REF] reported some questions about the accuracy of the measurements of the slip length using a rheometer. They showed that the uncertainty in the measurements can be of the same order of magnitude as the slip length.

(b) Particle image velocimetry experiments

In order to measure the flow in the vicinity of such surfaces, particle image velocimetry has been used by several authors [START_REF] Ou | Direct velocity measurements of the flow past drag-reducing ultrahydrophobic surfaces[END_REF][START_REF] Joseph | Slippage of water past superhydrophobic carbon nanotube forests in microchannels[END_REF][START_REF] Truesdell | Drag reduction on a patterned superhydrophobic surface[END_REF]. In such techniques, the velocity profile of the liquid can be measured with high accuracy, and the slip length can be extracted from the equation [24] b = v s /vv s /vz, where v s and vv s /vz are the velocity and the shear rate close to the surface.

The common set-up used in particle image velocimetry is based on the measurement of the velocity of fluorescent particles embedded in water. An example of these set-ups is used by Joseph et al. [START_REF] Joseph | Slippage of water past superhydrophobic carbon nanotube forests in microchannels[END_REF] to study water flow on a superhydrophobic surface made of a carbon nanotube (CNT) forest. These CNTs were embedded in a microfluidic set-up with a microchannel, as shown in figure 7. A pressure-driven flow is induced by controlling inlet and outlet pressure. Usually, the observation was made using an objective with a large numerical aperture (NA = 1.3) to achieve a narrow depth of field. A scan in the vertical direction allows a change of the position of the focal plane and measurement of the velocity profile v(z) in the whole microchannel. The focal plane was controlled with a piezo. [START_REF] Joseph | Slippage of water past superhydrophobic carbon nanotube forests in microchannels[END_REF]. A pressuredriven flow is conducted in a microchannel made from glass/carbon nanotube forest surfaces. The focal plane Z is controlled with a piezo. Velocimetry is measured using particle image velocimetry over the whole channel depth. Reproduced with permission.

On CNT forest surfaces, Joseph et al. [START_REF] Joseph | Slippage of water past superhydrophobic carbon nanotube forests in microchannels[END_REF] measured the velocity profile and extracted the slip length by fitting the data using the predicted Poiseuille profile. They studied flow on CNT surfaces in the Cassie state (trapped air pocket) and also surfaces in the Wenzel state where the CNT surfaces are impregnated. (Descriptions of various wetting regimes can be found in Nosonovsky & Bhushan [START_REF] Nosonovsky | Multiscale dissipative mechanisms and hierarchical surfaces: friction, superhydrophobicity, and biomimetics[END_REF].) In the latter surface, the effective solid-liquid interface is fixed at the mean position between the top and the bottom of the CNT. The characteristic lateral scale of these surfaces is related to the nanotube interdistance, which can be tuned during the deposition process of the nanotube on the substrate. The slip length was measured for a surface having a lateral characteristic length between 1 and 6 mm. For the surfaces prepared in the Cassie state, the slip length increases with increasing length scale L, as shown in figure 8. However, for the surfaces that are in the Wenzel state, the slippage is strongly reduced owing to the absence of trapped air on the surface and thus increases of liquid-solid friction. The error bars on figure 8 correspond to dispersion in the measurements.

Ou & Rothstein [START_REF] Ou | Direct velocity measurements of the flow past drag-reducing ultrahydrophobic surfaces[END_REF] used particle imaging on superhydrophobic surfaces fabricated from silicon wafers using photolithography and designed to incorporate precise patterns of micrometre-sized ridges aligned in the flow direction. A schematic of the surface is presented in figure 9a. In figure 9b, the velocity profile across two different microchannels is shown, both with H = 80 mm height. For one of them, the superhydrophobic surface consists of w = 30 mm wide microridges spaced at d = 30 mm, and for the second surface, w = 30 mm and d = 60 mm. As given by the boundary conditions, the velocity profile on the microridges is zero (no-slip). Above the microridges, the velocity increases to a maximum near the centre line and returns to zero at the opposite surface, which is a smooth surface with no-slip boundary conditions. However, when the measurement was on the shear-free air-water interface (between the ridges), the velocity is not equal to zero, and it is maximum for measurements at the centre of the shear-free interface. For the surface having microridges (w = 30 mm and d = 60 mm), the slip length extracted from the measurement is found to be equal to b = 7.5 mm. The measurements made are in agreement with their simulations, as shown in figure 9b.

In summary, particle image velocimetry is a direct technique that studies liquid flow close to solid surfaces. It allows direct access to the velocity profile and extraction of the slip length with high accuracy. It was used to measure the effective slip length on different superhydrophobic surfaces.

(c) Pressure-drop experiments

In such experiments, the slip length is extracted from the measurements of pressure drop as a function of the flow rate in a microchannel [START_REF] Ou | Direct velocity measurements of the flow past drag-reducing ultrahydrophobic surfaces[END_REF][START_REF] Ou | Laminar drag reduction in microchannels using ultrahydrophobic surfaces[END_REF]. Using the Navier boundary conditions, the volume flow rate per unit length q of the fluid between two infinite parallel plates separated by a height H is given by

q = H 3 4h - dp dz 1 3 + b H + b , ( 3.2)
where (dp/dz) is the pressure gradient and b is the effective slip length (in the above expression, one of the surfaces is supposed to be superhydrophobic).

Ou & Rothstein [START_REF] Ou | Direct velocity measurements of the flow past drag-reducing ultrahydrophobic surfaces[END_REF] showed that the presence of an air-water interface on superhydrophobic surfaces induces less resistance on the flow compared with noslip surfaces. They measured the drag reduction coefficient, where Dp SH is the experimentally measured pressure drop on the superhydrophobic surfaces and Dp no-slip is the pressure drop for flow over a no-slip surface at the same flow rate.

D R = Dp no-slip -Dp SH Dp no-slip , (3.3)
An important consequence of the surface slip is the reduction of drag or pressure drop, which is directly measured by the flow-rate measurement in the experiment. They achieved, in some cases, drag reduction D R >40%, which corresponds to a slip length larger than 20 mm. They have also shown that the drag reduction increased with increasing shear-free area. Figure 10 shows the drag reduction versus the fraction of the free-shear area.

(d) Surface force apparatus experiments

Steinberger et al. [START_REF] Steinberger | High friction on a bubble mattress[END_REF] used a homemade SFA to study the liquid flow on superhydrophobic surfaces. They used a microstructured surface embedded with a square lattice of calibrated cylindrical holes. The plain microstructured surface was a hydrophilic surface wetted by water in the Wenzel regime where the holes were filled by water. Silanized microstructured surfaces are superhydrophobic surfaces, and the wetting state is the Cassie regime with air bubbles trapped in the cylindrical holes. They studied the flow of a mixture water and glycerol (viscosity h = 39 ± 2 mPa s) using the dynamic SFA. The dynamic SFA measures the real G (u) and the imaginary part G (u) of the complex force acting on the surfaces when the sphere of the SFA is vibrated in the normal direction with a frequency u/2p.

The imaginary part of the complex force response G (u) is the viscous damping of the liquid flow, and it depends linearly on the gap for a distance larger than the slip length. The linear extrapolation of the inverse of the viscous damping G -1 (u) intercepts the distance axis at a position given by the effective slip length. The real part G (u) is related to the elastic interaction owing to the confined liquid.

For the surface wetted in the Wenzel state, some of the liquid inside the holes contributes to the flow. The real no-slip plane lies underneath the top of the wall, and thus the measured G -1 (u) gives an effective slip length of b = 105 nm (figure 11a). The real part G (u) is equal to zero for all distances, as there is no elastic deformation of the surface (figure 11b).

For the surface wetted in the Cassie state, the presence of air bubbles trapped in the holes reduces the value of the effective slip length, as shown in figure 11c. The asymptote of the inverse of the viscous damping intercepts the distance axis at a distance of b = 20 nm. The real part G (u) of the response is not zero, and it shows the contribution of the bubbles that deform under the pressure induced by the flow (figure 11d). This elastic response confirms the presence of the bubbles in the holes.

Steinberger et al. [START_REF] Steinberger | High friction on a bubble mattress[END_REF][START_REF] Steinberger | Nanoscale flow on a bubble mattress: effect of surface elasticity[END_REF] and Hyväluoma & Harting [START_REF] Hyväluoma | Slip flow over structured surfaces with entrapped microbubbles[END_REF] studied the influence of the shape of the meniscus (contact angle of the bubbles) on the slip length. Their numerical calculations showed a decrease of the effective slip length owing to the meniscus curvature (figure 12). A huge decrease of the effective slip length is obtained for menisci that correspond to a contact angle q > 45 • . The obtained results are in agreement with the predictions of Richardson [START_REF] Richardson | On the no-slip boundary condition[END_REF], which shows that the boundary conditions for a shear-free interface become no-slip if the surface is sufficiently rough.

In summary, this SFA experiment pointed out the role of the liquid-gas menisci on the boundary conditions of the liquid flow. The gas trapped at the solid surface can act as an anti-lubricant and promote high friction. These results cast some doubts on whether the presence of nanobubbles can be responsible for the large slip length on a smooth hydrophobic surface. Bhushan et al. [START_REF] Bhushan | Boundary slip study on hydrophilic, hydrophobic, and superhydrophobic surfaces with dynamic atomic force microscopy[END_REF] have used a dynamic mode AFM experiment to measure slip length on hydrophobic and superhydrophobic surfaces (also see Wang & Bhushan [START_REF] Wang | Boundary slip and nanobubble study in micro/nanofluidics using atomic force microscopy[END_REF]). Their experiments were calibrated using a hydrophilic surface (mica) where the no-slip boundary conditions were well known. In their experiment, a glass sphere of 42.4 ± 0.8 mm diameter was glued to the end of a cantilever (figure 13a). The substrate sample was approached to the sphere at very low velocity, and the deflection signal was monitored as an indicator to determine a hard contact position. They measured, using a lock-in-amplifier, the amplitude and phase shift data as a function of piezo displacement during the cantilever approach to the hydrophilic, hydrophobic and superhydrophobic surfaces.

From the amplitude and phase, they extracted the total damping coefficient as a function of separation distance between the sphere and the sample surfaces [START_REF] Jai | Analytical description of the motion of an acoustic-driven atomic force microscope cantilever in liquid[END_REF],

1 + g H g 0 = - A 0 A sin(4) Q 0 1 + Q 2 0 + 2(A 0 /A) cos(4) + A 2 0 /A 2 1 + Q 2 0 , (3.4)
where A 0 is the free oscillation amplitude of the cantilever far from the surface, while A is the amplitude of the cantilever at a given position, Q 0 is the quality factor of the cantilever far from the surface, 4 is the phase shift at a given position and g 0 is the bulk viscous damping coefficient far from the substrate surfaces. The damping coefficients versus distance are shown in figure 13b. For a given separation distance, the total damping coefficient decreases with the sequence of hydrophilic, hydrophobic and superhydrophobic surfaces, indicating an increasing degree of boundary slip with increasing hydrophobicity. For each surface, the calculated damping coefficient increases with decreasing separation distance.

In the analysis of the experiment of Bhushan et al. [START_REF] Bhushan | Boundary slip study on hydrophilic, hydrophobic, and superhydrophobic surfaces with dynamic atomic force microscopy[END_REF] to eliminate the influence of surface roughness of hydrophobic and superhydrophobic surface on boundary slip, they took the mean surface as a virtual plane where the solid-liquid interface is located.

To evaluate the degree of the slip length at the solid-liquid interfaces of the three surfaces, the measured hydrodynamic damping coefficients were fitted using

g H = 6pR 2 h D f * , (3.5a)
where f * is the correction function that takes into account the boundary slip on the surfaces, given by Vinogradova [START_REF] Vinogradova | Drainage of a thin liquid-film confined between hydrophobic surfaces[END_REF],

f * = 1 4 1 + 6D 4b 1 + D 4b ln 1 + 4b D -1 , (3.5b)
where D is the closest separation distance between the sphere and the solid surface, h is the viscosity of the liquid and g H is the hydrodynamic damping coefficient. Slip lengths of about 43 and 236 nm are obtained on hydrophobic and superhydrophobic surfaces, which indicate that the slip length increases with increasing hydrophobicity. In this experiment [START_REF] Bhushan | Boundary slip study on hydrophilic, hydrophobic, and superhydrophobic surfaces with dynamic atomic force microscopy[END_REF], the authors have assumed that the boundary slip is described by the Vinogradova equation shown above. More recently, Asmolov et al. [START_REF] Asmolov | Drag force on a sphere moving toward an anisotropic superhydrophobic plane[END_REF] showed that on textured superhydrophobic surfaces, the effective slip length is not independent of the distance in the drainage experiments. Asmolov et al. [START_REF] Asmolov | Drag force on a sphere moving toward an anisotropic superhydrophobic plane[END_REF] suggested a new expression to evaluate the slip length on textured superhydrophobic surfaces, and it will be useful to check such expressions in forthcoming AFM experiments.

Summary and outlook

This paper presents a review of different techniques for measurement of the large slip length on superhydrophobic surfaces. First, the theoretical models used to calculate the effective slip length on superhydrophobic surfaces in different configurations of liquid flow are presented. Then, details of different experimental techniques are presented that are used to measure the slip on these superhydrophobic surfaces.

Particle image velocimetry is a direct measurement of the slip length since it measures the velocity profile of the fluid in the vicinity of the superhydrophobic surfaces. The measurements in this case are also more local and can probe the liquid flow in different parts of the surface; on the shear-free part (liquidgas interface) and also on the low-slip part (solid-liquid interface). The other techniques described in this review are indirect measurements. In these cases, the slip length is extracted from torque measurements (rheometry experiments), flow rate versus pressure-drop measurements (pressure-drop experiments) and from damping versus distance (SFA and AFM experiments). The rheometry and pressure-drop experiments are useful only for a large slip length due to the low resolution of such techniques. Furthermore, the measured flow is an average flow on the whole surface and on the whole gap of the microchannel. The SFA also provides an average value of the flow; however, it has the benefit of probing the flow at the nanoscale by changing the gap between the confining surfaces. The AFM experiments can provide local measurements of the boundary conditions by changing the probed position on the surfaces. They also have the benefit of controlling the gap of the flow (from a few of micrometres to nanometres), and produce measurements with high accuracy down to a few nanometres.

Although the field of slip of the liquid flow on a superhydrophobic surface has been extensively studied, as we have shown in this review; however, many challenges remain. It has been shown by the SFA experiments that the liquid-gas menisci on superhydrophobic surfaces play an important role on the boundary conditions of the liquid flow. The gas trapped at solid surfaces can act as an anti-lubricant and promote high friction. Therefore, controlling the shape of the meniscus will allow the tuning of the boundary conditions at the surface. This may be done by use of the electrowetting process, which is a promising technique [START_REF] Bhushan | Integrating electrowetting into micromanipulation of liquid droplets[END_REF][START_REF] Bhushan | Role of electric field on surface wetting of polystyrene surface[END_REF].

Liquid flow at turbulent flow (high Reynolds numbers) is another interesting topic to be examined. For laminar flow, the reported slip was independent of flow rate; however, for turbulent flow, recent results demonstrated that the slip increased with increasing Reynolds number, and substantial reduction of friction coefficient was also reported [START_REF] Daniello | Drag reduction in turbulent flows over superhydrophobic surfaces[END_REF]. The upper value of the slip and the maximum drag reduction is still unknown in the high Reynolds number regime. The role of the curvature of the shear-free interface on the slip in turbulent flow needs to be studied.

The authors thank Yunlu Pan for critical reading of the paper.

Figure 1 .

 1 Figure 1. Schematic of the microscopic hydrodynamic boundary condition, alternating stripes of slip length b 0 and slip length b 1 . Adapted from [12]. (Online version in colour.)

Figure 2 .

 2 Figure 2. The effective slip length b as a function of the gas fraction area f g = (La)/L (fraction of the surface where the slip is infinity). The curves represent the analytical calculations and the circles the semi-analytical calculation: (a) flow perpendicular to the stripes and (b) flow parallel to the stripes. Reproduced with permission from [12].

Figure 3 .

 3 Figure 3. (a) Effective slip length b on a surface composed of alternating stripes with slip length b 0 and vanishing slip length (b 1 = 0) as a function of the slip length b 0 (the gas fraction area is f g = 97%) and (b) effective slip length b on a surface composed of alternating stripes with slip length b 1 and infinite slip length b 0 → ∞. The semi-analytical calculation was carried out for different values of gas fraction area f g . The short-dashed lines are a guide to the eye. In the inset, the same points are plotted on a linear scale. The dashed lines are linear fits to the measured values of b. The shear flow in both cases is parallel to the stripes. Reproduced with permission from [12]. Circles, f g = 0.833; squares, f g = 0.5; triangles, f g = 0.25; inverted triangles, f g = 0.1. (Online version in colour.)

  ,b shows the evolution of the macroscopic slip length b as a function of the microscopic slip length b 0 and b 1 , respectively.

Figure 4 .

 4 Figure 4. (a) Schematic of the liquid interface in a surface with two-dimensional patterning with posts organized on a square lattice. (b) Normalized effective slip length b eff /L for a composite surface plotted against 1/ √ f s . The solid line represents the linear regression (b eff /L = 0.325/ √ f s -0.44). Reproduced with permission from [13]. (Online version in colour.)

Figure 5 .

 5 Figure 5. Schematic of a cone-and-plate rheometer, showing geometric parameters and reference frames. Reproduced with permission from [19].

Figure 6

 6 Figure 6. The effect of the pitch on the slip length with a gas fraction fixed at 98%[START_REF] Lee | Structured surfaces for a giant liquid slip[END_REF]. The solid line represents theory for grates[START_REF] Lauga | Effective slip in pressure-driven stokes flow[END_REF]; dashed line, theory for posts[START_REF] Ybert | Achieving large slip with superhydrophobic surfaces: scaling laws for generic geometries[END_REF]; triangles, experimental data for grates[START_REF] Lee | Structured surfaces for a giant liquid slip[END_REF]; circles, experimental data for posts[START_REF] Lee | Structured surfaces for a giant liquid slip[END_REF]. Reproduced with permission. (Online version in colour.)

Figure 7 .

 7 Figure 7. Experimental particle image velocimetry set-up used by Joseph et al.[START_REF] Joseph | Slippage of water past superhydrophobic carbon nanotube forests in microchannels[END_REF]. A pressuredriven flow is conducted in a microchannel made from glass/carbon nanotube forest surfaces. The focal plane Z is controlled with a piezo. Velocimetry is measured using particle image velocimetry over the whole channel depth. Reproduced with permission.

Figure 8 .

 8 Figure 8. Measured slip length on the CNT surface as a function of the characteristic lateral length scale L of the underlying CNT forests. Circles, surfaces in the Cassie state; squares, surfaces in the Wenzel state. The dotted line is the theoretical prediction for the slip length [22]. Reproduced with permission. (Online version in colour.)

g la s s to p s u rf a c e fl o w d ir e c ti oFigure 9 .

 9 Figure 9. (a) Schematic of the superhydrophobic surface to perform particle image velocimetry for measuring the slip length of water on these surfaces and (b) measurements of the velocity profile and comparison with numerical simulations for the flow on the microchannel. The data are for measurements for the vertical profile above the centre of the microridges (triangles), above the centre of the free-shear interface for the surface having microridges spaced by d = 30 mm (filled squares) and above the centre of the shear-free interface for the surface d = 60 mm (open squares).Reproduced with permission from[START_REF] Ou | Direct velocity measurements of the flow past drag-reducing ultrahydrophobic surfaces[END_REF].

Figure 10 .

 10 Figure 10. The drag reduction as a function of the fraction of the shear-free interface area. The solid line is a prediction of fluid dynamics simulations for a surface of 20 mm microridges for various space widths [8]. Reproduced with permission. (Online version in colour.)

1 )Figure 11 .

 111 Figure 11. Flow between a Pyrex sphere and a textured microstructured plane (see insets). In (a,b), the surface is hydrophilic; the liquid fills the holes. In (c,d), the surface is superhydrophobic and the liquid does not fill the holes. (a,c) Show the evolution of the inverse of the viscous damping G -1 (u) as a function of the distance between the surfaces. The far-field asymptote (dashed line) allows the determination of the slip length. (b,d) Show the evolution of the elastic part G (u) as a function of the distance. Reproduced with permission from [27]. (Online version in colour.)

Figure 12 .

 12 Figure 12. Evolution of the slip length with a menisci shape. The dashed line corresponds to a numerical value in the hydrophilic case when the liquid fills the holes. Reproduced with permission from[START_REF] Steinberger | High friction on a bubble mattress[END_REF]. (Online version in colour.)

Figure 13 .

 13 Figure 13. (a) A schematic of an oscillating sphere approaching a surface and profiles of the velocity amplitude of fluid flow with and without boundary slip. The arrows above and below the solid surface represent directions for fluid flow and relative movement of the sample surface to the sphere, respectively. (b) Total damping coefficient as a function of separation distance between the sphere and sample surfaces and corresponding fitting curves on the hydrophilic, hydrophobic and superhydrophobic surfaces. Adapted from [31]. (Online version in colour.)


Ahmed Boubrima

Azzedine Boukerche

Walid Bechkit

Hervé Rivano

WSN Scheduling for Energy-Efficient Correction of Environmental Modelling

Keywords: Wireless sensor networks (WSN), activity scheduling, lifetime maximization, environmental modelling, data assimilation

Wireless sensor networks (WSN) are widely used in environmental applications where the aim is to sense a physical parameter such as temperature, humidity, air pollution, etc. Most existing WSN-based environmental monitoring systems use data interpolation based on sensor measurements in order to construct the spatiotemporal field of physical parameters. However, these fields can be also approximated using physical models which simulate the dynamics of physical phenomena. In this paper, we focus on the use of wireless sensor networks for the aim of correcting the physical model errors rather than interpolating sensor measurements. We tackle the activity scheduling problem and design an optimization model and a heuristic algorithm in order to select the sensor nodes that should be turned off to extend the lifetime of the network. Our approach is based on data assimilation which allows us to use both measurements and the physical model outputs in the estimation of the spatiotemporal field. We evaluate our approach in the context of air pollution monitoring while using a dataset from the Lyon city, France and considering the characteristics of a monitoring system developed in our lab. We analyze the impact of the nodes' characteristics on the network lifetime and derive guidelines on the optimal scheduling of air pollution sensors.

I. INTRODUCTION

Wireless sensor networks (WSN) play an important role in the infrastructure of the Internet of Things (IoT) by connecting the physical world to the Internet through real time measurements which are carried out by a set of sensors placed in a deployment field and connected to collector nodes called sinks [START_REF] Yick | Wireless sensor network survey[END_REF]. The applications of wireless sensor networks can be grouped into two main classes: 1) target detection applications where sensors have a sensing range allowing them to capture the variations of target objects within a certain distance (like acoustic and video sensors); And 2) environmental applications (temperature or air pollution monitoring for instance [START_REF] Boubrima | Error-bounded air quality mapping using wireless sensor networks[END_REF] [START_REF]Optimal wsn deployment models for air pollution monitoring[END_REF]) where data is available only at the points where sensors are deployed and the objective is to construct a spatiotemporal field (or set of maps) with a good estimation precision [START_REF] Du | Sensor placement and measurement of wind for water quality studies in urban reservoirs[END_REF].

Most existing WSN-based environmental monitoring applications use interpolation or regression techniques to construct the spatiotemporal field of the physical parameters while focusing only on the measurements gathered by sensors [START_REF] Roy | Spatio-temporal sensor management for environmental field estimation[END_REF]. However, these fields can be also constructed using physical models which simulate the phenomenon of physical parameters. For example, air pollution maps can be obtained using atmospheric dispersion models based on traffic and weather data [START_REF] Soulhac | The model sirane for atmospheric urban pollutant dispersion; part II, validation of the model on a real case study[END_REF]. Physical models like atmospheric dispersion models allow us to get fields with much higher spatial and temporal granularity than interpolation-based techniques. However, the estimations are less accurate than the measurements performed by the sensor nodes due to the uncertainty in the inputs of physical models. In this context, the measurements of sensors can be used to correct the simulations of physical models rather than directly constructing the spatiotemporal field by interpolation or regression. This approach is called data assimilation and its efficiency has been already proven in the literature [START_REF] Tilloy | Blue-based no2 data assimilation at urban scale[END_REF].

In this paper, we focus on environmental wireless sensor networks that are used to correct the errors of physical models and we tackle the problem of sensor scheduling in order to maximize the lifetime of the network. Maximizing the lifetime of the network is a major issue in wireless sensor networks which usually operate using batteries. The most used definition of network lifetime is the time period during which the network is operational; this means that coverage is ensured (the application sensing requirements are verified) and the network is connected (every sensor is capable of sending its data to at least one sink node) [START_REF] Yetgin | A survey of network lifetime maximization techniques in wireless sensor networks[END_REF]. Several works in the literature have targeted the problem of network lifetime maximization and at different design levels: deployment, sensor scheduling, communication load balancing, transmission rate selection, transmission power selection, routing, etc.

Existing works on sensor activity scheduling usually assume sensors to have two operation modes: active mode where sensing, communication and computation can be performed; and sleep mode where the sensor consumes a very small amount of energy [START_REF] Biswas | Energy-efficient connected target coverage in multi-hop wireless sensor networks[END_REF]. Activity scheduling consists of keeping only a subset of sensors in active mode and can be performed in two ways: 1) in a distributed way where a sensor communicates with its neighbors to decide whether it should turn off or not; or 2) a central way where an optimal sensor schedule is determined by a central node (the sink node for instance) [START_REF] Tashtarian | On maximizing the lifetime of wireless sensor networks in event-driven applications with mobile sinks[END_REF]. Without loss of generality, we focus on the second case because the physical models are executed only on the sink nodes due to their high computation requirements.

We propose in this paper a scheduling approach for the application of correcting physical models simulations. We use data assimilation to design a mathematical formulation and then extract from this formulation a mixed integer programming model (MILP) using linearization techniques. We then analyze the complexity of the model and propose a heuristic algorithm to solve large instances of the scheduling problem. We evaluate our approach while considering the application case of air pollution monitoring and using a data set from the Lyon city, France.

To the best of our knowledge, this is the first work targeting the scheduling problem of WSN for the application of correcting physical models' simulations. The main contributions of this paper are as follows:

1) We provide a mathematical formulation of WSN scheduling based on data assimilation in order to correct the simulations of physical models. 2) We design an optimization model using linearization techniques. 3) We design a heuristic algorithm using linear relaxation. 4) We evaluate our proposal on a real application case, that is air pollution monitoring, in order to derive engineering insights on the effective scheduling of pollution sensors for physical models correction. The remaining of this paper is organized as follows: we first review the existing works on WSN scheduling in section II while highlighting the differences with respect to the application of physical models correction. Then, we present the problem formulation in section III and our method to physical models correction in section IV. Next, we present the design of our MILP optimization model in section V and the heuristic algorithm in section VI. After that, we present the simulation results in section VII. Finally, we conclude the paper and discuss the open perspectives in section VIII.

II. RELATED WORKS: WSN ACTIVITY SCHEDULING FOR LIFETIME MAXIMIZATION

In this section, we review the related prior works which target the scheduling problem of sensor networks. We focus only on the works providing optimization techniques for the centralized sensing activity scheduling problem as in our proposal.

In [START_REF] Hu | Hybrid genetic algorithm using a forward encoding scheme for lifetime maximization of wireless sensor networks[END_REF], authors consider the case where sensors are deployed randomly for target detection. They focus on the problem of finding the maximum number of disjoint sets where each sensor set is capable of ensuring coverage. They propose to solve this problem, known as the disjoint set covers problem, using a hybrid genetic algorithm. In this work, authors don't take into account the connectivity of the network and assume that sensors have the same sensing range and the same initial amount of energy. They justify their choice of not taking into account the network connectivity by the fact that in most target detection applications, a relationship between the sensing range and the communication range of nodes allows them to get a connected network by ensuring only the coverage constraints.

In [START_REF] Liaskovitis | Leveraging redundancy in samplinginterpolation applications for sensor networks: A spectral approach[END_REF], authors focus on field estimation applications where there is no sensing range and sensors are deployed in order to construct a spatiotemporal field of a given physical parameter. In order to maximize the lifetime of the network, they determine the maximum number of disjoint sensor sets as in the work of [START_REF] Hu | Hybrid genetic algorithm using a forward encoding scheme for lifetime maximization of wireless sensor networks[END_REF]. The sensor selection in the scheduling process is based on the data that is gathered by sensors right after the deployment. They propose to first learn the characteristics of the physical stochastic process and then use these characteristics as a reference in the definition of the sensing schedule. This work makes a strong assumption, that is the stochastic process describing the physical phenomenon is stationary. By minimizing the number of active sensors, less packets are sent to the sink node, which minimizes the energy consumption due to connectivity. However, authors don't take into account the connectivity constraints between active nodes and the sensing capabilities of nodes are assumed to be homogeneous.

In [START_REF] Lin | An ant colony optimization approach for maximizing the lifetime of heterogeneous wireless sensor networks[END_REF], authors consider the application case of target detection using heterogeneous sensor nodes. They propose to use a heuristic based on ant colony optimization to maximize the number of disjoint connected sensor cover sets and hence get a sensor schedule which maximizes the lifetime of the network. Compared to the work of [START_REF] Hu | Hybrid genetic algorithm using a forward encoding scheme for lifetime maximization of wireless sensor networks[END_REF], they take into account both connectivity and nodes' heterogeneity while considering two types of nodes: sensors and relay nodes.

In [START_REF] Du | Maximizing the lifetime of kdiscrete barrier coverage using mobile sensors[END_REF], authors focus on a particular case of target detection which is the barrier cover for border surveillance. They consider in their work a set of already deployed mobile sensor nodes and propose to maximize the lifetime of the network by minimizing the energy consumption due to both sensing and redeployment mobility. Authors assume in this work that sensor nodes are identical. In addition, they don't take into account connectivity constraints.

In [START_REF] Keskin | Wireless sensor network lifetime maximization by optimal sensor deployment, activity scheduling, data routing and sink mobility[END_REF], authors focus once again on the application of target detection while providing a joint deployment and scheduling optimization approach. They consider the case of heterogeneous WSN and provide a linear mathematical model and heuristic algorithms. Their approach is to optimize the network lifetime given a deployment budget by minimizing the energy that is due to both sensing and data routing. Unfortunately, their design is adapted only to target detection applications. As in [START_REF] Keskin | Wireless sensor network lifetime maximization by optimal sensor deployment, activity scheduling, data routing and sink mobility[END_REF], authors in [START_REF] Mini | Sensor deployment and scheduling for target coverage problem in wireless sensor networks[END_REF] focus on both deployment and scheduling for the application case of target detection. They consider a fixed set of sensors to be deployed and scheduled with different sensing capabilities. However, they don't take into account the connectivity of the network.

Recent works investigate more constraints in their scheduling approaches. In [START_REF] Chen | Efficient coverage and connectivity preservation with load balance for wireless sensor networks[END_REF], Authors propose a scheduling technique to select not only sensing nodes to turn on but also they schedule the communication paths during data collection while balancing the load due to communication. In [START_REF] Deng | Sensor scheduling for multi-modal confident information coverage in sensor networks[END_REF], authors focus on using WSN for monitoring multiple physical parameters at the same time. In [START_REF] Lu | Maximum lifetime scheduling for target coverage and data collection in wireless sensor networks[END_REF], authors focus on both scheduling and data collection while assuming that nodes have the same sensing range and communication range.

A. Discussion

Even if recent works are considering connectivity and nodes' heterogeneity constraints, they focus on mainly two applications: target detection and field estimation using interpolation. Target detection applications consider the use of sensors having a given sensing range like presence sensors. Field construction applications consider the use of environmental sensors' data to construct spatiotemporal fields of a given physical parameter. Our aim in this paper is different in the sense that we consider a different and a new challenging application which is the correction of physical models which are used to simulate the variations of physical parameters but are not accurate and need to be corrected using measurements. In this context, our objective is to turn off sensors where data is not needed for the correction of the physical model in question. To that end, we propose in what follows to use the data assimilation technique which allows us to combine the physical model outputs and sensor measurements in the estimation of environmental fields.

III. PROBLEM FORMULATION

We consider as input the map of a given urban area that we call the region of interest where the estimation of a given physical parameter can performed with a physical model. Let P (with |P | = N ) be a set of discrete points within the region of interest where sensors are already deployed. Sensor measurements are used in order to correct the physical model outputs at each point of the set P.

Let T denote the set of time periods over which we would like the network to be operational. We assume that at the end of each time period, sensors send their data to the sink node located at a given point p ∈ P. Based on the collected data, the sink node estimates the future errors of the physical model at each point of the region of interest. In order to save the energy of sensor nodes, we aim at turning off the nodes where the physical model future estimated errors do not exceed a critical level. We use binary decision variables x t p to specify if sensor p is scheduled to be active during the time period t or not. Our objective is therefore to determine for each future time period t ∈ T , the state of sensors (ON for active mode and OFF for sleep mode) based on the already collected data and the threshold of the estimation error that should not be exceeded after the correction of the physical model. This means that when sensor at point p is turned off, the data of the neighboring active sensors should be sufficient to correct the physical model output at point p. In addition, the network should remain connected during each time period in order to ensure the communication of sensor data to the sink node. The main notations used in this paper are presented in TABLE I.

Sets and parameters

IV. FORMULATION OF PHYSICAL MODELS CORRECTION

We propose to use data assimilation in order to correct the simulations of physical models: the estimated value Z t p at a given location p ∈ P where the sensor is OFF is formulated as the sum of M t p , which is the physical model simulation value at p, and a weighted combination of the difference between the physical model values M t q and the measured values at neighboring sensor nodes Z t q , q ∈ P where x t q = 1 [START_REF] Asch | Data assimilation: methods, algorithms, and applications[END_REF]. The weights used for the estimation are called correlation coefficients and can be evaluated in a deterministic way based on the distance between the location of the measured value and the location of the estimated value. These coefficients can be also evaluated in a stochastic way, but, without loss of generality, we focus in this paper on the case of deterministic data assimilation. In this case, Z t p is calculated using formula 1 where W pq denote the correlation coefficients [START_REF] Asch | Data assimilation: methods, algorithms, and applications[END_REF].

Z t p = M t p + q ∈ P W pq • x t q • (Z t q -M t q ) q ∈ P W pq • x t q (1)
Let G t p denote the ground truth (or real) value at point p. We denote by m t p (respectively s t p ) the physical model error (respectively the sensing error of nodes) which is defined as the difference between M t p and G t p (respectively the difference between Z t p and G t p ). With these definitions, formula 1 can be transformed into formula 2.

Z t p = M t p - q ∈ P W pq • x t q • (m t q -s t q ) q ∈ P W pq • x t q (2)
The data assimilation equation in formula 2 is constrained by formula 3, which ensures that the denominator is never equal to 0. B pq parameters define whether there is a correlation between points p and q or not; that is, B pq = 1 when W pq > 0.

q ∈ P B pq • x t q ≥ 1 (3)
Given the formula of the assimilation estimated value Z t p , the assimilation error with respect to the ground truth value (the difference between Z t p and G t p ) can be derived as in formula 4.

E t p = m t p - q ∈ P W pq • x t q • (m t q -s t q ) q ∈ P W pq • x t q (4)
Note that both physical model simulation errors (m t p and m t q ) and sensing errors (s t q ) are unknown values because the t index corresponds to a future time period. Therefore, we propose in this paper to consider these errors as random variables where only the the variance and the expectation are known by means of empirical analysis of the already collected data. We assume that the expectation of the errors is equal to 0. This is not a strong assumption since both the physical model and sensors can be calibrated to get an error expectation equal to 0 by adding or subtracting the real expectation. That is, the variance defines how much the model (or the sensors) are incorrect at a given point. Based on these assumptions, we define the coverage quality at a given point p and during a future time period t with respect to the set of active sensors as the variance of the assimilation error. To get this formulation, we apply the variance function to formula 4 while assuming that sensing errors are independent between them and are also independent with respect to the physical model errors. Hence, we get formula 5 where V ar (respectively Cov) denotes the variance (respectively covariance) function.

V ar

(E t p ) = V ar (m t p ) + q∈P W 2 p q •x t q •(V ar (m t q )+V ar (s t q )) ( q∈P W p q •x t q ) 2 -2 • q∈P W p q •x t q •Cov(m t p ,m t q ) q∈P W p q •x t q + q 1 p q 2 p, q 1 W p q 1 •W p q 2 •x t q 1 •x t q 2 •Cov(m t q 1 ,m t q 2 ) ( q∈P W p q •x t q ) 2 (5)
Note that the covariance Cov(m t p , m t q ) is mathematically a function of correlations W pq and variances V ar (m t p ) and V ar (m t q ) as in formula 6 [START_REF] Revelle | An introduction to psychometric theory with applications in r[END_REF].

COV (m t p , m t q ) = W pq • V AR(m t p ) • V AR(m t q ) (6)
In order to ensure that the application requirements are met during the network operation, the coverage quality presented in formula 5 should not exceed a threshold (or required) value E during each time period t where a sensor p is turned off. This constraint can be formulated as follows:

V ar (m t p ) + q∈P W 2 p q •x t q •(V ar (m t q )+V ar (s t q )) ( q∈P W p q •x t q ) 2 -2 • q∈P W p q •x t q •Cov(m t p ,m t q ) q∈P W p q •x t q + q 1 p q 2 p, q 1 W p q 1 •W p q 2 •x t q 1 •x t q 2 •Cov(m t q 1 ,m t q 2 ) ( q∈P W p q •x t q ) 2 ≤ E, p ∈ P (7)
We seek later in this paper a linear optimization model. Therefore, we need to linearize constraint 7 by eliminating the fraction and the multiplications between the decision variables.

We first multiply both sides of formula 7 by the denominator of the fraction. Next, we simplify the parts where the square function is applied to variables x t q . Hence, we obtain the linear form of our coverage formulation in formulas 8 and 9 where expressions expr 1 and expr 2 are detailed in formulas 10 and 11 respectively. Finally, real variables v t q 1 q 2 correspond to the linear form of the product of decision variables x t q 1 and x t q 2 thanks to constraints 12.

(V ar (m t p ) -E) • expr 1 + q ∈ P W 2 pq • x t q • (V ar (m t q ) + V ar (s t q )) -2 • expr 2 + q 1 p q 2 p,q 1 W pq 1 • W pq 2 • v t q 1 q 2 • Cov(m t q 1 , m t q 2 ) = C t p , p ∈ P (8)
C t p ≤ 0, p ∈ P (9)

expr 1 = q 1 ∈ P q 2 ∈ P W pq 1 • W pq 2 • v t q 1 q 2 ( 10
)
expr 2 = q 1 ∈ P q 2 ∈ P W pq 1 W pq 2 v t q 1 q 2 Cov(m t p , m t q 1 ) (11) v t q 1 q 2 ≤ x t q 1 , q 1 , q 2 ∈ P v t q 1 q 2 ≤ x t q 2 , q 1 , q 2 ∈ P v t q 1 q 2 ≥ x t q 1 + x t q 2 -1, q 1 , q 2 ∈ P (12)
V. OPTIMIZATION MODEL

In this section, we present our MILP optimization model based on the formulation of the physical models correction presented in the previous section. Our objective is to ensure that the network is operational during a maximum number of time periods while ensuring the coverage requirements when sensors are turned off and while ensuring also the network connectivity. We denote by α t = 1 the fact that the network is operational during time period t. a) Coverage requirements: based on the coverage formulation presented in formulas 8 and 9, we get the following constraint while integrating α t variables. Here, M is a sufficiently big number used to ensure that there is always a feasible solution in the case where the network cannot be operational during the whole number of time periods.

C t p ≤ (1 -α t ) • M (13)
b) Connectivity constraints: we formulate the connectivity problem as a graph flow. We first denote by Γ(p), p ∈ P, the set of neighbors of the sensor located at point p. This set can be determined using sophisticated path loss models. It can also be determined using the binary disc model, in which case Γ(p) = {q ∈ P where q ∈ Disc(p, R)} where R is the communication range of sensors. Then, we define the decision variables g t pq as the number of flow units transmitted from a sensor located at point p to another sensor at point q. We also use variables f t pq to denote the flow quantity transmitted from a sensor located at point p to a sink node located at point q. We suppose that each active sensor of the resulting WSN generates a flow unit in the network during each time period, and then verify if these units can be recovered by sinks.

The following constraints ensure that the active sensors form together with sink nodes a connected wireless sensor network; i.e. each sensor can communicate with at least one sink. Here y p denotes whether there is a sink node located at p or not. Since we are focusing on the scheduling of the sensing activity, we assume that there are no energy constraints on sink nodes.

q ∈Γ(p)

g t pq + q ∈Γ(p) f t pq - q ∈Γ(p) g t qp = x t p ( 14
)
q ∈Γ(p)

g t pq + f t pq ≤ N • x t p ( 15
)
q ∈Γ(p)

f t pq ≤ N • y q ( 16
)
p ∈ P q ∈Γ(p)

f t pq = p ∈ P x t p ( 17
)
Constraints 14 ensure that each active sensor node, i.e. such that x t p = 1, generates a flow unit in the network during each time period. In addition, constraints 15 ensure that non active nodes, i.e. x t p = 0, do not participate in the communication. Thanks to constraints 16, only the deployed sinks, i.e. y p = 1, are taken into account. Finally, constraints 17 ensure that the overall flow is conservative, i.e. the flow sent by the active sensor nodes has to be received by sinks. c) Energy consumption constraints: first, let EI p , ES t p , ET t pq and E R t pq denote respectively the initial amount of energy (battery capacity) of sensor p, the energy consumption due to sensing during time period t, the energy consumption due to the transmission to a neighbor during time period t and the energy consumption due to the reception from a neighbor during time period t. The following constraints ensure that sensor nodes cannot consume more than their initial amount of energy.

∀p ∈ P t ES t p • x t p + t,q ET t pq • g t pq + t,q E R t pq • g t qp ≤ EI p (18)
d) Lifetime of the network: finally the network lifetime to maximize corresponds to the number of time periods during which the network is operational as follows:

Maximize t ∈ T α t (19) α 0 ≥ α 1 ≥ α 2 ≥ ... ≥ α | T | (20)

VI. THEORETICAL COMPLEXITY AND HEURISTIC ALGORITHM

A. Complexity of the scheduling model

The proposed optimization model is based on integer linear programming that can be solved using exact MILP solvers. In this paper, we use the IBM Cplex solver later in the simulation part. In terms of complexity, the execution time of the MILP solvers increases exponentially with the size of the problem. In fact, what makes our MILP model difficult to solve is the number of binary variables which causes an exponential increase in the number of iterations when using the exact MILP solvers. That is, the complexity of the model is mainly due to the number of sensors and time periods. In order to alleviate the resolution process of the proposed model while being able to get solutions with a sufficiently good quality, the exact MILP solver can be used with an input integrality gap value. The integrality gap defines the quality gap between the theoretical optimal solution and the current solution of the MILP solver during its execution time.

B. Scheduling heuristic

In order to solve our optimization model on large instances in a reasonable time while getting good solutions, we propose to use the concept of linear relaxation. We first define the linear programming model LP while considering the same objective function and constraints as our initial model and relaxing binary variables x t p ; i.e. binary variables are considered in the range of [0, 1], this means that the solutions of the LP model are not necessarily binary. Note that in a given solution of LP where variables x t p are fractional, the variable having the maximum value (i.e. the closest variable to 1) corresponds to the most important variable in the satisfaction of coverage and connectivity constraints. Based on this fact, we propose in each iteration of our heuristic algorithm to activate sensor p during time period t where x t p is the closet variable to 1. The loop performs iterative rounding and stops once the scheduling variables are equal to either 0 or 1, all the coverage and connectivity constraints are ensured and the network lifetime can no longer be extended. In this section, we present the simulations that we have performed in order to evaluate our scheduling approach. We first present the data set that we used and the common simulation parameters. Then, we provide a proof-of-concept to show how we execute our approach on a real dataset. Next, we evaluate the physical models coverage results. Finally, we assess the impact of sensing frequency and transmission power on the the network lifetime.

A. Dataset

We perform the evaluation of our proposal on monthly pollution data corresponding to the 2008 Nitrogen Dioxide (NO 2 ) concentrations in the Lyon district of La-Part-Dieu, which is the heart of the Lyon City. This pollution data set has been generated by an enhanced atmospheric dispersion simulator called SIRANE [START_REF] Soulhac | The model sirane for atmospheric urban pollutant dispersion; part II, validation of the model on a real case study[END_REF], which is designed for urban areas and takes into account the impact of street canyons on pollution dispersion. The dataset has been provided by LMFA, which is a research lab specialized in fluid mechanics in the Lyon city, France.

The region of interest has a spatial resolution of 100 meters and is depicted in Fig. 1a with a set of 25 pollution sensors. We calculate the correlation coefficients W pq using an exponential decay function. That is, the correlation between points decreases exponentially with the euclidean distance.

We recall that the main input of our scheduling approach is the variance of the errors of the physical model. In this evaluation part, we assume that the errors of the model are linearly correlated with its concentrations. Let γ express the linear relationship between the model concentrations and the model errors. Thus, we first calculate the variance of the concentrations of the physical model based on the pollution maps and then we multiply these variances by γ 2 to get the variance of the physical model errors. We calculate the γ parameter by evaluating the linear regression between the concentrations of the dataset of the physical model and the real data of the few monitoring stations which are already deployed in the Lyon city. For illustration purposes, we depict in Fig. 1b the variance map of the physical model errors corresponding to the month of June.

In order to consider a real scheduling scenario, we use the characteristics of an air pollution monitoring system that we developed in our lab (see Fig. 1c). The energy consumption required for sensing, transmission and reception per month as well as the default simulation parameters are depicted in TABLE II. In all the results presented in this section, we used the Cplex MILP solver to optimally solve our scheduling model while using a PC with Intel Xeon E5649 processor under Linux.

B. Proof of concept

In order to provide a proof of concept of our scheduling approach, we run our model while considering the default simulation parameters and we set coverage requirements to 1(µg) 2 /m 6 . We get an optimal network lifetime equal to 3 months as depicted in Fig. 2 (from left to right: January, February and March). We notice that the number of active sensors is different depending on the time period. In fact, the correction of the air pollution physical model requires more nodes depending on the variability of air pollution within each time period. We also notice that some nodes are active in more than one time period. Indeed, according to the default simulation values, the battery capacity allows sensor nodes to run up to two months without stop.

We also evaluate the assimilation error provided by the sensor nodes that are active during the months of January, February and March. Results are depicted in the same Fig. 2. We notice that the assimilation error never exceeds the coverage requirement set to 1(µg) 2 /m 6 thanks to the coverage constraints of our optimization model.

C. Impact of coverage requirements on the network lifetime

In this simulation scenario, we investigate the impact of coverage requirements on the maximization of the network lifetime. We consider 3 different values of the battery capacity and we vary the the coverage requirements from 0.5(µg) 2 /m 6 Fig. 2: Proof-of-concept: optimal WSN activity scheduling and the corresponding estimation errors' variance ((µg) 2 /m 6 ) while considering coverage requirements equal to 1(µg) 2 /m 6 . Active sensors are depicted in blue squares. to 2.5(µg) 2 /m 6 . The obtained results are depicted in Fig. 3. We notice that the higher the assimilation error threshold, the higher the network lifetime. This is expected since with less coverage requirements, less sensors are activated during each time period allowing us to get multiple subsets of sensor nodes capable, each of which, of ensuring coverage and connectivity constraints.

In terms of battery capacity impact, as expected the higher the initial amount of energy per node, the largest the network lifetime. However, doubling the battery capacity does not always involve twice lifetime as shown in the curves when the assimilation error threshold is equal to 1.4 (µg) 2 /m 6 . Indeed, this depends on the variability of pollution within each time period and which may require different number of active nodes.

D. Impact of sensing frequency on the network lifetime

We now evaluate the impact of sensing quality on the network lifetime. We consider two scenarios varying the sensing frequency of nodes: in the first case, sensors analyze the quality of the air without stop whereas in the second case, sensing is performed during only 30 seconds every minute. It is worth mentioning that in the latter, the node is kept active during the non sensing 30 seconds which is necessary for the air pollution electrochemical sensing probes. Reducing the sensing frequency impacts both energy consumption and the correction of the physical model. In order to understand how the the sensing frequency impacts the network lifetime, we depict in Fig. 4 the obtained results while varying the coverage requirements. Results show that first, coverage requirements cannot be met with low sensing frequency for an assimilation error threshold that is less than 1(µg) 2 /m 6 . Indeed, with low sensing frequency, the sensing error of nodes goes up from 0.5(µg) 2 /m 6 to 1(µg) 2 /m 6 . Moreover, we still get lower network lifetime when using low sensing frequency. However, as the assimilation error threshold of the physical model goes up, we tend to get the same results for both sensing techniques because the sensing errors become tolerable with respect to the assimilation process.

E. Impact of transmission power on the network lifetime

In the last simulation scenario, we analyze the impact of transmission power on the network lifetime. Increasing the transmission power allows us to get larger transmission range but at the cost of energy consumption. We consider three different cases: i) high transmission power (our LoRa powered nodes with 20dbm transmission power) allowing us to get onehop communication to the sink node (communication range equal to 500m); ii) low transmission power (14 dbm) with a communication range equal to 200 in non line of sight while considering the sink node in the corner of the map; and finally iii) low transmission power with the sink node in the center. Results are depicted in Fig. 4 and show that despite the multi-hop communication, low transmission power leads to better network lifetime compared to high transmission power. However, this is not always the case as the results also show that one-hop communication is preferable over multihop communication if the sink node is not well positioned (in the corner rather than in the center of the map for instance).

VIII. CONCLUSION

In this paper, we focus on extending the lifetime of the network by scheduling the sensing activity of sensors. Our objective is to decide which sensors should be turned off in order to let the network operate longer while ensuring the application requirements. We design a mathematical model and a heuristic algorithm that we use to determine which sensors should be turned off based on the collected data. Related works in the literature focused only on two main applications: target detection and spatiotemporal field construction. Our work is different in the sense that we consider a different and a new challenging application which is the correction of physical models. We use the data assimilation technique to provide a mathematical formulation to the locations where data is not needed for the correction of the physical model. We apply linearization techniques and linear relaxation to get solvable optimization models. Finally, we evaluate our proposal on an air pollution dataset from the Lyon city, France and provide engineering insights on the optimal scheduling of air pollution sensors.

As future work, we plan to perform simulations on different environmental applications while studying the quality of our approach depending on the nature of the physical model. Another perspective of our work would be to analyze the impact of weather conditions on the longevity of the electronic components of sensor nodes.

Algorithm 1

 1 Scheduling heuristic Inputs: P Outputs: {x t p , α t } repeat Solve the LP scheduling model Let f be the maximum fractional variable among x t p Add constraint f = 1 to the LP model until all the variables are binary VII. SIMULATION RESULTS: APPLICATION TO AIR POLLUTION MONITORING

Fig. 1 :

 1 Fig. 1: WSN nodes locations, simulation errors corresponding to the district of La-part-dieu, Lyon, France during June, 2008 and an internal view of our sensor nodes.

Fig. 3 :

 3 Fig. 3: Impact of battery capacity.

Fig. 4 :

 4 Fig. 4: Impact of sensing quality.

Fig. 5 :

 5 Fig. 5: Impact of transmission power.

TABLE I

 I

: Main notations used in our approach.

TABLE II :

 II Default values of main simulation parameters.

ACKNOWLEDGMENT

We would like to thank the anonymous reviewers for their comments and feedback. This work has been supported by Mitacs, Canada and the "LABEX IMU" (ANR-10-LABX-0088) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR).