Christopher Yeates

Souhail Youssef

Elise Lorenceau

Accessing preferential foam flow paths in 2D micromodel using a graphbased 2-parameter model

Keywords: Foams, Micromodels, Preferential paths, Network analysis

published or not. The documents may come

Accessing preferential foam flow paths in 2D micromodel using a graph-based 2-parameter model

Introduction

Foams constitute an attractive method to enhance oil recovery due to their high viscosity, high potential for conformance control, low cost and low environmental impact.

To implement a large-scale foam operation, prior validation and injection optimization must be performed using a reservoir modelling software. While some attempts have been made to integrate microstructural parameters such as average pore size into upscaled foam models [START_REF] Ettinger | Influence of Texture on Steady Foam Flow in Berea Sandstone[END_REF][START_REF] Gassara | Equivalence Between Semi-empirical and Population-Balance Foam Models[END_REF]) a great deal of foam behavior laws are derived from Darcy-scale laboratory experiments, and require significant a posteriori parameter fitting to match observations. The numerical load would be greatly alleviated by the use, or at least constraint of, fitting parameters with a value obtained through a purely structural characterization of the porous medium.

Recent works on 2D micromodels [START_REF] Géraud | The flow of a foam in a two-dimensional porous medium[END_REF][START_REF] Yeates | New insights of foam flow dynamics in high-complexity 2D micromodels[END_REF]) have shown that a range of foam behaviors exist that are highly likely to contribute to Darcy-scale flow properties. Indeed, trapped foams and high-velocity preferential paths have been repeatedly observed, displaying ranges of velocity much lower and higher, respectively, than those seen in Newtonian flow Furthermore, it has been shown that high-velocity zones are accessed primarily by larger bubbles in a heterogeneous bubble distribution flow, containing the majority of flow in the model. The interest of characterizing the specific paths that these bubbles will take, as well as the number and rank of equivalent or alternative paths, becomes obvious.

A purely local description of flow has shown to be insufficient to difficultly access larger scale phenomena such as preferential paths and trapped zones For this reason, here we attempt a novel characterization approach, by integrating the notion of overall path (spanning inlet to outlet of the model) and ranking paths according to a numerical value derived from the components that make up it.

Modelling porous media with graph models has been recently given new life. Graph-based models have been used to find least resistance paths and predict breakthrough points successfully in heterogeneous porous media [START_REF] Rizzo | Minimum Hydraulic Resistance and Least Resistance Path in Heterogeneous Porous Media[END_REF].Graph models have also found use in discrete fracture modelling to rapidly characterize, query, and interrogate fracture network connectivity [START_REF] Viswanathan | Advancing Graph-Based Algorithms for Predicting Flow and Transport in Fractured Rock[END_REF] The path-proposing is done via an algorithm that makes use of a graph representation of the porous network. Nodes of the graph represent the pores, while the graph edge weights represent throats. This gives a labelled graph of connected objects, with edge weights given by a throat property of our choosing. The sum of edge weights of a given path between the source node and the target node p(s → t), otherwise known as the path length, is then given by ∑ w(i, j)

i,j∈p (s→t) . The shortest path is then path that minimizes this sum. Many algorithmic solutions exist for this problem including Dijkstra's shortest path algorithm [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF])An extension of this problem is interested with the first K paths, ranked by increasing value of path length (K-shortest paths), and can be solved with Yen's algorithm for example [START_REF] Yen | An algorithm for finding shortest routes from all source nodes to a given destination in general networks[END_REF].When applied to our porous network graph, solutions to the K-shortest paths problem often are small variants of the best path, with one or two elements of the paths changed, even for large values of K (>5000). When combined with our observations of preferential path flow in the micromodels, in which a multitude of somewhat distinct paths serve as preferential paths at the same time, we are motivated to assimilate our problem to an extension of the K-shortest paths problem, finding the K-Shortest Paths with Limited Overlap (k-SPwLO). The path solutions to this problem are the K-shortest paths respecting a condition of maximal path overlap between one another. We make use of the recently published OnePass algorithm, a solution algorithm that efficiently proposes candidate paths with user-chosen maximal path overlap [START_REF] Chondrogiannis | Exact and Approximate Algorithms for Finding k-Shortest Paths with Limited Overlap[END_REF]Experimentally observed flow maps serve as the basis for comparison of proposed paths. Experimental images showing long preferential paths of flow in a fully characterized micromodel are analyzed according to a well-defined statistical workflow to quantify active and dead zones in the network. Paths are then proposed using a graph representation of the porous model. The path fit to experimental data is evaluated by counting the proportion of elements in the proposed paths that are active in experimental flow map, corresponding to a match value.

Materials and methods

Foam data acquisition procedure and injection conditions

Our experimental data acquisition setup and procedure is described in detail in [START_REF] Yeates | New insights of foam flow dynamics in high-complexity 2D micromodels[END_REF]. We create velocity heatmaps by averaging local tracked bubble velocities over multiple sets of 200 highfrequency images. Fixed injection rates and gas fraction are injected into a foaming device upstream of the micromodel, then connected in series to the model. Pressure gradient is measured, and images are acquired at steady state (constant pressure gradient over the model). The image processing procedure leading to bubble segmentation as well as the tracking procedure is also described in [START_REF] Yeates | New insights of foam flow dynamics in high-complexity 2D micromodels[END_REF].

Injections were done with purified N2 gas. The surfactant solution was composed of purified water with 30 g/l of salt and an AOS surfactant at 10 times the CMC, supplemented with betaine-based foam booster.

Two experiments are described in this study. We compare injections of same gas fraction but different injection rates. The two injections were performed in succession. Experiment 1 was at a gas fraction of 79% at an injection rate of 2.53 × 10 -2 cm 3 /min. For Experiment 2, as lower injection rates were not obtainable with the gas flow controller, a biphasic reservoir, essentially a long coiled up tube, was used to enclose the gas and surfactant solution at the chosen gas fraction. The mixture was then pushed using the liquid pump, reaching much lower injection rates, into the foaming device. The injection rate for Experiment 2 was 8.9 × 10 -3 cm 3 /min. Despite a reduction in injection rate of 65%, the measured pressure drop over the model decreased from 275 ± 4 mBar for the high rate experiment to 189 ± 2 mBar for the low rate experiment, corresponding to a 32% decrease in pressure gradient, showing a highly shear-thinning apparent viscosity profile coherent with other foam studies [START_REF] Alvarez | Unified Model for Steady-State Foam Behavior at High and Low Foam Qualities[END_REF].

As previously shown [START_REF] Géraud | The flow of a foam in a two-dimensional porous medium[END_REF][START_REF] Yeates | New insights of foam flow dynamics in high-complexity 2D micromodels[END_REF], the majority of flow in a polydisperse foam injection into a heterogeneous medium is carried by the largest bubbles, in a series of preferential paths. Furthermore these paths show higher relative flow intensity than for Newtonian fluids [START_REF] Yeates | New insights of foam flow dynamics in high-complexity 2D micromodels[END_REF]. Correct characterization and modelling of the flow within these preferential paths will constitute a significant part of macroscopic foam flow modelling due to the large contribution of flow solely in these paths. For the experimental comparison, we therefore only use the flow maps of only the right-hand side of the bubble size distributions.

We show in Figure 1 the bubble size spatial probability distributions for each experiment, and highlight the bubbles used contributing to the flow maps used in this study. The spatial probabaility distributions are weighted histograms of bubble frequency histograms, with each histogram frequency peak weighted by the average bin bubble size (taken as bin center value). Such distributions give a more accurate depiction of the bubble size distribution inside the medium. The bubble size distributions are established by counting bubbles measured in the entire model. The bubble sizes contributing to the flow maps in this study are highlighted in red.

Figure 1: Bubble frequency histograms and bubble sizes contributing to flowmaps in each case highlighted in red. We observe that the lower minimum of bubble sizes is approximately equivalent in each experiment, but the experiment 2 distribution shows a larger right-hand tail. It is to be noted that although the higher bound shown here was cut for visibility, larger bubbles were observed (especially in experiment 2) and took up non-negligible areas of the porous volume.

Micromodel structure and decomposition

The 2D micromodel used has previously been studied elsewhere [START_REF] Yeates | New insights of foam flow dynamics in high-complexity 2D micromodels[END_REF]. It is a 69.7% porosity glass micromodel with an wet-etched depth of 40 µm. The permeability of the total model with the flow spreading system before and after the model is 4.7 Darcy. A binarized version of the model was decomposed into pores and throats using an adjustable watershed algorithm [START_REF] Soille | <title>Determining watersheds in digital pictures via flooding simulations</title>[END_REF] with the tunable sensitivity parameter set to 2. The model decomposition gives a total of 3700 pores and 6284 throats. The throat radii show a unimodal distribution with sharp peak at 51 µm. A network extraction algorithm was developed and applied to the 2D decomposed porous medium to find neighborhood properties of each pore and a graph was created. We show the different steps of the graph extraction process in Figure 2: Note here that this graph simply displays graph topology by positioning nodes on the centers of mass of the pores and edge lengths or widths are not representative of edge weights.

Due to the quantized pixel nature of the model image, and the design of the watershed algorithm, some of the throats were found to not be perfectly straight and resemble S-shapes. As we make substantial use of the throat size in our study, use of the throat length as given by the watershed was inadequate. Instead, we use the maximal Feret diameter of the watersheds, as measured with ImageJ analyze particles tool, to access the exact size of the throats. Half of the measured Feret diameter was then approximated to represent the throat radius, referred to in the rest of the article simply as throat size.

We show the model and associated graph in their entirety in Appendix A, with the added artificial inlet and outlet pores (described below) visible at each extremity.

Path-proposing algorithm

We make us of the recently published One-Pass graph algorithm [START_REF] Chondrogiannis | Exact and Approximate Algorithms for Finding k-Shortest Paths with Limited Overlap[END_REF], solution to the K-shortest paths with limited overlap problem (K-sPwLO), to provide a series of paths that resemble the preferential paths visible in our experiments.

In our context, to capture the specificity of the flow in the porous network, edge weights are represented by a function of the throat properties, initially the throat radii, located at the distance map minima, that will contribute the most to variations in hydraulic conductivity. In this sense, we try to minimize overall sum of contributions of throat radii, rather than distance. Note, that this problem is distinct from the shortest bottleneck, or widest path problem, which yields the shortest source-to-target path with the largest minimum throat size possible.

Naturally, as we desire to find a path that spans the entirety of the model, we construct an artificial source node and target node representing the model inlet and outlet, sharing connections with all the nodes at each extremity of the model, as to not restrict the path-proposing procedure to a specific point in either side of the model.

As we want the largest throats to contribute the less to the overall sum, and hence be chosen in the paths, we use functions of the inverse of the throat size for edge weights. However, it is unclear which (positive) exponent α should be used to integrate this throat size as an edge weight. The shortest path will then

minimize the function ∑ ( 1 r(i,j) ) α i,j∈p(s→t)
in which r(i, j) is the radius of the throat connecting pores i to j, and s and t represent the source or target nodes, or vice-versa, as the graph is undirected.

While individual throats radii will preserve their rank if they are raised to a power α = 1 or α = 2, the rank of a sum of throat radii forming a path will not be conserved for each exponent. The larger the value of α, the more a small, difficult to pass throat will contribute towards the sum. The proposed paths will then avoid smaller throats, at the price of increasing the overall path length and tortuosity. Inversely, a lower value of α will render the algorithm less sensitive to the actual value of throat size, as the virtual length of the path (or number of throats that compose it) becomes more significant. At the extreme of α = 0 , all weights are the same and we retrieve the shortest path based on number of elements alone. We give an example of different proposed paths for different values of α in Figure 3. We see that different throat radii exponents produce different shortest paths when integrated as edge weights. The shortest paths connecting the two nodes marked by red dots are shown for values of α = 0,1,2. For α = 0, all the weights are equal to 1 and the optimal path is simply the one with the less throats, shown as a full black line passing a throat of size 2. For α = 1 the thinly dashed path is chosen with two throats of size value 5. Indeed, the sum edge of weights given by ∑ (

1 r(i,j) ) i,j∈p(s→t)
is the smallest for the thinly dashed path, such as

1 5 + 1 5 < 1 7 + 1 7 + 1 7 <
1 2 or 0.2 < 0.428 < 0.5. For α = 2, the thickly dashed path is chosen, composed of three throats of size 7. Indeed, the sum edge of weights given by ∑ (

1 r(i,j) ) 2 i,j∈p(s→t)
is smallest for the thickly dashed path, such as

3 7² < 2 5² < 1 2²
or 0.061 < 0.08 < 0.25, despite being the path with the most throats.

Experimental path match

To evaluate a candidate path, we make use of a system of classification of active/inactive pores of the experimental data, as previously described in [START_REF] Yeates | New insights of foam flow dynamics in high-complexity 2D micromodels[END_REF], Appendix B. Average pore velocity intensity values are compared to values for pores of similar sizes (split between 20 equidistant intervals of area). If a threshold is passed, then the pore is considered active. The threshold values are chosen as fractions of the average pore intensities for each of the pore sizes categories. In this way we extract ourselves from making any hypotheses on the way velocity intensity is spread out over the different pore sizes and simply recalculate a threshold value for each size category. We chose 19 size categories. We define two classification thresholds: half of average (low threshold) and the average intensity (high threshold). The pores contained the proposed paths can then be tallied up and the proportion of active pores within the proposed path (situated between 0-1, i.e. the fraction of active pores within the candidate path) assesses quantitatively the match of the candidate path to the preferential path in the experimental image. This binary classification (i.e. active or inactive) method is preferred over a more direct method of taking average pixel intensity values of the whole path for two reasons. Firstly, a direct comparison of average intensity values for pores of different sizes seems prejudicial to larger pores, in which flow is more spread out, creating lower average intensity values. Secondly, due to the high spread of average pore intensities, taking the average pixel value over the path may create the illusion of a satisfactory fit when only a restricted number of high-intensity pores are contributing, and the rest of the candidate path is incorrect.

In Figure 4 we show the active/inactive pores for two different threshold values overlaid on the velocity map for both experiments studied. The active pores are shown in green and the inactive pores in red. We notice the path-like nature of the active pores and the low-flow inactive pores, in red, situated in between these paths. In these experimental conditions, we observe an important number of active paths and flow distributed equally throughout all sections of the model. In Table 1 we display the proportion of active pores for each experiment and threshold value, that serve as baseline values for assessment of the match of proposed paths to the experimental data. We note that the ratio of low to high threshold pores is equal for both experiments, demonstrating the equivalence between the two thresholds for both experiments.

Results

In the results section we proceed by first describing the model used, and give associated statistics of the proposed paths in terms of geometrical and network properties only. Secondly we show the match with experimental data of the proposed paths. This sequence is repeated twice: for the 1-parameter model, and for the 2-parameter model.

1-parameter model -description and path properties

To propose a series of realistic paths, the parameters we can input into the proposing algorithm are: the graph, that includes a function of throat sizes as edge weights; the number of desired paths, K; and finally the maximal allowed overlap value θ, given between 0, for no overlap, and 1, if one path is fully contained in the other. We share here some statistics for different functions of throat sizes, with varying maximal path overlap value. The number of output paths, K, is maintained at a fixed value of 7. The throat size exponent spans 31 values, varying from 0 to a large value of 6, even if physically unrealistic, for the sake of assessing limiting behavior. The maximal overlap is restricted to 3 values for comparison: 0, 0.1, 0.5. The θ = 0 case describes a situation in which no pore can be shared between any of the proposed paths, effectively slicing the model in two separate parts for each new path found. The θ = 0.1 case allows limited overlap, in which some short passages that could be structurally attractive for flow, can be proposed in all the paths output by the algorithm. The θ = 0.5 case represents a higher boundary, in which the number of shared pores between proposed paths can be up to half of the pores in either of the paths. See [START_REF] Chondrogiannis | Exact and Approximate Algorithms for Finding k-Shortest Paths with Limited Overlap[END_REF] for a detailed description of θ calculation.

Note that for some cases, for θ = 0.1 notably, when calculation time was too long for 7 paths, the algorithm was restarted, aiming for 6 paths, and so on, until it reached a solution in a practical amount of time (here 2 000 seconds). Figure 5 shows some characteristics of the proposed paths for each overlap value and throat size exponent. Statistics for each point are calculated over the average for all the 7 paths output by the algorithm for each throat size exponent α and overlap value θ. We observe that as the throat size exponent increases, the output path lengths increase, as described previously. We observe that for α = 0, the average throat size for the paths is close to the average model, throat size, providing confirmation that these paths are independent of throat size. However, the average throat size seems to reach a plateau around a value of approximately α = 2.5. This is understood as follows: for higher exponent values, the algorithm outputs paths that will avoid smaller and smaller throats at the expense of creating longer paths. However there exists a point in which the smaller, avoided throat contributes a negligible amount to the average, as the number of throats contributing to the average is ever increasing.

Another observation is the difference between the characteristics of the θ = 0 paths on one side, and the θ = 0.1 and θ = 0.5 paths on the other side. The number of throats for the θ = 0 paths at the lowest throat size exponent value is smaller and the value of the plateau of the throat size average is smaller. This translates the fact that each proposed path for the θ = 0 output effectively separates the model into two parts and makes it harder to find short paths with small throat sizes in the created sections of the model.

Experimental match: 1-parameter model

We now perform a scan over the proposed paths for each setting and attempt to retrieve the highest fit to the experimental data. We display in Figure 6 the active pore match values for both threshold values, for different values of θ and the 31 throat size exponents for both experiments. Each exponent and θ pair form a algorithm input setting, from which all the pores in all the proposed paths are noted and compared against a dictionary of active pores for both thresholds and both experiments. One input setting is then evaluated in four different ways, as shown by the four distinct plots in Figure 6. The match value is the ratio of active pores to all pores in the proposed paths for the given threshold value, and is referenced interchangeably as a match percentage. The observed match is globally very high. For Exp. 1, most input settings propose paths in which 90% of pores are active as defined by the low threshold and 65% as defined by the high threshold. In either threshold such ratios are more than 25% above the model average (see Table 1). For Exp. 2, in most input settings 80% of the proposed pores are active as defined by the low threshold and 60% as defined by the high threshold. The best fits occur for throat size exponents between 1-3 for the low threshold pores and 1-2 for the high threshold pores. This still gives a large degree of uncertainty as to how to best describe foam preferential paths in a characterized porous medium. The match values for θ = 0.1, 0.5 are not significantly larger than the values obtained for θ = 0. The match values for the high threshold for θ = 0.5 however are significantly lower than the other values. For this reason, and for the reasons we previously mentioned issues occurring with the θ = 0 value, we shall use the value of θ = 0.1 for the rest of the analysis.

We first give some examples of relatively satisfying results. In Figure 7 we give three images of an overlay of the proposed paths calculated using α = 2, θ = 0.1 and displayed on the flow data of Experiment 2.

The active pores represented are from the high threshold activity classification. This is the image with the lowest number of active pores, so statistically hardest to match.

Figure 7: Examples of successful matches between model and experiment, flow is from left to right. Experimental velocity fields are shown in greyscale inside the pores (darker is higher velocity). The grains are in solid gray. The color code is the following: green pores are active in both experimental and proposed cases, red pores are proposed active but experimentally are inactive, purple pores are not proposed active but are experimentally active.

We observe long sections high intensity flow are identified by the proposed paths (green overlay). The top image of Figure 7 shows a chain of 30 pores in which 28 are correctly identified by a proposed path.

The bottom right image of Figure 7 shows a chain of pores in which 20 out of 21 are colored in green.

While these results are very encouraging, we still note some parts of the proposed paths that don't adhere well to the flow reality. One telling example is shown in the bottom right image of Figure 7 in which the proposed path takes makes a sharp curve around an obstacle while the flow takes a straighter alternative path. Here we can start to see the limitations of using a simplified model of the network that only takes the throat size into account and disregards any notion of path distance or tortuosity. In Figure 8 we give three more examples of proposed path failures occurring for the paths calculated using a throat size exponent of 2, θ = 0.1 and displayed on the flow data of Experiment 1. The active pores represented are from the low threshold activity classification. This is the image showing the largest number of active pores, in which a failed path match should statistically happen the less often. It's worth noting that although many purple pores are visible, these do not describe a model failure as only a limited number of paths were calculated, and it is possible that a larger number of output paths could have included these. In the three cases shown here we identify a common error: proposed paths go off the experimental track by moving into pores that are perpendicular to the flow direction and aren't visited by the foam, despite not showing significantly small throat entrances.

2-parameter model -description and path properties

We make an addition to the 1-parameter model to offset the proposed path failures occurring for pores branching outwards to the flow direction. We introduce a second local parameter in the description of the flow, the angle of the pore centers surrounding a throat in relation to the flow direction: noted γ. The angle is formed by the vector the connecting the two pores (oriented in the flow direction) and direction of pressure gradient. An example of the angle definition is shown in Figure 9. We refer the reader to the Appendix B for some statistics relating to this parameter and the independence of the 2 parameters. The value of γ has therefore a maximal value of 90° when two neighboring pores are perpendicular with respect to flow direction, and a minimal value of 0° for pores aligned in parallel to flow. Inclusion of the γ parameter in the edge weight of the throat connecting pores i to j in the graph, as noted by w(i, j), is done by taking a function the cosine of γ and multiplying it by the throat size term. A similar function to the throat size term is used, with a simple exponent β that is left free to vary. We then write the edge weight function, with both throat size and γ angle terms as:

w(i, j) = ( 1 r(i, j) ) α ( 1 cos(γ(i, j)) ) β
For simplicity, the allowed values of the β parameter vary in the same way as the α parameter, in 31 equally spaced values from 0 to 6. We show statistics associated to the proposed paths in the 2D parameter space in Figure 10. In accordance with the explanation of the 1parameter model results, path length increases as the algorithm is fed more and more widely spread edge weight values in the graph (i.e. some throat weights become so large that it chooses to take a series of smaller values with a lesser contribution to the overall path sum). The (α, β) = (6,0) path length maximum is then understood as occurring because of a widely spread throat size term only, with no contribution from the angle term, explaining the average throat size rise in this zone. The (α, β) = (0,6) path length maximum is then understood as the equivalent phenomenon, in which a longer path is chosen purely on the basis on avoiding unfavorable pore-to-pore angles. The average throat sizes in (α, β) = (0,6) (Figure 10b) and angle cosine in (α, β) = (6,0) (Figure 10c) confirm this observation. For these parameter values, they respectively take the model average throat size and the model average cosine value (for the full model: r̅ = 51 µm; cos(γ) ̅̅̅̅̅̅̅̅ = 0.63). What is not easily understood is the decrease in path length across the α = 6 parameter series, before the rise to the (α, β) = (6,6) maximum.

Increasing β with a fixed α should create a larger weight spread and lead to longer paths. We refer the reader to Appendix C for this analysis.

Experimental match: 2-parameter model

We show in Figure 11 colormaps of experiment match values for each edge weight parameter combination for both experiments and both thresholds. We observe that many parameter combinations match well the experimental paths. Low values of β give the best matches.

In fact, for experiment 2, the best matches are observed for β = 0, i.e. the one-parameter model. For experiment 1, while a large zone of good experimental matches is observed for the low intensity threshold for values of α < 3 and values of β < 1.6, the best values for the high intensity threshold are observed for (α, β) = (1.8 -2.4 , 0.3 -0.9), i.e. for a non-zero β, with the angle parameter adding capacity to match the data more closely. The best fit for experiment 1, high threshold (bottom left of Figure 11) using a 1-parameter model is obtained at (α, β) = (1 , 0) obtaining a match score of 0.749, while the best fit using a 2-parameter model : (α, β) = (2.2 , 0.6), achieves 0.7573, a small but non-negligible 1.1% increase.

To interpret this result, we consider the varying flow rates in both experiments. The lower flow rate of experiment 1 could shows no angle dependence. Indeed, at lower flow rates, we observe that preferential paths are less chosen on the basis of immediate path straightness or pore-to-pore alignment with flow, as described by the angle parameter, and instead are mainly chosen as tradeoffs between path shortness and throat size, as described by the one-parameter model. With faster flow rates, while the paths seen in the lower velocity experiment are also accessed, we observe more paths are accessed, which in turn creates a more parallel flow. In this way we observe less paths don't branch off perpendicularly to the flow direction.

The colored imaged showing the results from the combination displaying the experiment 1, high threshold best fit, (α, β) = (2.2 , 0.6), shown by a small marker in Figure 11, is shown in its entirety in Appendix D. We observe large sections of the flow-carrying paths are correctly matched by the algorithm. While a large number of incorrectly proposed pores, colored in red, do appear, upon inspection they seem to still carry flow and aren't entirely incorrect such as in the examples given in Figure 8. Furthermore, as active/inactive pores are determined on a local basis and not included in a path, it is not even guaranteed that a modelspanning path of active pores will truly exist. Rather, the path should be evaluated in its capacity to capture the largest number of active pores in the model.

Discussion

Path-based flow characterization viability

Here we stress the effectiveness of the path-based formalism by showing the results for the best model fit for experiment 1, (α, β) = (2.2 , 0.6). In Figure 12 we display how different models successfully capture high-flow pores, even with the simplest model explored (α, β) = (0 , 0), in which no structural parameters are included into the edge weights, and the paths are simply the shortest (in graph representation). However, the (α, β) = (0 , 0) model does implicitly take the porous structure into account as the decomposition into pores and throats itself was performed using a distance-map based watershed. A comparison is also made to a model proposing paths that minimize physical path distance.

The paths minimizing the physical distance were obtained by taking throat weight the sum of the distance between the pore center and throat center for both pores associated with the throat. This distance is named ∆𝑃𝑡𝑃. The sum of the weights for all throats in the path then was used as physical path distance.

In Figure 12, left we show for varying values of N, the percentage of N-highest average velocity pores in the entire model that are captured in the paths. For emphasis, the model-proposed paths are also compared to a random choice of pores (right). The success probability for the random choice of pores was simply the number of chosen pores (same number as in the proposed paths) over the total number of pores.

We observe that the path-based models perform better than a random choice of pores (Figure 12, right). Even the simplest model, based purely on minimizing the number of pores, successfully predicts a large degree of high flow pores, achieving similar results to the physically shortest paths. The similarity between the results of the shortest distance model and the smallest number of pores model seems to indicate that both descriptions are similar and could originate from a spiked distribution of distances between first neighbors. However, upon inspection of the proposed paths, this appears to be untrue and the proposed paths are quite distinct. The paths are shown in the Appendix E. Finally, the best model fully predicts the 3 highest flow-carrying pores and almost 50% of the top 20. This demonstrates the viability of a pathbased graph model for characterizing high flow zones for foam injections in a heterogeneous porous medium.

Difference between experiments

We note that the best models for each experiment take different parameter values. The active pores in Experiment 1, the high rate experiment, are best predicted by the 2-parameter model with (α, β) = (2.2 , 0.6). However, Experiment 2, the low rate experiment, is best predicted by the 1-parameter model with (α, β) = (1.2, 0). The problem this study attempts to solve, i.e. predicting preferential paths of a polydisperse foam in a heterogeneous porous medium, is akin to finding the paths of least resistance in a medium, evaluated for combinations of local structural parameters taken as hydraulic resistance. The path of least resistance then minimizes the sum of resistances between both ends of the model. For a rectangular channel, of width W and of height H, length L and of fluid viscosity μ, the hydraulic resistance R h can be approximated as [START_REF] Tanyeri | A microfluidicbased hydrodynamic trap. Design and implementation[END_REF], for H W ⁄ < 1:

R h = 12μL WH 3 (1 -0.63 H W ⁄ )
The path of least resistance then minimizes the sum R tot = ∑ R h for all the channels in the model. The validity of this equation depends on the exact ratio H W ⁄ . For values of H/W close to 0.5, the error is approximately 0.15% and for H/W = 0.1 the error drops to 0.003%. We show a histogram of H/W values in Figure 13 showing that most of the throats in our model respect this condition. The calculated hydraulic resistance from the equation above is also shown in a Log-Log plot making the throat width dependencies apparent (right). For simplicity, the throat length is taken a constant value here, given by half the mean ∆𝑃𝑡𝑃 value for all first neighbors, and water viscosity is used. The R h values are shown as red dots. To guide the reader, we give two lines displaying on one side the asymptotic dependency of R h as W -1 at large throat widths and provide a W -2 line tangent to R h values at lower widths. We supplement the plot with a histogram showing the distribution of widths in logarithmic space on top of the Figure 13. This plot shows that most throats are situated in an intermediate zone of mixed dependence between W -1 and W -2 , demonstrating that a series of models with values of α between 1 and 2 may in fact be physically appropriate. For this reason, and with the understanding that the average throat size in both models does not change significantly, (see Figure 10, center), we will not justify the higher α exponent in the best fit for experiment 1 as its physical justification is not obvious.

The best model fit for Experiment 1 however, includes a non-zero contribution from the angle parameter, with (α, β) = (2.2 , 0.6). This is a demonstration of difference in behavior for higher velocity foams, in which, while showing a strong dependence to throat size, preferential paths are also chosen based on pore-to-pore alignment with pressure gradient, despite their capillary ease of entry.

In supplement, Figure 14 displays the average velocity in the pores, rather than the number of active pores, to reinforce the validity of the model as best fit in both experimental matching frameworks. We observe a maximum located in the same parameter pair as found using the active/inactive pore classification. The presence of a non-zero β parameter in the best fit is conformed microscopically by looking at the bubble displacement distributions. Access to individual bubble tracks can provide a displacement term in the Y-direction (perpendicular to pressure gradient) compared to total displacement in all directions. Ydisplacement and total displacement were taken simply as Euclidean distances between start and end point of the bubble tracks. Taking the inverse sine of the ratio of Y to total displacement yields the angle of deviation of the bubble from the pressure gradient direction, named 𝜙. We take care here to distinguish the displacement of flowing bubbles from the quasi-zero displacement of static bubbles, which adds noise to the final probability distribution by filling it with the displacement ratios of static bubbles. We show in Figure 15 the different bubble displacement distributions that are bimodal in log-space, allowing us to segment the flowing bubbles with ease as the right side of the bimodal local minimum for each distribution. We offer the explanation that the higher flow rate and pressure gradient enable a larger number of paths to be accessed by the foam and remain aligned with pressure gradient to minimize the number of flow path intersections in which a large amount of energy is dissipated through viscous shearing. Inactive paths are observed more frequently in the lower velocity case. They can originate from either capillary difficulty of entry or insufficient pressure gradient required to dislodge trapped foams. The larger pressure gradient of experiment 2 facilitates foam mobilization creating a larger number of active paths. In turn, a flow path to crossing the medium in the Y-direction would imply numerous intersections with other active flow paths in which crossroads situations arise. Flow paths then can't freely traverse the medium in the Ydirection due to pre-existing flow paths, staying aligned with pressure gradient. In this sense, the initial motivation of overlapping preferential flow paths requires some refining. While we observe that some flow paths join and create higher velocity sections for a given distance, or otherwise split into two branching downstream paths, we rarely observe point-like intersections between two continuously flowing paths from different directions. More exactly, these kinds of intersections create either flow paths that deflect off each other, if the geometry permits, or otherwise creates intermittent traffic-like flow. We give two examples of flow paths that show large vertical displacement in lower velocity experiment, while being contained by surrounding active paths in the high-velocity experiment. In Figure 16 we can observe two distinct cases in which pre-existing flow paths in the higher velocity experiment inhibits the foam access to these zones. Following the flow from left to right, we see high intensity paths describing significant movement in in the Y-direction only in experiment 2 in two key areas circled in green.

As well as the different injection rates, the experiments also show variable bubble size distributions within the model. Larger bubbles contribute to the flow map in Experiment 2 as shown in Figure 1. This could also be a source of dissimilarity in the fitted models. Indeed, while the lower bound of bubble sizes is roughly equivalent, larger bubbles could in fact guide the flow within the experiment 2. For a lower injection rate and lower induced pressure gradient, the necessary Laplace pressure gradient needed to enter certain smaller throats may not be achieved, and large bubbles will in fact test multiple downstream throats before entering the widest. The smaller bubbles shown in Experiment 1 will mostly be of the order of the throat size, therefore not probing multiple throats for the widest possible entrance. Instead, in this situation, straighter paths therefore used.

The presence of the higher throat size exponent α = 2.2 in the Experiment 1 best fit is more intriguing. This exponent could either reflect a physical phenomenon for faster flowing foams, or could be a consequence of the intermediate dependence (between α = 1 and α =2) of the hydraulic resistance for our model throats, as shown in Figure 13b. Finally, this effect could more probably be a product of the interaction between the two parameters within the model. Indeed, we can remark that the good model fit zone, shown in bright colors in Figure 11, has a shape that evolves diagonally in the positive direction of both α and β. This behavior could show that the two components of the graph weight function interact non-trivially, and is confirmed by the lower skew in edge weight distribution in this area, as shown in Appendix C. In this sense, the higher exponent of the best fit for Exp. 1 could simply be a consequence of the best model necessitating a non-zero value of β.

Conclusion

In this study we show how a simple graph-based model can successfully capture zones of high-flow in experimental foam data in a porous medium. From experimentally established flow maps for different experiments, we attempt to characterize the paths chosen by the largest bubbles, that have previously been shown to contain a significant portion of the flow in a series of distinct paths. The two experiments investigated here are performed with varying injection rates and display slightly different bubble size distributions. A graph representation of the porous network is used in which edge weights are chosen as functions of two separate throat properties. A model parameter optimization is done to establish the best model fits for each experiment individually. We show how a simple 1-parameter model based on throat size retrieves a large amount of the high flow areas in a complex medium. Furthermore, we include a 2 nd local structural parameter that describes the alignment of pores with respect to pressure gradient. We observe that a larger sensitivity to the pore alignment parameter exists for the experiment with a higher injection rate, with smaller sized bubble sizes, whereas the 1-parameter model best describes slower injections with the presence of large bubbles. While this path-based approach shows promise by identifying high-flow zones in a complex medium using a simple model, further investigation into the effect of varying model structure and injection conditions in necessary to generalize the behavior seen here.

Figure 2 :

 2 Figure 2: Image and network analysis process. Porous area is shown in black. Throats are shown in light green. Extracted graph nodes are shown as red circles with edges represented as grey rectangles joining them. The displayed graph is superimposed on the initial network for clarity.

Figure 3 :

 3 Figure 3: Motivational example. Optimal paths for different values of 𝛼 in a schematic example. Porous space is shown in white between gray solid obstacles. Throats are shown in blue. Multiple paths link the sources and destination pores, shown in red. Throat size values are given in white, next to the relevant throats. For 𝛼 = 0 the path shown by full black line is optimal, for 𝛼 = 1 the path shown by the thinly dashed path is optimal, for 𝛼 = 2 the path shown by thickly dashed line is optimal.

Figure 4 :

 4 Figure 4: Velocity maps with overlaid active/inactive pore classification. Measured velocity in the pores is shown as dark shade in the model. Each pore is colored and overlaid according to its activity. Active pores are shown in green while inactive pores are red. White areas of the image correspond to solid grains. Flow is from left to right.

Figure 5 :

 5 Figure 5: Characteristics of proposed paths for different throat size exponent values.

Figure 6 :

 6 Figure 6: Path matches for high and low threshold values of both experiments for different overlap and throat size exponent settings,

Figure 8 :

 8 Figure 8: Examples of path match failure. Flow is from left to right. See Figure 7 for color code.

Figure 9 :

 9 Figure 9: Simple example of γ. Two pores in black are connected by their centers (red dots) and projected onto the axis of flow, here shown as left to right.

Figure 10 :

 10 Figure 10: Statistics of paths proposed in a 2D parameter spaceInteresting 2-dimensional patterns appear. We note that the path lengths (Figure10, left) are of the same order as for the 1-parameter model. While the path lengths heatmap seems to have a global minimum in the bottom left corner for (α, β) = (0,0) it has two local maxima in the bottom right and top right corners for values of (α, β) = (6,0) and (α, β) = (6,6) respectively. In accordance with the explanation of the 1parameter model results, path length increases as the algorithm is fed more and more widely spread edge weight values in the graph (i.e. some throat weights become so large that it chooses to take a series of smaller values with a lesser contribution to the overall path sum). The (α, β) = (6,0) path length maximum is then understood as occurring because of a widely spread throat size term only, with no contribution from the angle term, explaining the average throat size rise in this zone. The (α, β) = (0,6) path length maximum is then understood as the equivalent phenomenon, in which a longer path is chosen purely on the basis on avoiding unfavorable pore-to-pore angles. The average throat sizes in (α, β) = (0,6) (Figure10b) and angle cosine in (α, β) = (6,0) (Figure10c) confirm this observation. For these parameter values, they respectively take the model average throat size and the model average cosine value (for the full model: r̅ = 51 µm; cos(γ) ̅̅̅̅̅̅̅̅ = 0.63). What is not easily understood is the decrease in path length across the α = 6 parameter series, before the rise to the (α, β) = (6,6) maximum. Increasing β with a fixed α should create a larger weight spread and lead to longer paths. We refer the reader to Appendix C for this analysis.

Figure 11 :

 11 Figure 11: Path matches for 2-parameter model combinations for both experiments and thresholds.

Figure 12 :

 12 Figure 12: Comparison of path prediction performance for experiment 1 for the best fit model, the shortest number of pores model, and the shortest distance model.

Figure 13 :

 13 Figure 13: Aspect ratio H/W distribution (left) and distribution of hydraulic resistance (right).

Figure 14 :

 14 Figure 14: Average velocity intensity for all paths in each setting. The best model is the same as shown in the active/inactive high threshold pore framework.

Figure 15 .

 15 Figure 15. Bubble displacement distribution (top) and bubble track deviation angle from pressure gradient (bottom) for both experiments.The deviation angle distribution confirms microscopically the different model fits for each experiment. While both angle probability distributions are symmetric and centered on zero, the faster flowing foam (experiment 2) is more spiked in zero, showing a larger tendency for the bubbles to remain in parallel tracks to flow.

Figure 16 :

 16 Figure 16: Comparison of flow maps for two areas of interest for each experiment. We observe flow path deviation perpendicular to pressure gradient in experiment 2, while access to these zones is denied by the presence of other flow paths in Experiment 1. We circle in green the exact access links that are crucial in keeping the flow in the straighter paths for the high-velocity experiment. Flow is from left to right in both examples. Note: colormap contrast has been forced in each example for clarity.

Table 1 :

 1 Proportion of active pores for each experiment and threshold value

	Exp. 1	Exp. 2

Appendix A: Porous network and extracted graph

Regarding the structural parameters used for the weights of the edges connecting the inlet/outlet nodes to the rest of the network; the throat size for each edge was trivially the throat that connected the network to the inlet/outlet zone. However, the angle γ that shows the pore-to-pore alignment with pressure gradient was set equal to 0 (perfect alignment with gradient) for all nodes connected to either to inlet or outlet nodes. Indeed, as the center of mass of the inlet and outlet nodes are positioned in the central axis, taking the alignment angle with this position would create erroneous bias against pores outwards from the central axis.

Appendix B: 𝐜𝐨𝐬(𝛄) properties and parameter independence

As the parameters are more and more spread out with the increasing exponents α or β, the path lengths along the α = 6 parameter series should therefore monotonically rise to the (α, β) = (6,6) maximum. This discrepancy reveals either a hidden relationship between the two supposed independent parameters r and cos(γ), or otherwise an oversight on the authors behalf in the way the weight spread is interpreted to contribute to path length.

We first explore the parameter independence. We display in Figure B1 The histogram shows one-sided distribution for the cos(γ) values. This is a natural consequence of the cosine function, and a similar distribution occurs when plotting the cosine of a random sequence of angles from 0-90°. However, a noticeable spike occurs around the 0.7 value. The boxplot shows a median throat value that raises slightly in the central quartile and shows a smaller data spread towards the central values of cos(γ). We believe this behavior, as well the spike occurring around 0.7 for the cosine histogram are consequences of the way in which the watershed algorithm establishes throats and pores from the binarized model but does not represent a significant dependence between the two structural parameters used.

Appendix C: Edge weight characterization with respect to parameter combinations in the 2-parameter model

Different combinations of parameters produce different edge weight distributions for the graph input to the path-proposing algorithm. Properties of this distribution are believed to be linked to the path length, we explore here the edge weight distributions and how they affect the proposed paths. Values of the throat size term, r take values all above 1, displaying a roughly gaussian distribution [START_REF] Yeates | New insights of foam flow dynamics in high-complexity 2D micromodels[END_REF]). The values of ( We observe in Figure C1a that increasing the β parameter spreads the distribution in the positive direction, while maintaining many values at w = 1. Raising the α parameter as shown in Figure C1b spreads the distribution in both direction and shifts it towards smaller values of w. As the parameters are independent, increasing both in combination therefore combines both effects without any further interaction; it spreads the distribution both ways and shifts it to the left, but creates a more significant spread on the right-hand tail due to the raised β parameter. As we can see, both the variance σ² and mean μ of edge weight distributions change with each parameter combination. We show the parameter distribution skew in Figure C2. We suspect that the non-monotony of the path lengths function visible on the α = 6 parameter series is related to the skewness of the distribution. Indeed, a decrease in the calculated skew of the edge weight distribution occurs where the path lengths dips along the α = 6 series. This indicates an interaction between the parameters in the model that explains why the proposed paths match experimental data relatively well in a large zone diagonally situated in the positive α and β directions. The edge weight distribution is significantly less skewed in this region, which leads the path proposing algorithm to output shorter paths that closer matches to the preferential paths observed in the foam experiments.


Elise Lorenceau

Florence Rouyer
email: florence.rouyer@u-pem.fr

Lifetime of a single bubble on the surface of a water and ethanol bath

The lifetime of a single bubble on the surface of an alcohol-water liquid bath is measured at intermediate humidity rates. Thanks to an automatic bubble generation device associated with a lifetime measuring device, we carry out a statistical study of the bubble lifetime, which extends over three decades, from 0.1 to 100 s. Up to alcohol concentrations of 12% (g/g), the average bubble lifetime increases with the alcohol concentration. In addition, normalized lifetime probability density functions are independent of the alcohol concentration and relatively well described by a Weibull distribution with an exponent of 4/3. We compare our experimental data with the predictions of a recent model highlighting the importance of Marangoni stresses induced by water contaminants to the fate of a bubble. In our case, these Marangoni stresses are probably induced by heterogeneities in the alcohol concentration maintained by evaporation.

I. INTRODUCTION

A bubble that bursts at the air/liquid interface typically ejects hundreds of air drops due to the disintegration of the thin hemispherical liquid film separating the interior of the bubble from the atmosphere. This phenomenon has important consequences in volcanology, for example. Since the magmatic processes that govern explosive volcanic eruptions are impossible to observe directly, indirect observations-of the average size or size distribution of the ejected droplets, for exampleprovide a better understanding of the underground magmatic activity [START_REF] Nguyen | Film drainage and lifetime of bubbles[END_REF]. In addition, the droplets ejected when bubbles burst at the sea surface are also a key factor in the study of air-sea interactions. Indeed, they contribute significantly to the evaporation of seawater [START_REF] Andreas | The spray contribution to net evaporation from the sea: A review of recent progress[END_REF] and are among the main sources of sea salt aerosol. By efficiently scattering solar radiation and serving as nuclei for cloud formation, these aerosols are particularly relevant to understanding radiative forcing of the climate. [START_REF] Murphy | Influence of sea-salt on aerosol radiative properties in the Southern Ocean marine boundary layer[END_REF][START_REF] Schwartz | The Whitehouse effect-Shortwave radiative forcing of climate by anthropogenic aerosols: An overview[END_REF].

In this context, the lifetime of a bubble on the surface of a liquid is an important parameter: the longer the lifetime, the more the film is drained and the finer the aerosol formed. Overall, the life of a bubble at an interface is the sum of the characteristic times of drainage, evaporation, hole initiation, and bursting [START_REF] Lhuissier | Bursting bubble aerosols[END_REF][START_REF] Rio | Thermodynamic and mechanical timescales involved in foam film rupture and liquid foam coalescence[END_REF][START_REF] Poulain | Ageing and burst of surface bubbles[END_REF][START_REF] Miguet | Stability of big surface bubbles: Impact of evaporation and bubbles sizes[END_REF][START_REF] Poulain | Biosurfactants Change the Thinning of Contaminated Bubbles at Bacteria-Laden Water Interfaces[END_REF]. It therefore subtly depends on various parameters such as the size of the bubble, the viscosity of the liquid, and the presence of surfactant molecules [START_REF] Nguyen | Film drainage and lifetime of bubbles[END_REF][START_REF] Lhuissier | Bursting bubble aerosols[END_REF][START_REF] Poulain | Ageing and burst of surface bubbles[END_REF][START_REF] Poulain | Biosurfactants Change the Thinning of Contaminated Bubbles at Bacteria-Laden Water Interfaces[END_REF][START_REF] Kočárková | Film drainage of viscous liquid on top of bare bubble: Influence of the Bond number[END_REF][START_REF] Debregeas | The life and death of bare viscous bubbles[END_REF][START_REF] Modini | Effect of soluble surfactant on bubble persistence and bubble-produced aerosol particles[END_REF][START_REF] Champougny | Life and death of not so bare bubbles[END_REF]. Indeed, by modifying the structure of the interface, the presence of surfactant molecules can increase the lifetime of air bubbles by several orders of magnitude [START_REF] Lhuissier | Bursting bubble aerosols[END_REF][START_REF] Poulain | Biosurfactants Change the Thinning of Contaminated Bubbles at Bacteria-Laden Water Interfaces[END_REF][START_REF] Champougny | Life and death of not so bare bubbles[END_REF][START_REF] Ross | Foaminess of binary and ternary solutions sydney[END_REF]. The lifetime of bubbles generated from a liquid consisting of a binary solution can also be very long as observed for saline solutions [START_REF] Bikerman | Persistence of bubbles on inorganic salt solutions[END_REF][START_REF] Burger | The persistence of air bubbles at a seawater surface[END_REF][START_REF] Struthwolf | The residence time of air bubbles less than 400 μm diameter at the surface of distilled water and seawater[END_REF], water-alcohol mixtures [START_REF] Poulain | Ageing and burst of surface bubbles[END_REF][START_REF] Samanta | Coalescence of air bubbles in aqueous solutions of alcohols and nonionic surfactants[END_REF][START_REF] Oolman | Bubble coalescence in stragnant liquids[END_REF], oil [START_REF] Chandran Suja | Evaporation-induced foam stabilization in lubricating oils[END_REF][START_REF] Arangalage | Asphaltenes at oil/gas interfaces: Foamability even with no significant surface activity[END_REF], and molten glass [START_REF] Kočárková | Film drainage of viscous liquid on top of bare bubble: Influence of the Bond number[END_REF][START_REF] Pigeonneau | Stability of vertical films of molten glass due to evaporation[END_REF]. For these binary systems, the bubble stabilization mechanism at the liquid/air interface is less well known. In the case of water/alcohol mixtures, the increase in stability may be due to adsorption of the alcohol at the interface. The high concentration of -OH dipoles at the air/water interface would induce repulsive van der Waals forces between the interfaces, leading to a disjoining pressure able to balance capillary and gravity drainage of the film [START_REF] Samanta | Coalescence of air bubbles in aqueous solutions of alcohols and nonionic surfactants[END_REF][START_REF] Oolman | Bubble coalescence in stragnant liquids[END_REF][START_REF] Drogaris | Coalescence behavior of gas bubbles in aqueous solutions of n-alcohols and fatty acids[END_REF][START_REF] Wang | Effects of surface forces and film elasticity on foam stability[END_REF]. This positive adsorption of alcohol has been demonstrated indirectly in experiments on adsorption isotherms and directly in molecular dynamics simulations for methanol and ethanol [START_REF] Kipling | Composition of the surface film on aqueous solutions of aliphatic alcohols[END_REF][START_REF] Wilson | Adsorption and solvation of ethanol at the water liquid-vapor interface: A molecular dynamics study[END_REF]. One could also think of a modification of the surface tension or surface rheology due to the presence of alcohol: Since regions of high surface tension exert higher tangential stresses on the air/liquid interface than regions of low surface tension, a surface tension gradient induces a flow of liquid away from the regions of low surface tension. These Marangoni stresses can even partially counteract gravity or capillary stress, thus stabilizing liquid film [START_REF] Pigeonneau | Stability of vertical films of molten glass due to evaporation[END_REF], single bubbles [START_REF] Ghabache | Surface libre hors équilibre : de l'effondrement de cavité aux jets étirés[END_REF][START_REF] Yaminsky | Bubble vortex at surfaces of evaporating liquids[END_REF], or foam [START_REF] Tuinier | Transient foaming behavior of aqueous alcohol solutions as related to their dilational surface properties[END_REF]. In this very comprehensive framework, one question remains open: Through which mechanisms common to all these systems are these Marangoni stresses maintained? The answer probably lies in the difference in volatility between the two compounds of the binary system, as recently shown for multicomponent lubricating oils [START_REF] Chandran Suja | Evaporation-induced foam stabilization in lubricating oils[END_REF]. The overall mechanism can be easily understood within the geometry of a single bubble on the surface of an alcohol/water mixture. The alcohol evaporates everywhere at the same rate. At the foot of the bubble, this evaporation is counteracted by the proximity of the reservoir that replenishes the liquid film, thus ensuring that the temperature and alcohol concentration in the film are equal to those in volume. However, at the top of the bubble, the liquid reservoir is not directly available, which induces a lower alcohol concentration and possibly a lower temperature than in volume. As the surface tension of the alcohol-water mixture decreases with the alcohol concentration or temperature, this creates a surface tension gradient that pulls the liquid upward, partially neutralizing the thinning effects of gravity and capillarity [START_REF] Poulain | Ageing and burst of surface bubbles[END_REF][START_REF] Chandran Suja | Evaporation-induced foam stabilization in lubricating oils[END_REF]. Despite this very convincing mechanism, correlating the importance of these Marangoni flows with bubble life is a delicate task. Indeed, for high-energy surfaces such as water/alcohol solutions, bubble life is the result of both deterministic processes such as gravity, capillarity, and Marangoni stresses and stochastic processes such as the presence of impurities in the air or the solution, which can trigger film rupture [START_REF] Poulain | Ageing and burst of surface bubbles[END_REF]. To test the validity of this deterministic/stochastic coupling, it is necessary to work on large assemblies of bubbles and measure their lifetime distribution. This has been done recently for oils [START_REF] Chandran Suja | Evaporation-induced foam stabilization in lubricating oils[END_REF][START_REF] Arangalage | Asphaltenes at oil/gas interfaces: Foamability even with no significant surface activity[END_REF], hot or cold clean water-for which thermal rather than solutal Marangoni flows may appear-and salt solutions [START_REF] Poulain | Ageing and burst of surface bubbles[END_REF]. However, the case of alcoholic binary solutions has been less studied [START_REF] Poulain | Ageing and burst of surface bubbles[END_REF][START_REF] Samanta | Coalescence of air bubbles in aqueous solutions of alcohols and nonionic surfactants[END_REF][START_REF] Ghabache | Surface libre hors équilibre : de l'effondrement de cavité aux jets étirés[END_REF]. In particular, the evolution of the bubble lifetime distribution with the alcohol concentration has never been considered.

To fill this gap and to try to quantify the relative importance between stabilization induced by alcoholic adsorption at the interface or by Marangoni flow, we measure the lifetime of individual bubbles on the surface of an ethanol/water liquid bath in an atmosphere of controlled composition. To obtain sufficient statistics, we use a device that allows us to repeat the same experiment more than a hundred times. The average lifetime and bubble lifetime distribution are presented and analyzed for various concentrations of ethanol in the reservoir.

II. MATERIALS AND METHODS

A. Materials

The container is filled with a binary mixture of purified water (conductivity, 18.2 M /cm) and ethanol (Sigma Aldrich; 96% pure). The ethanol concentration varies from 0 to 12% (g/g). Within this range of ethanol concentrations, the physical properties of the mixture can vary significantly. Therefore, the density ρ, the surface tension γ , the dynamic viscosity η, and the diffusion coefficient D of the alcohol-water mixture are listed in Table I.

To form a bubble, ambient air is pumped by an aquarium pump operated by an automated electronic system, then blown through a vertical capillary tube with a diameter of 300 μm immersed a few centimeters below the air/liquid interface (see Fig. 1). To ensure proper alignment, the 063603-2 TABLE I. Liquid density, liquid/air surface tension, dynamical viscosity, diffusion coefficient, radius of the bubble cap, and number of bubbles for different percentages of ethanol. The liquid density, surface tension, and viscosity are derived from [START_REF] Khattab | Density, viscosity, and surface tension of water plus ethanol mixtures from 293 to 323K[END_REF], and the diffusion coefficient from [START_REF] Harris | Mutual diffusion-coefficients for the systems water-ethanol and water propan-1-op at 25 °C[END_REF] capillary tube is placed along the axis of symmetry of the 4.5-cm-wide cylindrical vessel. Once detached from the cylindrical vessel, the bubble rises up and stabilizes at the center of the air/liquid interface. Indeed, the surface is slightly convex thanks to a hydrophobic Teflon ring placed at the top of the container, which imposes a contact angle greater than 90 • along the walls. Unless otherwise specified, the experiments are performed at 24 • C ± 2 • C in an atmosphere with an ethanol partial pressure and a water partial pressure of 1365 ± 45 Pa, corresponding to a humidity level of 45.5% ± 1.5% (measured with a TESTO 625 probe).

FIG. 1. Setup used to form bubbles and to measure their lifetimes on the surface of the liquid bath. The dashed red arrow corresponds to the reflection of the laser beam when a bubble is at the interface, while the solid red arrow corresponds to the reflection of the laser beam when there is no bubble at the interface.

063603-3

B. Automatic bubble generation and lifetime measurement

Bubble lifetime on the surface of the liquid bath is measured by means of a laser beam inclined approximately 45 • and directed towards the center of the cylindrical container. If the liquid/air interface is flat and undistorted, the laser beam is reflected by the liquid interface and its intensity is recorded using a photodiode (see Fig. 1). If the liquid/air interface is distorted by the presence of a bubble, then the laser reflection is deflected and the photodiode does not record any light intensity. Once the bubble has burst, the liquid/air interface returns to its equilibrium position and the laser beam is reflected back to the photodiode, which stops recording the time. Thus, the lifetime of a bubble corresponds to the time interval during which there is no signal at the output of the photodiode. It is measured using a universal counter triggered by the photodiode signal switch. Five seconds after the bubble burst, the universal counter triggers the formation of a new bubble and the process is repeated. Measured in this way, the typical lifetime of a bubble t is determined with an accuracy of ±0.2 s, corresponding to the damping of the oscillations induced by the bursting of the bubble. It typically varies from a few milliseconds to 100 s.

The radius at the apex of the bubbles obtained with this setup varies from 2.8 to 3.5 mm. As the shape of the bubble results from an equilibrium between gravity and capillary forces [START_REF] Princen | Shape of a fluid drop at a liquid-liquid interface[END_REF], we introduce the capillary length l c = √ γ /(ρg), which varies from 2.6 to 2.2 mm. The dimensionless parameter R/l c ≈ 1.31 ± 0.04 as reported in Table II. This corresponds to a Bond number Bo = R 2 /l 2 c = ρgR 2 /γ ∼ 1.6 of the order of unity as reported in Tables I andII. It is reproducible within 5% tolerance due to slight variations in the ethanol concentration induced by evaporation.

This device can therefore generate a large number of data as long as the external experimental conditions (such as the liquid bath level, ambient atmosphere, and temperature or humidity level) do not change. The number of bubbles generated for each alcohol concentration, listed in Table I, varies from 250 to 1250.

III. RESULTS AND DISCUSSION

Figure 2 presents typical raw data on the bubble life in a 12% (g/g) ethanol-water bath. For this assembly of 340 bubbles, the lifetime varies over three decades, from 0.1 to 100 s. To quantify these data, we successively consider the average bubble lifetime, t b , as well as the probability density function (PDF) of observing a bubble with a given lifetime.

A. Average lifetime

First, we check that the average lifetime of the bubble t b does not change over time. To do this, we average over dozens of bubbles rather than over the entire range: plotted against the number of bubbles, the values obtained do not show any significant variation. Then we calculate t b over the entire range. For the data presented in Fig. 2,t The average life of the bubbles, plotted as a function of the ethanol concentration in Fig. 3, increases with the concentration of ethanol and is well described by a linear (empirical) adjustment passing through 0. This is not surprising: at 0% ethanol, the film of pure water at the top of the bubble being in no way stabilized, the bubble explodes in a few milliseconds, which is much too short a duration for our experimental automatic device and negligible compared to the average lifetime measured at higher ethanol concentrations. We were unable to work at an ethanol concentration greater than 12%: above this value, the contact angle of the binary mixture on the vertical wall is smaller, thus reducing the amplitude of the meniscus and its ability to stabilize the bubble in the center of the tank. We also carry out the same experiment but in the presence of a vapor saturated with ethanol and water. In the absence of evaporation, the bursting is quasi-immediate and the response time in our setup (of the order of 0.2 s) is too large to quantitatively probe this effect. All these results are consistent with a stabilization mechanism of bubbles induced by Marangoni effects and by the evaporation of ethanol as presented in Sec. I. Indeed, if the vapor is not saturated with ethanol, evaporation generates concentration heterogeneities at the bubble scale. As the surface tension of water/ethanol mixtures decreases with increasing ethanol concentration, these concentration heterogeneities generate surface tension gradients that can drag the liquid upwards. These Marangoni stresses can even counteract the gravitational and capillary thinning of the film, making it less likely to explode. However, if the vapor is saturated with ethanol, the ethanol cannot evaporate and the bubble bursts rapidly. This stabilizing effect observed over a range of ethanol concentrations from 0 to 12% (g/g) is also observed in [START_REF] Tuinier | Transient foaming behavior of aqueous alcohol solutions as related to their dilational surface properties[END_REF]. However, while we observe increasing stability with ethanol concentration up to 12% (g/g), Tuinier et al. observe maximum foamability for a concentration of 5%-7% (g/g) ethanol in water. However, comparison between the two experimental works cannot be made easily. In Tuinier and coworkers' experiment, the most fragile bubbles, located at the top of a foam column, are far away from the ethanol tank. Moreover, as the system is closed, unlike here, it is possible that the vapor quickly saturates with ethanol, thus reducing evaporation. The two main ingredients of the stabilization mechanism we describe in this paper, namely, evaporation and ethanol concentration gradients, are therefore probably limited in the geometry of the Tuinier et al. experiment.

B. Probability density function (PDF) of the bubble lifetime

In Fig. 4, the PDFs for the minimum and maximum ethanol concentrations are plotted on a lin-lin scale. The PDF of a random positive variable x is often described using a Weibull distribution, p(x) = β/x 0 (x/x 0 ) β-1 exp(-(x/x 0 ) β ), where β > 0 is the shape parameter and x 0 > 0 is the scale parameter. In the case where x is a value to failure, the failure rate and β are related: β = 1 corresponds to a constant failure rate for which the distribution is a monotonously decreasing exponential function, while β = 2 corresponds to an increasing failure rate for which the distribution is a Rayleigh distribution, also often used to describe statistics of a positive random variable based on diffusion, p(t ) = π/2(t/t 2 b ) exp(-π/4(t/t b ) 2 ). In the problem considered here, the bursting of the bubbles is triggered by the appearance of thinner areas of the film. By carefully monitoring the rate of film thinning and assuming that the PDF of the bubble lifetime is related to the frequency of appearance of these independent fragile parts, Lhuissier showed that the distribution of the lifetime of unpurified water bubbles is well described by a Weibull distribution with β = 4/3, given by q(t ) = 4 3

t 1/3 t 4/3 0 exp - t t 0 4/3 , ( 1
)
where t 0 is proportional to t b , t 0 = t b / (1 + 3/4) ∼ t b /0.92, with the gamma function [START_REF] Gilet | How long a bubble will be[END_REF]. And, indeed, we observe a reasonable agreement between Eq. ( 1) and the two sets of data in Fig 4 .  In Fig. 5, the normalized PDFs obtained for all ethanol concentrations are plotted on a lin/log scale as a function of t/t b , using for t b the previously determined average lifetime. We observe a collapse of the experimental PDFs. Yet for t/t b > 3, the data are very scattered. However, for this normalized time range, the number of data is low. To qualify the importance of the least statistically significant data, two types of symbols of the same color are used for each alcohol concentration. Filled symbols are data for which more than 10 events are counted and open symbols are those 063603-6 FIG. 4. Probability density function (PDF) of experimental bubble lifetimes for 2.5% (red diamonds; over 1251 bubbles) and 12% (blue circles; over 312 bubbles) alcohol in water (g/g). Dotted lines are Weibull distributions [START_REF] Lhuissier | Bursting bubble aerosols[END_REF][START_REF] Gilet | How long a bubble will be[END_REF], Eq. ( 1). Inset: Same data plotted on log-lin axis.

for which fewer than 10 events are counted. For all ethanol concentrations, we lack statistics as soon as t/t b > 3. These normalized experimental PDFs are now compared with the Weibull PDFs described previously. First, the exponential function (i.e., the Weibull distribution with β = 1) which corresponds to random external events causing a break and constant over time does not allow us to capture the increase in probability at low lifetimes. We also compare in Fig. 5 normalized PDFs along with Rayleigh distributions as done in [START_REF] Chandran Suja | Evaporation-induced foam stabilization in lubricating oils[END_REF] and [START_REF] Zheng | Laboratory measurements of water surface bubble life time[END_REF]. A Rayleigh distribution (i.e., Weibull distribution with β = 2) also fits our data correctly for t/t b < 2 but seems to fail to catch the tail of the distribution. Finally, a good agreement with Weibull's law with β = 4/3 [see Eq. ( 1)] is observed without any other adjustable parameter. This fit is particularly good up to t/t b <3, which is precisely the range of data for which we have a sufficient statistic. In addition, since PDFs seem to be independent of the alcohol concentration, we merge all the data in the same graph (see inset in Fig. 5). This increases the statistics for the tail of the distribution and confirms that the longterm distribution is better described by a Weibull function than by a Rayleigh function. Our results therefore suggest that the function's distribution of the bubble lifetime is more stretched than the Rayleigh function based solely on diffusive processes and narrower than the exponential function based solely on the time-independent failure rate.

C. Discussion

During the life of a bubble, liquid film thins inexorably until bursting. This thinning dynamics can be controlled either by the evaporation rate of the liquid (water being the main component of these films) or by the drainage of the liquid due to various mechanical stresses (such as gravity or capillary stresses).

It is presumably the composition of the liquid constituting the films and of the atmosphere as well as the bubble size that favor one or the other of the scenarios. For example, films of water bubbles at temperatures from 15 • C to 77 • C in unsaturated atmospheres-and thus subject to evaporation-are always thicker than those of bubbles in saturated atmospheres. This counterintuitive observation 063603-7 FIG. 5. Probability density function (PDF) of bubble lifetimes normalized by average lifetimes for five concentrations of ethanol in water, varying from 2.5% to 12% (g/g) (see legend). Two types of symbols of the same color are used for each concentration. Filled symbols are data for which more than 10 events are counted and open symbols are those for which fewer than 10 events are counted. The solid line represents the Weibull function [START_REF] Lhuissier | Bursting bubble aerosols[END_REF][START_REF] Gilet | How long a bubble will be[END_REF], and the dashed line the Rayleigh function. Inset: PDF for all normalized lifetimes-t/t bwhatever the ethanol concentration.

shows that for this system the mass loss induced by evaporation is negligible on the time scale of the bubble life [START_REF] Poulain | Ageing and burst of surface bubbles[END_REF]. On the contrary, films of bubbles stabilized by amphiphilic molecules adsorbed at the film interface are barely subject to drainage since, in this case, Marangoni stresses balance the gravitational or capillary stresses. This induces much longer lifetimes, for which mass loss induced by evaporation is no longer negligible [START_REF] Miguet | Stability of big surface bubbles: Impact of evaporation and bubbles sizes[END_REF][START_REF] Poulain | Biosurfactants Change the Thinning of Contaminated Bubbles at Bacteria-Laden Water Interfaces[END_REF]. Using the framework developed in [START_REF] Miguet | Stability of big surface bubbles: Impact of evaporation and bubbles sizes[END_REF] and [START_REF] Poulain | Biosurfactants Change the Thinning of Contaminated Bubbles at Bacteria-Laden Water Interfaces[END_REF], assuming that this is the liquid evaporation that limits the maximum lifetime accessible to very stable bubbles [i.e., evaporation of a film of critical thickness for which both film thinning phenomena induced by capillary drainage (lower and lower over time) and evaporation are of the same order (constant over time)], one can estimate the bubble lifetime

t evap b ≈ R(( ρ J ) 3 ( ηR γ l c ) 2 ) 1/5
, where J is the liquid mass evaporation rate.

Taking pure liquid ethanol evaporation rate values for ranges of intermediate humidity [START_REF] Spedding | Abnormal evaporation rate of ethanol from low concentration aqueous solutions[END_REF], as well as the bubble size R, surface tension γ , viscosity η, and liquid density ρ as listed in Table I, we find t evap b at a maximum of 10 s of the same order of magnitude as our observations. Yet, this value decreases 20% as the ethanol concentration increases from 1% to 12% (g/g). Film thinning induced only by evaporation is therefore in disagreement with our data.

Thus, for the bubbles we are considering, film thinning is likely to be controlled by liquid drainage. In this case, the addition of ethanol to the aqueous solution constituting the film can modify this drainage dynamic. Indeed, ethanol can adsorb at the interface and modify various surface parameters such as the surface viscosity, which would slow down the drainage of the film [START_REF] Seiwert | Velocity Field in a Vertical Foam Film[END_REF]. It can also induce solutal or thermal Marangoni stresses that would cause the liquid in the bath to rise up, thus slowing down the thinning rate due to capillary or gravity stresses [START_REF] Poulain | Ageing and burst of surface bubbles[END_REF]. To stabilize 063603-8 the bubble for a long enough time, these stresses must be sustained either by evaporation or by a temperature difference between the atmosphere and the bulk.

Stabilization induced by adsorption of surface active agents should not be sensitive to the composition of the vapor phase. However, when the experiment is carried out in an atmosphere saturated with alcohol-which should allow a maximum concentration of molecules adsorbed on the surface-the stabilizing effect disappears and all the bubbles burst in less than a second. This mechanism is therefore probably not dominant here, and in the rest of this paper, we assume that the thinning dynamics of films is governed by the generalized drainage model proposed in [START_REF] Lhuissier | Bursting bubble aerosols[END_REF] and [START_REF] Poulain | Ageing and burst of surface bubbles[END_REF].

To go one step further and understand what governs the death of the bubbles, we now recall the different pictures proposed in the literature. A classic deterministic scenario suggests that films thin out until they are thin enough for van der Waals forces or temperature fluctuations to break down. However, at the time of rupture, experimental measurements of film thicknesses are found to be highly dispersed, ranging from a few nanometers to a micrometer, suggesting a stochastic scenario. Recently, a very comprehensive framework has been proposed to understand what governs the death of a bubble [START_REF] Lhuissier | Bursting bubble aerosols[END_REF][START_REF] Poulain | Ageing and burst of surface bubbles[END_REF][START_REF] Poulain | Biosurfactants Change the Thinning of Contaminated Bubbles at Bacteria-Laden Water Interfaces[END_REF]. Stochastic local perturbations arising from intrinsic perturbations such as impurities from the solution or extrinsic ones such as impacts by mist droplets are responsible for the bubble's death. The authors show that the ability of these perturbations to kill the bubble are not all equivalent: they are highly efficient bubble killers upon reaching the interface at the bubble foot and when coupled with contaminants inducing Marangoni stresses large enough to trigger self-sustained lethal thinning.

Thus, by coupling a deterministic generalized bubble cap drainage model accounting for both curvature-pressure-induced drainage and Marangoni flows and the time scales involved in the sideto-side diffusion with the time of residence of the stochastic contaminant on the bubble, the authors of [START_REF] Lhuissier | Bursting bubble aerosols[END_REF] and [START_REF] Poulain | Ageing and burst of surface bubbles[END_REF] show (i) that the adequacy between the measured PDF and a Weibull distribution is theoretically justified and (ii) that the bubble lifetime reads as

t b = ατ Sc R l c 1/2 , ( 2
)
where τ = ηl c γ is the viscocapillary time and Sc the Schmidt number, Sc = η ρD , with D the diffusion coefficient of the impurities in the aqueous solution. In the following we assume that this diffusion coefficient is equal to the mutual diffusion of the ethanol-water mixture (cf. Table I) because in the present case a local variation of the ethanol concentration would affect the surface tension of ethanol-water, similarly to impurities at the surface of the water. Typical values of τ and Sc are listed in Table I. The prefactor α is a numerical coefficient accounting for the efficiency of Marangoni stresses to counteract gravitocapillary drainage, which reads

α = 1 1/3 -kBo , ( 3
)
where k is a numerical coefficient accounting for the reduction in the thinning rate of the films due to Marangoni stresses. For a hot film in cold air, k has been found to be of the order of its maximal value, k ∼ 0.05 < 1/3Bo, whatever the value of the temperature difference [START_REF] Poulain | Ageing and burst of surface bubbles[END_REF] for Bo ∼ 4.

To determine whether water/ethanol thin film stability can be described using this picture, we plot t b as a function of τ Sc(R/l c ) 1/2 in Fig. 6. Although Bo and ρ do not vary significantly, γ , η, D, and thus τ and Sc (whose values are listed in Tables I andII) all vary significantly with the ethanol contents in the probed range. The data in Fig. 6 are reasonably well described by a linear variation with a slope α of the order of 350. Note that a plot of t b versus τ shows a strong nonlinear variation; instead, t b ∝ τ 3.5 . The difference between the linear law observed between t b and τ Sc(R/l c ) 1/2 and a power law between t b and τ with an exponent greater than 3 may seem surprising. However, we recall that D ∝ 1/η, thus Sc ∝ η 2 and τ Sc(R/l c ) 1/2 ∝ η 3 ∼ τ n , with n ∼ 3. Using Eqs. consistent with the maximal value predicted in [START_REF] Poulain | Ageing and burst of surface bubbles[END_REF] and k ≈ 0.2. Our work thus reaffirms that even if eventually undergoing stochastic rupture, the lifetime of bubbles stabilized by upward Marangoni stresses (solutal due to the evaporation of ethanol in this case) rather than by surfactants can be modeled by an equilibrium between the diffusion time and the residence time of the contaminant at the foot of the bubble. The film at the foot of the bubble is very thin and subject to violent evaporation of ethanol in the present case. This very particular geometric configuration (confined film and evaporation) could be responsible for the stochastic appearance of large local concentration gradients.

The normalized PDF shown in the inset in Fig. 5 also agrees very well with the prediction from [START_REF] Lhuissier | Bursting bubble aerosols[END_REF] and [START_REF] Poulain | Ageing and burst of surface bubbles[END_REF]. Indeed, when the bubble lifetime is normalized by the mean bubble lifetime, all the data are observed to collapse on a master curve perfectly fitted by Eq. (4). A small deviation occurs for the data with the longest lifetime (t/t b > 4), which are a bit overrepresented. This deviation remains unexplained.

IV. CONCLUSION

A statistical study of thousands of bubbles shows that the average lifetime of bubbles on the surface of an ethanol-water bath increases considerably with the concentration of ethanol. Our data are consistent with the framework in [START_REF] Lhuissier | Bursting bubble aerosols[END_REF] and [START_REF] Poulain | Ageing and burst of surface bubbles[END_REF]: This increase can be understood by considering the ethanol concentration heterogeneities induced by evaporation, which in turn create upward Marangoni stresses. In order to validate this picture, several hypotheses remain to be validated, notably that significant ethanol concentration gradients may appear on the bottom of the bubble and that the stochastic appearance of these gradients are due to the evaporation of the very thin film of liquid. In addition, several questions remain unanswered, in particular, those concerning the 063603-10 thickening of the film on top of the bubble due to the evaporation of the surface-active species as observed on vertical films [START_REF] Kočárková | Film drainage of viscous liquid on top of bare bubble: Influence of the Bond number[END_REF][START_REF] Chandran Suja | Evaporation-induced foam stabilization in lubricating oils[END_REF][START_REF] Fournier | Tears of wine[END_REF] and on the surface of the bubbles [START_REF] Poulain | Ageing and burst of surface bubbles[END_REF].

  FIG.3. Average bubble lifetime, t b , as a function of the percentage (g/g) ethanol in water. The error bar corresponds to the standard deviation. The linear fit is a guide for the eyes.

9 FIG. 6 .

 96 FIG. 6. Average bubble lifetime as a function of τSc( R lc ) 1/2 .

  .

	% ethanol	ρ	γ	η	D	R	N bubbles
	(g/g)	(kg/L)	(mN/m)	(mPa s)	(m 2 /s)	(mm)
	0	0.997	71.1	0.89	1.22 × 10 -9	-	-
	1	0.996	65.9	0.92	1.19 ×10 -9	3.5	250
	2.5	0.993	60.9	0.97	1.14 × 10 -9	3.3	1250
	5	0.989	56.1	1.04	1.07 × 10 -9	3.1	1000
	7	0.986	52.7	1.11	1.02 × 10 -9	3.0	250
	10	0.986	49.1	1.21	0.94 × 10 -9	3.0	250
	12	0.979	46.5	1.28	0.89 × 10 -9	2.8	500

TABLE II .

 II Capillary length l c , viscous time τ , and Schmidt (Sc) and Bond (Bo) numbers calculated from the values in TableIfor different percentages of ethanol.

	% ethanol	l c	τ	Bo	Sc
	(g/g)	(mm)	(10 -6 s)
	0	2.70	33.8	1.68	735
	1	2.60	36.3	1.82	779
	2.5	2.50	39.7	1.79	850
	5	2.40	44.7	1.65	982
	7	2.33	49.1	1.65	1102
	10	2.25	55.6	1.76	1306
	12	2.20	60.5	1.62	1459

ACKNOWLEDGMENTS

We acknowledge the technical team of Laboratoire Navier for realizing the setup. We thank P. Willay and D. Dessalles for conducting experiments. We thank E. Rio, J. Miguet, Q. Magdelaine, A. Antkowiak, O. Atasi, and B. Scheid for fruitful and stimulating discussions. We gratefully acknowledge financial support from IFSTTAR and ENPC for internship remuneration.