Annika Vilmi

Gibert Corentin

Gilles Bret

Konsta Escarguel

Jani Happonen

Aurélien Heino

Sophia I Jamoneau

Félix Passy

Janne Picazo

Juliette Soininen

Tison-Rosebery

Corentin Gibert

Gilles Escarguel

Konsta Happonen

Jani Heino

Aurélien Jamoneau

Sophia I Passy

Félix Picazo

Janne Soininen

Juliette Tison-Rosebery

Jianjun Wang
email: jjwang@niglas.ac.cn

J Soininen

Dispersal-niche continuum index: a new quantitative metric for assessing the relative importance of dispersal versus niche processes in community assembly

Keywords: assembly processes, determinism, dispersal, environmental filtering, niche processes, null models, PER-SIMPER, stochasticity Dispersal-niche continuum index:

scientifiques de 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

Understanding how biological communities are assembled is one of the main goals of community ecology, macroecology and biogeography. The importance of correctly identifying the key mechanisms of community assembly increases as global change proceeds and relevant ecological information is urgently needed to help designing appropriate management and conservation programs [START_REF] Mittelbach | Ecological and evolutionary perspectives on community assembly[END_REF]. Ecological determinism and stochasticity are considered two main forces underlying community assembly [START_REF] Leibold | The metacommunity concept: a framework for multi-scale community ecology[END_REF][START_REF] Chase | Drought mediates the importance of stochastic community assembly[END_REF][START_REF] Vellend | Assessing the relative importance of neutral stochasticity in ecological communities[END_REF]. The deterministic paradigm postulates that biological communities are structured by niche properties, stemming from species' physiological responses, resource use and biotic interactions, such as competition [START_REF] Vellend | Conceptual synthesis in community ecology[END_REF]. The stochastic paradigm, by contrast, considers biological communities as products of random processes, such as dispersal [START_REF] Hubbell | The unified neutral theory of biodiversity and biogeography[END_REF]. However, it is likely that niche and dispersal processes are not mutually exclusive but constitute the end points of a continuum, similar to the continuum of niche and neutral models proposed by [START_REF] Gravel | Reconciling niche and neutrality: the continuum hypothesis[END_REF]. Within the dispersal-niche continuum, compositional variation of biological communities can be driven to a greater or lesser extent by each process [START_REF] Holyoak | A framework for large scale community ecology[END_REF][START_REF] Mouillot | Niche-assembly vs dispersal-assembly rules in coastal fish metacommunities: implications for management of biodiversity in brackish lagoons[END_REF][START_REF] Vellend | Conceptual synthesis in community ecology[END_REF][START_REF] Gibert | PER-SIMPER -a new tool for inferring community assembly processes from taxon occurrences[END_REF].

As patterns in community composition are scale-dependent [START_REF] Chase | Embracing scale-dependence to achieve a deeper understanding of biodiversity and its change across communities[END_REF]) and underlying drivers are generally difficult to distinguish [START_REF] Gotelli | Quantifying biodiversity: procedures and pitfalls in the measurement and comparison of species richness[END_REF], quantifying the main assembly processes in various systems and across different datasets has remained challenging [START_REF] Viana | Spatial scale modulates the inference of metacommunity assembly processes[END_REF]. Null model approaches are commonly used for studying whether communities differ from null-expectations and to what extent deterministic (i.e. niche-related) processes affect this deviation [START_REF] Chase | Using null models to disentangle variation in community dissimilarity from variation in a-diversity[END_REF]). Variation partitioning in constrained ordination [START_REF] Borcard | Partialling out the spatial component of ecological variation[END_REF]) is another widely implemented method. It assesses the effects of environmental conditions (a proxy for niche processes) and space (a proxy for dispersal) on community assembly, and it relies on detailed data on community composition, environmental factors and geographic locations. However, traditional variation partitioning has been criticized for statistical and ecological reasons in the context of community assembly research, emphasizing the need for further methodological innovations. Specifically, in addition to capturing the relative importance of environmental conditions and space, the results of variation partitioning can also reflect the impact of spatial arrangement of environmental conditions, making interpretations of underlying processes difficult [START_REF] Smith | Variation partitioning as a tool to distinguish between niche and neutral processes[END_REF]Lundholm 2010, Clappe et al. 2018). Simulations have suggested that variation partitioning fails to correctly represent environmental and spatial components of community variation, sometimes resulting in poor environmental models, poor spatial models or both [START_REF] Gilbert | Partitioning variation in ecological communities: do the numbers add up?[END_REF]Bennett 2010, Clappe et al. 2018).

To address some of the data and methodological limitations of previous methods, [START_REF] Gibert | PER-SIMPER -a new tool for inferring community assembly processes from taxon occurrences[END_REF] recently developed a method called PER-SIMPER. PER-SIMPER allows the qualitative identification of the dominant assembly mechanism based solely on a matrix of species occurrences across sites. Building on SIMPER [START_REF] Clarke | Non-parametric multivariate analyses of changes in community structure[END_REF], this permutation-based (hence, 'PER-') null model approach identifies the dominant assembly process accounting for compositional similarity percentages between site groups consisting of local communities (hence, 'SIMPER') within the same regional species pool [START_REF] Gibert | PER-SIMPER -a new tool for inferring community assembly processes from taxon occurrences[END_REF]. PER-SIMPER utilizes sites-by-taxa community matrices, generating three distinct null models, constraining rows ('niche assembly'), columns ('dispersal assembly') or both during permutations. It uses the SIMPER [START_REF] Clarke | Non-parametric multivariate analyses of changes in community structure[END_REF]) method to model the original community matrix compositional similarity pattern and to detect the empirical model profile, which is then compared to the three null distribution profiles (hereafter called 'null PER-SIMPER profiles', Fig. 1). PER-SIMPER makes a qualitative assessment by identifying which null profile matches best the empirical profile. However, the sensitivity of a qualitative approach to assembly mechanisms can be limited by the fact that most communities are structured by both niche and dispersal processes [START_REF] Mouillot | Niche-assembly vs dispersal-assembly rules in coastal fish metacommunities: implications for management of biodiversity in brackish lagoons[END_REF]. In addition, the qualitative nature of the PER-SIMPER decision-making process prevents precise comparisons of the assembly process in different metacommunities.

As a solution to these limitations, we here propose a new metric, the dispersal-niche continuum index (DNCI) with broad ecological application. Specifically, the PER-SIMPER analysis returns three E-metric distributions, which correspond to the deviation between the empirical SIMPER profile and the three null PER-SIMPER profiles generated from the three permutation modes of the original matrix, where each permutation mode is linked to one of the three main assembly perspectives (i.e. dispersal assembly, niche assembly or joint dispersal-niche assembly). The new DNCI proposed here is derived from these computed E-values. Quantitative identification of the assembly process is based on the subtraction of the standard effect size (SES) of E n (i.e. E-metric distribution from the 'niche' model) from the standard effect size of E d (i.e. E-metric distribution from the 'dispersal' model). Thus, DNCI offers a way to both quantify and compare the strengths of the main assembly processes across datasets (Fig. 1). Positive or negative DNCI values indicate that niche or dispersal assembly is the main process structuring the studied site groups consisting of biological communities, respectively. Higher absolute values of the index represent greater potential strength of the dominant assembly process.

To test the robustness of the DNCI, we performed robustness tests with simulations. We then calculated the DNCI for two example datasets to test the applicability of the DNCI with empirical data. Using stream monitoring data from the USA, we compared whether the DNCI values differ among metacommunities of organisms varying in trophic level and dispersal capacity, including diatoms, macroinvertebrates and fish. Using a field study dataset from China, we addressed whether the DNCI values of stream bacteria and macroinvertebrates differ from each other and whether they vary with Negative DNCI values indicate that dispersal assembly is the dominating process, while positive DNCI values indicate that niche assembly is the dominating process behind community dissimilarities among site groups. With increasing absolute values, the strength of the dominating process increases. DNCI values close to zero indicate that both niche and dispersal assembly processes affect distributions. In dataset A, the DNCI does not significantly differ from zero, indicating that both dispersal and niche processes control the empirical species distributions. In dataset B, DNCI is significantly greater than zero, indicating that niche assembly is the main driver of the empirical species distribution. In datasets C and D and the illustrated dataset, DNCI is significantly less than zero, indicating that dispersal assembly is the main driver of the empirical species distributions. Dataset D and the illustrated dataset do not have significantly different DNCI values, i.e. the control of the empirical species distributions by dispersal assembly is similar. Dataset C has a DNCI significantly higher than dataset D and the illustrated dataset, i.e. dispersal assembly control in dataset C is significantly weaker than in dataset D and the illustrated dataset. SES d = standard effect size of E d . SES n = standard effect size of E n . Modified from Gibert and Escarguel's (2019) Supporting information. spatial distance. Finally, we discussed the strengths and weaknesses of the DNCI in inferring community assembly processes and suggested potential research avenues.

Material and methods

Description of DNCI

The initial approach of PER-SIMPER [START_REF] Gibert | PER-SIMPER -a new tool for inferring community assembly processes from taxon occurrences[END_REF] compares an empirical SIMPER profile with three null PER-SIMPER profiles and qualitatively identifies which null profile and thus which assembly mode (i.e. dispersal assembly, niche assembly or both) is responsible for the composition of the analyzed communities (Fig. 1). To this end, the empirical SIMPER matrix is permuted either by maintaining the richness of localities (the row sums) in order to retain information on the potential number of ecological niches per sample locality, or by maintaining the number of localities where the taxa may have dispersed (the column sums) in order to retain information on the potential dispersal capacity of taxa. The third permutation mode maintains both row and column sums, retaining a larger part of the sampled information. While the number of localities where a species was sampled can provide indication for dispersal capacity, the taxonomic richness of a community can reflect the number of ecological niches, with richer communities having more niches [START_REF] Granot | Niche breadth and species richness: correlation strength, scale and mechanisms[END_REF]. The identification of the assembly mode relies on the E-metric, corresponding to the square of the deviations between the empirical SIMPER profile and multiple null PER-SIMPER profiles produced with the three permutation modes (by default, 1000 permutations for each assembly mode). This approach can be problematic for two main reasons. First, its sensitivity is low because most communities are structured by both niche and dispersal processes [START_REF] Mouillot | Niche-assembly vs dispersal-assembly rules in coastal fish metacommunities: implications for management of biodiversity in brackish lagoons[END_REF]. Second, being qualitative, this approach prevents precise comparisons (in space and time) of the assembly processes in different metacommunities. To address these weaknesses, we developed the DNCI based on the three E-metric values:

DNCI SES SES = - = - ae è ç ç ö ø ÷ ÷ - - = ( ) = ( ) å å d n i n d i dn dn i n n i n E E E n E E 1 1 1 1 s d dn dn E s ae è ç ç ö ø ÷ ÷
where SES d and SES n are the standard effect sizes of E d and E n , respectively, E d , E n and E dn are the E-metric values (i.e. deviation from the empirical SIMPER profile) for the three PER-SIMPER null models (i.e. 'dispersal assembly' permutation mode, 'niche assembly' permutation mode and 'dispersal and niche assembly' permutation mode, respectively), and the standalone n is the number of iterations selected for the PER-SIMPER method.

The standard deviation related to DNCI is obtained by simple error propagation from the standard deviations associated with SES d and SES n , such as:

s s s DNCI SES S ES = + d n 2 2
If the DNCI is not significantly different from 0, the dispersal and niche processes could be assumed to contribute equally to variations in community composition. If the DNCI is significantly lower than 0, dispersal processes are the dominant drivers of community composition. Conversely, if the DNCI is significantly higher than 0, niche processes are the primary determinants of community composition (Fig. 1).

It is noteworthy that negative DNCI values indicating the predominance of dispersal processes do not give information on actual dispersal rates.

DNCI robustness test with simulated communities

PER-SIMPER's original qualitative approach has shown accurate identification of the main assembly process underlying the composition of simulated communities, even when the analysed groups have different number of taxa or result from various sampling intensities. Furthermore, PER-SIMPER's qualitative identifications are particularly robust when one assembly process plays a role significantly larger than the other one. Conversely, when both processes are simulated, PER-SIMPER is more sensitive to sampling biases. In particular, the number of sampled sites (less than five sites in a site group) within the analysed dataset may bias the PER-SIMPER inference toward the niche-based expectation by underestimating the effect of dispersal processes [START_REF] Gibert | PER-SIMPER -a new tool for inferring community assembly processes from taxon occurrences[END_REF].

Here, the effects of number of taxa and number of sites are tested and measured on simulated datasets prior to using the new DNCI on empirical datasets. To ensure consistency between the quantitative and qualitative robustness tests, the simulated communities are obtained from cellular automata detailed by [START_REF] Gibert | PER-SIMPER -a new tool for inferring community assembly processes from taxon occurrences[END_REF]. The automata are designed for simulating the different displacement potential of various taxa (dispersal) as well as the influence of local environmental conditions on taxonomic richness (niche). The dispersal potential of taxa and the taxonomic richness limit inside cells can be tuned by, respectively, changing the number of individuals moving along the automaton grid (i.e. the more individuals, the more dispersal) and by changing the pattern limiting maximum richness within the cells. Using these two parameters, it is possible to simulate a complete continuum of communities from those strongly constrained by niche processes to those strongly constrained by dispersal as well as a set of intermediate types of communities.

The final occurrence matrices are composed of 30 sites per group (i.e. 60 sites in total) and 100-150 taxa. Three simulated matrices, one assembled only by niches, one assembled only by dispersal and one assembled jointly by niches and dispersal, are used as references. Their measured DNCI values are then compared to the DNCI values computed for site groups with increasing proportions of removed taxa and sites. Specifically, 10-70% of the sites are randomly removed, and 10-70% of the simulated taxa are randomly removed.

The removals were done in one site group only or in both site groups simultaneously.

Observational field survey data on diatoms, macroinvertebrates and fish from the United States

Datasets including stream diatoms, macroinvertebrates and fish were obtained from the National Water-Quality Assessment (NAWQA) Program of the U.S. Geological Survey. We selected four subregions of the Mid-Atlantic Region (hydrologic unit HUC02), situated in the northeastern part of the United States, including the Upper Hudson (HUC0202), Delaware (HUC0204), Susquehanna (HUC0205) and the Potomac (HUC0207) (map in Supporting information). In total, there were 131 sites sampled for diatoms, 128 sites sampled for macroinvertebrates and 43 sites sampled for fish across the four subregions (i.e. four site groups). Each of the subregions comprised a sufficient number of sites for the analyses (a minimum of eight sites in Delaware for fish and a maximum of 66 sites in the Potomac for macroinvertebrates), given that at least five sites per group is recommended by [START_REF] Gibert | PER-SIMPER -a new tool for inferring community assembly processes from taxon occurrences[END_REF]. Importantly, the subregions likely shared the same regional species pools as they belong to the same hydrologic region. The streams were sampled between 1993 and 2011. We included only summer months, i.e. June, July and August, to eliminate seasonal effects on community compositions and assembly processes, which can make interpretations more complex [START_REF] Li | Seasonal changes in metacommunity assembly mechanisms of benthic macroinvertebrates in a subtropical river basin[END_REF]). We then selected taxa identified at the species level for diatoms (444 species) and fish (93 species), and at species or genus levels for macroinvertebrates (429 taxa).

Observational field survey data on bacteria and macroinvertebrates from China

In October 2014, 89 individual streams were sampled for bacteria and macroinvertebrates in the southeastern part of the Tibetan Plateau, Yunnan province, China. The sampled streams include 52 tributaries of the Salween River and 37 tributaries of the Mekong River (see map in Supporting information). The elevational ranges were pronounced, spanning from 570 to 1592 m a.s.l. in the Salween catchment and from 1270 to 2100 m a.s.l. in the Mekong catchment. Macroinvertebrate samples were identified to genus level when possible. Bacterial OTUs (i.e. operational taxonomic units) were determined based on the 16S rRNA genes using bacterial universal primers 515F and 806R that target the V4 region. The bacterial dataset was rarefied to 10 000 sequences. Further details on field sampling and laboratory procedures are presented in [START_REF] Wang | Regional and global elevational patterns of microbial species richness and evenness[END_REF] and [START_REF] Vilmi | Ecological processes underlying the assembly of aquatic bacterial and macroinvertebrate communities under contrasting climates on the Tibetan Plateau[END_REF].

Statistical methods

We developed an R function ('DNCI_multigroup') for the calculation of DNCI values. The function first performs the PER-SIMPER procedure [START_REF] Gibert | PER-SIMPER -a new tool for inferring community assembly processes from taxon occurrences[END_REF] and, then, as described earlier, calculates pairwise DNCI values among site groups based on the E-values from the PER-SIMPER analysis. For the US data, we used the four subregions as site groups for all three presence/absence community matrices (i.e. diatoms, macroinvertebrates and fish). Due to differing numbers of sites per group, we selected the size of the smallest group (Delaware with eight sites in the fish dataset) and randomly resampled sites from all other groups to eliminate a potential group size bias. We performed the random resampling of eight sites from each of the four groups 200 times. For the Chinese data, we used elevational bins, each containing 12 or 13 sites, as site groups in our analyses (Supporting information). In both datasets, we removed singletons (i.e. taxa occurring at one site only) prior to calculation of the DNCI values.

We employed the Kruskal-Wallis test followed by posthoc Dunn tests with p-value adjustments (method = 'bh') to study whether the DNCI differed across diatoms, macroinvertebrates and fish in the US Mid-Atlantic hydrologic subregions, and the Wilcoxon rank sum test to study whether the DNCI of bacteria and macroinvertebrates differed from each other in the Chinese catchments. Separately for bacteria and macroinvertebrates, the Kruskal-Wallis test followed by posthoc Dunn tests with p-value adjustments (method = 'bh') was used to study whether within-and across-catchment DNCIs differed from each other in the Chinese dataset. With Mantel tests, we studied whether the DNCI of Chinese bacteria and macroinvertebrates varied with distance using pairwise DNCI values and pairwise spatial distances of the site groups. As the computed DNCI values were all negative, we performed Mantel tests with absolute DNCI values (Mantel test does not process negative values). The absolute DNCI values simply indicate the distances to zero. After analyses with absolute values, the original DNCI values (i.e. negative, in our case) were plotted against the pairwise spatial Euclidean distances of the site groups. As the relationships were apparently negative, the Mantel r statistic was changed to negative. It is important to note that Mantel tests can only be used for cases where all DNCI values are either positive or negative, the latter after the transformation described above. We made network plots based on the results from China to graphically illustrate the variation of DNCI in space.

All statistical analyses and plots were done using packages DNCImper [START_REF] Gibert | DNCImper: assembly process identification based on SIMPER analysis[END_REF], vegan [START_REF] Oksanen | Package 'vegan[END_REF], FSA [START_REF] Ogle | FSA: fisheries stock analysis[END_REF]), ggplot2 [START_REF] Wickham | ggplot2: elegant graphics for data analysis[END_REF], igraph [START_REF] Csardi | The igraph software package for complex network research[END_REF] and dichromat [START_REF] Lumley | dichromat: color schemes for dichromats[END_REF] with the R statistical software (<www.r-project.org>). The functions and codes for DNCI calculation and examples for preparing plots are available in GitHub (<https://github. com/Corentin-Gibert-Paleontology/DNCImper>) and as Supporting information.

Results

Robustness tests

If only niche effects are simulated, removal of 40% of the taxa or 30% of the sites in one site group of the dataset can lead the DNCI analysis to fail in identifying the correct assembly mechanism (Fig. 2a, 3a, blue boxplots). For assemblages simulated solely via dispersal effects, the removal of 50% of the sites in one site group can also lead to the identification of the wrong assembly mechanisms (Fig. 3a, orange boxplot). Consequently, the DNCI is more robust to variation in number of taxa than in number of sites but can also be highly sensitive to specific variation in number of sites (Fig. 2,3). These results are consistent with preliminary robustness tests conducted on the qualitative PER-SIMPER method [START_REF] Gibert | PER-SIMPER -a new tool for inferring community assembly processes from taxon occurrences[END_REF]. However, the computation of robustness tests on a quantitative index provides a more precise understanding of the consequences of sampling effort variation on the method's ability to detect assembly mechanisms and the intensity of their respective contributions. First, the computation of a single, robust DNCI value requires groups with a limited difference in number of taxa and sites (Fig. 2a,3b). Specifically, comparing groups with deviation in number of taxa exceeding 40% is not recommended as it leads to underestimation of the niche processes (Fig. 2a). In addition, to discriminate the respective contribution of niche and dispersal processes, the DNCI must be computed for site groups with less than 30% differences in site numbers, which is especially important for communities assembled by both niche and dispersal processes (Fig. 3a). The DNCI is, therefore, to some extent sensitive to the symmetry of the occurrence matrix of interest. Second, DNCI values can be compared accurately across datasets of different size both in terms of taxa and site numbers, as DNCI values are insensitive to symmetric removal of sites or taxa (at least up to 70% differences; Fig. 2b,3b). Consequently, random re-sampling [START_REF] Sanders | Marine benthic diversity: a comparative study[END_REF][START_REF] Raup | Taxonomic diversity during the Phanerozoic[END_REF][START_REF] Gotelli | Quantifying biodiversity: procedures and pitfalls in the measurement and comparison of species richness[END_REF]) may be necessary when computing DNCI between site groups of dissimilar size within a dataset, but not when comparing DNCI values between distinct datasets of dissimilar size. This is a significant improvement in terms of precision compared to PER-SIMPER, which confounds the effect of matrix symmetry with the effect of unequal sample size.

Specific applications of DNCI

Identifying the main community assembly process in major stream organismal groups

The DNCI values of diatoms, macroinvertebrates and fish from the Mid-Atlantic hydrological subregions of the US were negative, suggesting that dispersal was the dominant assembly process for all three organismal groups (Fig. 4). The DNCI values of diatoms significantly differed from the DNCI values of macroinvertebrates and fishes (p < 0.001). Dispersal assembly was, in a relative sense, most intense for diatoms and less intense for macroinvertebrates and fishes (Fig. 4).

Identifying the main community assembly process for stream bacteria and macroinvertebrates along an elevational gradient

The DNCI values of bacteria and macroinvertebrates were negative, suggesting dominance of dispersal assembly in the two Chinese catchments. The DNCI differed significantly between the two organismal groups (p = 0.016), with bacteria being more strongly assembled by dispersal than macroinvertebrates. The DNCI of bacteria showed a decreasing trend with increasing spatial distances (Mantel r = -0.486, p = 0.043; Fig. 5a). For macroinvertebrates, a relationship between the DNCI and distance was not evident (Mantel r = -0.153, p = 0.251; Fig. 5b). The DNCI of bacteria was constant within and across the studied catchments (Fig. 5c; p = 0.326). The difference in DNCI of macroinvertebrates within the Mekong and Salween catchments was close to the p = 0.05 level of statistical significance (overall Kruskal-Wallis test p = 0.056, post-hoc Dunn test for comparisons of within-catchment DNCIs p = 0.050; Fig. 5d), with DNCI values being more negative within the Salween than the Mekong catchment. Network plots (Fig. 5e-f ) illustrated the spatial setting and relationships between site groups (here, elevational bins), and complement the information presented in Fig. 5a-d. For bacteria, dispersal assembly appeared to be strongest between the lower mid elevations (S-MID2) and high elevations (S-HIGH) within the Salween catchment, as well as between the lower mid elevations of the Salween catchment (S-MID2) and high elevations of the Mekong catchment (M-HIGH) (Fig. 5e). For macroinvertebrates, the strongest dispersal assembly signal was observed across the catchments, i.e. between the low elevations in the Salween catchment and mid elevations of the Mekong catchment (Fig. 5f ). Interestingly, for both bacteria and macroinvertebrates, dispersal assembly was weaker in the highest elevations across the two catchments (S-HIGH and M-HIGH). For macroinvertebrates, community dissimilarities were more weakly driven by dispersal assembly across the highest elevations in the Mekong catchment and all other elevations within and across catchments (M-HIGH versus all other elevational bins).

Discussion

Here, we developed a novel quantitative index, DNCI, to disentangle the importance of niche and dispersal processes in community assembly, with distinct advantages over existing methods, including earlier null model approaches [START_REF] Chase | Using null models to disentangle variation in community dissimilarity from variation in a-diversity[END_REF], variance partitioning in constrained ordination [START_REF] Borcard | Partialling out the spatial component of ecological variation[END_REF]) and PER-SIMPER's original qualitative inference procedure [START_REF] Gibert | PER-SIMPER -a new tool for inferring community assembly processes from taxon occurrences[END_REF]. The advantages of the DNCI include, 1) reliance on taxa distributions rather than summary statistics, with the potential to better assess assembly processes (cf. [START_REF] Chase | Using null models to disentangle variation in community dissimilarity from variation in a-diversity[END_REF], 2) relative simplicity, given that it does not require abundance, environmental or spatial data and, 3) high comparability across datasets, even those that differ greatly in numbers of sampled localities and taxa. Thus, the DNCI allows evaluating whether the dominant assembly mechanism and its strength vary across different datasets.

Based on observational data from the US and China, the DNCI-based results suggested that communities of distinct organismal groups, from microbes to fishes, were mainly assembled by dispersal processes. Differences in the DNCI were partly associated with organismal groups and geographical contexts. In the empirical datasets from the US and China, DNCI values varied among organismal groups, with dispersal assembly being slightly but significantly stronger for the metacommunities of smaller organisms. Furthermore, in the Chinese data, our results implied that the intensity of dispersal assembly may increase at greater distances for microbes. Although the DNCI does not inform us about actual dispersal rates, it can perhaps be reasoned that the distance-related DNCI pattern observed for the Chinese microbes may suggest that the occurrences are, to some degree, limited by dispersal in the two catchments. Overall, dispersal was the dominant assembly process in the studied regional datasets, but further research is needed to assess whether the intensity of dispersal assembly changes across larger spatial distances and environmental gradients in other datasets.

For the Chinese data, some pairs of site groups were very strongly structured by dispersal both within and across catchments, while some pairs of site groups were structured by a more balanced mix of dispersal and niche processes within and across catchments (Fig. 5e-f ). These two examples illustrate some of the limitations of the DNCI approach related to the use of occurrence matrices and pairwise comparisons. First, the use of a single occurrence matrix without additional environmental information makes the use of DNCI extremely straightforward, but it also implies a lack of information when it comes to the effect of niche mechanisms on the distribution of taxa. Only the assumed effect of niche mechanisms on the presence/absence of taxa and, therefore, on the richness of localities (i.e. row sums) are compared with the assumed effect of dispersal potential of taxa (i.e. column sums). This bias might be important especially for small organisms, such as bacteria, where large differences in abundance between sites rather than richness may be observed. Ecological differences between and within the two catchments were not strong enough compared to differences in dispersal potentials to explain the dissimilarity structure and thus the composition (in terms of occurrence) of assemblages. This was evidenced by the highly negative DNCI values obtained from the datasets collected from the two catchments. For macroinvertebrates, the three site groups within the Mekong catchment showed a slightly more balanced situation between dispersal and niche processes, which was indicated by higher DNCI values and, therefore, a comparatively increased role of niche conditions on the presence and absence of the sampled taxa. Second, the DNCI approach analyses pairs of assemblages. Therefore, a single analysis does not allow us to conclude on the nature of the structure of one community but on the nature of the assembly processes that differentiate groups of local communities. Here, network plots can be useful for interpreting results. By simultaneously observing the DNCI values of multiple pairs of site groups, we can identify the pairs intensely structured by dispersal (such as S-MID2-M-HIGH for bacteria and S-LOW-M-MID for macroinvertebrates) as well as those comparatively more marked by the effect of niche conditions (such as M-HIGH paired with all other site groups for macroinvertebrates). Our general results that communities from different spatial site groups were more differentiated by dispersal than niche processes can be partly seen as a consequence of the pairwise nature of the method and, in our example datasets, of a comparatively weak environmental effect, assumingly portrayed by the local taxonomic richness, on the presence or absence of taxa.

Although the DNCI values based on our observational example datasets were all negative, indicating that dispersal assembly was dominant for all tested aquatic metacommunities from two continents, it is likely that datasets from different biological, environmental, spatial and temporal settings may produce both positive and negative DNCI values. Exploring the sensitivity of DNCI can guide future research to assess how community assembly is determined by the following factors, for instance: 1) spatial extent of studied areas (i.e. site groups in the analysis); 2) spatial distances among studied areas (site groups) and within studied areas (sites); 3) physical connectivity within and across site groups; 4) organismal groups, taxonomic levels and functional traits (e.g. dispersal ability, stress tolerance or feeding behaviour), 5) environmental heterogeneity versus homogeneity; 6) successional phase; 7) seasonality; 8) temporal environmental trends, such as progressive eutrophication, land use changes and climate change. The more the proposed DNCI approach is applied, the more information we obtain on how accurate the pairwise calculations based on presence-absence matrices of site groups are for inferring the intensities of dispersal and niche processes in community assembly in various study settings.

There are a few caveats for applying the DNCI approach. First, an appropriate spatial scale must be defined, whereby all local communities share at least partly the same regional species pool (i.e. dispersal of organisms potentially occurs within the study context). Second, based on our simulations, the numbers of taxa and sites per group should not vary more than 40% and 30%, respectively. Third, since this approach uses occurrence data, the assembly processes identified by this approach are limited to those resulting from the presence or absence of taxa, and the driving mechanisms of abundance variation are left undetected, suggesting a further need for methodological development of the DNCI. Fourth, distances between sites and site groups and, more generally, the overall spatial setting, are not directly accounted for (but indirectly via the occurrence matrix structure), potentially complicating interpretation of results. Fifth, unlike methods which use numerous environmental variables as predictors (e.g. variation partitioning), the DNCI alone cannot show which types of environmental variables contribute to variation in community structure across sites and site groups.

Keeping in mind the requirements for shared regional species pools and limited size variations among site groups as well as the caveats described above, the proposed DNCI approach can shed light into the importance of niche and dispersal processes by measuring the strength of community assembly processes in a way that is simple, quantifiable and easily comparable across datasets. We believe that this method is a useful addition to the toolbox of biogeographers, ecologists and palaeontologists; if used for datasets from different environmental and spatial settings, it would increase understanding of community assembly in nature.

Figure 1 .

 1 Figure 1. Flowchart of PER-SIMPER derivation from the initial SIMPER method and the development and interpretation of the DNCI. Negative DNCI values indicate that dispersal assembly is the dominating process, while positive DNCI values indicate that niche assembly is the dominating process behind community dissimilarities among site groups. With increasing absolute values, the strength of the dominating process increases. DNCI values close to zero indicate that both niche and dispersal assembly processes affect distributions. In dataset A, the DNCI does not significantly differ from zero, indicating that both dispersal and niche processes control the empirical species distributions. In dataset B, DNCI is significantly greater than zero, indicating that niche assembly is the main driver of the empirical species distribution. In datasets C and D and the illustrated dataset, DNCI is significantly less than zero, indicating that dispersal assembly is the main driver of the empirical species distributions. Dataset D and the illustrated dataset do not have significantly different DNCI values, i.e. the control of the empirical species distributions by dispersal assembly is similar. Dataset C has a DNCI significantly higher than dataset D and the illustrated dataset, i.e. dispersal assembly control in dataset C is significantly weaker than in dataset D and the illustrated dataset. SES d = standard effect size of E d . SES n = standard effect size of E n . Modified from Gibert and Escarguel's (2019) Supporting information.

Figure 2 .

 2 Figure 2. Distribution of DNCI (1000 iterations) after degradation of three reference datasets: niche assembled (N, blue line), niche and dispersal assembled (N + D, grey line) and dispersal assembled (D, orange line). The number of taxa is randomly decreased by 10-70% within a single group of sites of the reference dataset (a). The number of taxa is randomly decreased by 10-70% in each of the two groups of sites in the reference dataset (b).

Figure 3 .

 3 Figure 3. Distribution of DNCI (1000 iterations) after degradation of three reference datasets: niche assembled (N, blue line), niche and dispersal assembled (N + D, grey line) and dispersal assembled (D, orange line). The number of sites is randomly decreased by 10-70% within a single group of sites in the reference dataset (a). The number of sites is randomly decreased by 10-70% in each of the two groups of sites in the reference dataset (b).

Figure 4 .

 4 Figure 4. DNCI values of diatoms (Dia), macroinvertebrates (Inv) and fish (Fish) from four Mid-Atlantic hydrological subregions in the US. Each dot represents the DNCI value between two groups computed in one iteration, i.e. the total number of DNCI values computed for each organismal group is 1200 (six pairwise comparisons across site groups, 200 iterations). Negative DNCI values indicate that dispersal is the dominant assembly process. Higher absolute values (darker colour) indicate greater strength of the dominant assembly process. Different letters represent statistically different DNCI values following the Kruskal-Wallis and post-hoc Dunn tests (p < 0.001). The black dots illustrate the mean values and the lines represent the standard deviation.

Figure 5 .

 5 Figure5. The DNCI values for bacteria (left column) and macroinvertebrates (right column) in the Chinese dataset. There was a significant negative relationship between DNCI of bacterial communities and distance (a), suggesting that intensity of dispersal assembly grew with increasing distances. There was no clear relationship between the DNCI of macroinvertebrate communities and distance (b). Based on the Kruskal-Wallis test, the DNCI of bacteria within and across catchments did not differ from each other (c; p = 0.326). Based on the Kruskal-Wallis and post-hoc Dunn tests, the DNCI of macroinvertebrates was almost statistically significantly different within the Mekong and the Salween catchments (d; overall Kruskal-Wallis p = 0.056, post-hoc Dunn p = 0.050). Network plots illustrate the DNCI, i.e. the strength of the dominating assembly process, within and across catchments for bacteria (e) and macroinvertebrates (f ). The nodes in the network plots illustrate site groups, each containing 12 or 13 sites, along elevational gradients in the Salween (S) and Mekong (M) catchments. The edges connecting the nodes in the network plots illustrate the intensity of the dominant assembly process; the darker the colour and wider the line, the more intense the process (i.e. the bigger the absolute value of the DNCI).

Acknowledgements -We acknowledge M. Luo and C. Qin for their help in sampling the field data from Yunnan and C. Yeh for help with processing the data. We thank Emeritus Professor Jari Oksanen for clarifying discussions on SIMPER and PER-SIMPER during the early phases of this work.

Funding -This study was supported by National Natural Science Foundation of China (91851117), the Second Tibetan Plateau Scientific Expedition and Research (STEP) Program (2019QZKK0503), the Chinese Academy of Sciences Key Research Program of Frontier Sciences (QYZDB-SSW-DQC043), National Natural Science Foundation of China (41871048, 41571058), the Chinese Academy of Sciences President's International Fellowship Initiative (2018PS0007), and the project FRESHABIT LIFE IP (LIFE14/IPE/FI/023). The field trip to Yunnan was partly supported by Air and Water Conservation Fund (GEFC12-14, National Geographic).

Author contributions -AV and CG contributed equally to this paper. AV and JW conceived the idea. AV, CG, GE and JW designed methodology; AJ, AV, CG, KH and JW generated the functions; CG ran the simulations; SP provided the US data; JW provided the Chinese data; AV analysed the empirical data; AV and CG led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

Data availability statement

The R functions, package and codes are available in GitHub (<https://github.com/Corentin-Gibert-Paleontology/ DNCImper>).


El Kamal

Eliane Omari

Teodor Younes

Cathy Burghelea

Yann Castelain

Yves Moguen

Le Guer

Omari

Eliane Younes
email: eliane.younes@univ-nantes.fr

Teodor Burghelea
email: teodor.burghelea@univ-nantes.fr

Cathy Castelain
email: cathy.castelain@univ-nantes.fr

Yann Moguen
email: yann.moguen@univ-pau.fr

Yves Le Guer
email: yves.leguer@univ-pau.fr

Active chaotic mixing in a channel with rotating arc-walls

I. INTRODUCTION

The efficient transfer of mass or heat in open channel flows of viscous fluids is of importance in many industrial processes: production of food pastes, waste water treatment, pulp and paper manufacturing, inline production of polymer blends with controlled poly-dispersity, etc. In open channel flow mixers, the active zone is continuously fed with the fluids to mix, which spend only a definite time in this zone, in contrast with closed flows mixers. Efficiency is then a crucial issue. To achieve efficient mixing, two categories of open flow mixers can be distinguished, depending on nature of the stirring protocol [START_REF] Cai | A review on micromixers[END_REF][START_REF] Hessel | Micromixers: a review on passive and active mixing principles[END_REF][START_REF] Stone | Engineering flows in small devices: microfluidics toward a Lab-on-a-Chip[END_REF]: passive or active.

Passive elements use the energy of the flow to generate flow structures enhancing mixing efficiency. This may be carried out by successive changes of direction of the channel walls, for example using twisted pipes [START_REF] Boesinger | Experimental study of reactive chaotic flows in tubular reactors[END_REF][START_REF] Castelain | Experimental study of chaotic advection regime in a twisted duct flow[END_REF][START_REF] Jones | Chaotic advection by laminar flow in a twisted pipe[END_REF] or grooved patterns [START_REF] Stroock | Chaotic mixer for microchannels[END_REF]. Passive elements can also move in the flow, such as flexible structures like flaps [START_REF] Ali | Heat transfer and mixing enhancement by free elastic flaps oscillation[END_REF][START_REF] Lambert | The role of elastic flap deformation on fluid mixing in a microchannel[END_REF]. Pursuing this idea further, viscoelastic polymers can be added in the fluids to generate a regime of elastic turbulence favorable for mixing [START_REF] Burghelea | Chaotic flow and efficient mixing in a microchannel with a polymer solution[END_REF][START_REF] Burghelea | Mixing by polymers: experimental test and decay regime of mixing[END_REF][START_REF] Groisman | Efficient mixing at low Reynolds numbers using polymer additives[END_REF]. More examples of passive elements are listed in Table I, left column. For a given geometry of an open flow mixer with only passive stirring elements, the bulk velocity is the sole parameter that controls the mixing efficiency. However, for highly viscous fluids, the bulk velocity may be difficult to tune and, consequently, the mixing efficiency is rather limited.

To efficiently mix highly viscous fluids, active elements are expected to perform better. Such elements require an energy source other than the energy associated to the mean flow. Various types of forcing schemes may be used, such as hydrodynamic forcing [START_REF] Bottausci | Mixing in the shear superposition micromixer: three-dimensional analysis[END_REF][START_REF] Okkels | Spatiotemporal resonances in mixing of open viscous fluids[END_REF] or structural-based forcing (egg-beater [START_REF] Gouillart | Open-flow mixing: experimental evidence for strange eigenmodes[END_REF], rotating walls [START_REF] Ober | Active mixing of complex fluids at the microscale[END_REF], artificial cilia [START_REF] Den Toonder | Artificial cilia for active micro-fluidic mixing[END_REF][START_REF] Khatavkar | Active micromixer based on artificial cilia[END_REF], etc.). More examples are listed in Table I, right column.

Passive stirring

Active stirring Converging-diverging channel [START_REF] Amon | Lagrangian chaos, Eulerian chaos, and mixing enhancement in converging-diverging channel flows[END_REF][START_REF] Gepner | Flow dynamics and enhanced mixing in a converging-diverging channel[END_REF] Hydrodynamics [START_REF] Bottausci | Mixing in the shear superposition micromixer: three-dimensional analysis[END_REF][START_REF] Okkels | Spatiotemporal resonances in mixing of open viscous fluids[END_REF] Rigid or flexible structures [START_REF] Ali | Heat transfer and mixing enhancement by free elastic flaps oscillation[END_REF][START_REF] Lambert | The role of elastic flap deformation on fluid mixing in a microchannel[END_REF] Rigid or flexible structures [START_REF] Den Toonder | Artificial cilia for active micro-fluidic mixing[END_REF][START_REF] Gouillart | Open-flow mixing: experimental evidence for strange eigenmodes[END_REF][START_REF] Khatavkar | Active micromixer based on artificial cilia[END_REF][START_REF] Ober | Active mixing of complex fluids at the microscale[END_REF]] Staggered herringbone channel [START_REF] Stroock | Chaotic mixer for microchannels[END_REF] Acoustics [START_REF] Collins | Continuous micro-vortex-based nanoparticle manipulation via focused surface acoustic waves[END_REF][START_REF] Huang | An acoustofluidic micromixer based on oscillating sharp-edges[END_REF] Partitioned pipe mixer [START_REF] Khakhar | A case study of chaotic mixing in deterministic flows: the partitioned-pipe mixer[END_REF] Electro-hydrodynamics [START_REF] Chang | Electrokinetic mixing in microfluidic systems[END_REF] Twisted pipes [START_REF] Boesinger | Experimental study of reactive chaotic flows in tubular reactors[END_REF][START_REF] Castelain | Experimental study of chaotic advection regime in a twisted duct flow[END_REF][START_REF] Jones | Chaotic advection by laminar flow in a twisted pipe[END_REF] Magneto-hydrodynamics [START_REF] Gopalakrishnan | Electromagnetically induced chaotic mixing in a pipe mixer[END_REF] Focusing channels [START_REF] Knight | Hydrodynamic focusing on a silicon chip: mixing nanoliters in microseconds[END_REF] Acoustically driven bubbles [START_REF] Rallabandi | Analysis of optimal mixing in open-flow mixers with time-modulated vortex arrays[END_REF] Elastic turbulence [START_REF] Burghelea | Chaotic flow and efficient mixing in a microchannel with a polymer solution[END_REF][START_REF] Burghelea | Mixing by polymers: experimental test and decay regime of mixing[END_REF][START_REF] Groisman | Efficient mixing at low Reynolds numbers using polymer additives[END_REF] Dielectrophoresis [START_REF] Lee | Optimizing micromixer design for enhancing dielectrophoretic microconcentrator performance[END_REF] Corrugated channel [START_REF] Ghanem | Mixing performances of swirl flow and corrugated channel reactors[END_REF] Electrokinetics [START_REF] Oddy | Electrokinetic instability micromixing[END_REF] TABLE I: Examples of configurations for passive stirring, and types of forcing for active stirring, for open flow mixing (not exhaustive list).

For many active open flow mixers, the mixing is carried out along a channel by moving impeller-like shapes. Yet the mixing length, defined as the characteristic length the fluids need to travel downstream to mix, may be actually smaller than the length of the entire processing line. This can lead to additional energetic costs which would ultimately increase the price of the final product. A second drawback relates to the quality of the final product, like its texture. Many highly viscous pasty materials are made of high molar mass molecular constituents which are mechanically fragile. "Over-mixing" such products may result in mechanical degradation or in a loss of the targeted texture. To circumvent these problems, which have hitherto never been quantitatively addressed [START_REF] Hessel | Micromixers: a review on passive and active mixing principles[END_REF], certain manufacturers are little inclined to use inline mixers. It is clear that for the latter the ability to easily control the level of mixing is desirable.

Following a preliminary work by two of the authors [START_REF] Omari | Chaotic mixing in a plane channel with rotating arc walls[END_REF], the central aim of the present paper is to study an active open flow mixer that allows both good mixing efficiency and easy control of the mixing level for highly viscous fluids.

By the large extent of the mixing devices considered above, one can infer that, even in the case of a Newtonian fluid studied in the present paper, the design of the mixer and the choice of the mixing protocol remain challenging tasks. One reason for that is the lack of understanding of the mixing mechanisms, which involve a broad range of spatial and temporal scales and complex interactions between processes at different scales. However, it is generally accepted that in a highly viscous fluid, the mixing of a passive scalar, hereafter referred to as concentration, requires two ordered steps [START_REF] Eckart | An analysis of the stirring and mixing processes in incompressible fluids[END_REF][START_REF] Villermaux | Mixing versus stirring[END_REF]. The first step is the stretching of the concentration blobs into striations (also called filaments, or lamellae in 2-D) with possibly enhanced diffusion by large shear [START_REF] Foures | Optimal mixing in two-dimensional plane Poiseuille flow at finite Péclet number[END_REF] (also called Taylor dispersion). The second step consists of the folding of the passive scalar filaments allowing for an efficient pure diffusion due to the generated high concentration gradients. One of the simplest mechanisms that permit to achieve the first step is the roll-up of concentration blobs in spirals [START_REF] Meunier | How vortices mix[END_REF][START_REF] Rhines | How rapidly is a passive scalar mixed within closed streamlines?[END_REF]. Thanks to its rotating arc-walls, the mixer we propose can generate such suitable spirals, provided that the transverse velocity is sufficiently large with respect to the bulk velocity. (What is meant by "sufficiently large" is studied in the present paper.) The step of folding is performed by reorientation of the page 2 transverse velocity field, namely by the rotational changes of direction of the arc-walls. For simplicity, the rotations of the arc-walls are chosen time periodic. Obviously, an infinite number of protocols could be used to achieve mixing by rotation of the arc-walls. There are practically no theoretical tools available for allowing an a priori choice of a particular form of protocol. Based on qualitative considerations, we thus have chosen a form of protocol that allows us to generate a flow in which hyperbolic and elliptic points coexist, which is recognised as a favorable condition for efficient mixing [START_REF] Rodrigo | On the optimization of mixing protocol in a certain class of three-dimensional Stokes flows[END_REF]. We will show that efficient mixing is achieved by a careful choice of the period of rotation and the transverse velocity relative to the bulk velocity. We do not consider here the issue of how to control the mixing level, as done, e.g., by feedback control by Aamo et al. [START_REF] Aamo | Control of mixing by boundary feedback in 2D channel flow[END_REF], by optimal control by Foures et al. [START_REF] Foures | Optimal mixing in two-dimensional plane Poiseuille flow at finite Péclet number[END_REF], or by entropy maximization by D'Alessandro et al. [START_REF] D'alessandro | Control of mixing in fluid flow: a maximum entropy approach[END_REF]. The control of mixing often requires a large computational effort. On the contrary, our approach here aims at obtaining choices of control parameters based on extremely simple arguments. For that purpose, a heuristic model based on dimensional and phenomenological arguments will be proposed, giving insights into suitable choices of the control parameters for good mixing.

In the present study a 2-D flow model is considered, as a first step towards understanding of the mixing mechanisms in a realistic 3-D channel mixer. It is recognized that for channel flows, 3-D effects may increase mixing beyond the 2-D mixing mechanisms. Moreover, in 3-D flows, 2-D and 3-D transport mechanisms may be related to each other [START_REF] Pouransari | Formation of coherent structures by fluid inertia in three-dimensional laminar flows[END_REF][START_REF] Ravu | Global organization of three-dimensional, volumepreserving flows: Constraints, degenerate points, and Lagrangian structure[END_REF]. This motivates the present 2-D study, which may thus be thought of as giving insights into mixing mechanisms in the channel, despite being only a first step towards understanding of the 3-D mechanisms.

The paper is organised as follows. In section II, the geometry of the proposed mixer is presented, along with the selected flow modulation protocol. In section III, heuristic conditions for mixing, based on dimensional and phenomenological arguments, are proposed. In section IV, the numerical setting-up for two-dimensional simulations is described. Results by the heuristic and numerical approaches are compared in section V. The presence of chaotic advection in the cases of good mixing is evidenced. The paper closes with a summary of the main findings and concluding remarks in section VI.

II. THE ROTATING ARC-WALL MIXER

With the objective of mixing efficiently with control of the mixing level, we propose a mixing unit, named Rotating Arc-Wall (or RAW) mixer, that combines the axial flow through a channel and the cross-flow generated by three rotating circular arc-walls.

A. Geometry and kinematics of the RAW mixer

For the two-dimensional configuration considered for simplicity in the present paper, the RAW mixer consists of a channel of length L and width W , with three circular cylinders of equal radii R having their axes of rotation perpendicular to the plane of the channel (see Fig. 1). The cylinders penetrate into the channel with the depth ∆. Two cylinders are located on the same channel wall, at a distance of 2R, and one other is located on the opposite wall, in the middle of the two facing cylinders. The three cylinders can rotate independently, with potentially different directions of rotation, with or without temporal modulations of their angular velocities. This flexibility allows a rich variety of flow patterns in the mixer.

B. A protocol for mixing

As a first step towards the choice of a suitable protocol for the angular velocities of the cylinders, we identify the steady flow topologies that may appear in the RAW mixer without modulation of the angular velocities. Combining the possibilities of rotation of the three cylinders in clockwise and counter-clockwise directions, we obtain 2 3 typical steady-state solutions whose streamlines are shown in Fig. 2. One, two or three elliptic regions [START_REF] Aref | Frontiers of chaotic advection[END_REF][START_REF] Haller | Lagrangian coherent structures[END_REF] appear (see the light-blue zones in Fig. 2), surrounded by separatrices that are transport barriers for the tracers from one elliptic region to another [START_REF] Biemond | Onset of chaotic advection in open flows[END_REF]. To break them, the flow is made time-dependent through the rotation of the three cylinders. The modulations of the angular velocities are designed to allow for the reorientation of the transverse velocity field [START_REF] Horner | Transport enhancement mechanisms in open cavities[END_REF]. This is simply achieved by an alternation of the rotation directions. Four couples of steady solutions can then be selected (see Fig. 2) : {(1), (2)}, {(3), (4)}, {(5), (6)} and {(7), (8)}, for which the streamlines shown in the figure evolve from one solution to another. At each change of direction, attachment points of the separatrices are either created, removed or destabilized, with a beneficial effect for material transport and thus for mixing. The case of a counter-rotating protocol, in which the cylinders located on the opposite walls rotate in opposite directions, allows the generation of hyperbolic points beneficial to the mixing efficiency. It is also desirable that these points are not page 3 FIG. 1: Sketch of the two-dimensional rotating arc-wall (RAW) mixer (not to scale). v bulk represents a given bulk velocity of the entering flow.

confined close to the fixed walls [START_REF] Gouillart | Open-flow mixing: experimental evidence for strange eigenmodes[END_REF], but on the contrary are well distributed in the channel, including in the middle area. This leads us to choose equal rotation frequencies and angular velocity amplitudes for the three cylinders. Finally, a smooth alternation between the directions of rotation of the cylinders is chosen so as to reduce the risk of wall slip, with a view to experimental studies of the mixer.

On the basis of many simulations (not shown), a scenario with a smooth alternation of the configurations ( 7) and (8) in Fig. 2 is selected, and we choose the protocol given by

Ω 1 (t) = Ω 2 (t) = Ω 0 sin 2πt T , Ω 3 (t) = -Ω 0 sin 2πt T , (1)
where t represents the time, T the forcing period, and Ω 0 the amplitude of the angular velocity.

We emphasise that only simple qualitative arguments have been used to end up with Eqs. [START_REF] Aamo | Control of mixing by boundary feedback in 2D channel flow[END_REF]. The form (1) of the protocol allows, depending on the parameter choices, to vary significantly the level of mixing, as will be shown in section V. This will give us the possibility of deducing conclusions of a certain generality, concerning the relationships between flow characteristics and the level of mixing.

III. HEURISTIC CONDITIONS FOR FOLDING AND MIXING

In this section, relationships between the dimensions of the RAW mixer and the operating conditions are derived using phenomenological arguments similar to those of Stroock et al. [START_REF] Stroock | Chaotic mixer for microchannels[END_REF] for a passive mixing device, and to those of Ober et al. [START_REF] Ober | Active mixing of complex fluids at the microscale[END_REF] for an active mixing device. This will provide heuristic conditions for folding and mixing. Here we took Re bulk = 1 and K 0 = 10; see section III A.

A. Non-dimensional numbers

With the given bulk velocity v bulk as reference velocity and the width channel W as reference length (see Fig. 1), a bulk Reynolds number is defined as

Re bulk = v bulk W ν , ( 2
)
where ν is the fluid viscosity, and a bulk Péclet number is defined as

Pe bulk = v bulk W D , (3)
where D is the molecular diffusivity. Next, introducing a characteristic cross-flow velocity denoted by v cross , a Reynolds number associated to the channel cross flows is defined as

Re cross = v cross W ν . (4)
The ratio

K = Re cross Re bulk = v cross v bulk , (5)
page 5

which represents the strength of the cross-flow velocity to the bulk velocity, will be used below, along with the bulk Strouhal number defined as

St bulk = W T v bulk . (6)
Our choice for the definition of v cross is given later.

B. Condition for folding

The stirring mechanism begins when the interface between the two different values of the concentration, assumed to be initially located at mid-width of the channel, is deformed toward the walls by the action of the cylinders rotation. Then, mixing efficiency depends on the relative strength of the bulk and cross-flow velocities. If the cross-flow velocity is too low with respect to the bulk velocity, the deformation of the interface will be almost suppressed by the next rotational change of direction. In such a case, no folding operation is realised. Otherwise, the pulled thread is long enough to be folded by the next rotation of the cylinders in the opposite direction. In terms of displacement, a minimal condition for achieving a folding of the interface is that the characteristic transverse displacement during one half period, v cross T /2, of a tracer initially located near a rotating arc-wall, exceeds the half-channel width W/2. In this case, however, the length of the generated striations can be negligibly small. To avoid this situation, we rather consider the condition of a characteristic displacement during a half period larger than the full channel width W . The length of the generated striations is then larger than the half-channel width. This is a reasonable condition to allow successive foldings of the concentration interface. With Eq. ( 6), this results in

K ≥ 2St bulk . (7)
To obtain from (7) a manageable condition, we have to clarify our choice for the value of v cross . Since the condition (7) concerns the transverse displacement of a tracer during one half period T /2, a definition of the characteristic cross-flow velocity v cross that is relevant here is by averaging the angular velocity on this duration. Eq. ( 5) and the condition (7) then give

K 0 ≥ πSt bulk , (8)
where K 0 is defined by

K 0 = RΩ 0 v bulk . (9)
The inequality [START_REF] Bottausci | Mixing in the shear superposition micromixer: three-dimensional analysis[END_REF] provides guidance on a minimal condition for folding needed to achieve mixing. However, as will be illustrated in section V A, this is only a necessary (but not sufficient) condition. Another condition that takes into account the diffusion process, thus complementary to the condition [START_REF] Bottausci | Mixing in the shear superposition micromixer: three-dimensional analysis[END_REF], is discussed in the next section.

C. Condition for mixing

If the cross-flow velocity is large enough with respect to the bulk velocity, then the successive changes of direction of the rotating cylinders can generate stretched and folded lamellae of the concentration values down to spatial scales

d = W 2 n . ( 10
)
Since the interface between the concentration values is initially located at mid-width of the channel, the number of cycles of stretching and folding operations is n -1, with n ≥ 1. The characteristic time needed for the diffusion to be effective after n -1 cycles is then

t diffusion = d 2 D = W 2 4 n D . ( 11
)
page 6

A phenomenological condition for the diffusion to be effective after n -1 cycles is

t residence ≥ t diffusion , (12)
with the characteristic residence time defined as

t residence = L mixing v bulk , (13)
where L mixing is the length of the active mixing zone (see Fig. 1). Using the bulk Péclet number of Eq. (3), the characteristic residence time can be written as

t residence = W L mixing D Pe bulk . ( 14
)
With the diffusion time given by Eq. ( 11) and the residence time given by Eq. ( 14), the mixing condition [START_REF] Castelain | Experimental study of chaotic advection regime in a twisted duct flow[END_REF] then becomes

n ≥ 1 ln(4) ln Pe bulk ζ , (15)
where ζ = L mixing /W is the non-dimensional extent of the mixing zone. Noting that under the folding condition of section III B,

ζK = t residence v cross W ≥ n, (16)
the condition [START_REF] D'alessandro | Control of mixing in fluid flow: a maximum entropy approach[END_REF] gives

K ≥ 1 ζ ln(4) ln Pe bulk ζ . (17)
Since in ( 16) t residence v cross is the transverse displacement of a tracer during the characteristic residence time t residence , a relevant definition of the characteristic cross-flow velocity v cross in ( 17) is now by averaging the angular velocity on t residence :

v cross = RΩ, with Ω = 1 t residence t residence 0 Ω 0 sin 2πt T dt. ( 18
)
We shall find it useful to write the characteristic residence time as

t residence = N T + δt, (19)
where 0 ≤ δt < T . Then,

Ω =        Ω 0 πζSt bulk 2N + sin 2 (πζSt bulk ) , if 0 ≤ δt ≤ T 2 , Ω 0 πζSt bulk 2(N + 1) -sin 2 (πζSt bulk ) , if T 2 < δt < T . (20)
Noting that 2(N + 1) -sin 2 (πζSt bulk ) ≥ 2N + sin 2 (πζSt bulk ), the inequality

Ω ≥ D RW ln(4) Pe bulk ζ ln Pe bulk ζ , (21)
which is equivalent to [START_REF] De Moura | Mixing and chaos in open flows[END_REF], is satisfied by both expressions in [START_REF] Omari | Chaotic mixing in a plane channel with rotating arc walls[END_REF] when the following mixing condition is fulfilled: The inequality ( 22) may be thought of as a guidance on a necessary condition for mixing, at any given value of the bulk Péclet number Pe bulk . An example of obtained separating curve in the parameter plane (St bulk , K 0 ) is shown in Fig. 3 (black line). Below the curve are values of St bulk and K 0 that can be considered unfavorable for mixing. Accounting for both conditions (8) and ( 22), values of St bulk and K 0 that would allow mixing may be determined (coloured zone in Fig. 3). However, at this point we cannot strictly guarantee mixing for such values. This issue is investigated through numerical experiments in the next sections, by comparing the results obtained from the heuristic conditions ( 8) and [START_REF] Omari | Thermal chaotic mixing of power law fluids in a mixer with alternately-rotating walls[END_REF] to those obtained by direct simulations of the mixing process.

K 0 ≥ πSt bulk ln(4) ln Pe bulk ζ 2ζSt bulk + sin 2 (πζSt bulk ) . ( 22

IV. NUMERICAL SETUP

In this section are presented the settings and the numerical method used for the direct simulations to be compared to the heuristic conditions ( 8) and ( 22) for folding and mixing.

A. Equations for the velocity and concentration fields

For the system we are interested in, the RAW mixer with a Newtonian fluid, we choose the channel width W as length scale, and the forcing period T as time scale. The amplitude of tangential velocity near the rotating arcs v 0 = RΩ 0 is chosen as velocity scale. A viscous scaling is taken for the pressure, with µv 0 /W . Assuming an incompressible flow, we thus need to solve, here in a two-dimensional domain, the continuity equation

∇ • v = 0, (23)
the momentum equation

Re 0 [St 0 ∂ t v + (v • ∇)v] = -∇p + ∇ 2 v, (24)
page 8

and the equation governing the advection and diffusion of the concentration c,

St 0 ∂ t c + v • ∇c = 1 Pe 0 ∆c. ( 25
)
In the dimensionless equations ( 23)-( 25), v is the velocity, p is the pressure, and

Re 0 = v 0 W ν , St 0 = W T v 0 , and Pe 0 = v 0 W D . ( 26
)
We conclude that for a fixed value of the Péclet number Pe 0 , two independent non-dimensional parameters in Eqs. ( 23)-( 25) are sufficient for describing the system: Re 0 and Re 0 St 0 . Turning to the inlet conditions, two fluid streams with two distinct concentration levels, c = 0 and c = 1, are evenly injected into the channel. The velocity field at the inlet is taken in the form of a Poiseuille profile with Re bulk = v bulk W/ν = 1. Thus Re 0 = Re 0 /Re bulk = K 0 and Re 0 St 0 = W 2 /(νT ) = Re bulk St bulk = St bulk . From the above non-dimensionalization, we conclude that K 0 and St bulk are two independent parameters able to describe the system. Their relevance as control parameters of the RAW mixer in the folding and mixing conditions obtained in sections III B and III C is thus confirmed.

At the initial time t = 0, the interface between the two fluid streams is placed at the middle of the channel, along the x-axis (see Fig. 4a). As boundary conditions we use no slip conditions at both the static and moving walls.

B. Dimensions, numerical method, and mixing indicator

With the channel width W as reference length scale, the non-dimensional geometrical parameters of the domain are as follows. The channel length is 5, the radius of the cylinders is 0.5, and their penetration depth is 0.167.

The used numerical method is based on a finite-volume formulation with co-located variable arrangement on unstructured meshes, implemented in the in-house code Tamaris. For technical details, the interested reader is referred to [START_REF] Bammou | A numerical study of the longitudinal thermoconvective rolls in a mixed convection flow in a horizontal channel with a free surface[END_REF][START_REF] Omari | Alternate rotating walls for thermal chaotic mixing[END_REF][START_REF] Omari | Thermal chaotic mixing of power law fluids in a mixer with alternately-rotating walls[END_REF][START_REF] Guer | Chaotic advection for thermal mixing[END_REF]. Pressure and velocity are coupled using the SIMPLE algorithm combined with a momentum interpolation method. The transient terms in Eqs. ( 24) and ( 25) are discretised using the second-order backward Euler scheme (BDF2). The diffusion, convection and pressure terms are discretised using second-order accurate schemes. For the discretisation of the convection fluxes, the nonlinear CUBISTA scheme is chosen for its low numerical diffusion allowing for accurate representation of the concentration transport.

Generally speaking, in fluid flow simulations the effective overall diffusion encompasses the molecular diffusion in the fluid, which is often prescribed as a numerical setting, and the numerical diffusion inherent to the chosen schemes, mesh size and CFL number, which is often difficult to assess. In the results shown hereafter, the molecular diffusion D is set such that the corresponding bulk Péclet number Pe bulk (see Eq. ( 3)) is equal to 10 6 . We use a 59 000-cell grid with CFL = 0.7. The effective Péclet number, which is thus less than 10 6 due to the numerical diffusion, is the same in all simulations. Thus, we do not focus on understanding the dependence of the mixing efficiency on the effective Péclet number.

The measure of the level of homogeneity of a mixture at finite Péclet number often requires complementary approaches and depends on the relative importance of advection and diffusion processes [START_REF] Danckwerts | The definition and measurement of some characteristics of mixtures[END_REF][START_REF] Finn | Mixing measures for a two-dimensional chaotic Stokes flow[END_REF][START_REF] Foures | Optimal mixing in two-dimensional plane Poiseuille flow at finite Péclet number[END_REF][START_REF] Lin | Optimal stirring strategies for passive scalar mixing[END_REF][START_REF] Mathew | A multiscale measure for mixing[END_REF][START_REF] Thiffeault | Using multiscale norms to quantify mixing and transport[END_REF]]. In the present work, we adopt for its simplicity and relevance at finite Péclet number the widely accepted measure by the standard deviation of the concentration c, σ = 1

I a I I a I (c I -c mean ) 2 , ( 27
)
with c I the concentration in cell I, a I the area of cell I, and c mean = a I c I /( I a I ). The standard deviation σ is evaluated on a rectangular area of non-dimensional width 0.285 distributed over the entire width of the canal downstream of the three cylinders; its location is shown in Fig. 4a.

V. RESULTS

In this section, the basic mixing mechanisms in the proposed RAW mixer are first described. Then, the mixing efficiency obtained by the phenomenological and numerical approaches are compared.

page 9

A. Typical patterns of the concentration A typical evolution of the two-dimensional field of the concentration in the RAW mixer is shown in Fig. 4, where St bulk = 1 and K 0 = 20. If we refer to the parameter plane in Fig. 3, these values of the parameters are clearly favorable for mixing from the heuristic conditions ( 8) and [START_REF] Omari | Thermal chaotic mixing of power law fluids in a mixer with alternately-rotating walls[END_REF]. Indeed, as can be observed in Fig. 4, the cross-flow velocity is sufficiently large with respect to the bulk velocity for complex crossing trajectories to be generated, along with backward flows. Thus, in this case, the heuristic conditions and the direct simulation are in agreement in terms of mixing efficiency.

As shown in Fig. 5, the concentration field is approximately periodic after about 10 periods of rotation of the cylinders. Then, the standard deviation of the concentration oscillates around a constant value (not shown). Recurrent patterns clearly appear after few periods of rotation of the cylinders (compare Figs. 4e and4f), as observed by Gouillart et al. [START_REF] Gouillart | Open-flow mixing: experimental evidence for strange eigenmodes[END_REF] for another open flow active mixer. Few non mixed regions do remain at the exit of the mixer, as shown in Figs. [START_REF] Aref | Frontiers of chaotic advection[END_REF] (see the right boundary of each frame) and 5b. Note that the RAW mixer with three cylinders shown in Fig. 1 can be considered as a single stirring cell unit. Thus, if the level of mixing is not satisfactory, it can be increased by lining up other RAW stirring units. FIG. 4: Evolution of the concentration c along the RAW mixer, for St bulk = 1 and K 0 = 20. Based on the bulk velocity v bulk , the duration needed for a tracer to cross the mixer is 5T . The solution may be considered as approximately periodic after about 10 periods of rotation of the cylinders (see Figs. 4e and4f). The colour code is: blue for c = 0, yellow for c = 0.5 and red for c = 1. In Fig. 4a, the white dashed lines indicate the rectangular area where the standard deviation σ of the concentration is evaluated.

We now describe the basic mixing mechanisms in the RAW mixer. Considering for clarity the first period of cylinder rotation, this period can be divided into two phases (see Figs. 4c and4d, to be compared to Figs. 2(7) and 2(8), respectively). During the first phase, which corresponds to the first half period, three uniform concentration regions are positioned opposite the rotating cylinders on the three elliptic regions visible in Fig. 2(7). The corresponding three uniform concentration regions take the form of anvils (red, blue, then red again), separated by narrower areas each consisting of two bands of different concentrations (see Fig. 4c). In subsequent periods, concentration regions are no longer monochrome, but made up of increasingly complex concentration gaps. However, the mechanism remains basically the same: During each first half period, the three elliptic regions serve as "waiting rooms" for concentration areas, before stretching them during the second half period. But before this second phase, the cylinders change direction of rotation at the half period. At that time, the elliptic regions present during the first half period disappear and the concentration areas that were on the elliptic regions are advected by the longitudinal flow. Then, as the angular velocity of the cylinders increases, three new elliptic regions emerge near the cylinders (see Fig. 2(8)); the concentration areas that were in the three elliptic regions during the first half period are now entrained between the rotating cylinders and the new elliptic regions. As these regions are very close to the cylinders, the concentration areas are then strongly stretched. In addition, the bands of concentrations separating the areas above mentioned are wound inside the new elliptic regions. Striations of different concentration levels are then generated. We thus observe page 10 that the combined action of elliptic regions close to the walls is a key condition for mixing in the RAW mixer, since they allow for the generation of concentration filaments and their folding. The RAW mixer then takes advantage of the possible beneficial role of elliptic regions for mixing in open flows [START_REF] De Moura | Mixing and chaos in open flows[END_REF]. FIG. 6: Evolution of the concentration along the RAW mixer for St bulk = 2 and K 0 = 5. This case is unfavorable for mixing if we refer to the parameter plane in Fig. 3. Based on the bulk velocity v bulk , the duration needed for a tracer to cross the mixer is 10T . The color code is: blue for c = 0, yellow for c = 0.5 and red for c = 1.

We now consider the parameters St bulk = 2 and K 0 = 5, which are, referring to Fig. 3, unfavorable for mixing. Indeed, no mixing is observed in Fig. 6, except in a thin area in the vicinity of the interface between the values 0 and 1 of the concentration. The concentration patterns in Figs. 4 and6 are completely different. First, the shapes of the interfaces are almost identical in Figs. 4d, 4e and4f (despite the fact that the values of the concentration fields in Fig. 4d are different from the values in Figs. 4e or 4f). However, the shape of the interface in Fig. 6b is different from the shapes in Figs. 6c or 6d. This means that, in such a case unfavorable for mixing, the duration needed to observe an approximately periodic regime is longer than in case favorable for mixing. Second, in the case favorable for mixing (Fig. 4), the interface is pushed to the walls thanks to the transversal velocity component which is sufficiently large with respect to the bulk velocity in that case. There is a noticeable difference from the results shown in Fig. 6, where the interface remains confined near the middle of the channel.

According to the parameter plane in Fig. 3, a parameter choice less favorable for mixing than (St bulk = 1, K 0 = 20) but more favorable than (St bulk = 2, K 0 = 5) is (St bulk = 2, K 0 = 20). It is apparent in Fig. 7 that the characteristic displacement of a tracer located near the rotating cylinders ( 2) and (3) during a half period is larger than the half channel width. Note that the requirement of this minimal displacement was used in section III C for the derivation page 11 of the mixing condition [START_REF] Omari | Thermal chaotic mixing of power law fluids in a mixer with alternately-rotating walls[END_REF]. However, it is clear from Figs. 7d and 8 that this condition is insufficient for mixing, since the mixing level is relatively low at the outlet of the channel. Although large transverse displacements of the tracers are observed, the deformations of the interface are somewhat canceled out by the sequence of rotations of the cylinders.

Finally, since the black line in Fig. 3, which represents the mixing condition [START_REF] Omari | Thermal chaotic mixing of power law fluids in a mixer with alternately-rotating walls[END_REF], has the vertical asymptote St bulk = 0, it is of interest to consider an example with St bulk 1. We thus take St bulk = 0.1 and K 0 = 20. The obtained results are quasi-identical to those without cylinder rotation (not shown). This means that the vertical asymptote St bulk = 0 of the phenomenological mixing model is meaningful, even if, for St bulk 1, it is difficult to infer from the model precise parameter values ensuring mixing efficiency.

B. Residence time distributions

We now turn our attention to more quantitative Lagrangian flow analysis. Insight into the mixing mechanisms in open flows results from analysis of the distribution of the residence times of tracers in the mixing zone. The residence time of a tracer in the mixing zone is the time spent by the tracer within the mixing channel [START_REF] Nauman | Residence time theory[END_REF]. For simplicity, we monitor the evolution of Lagrangian tracers along the entire flow channel. At the initial time t = 0, 10 000 tracers are uniformly distributed over the inlet line. Their trajectories are determined by solving:

ẋ(t) = v [x(t), t] , x(0) = x 0 ,
which describes the evolution of the positions x of the tracers released at x 0 , and represents the fact that the rate of change of position of each tracer is given by its velocity. Here, v is the velocity field obtained by the direct simulations.

We first consider the flow without rotation of the cylinders (K 0 = 0). This is approximately a Poiseuille flow, which can be seen as the worst of the cases considered in the present study in terms of mixing efficiency, since it exhibits no transverse displacements of the tracers. A log-log plot of the residence time distribution, shown in Fig. 9, suggests page 12 that the statistical distribution follows approximately a power law t -γ with an exponent γ between -2 and -2.5. Adapting the calculation proposed by Raynal and Carrière [START_REF] Raynal | The distribution of "time of flight" in three dimensional stationary chaotic advection[END_REF] for the 3-D Poiseuille flow to the 2-D case, we obtain that the exponent is equal to -3. Given that the cylinder protuberances in the mixer channel have the effect of slowing down the evacuation of the tracers downstream in comparison to the channel considered in the theoretical case, the value of the exponent obtained numerically is in agreement with the theoretical value. The approximate algebraic -2.5 decay rate is also evidenced in Fig. 10a by comparison with an exponential decay rate. In contrast, the cases with rotations of the cylinders shown in Figs. 10b to 10d exhibit decay rates becoming approximately exponential as the parameters St bulk and K 0 are favorable to the mixing efficiency (in agreement with the literature on chaotic advection in open flows, e.g. [START_REF] Gouillart | Open-flow mixing: experimental evidence for strange eigenmodes[END_REF][START_REF] Tél | Chaotic advection, diffusion, and reactions in open flows[END_REF]). This exponential decay rate is observed within time intervals in which a constant fraction of tracers, η, leave the mixing zone during each period. Indeed, at the end of the cycle n -1 of stretching and folding operations (as described in section III C ), the number of tracers present in the mixing zone is

N = N 0 (1 -η) n-1 .
Using t = (n -1)T and the estimate η ∼ v bulk T /L mixing , we have

N ∼ N 0 exp t T ln 1 - v bulk T L mixing ,
and for ζSt bulk sufficiently large with respect to unity:

N ∼ N 0 e -t /ζ ,
where t = t/(W/v bulk ). Here ζ = 5, so that -1/ζ = -0.2, which is approximately the slope that can be observed in Figs. 10c and 10d within time intervals in which the escape rates of the tracers are approximately constant.

Notice that, as the parameters St bulk and K 0 become favorable for a a high mixing efficiency according to the discussion of section V A, the mean and the standard deviation of the residence times of the 10 000 tracers decrease (see Table II), which can be of major practical interest, particularly in industrial applications.

C. Ageing of seeding tracers

In Fig. 11 are shown the trajectories of 100 tracers coloured by their ages for various values of the parameters (St bulk , K 0 ). At the initial time t = 0, the tracers are uniformly distributed over the inlet line, at the left of each frame. Considering first the case without cylinder rotation (K 0 = 0), it is shown in Fig. 11a that, to a large extent, the uniformity of the tracer distribution is preserved all along the flow channel. At the end of the simulation, the seeding tracers have left the channel, except for few tracers which are "trapped" near the walls. For the case with St bulk = 2 and K 0 = 5 (unfavorable for mixing, as shown in section V A), it is visible in Fig. 11b that the uniformity of the tracer distribution is rapidly lost, even at the left of the cylinder (1) (see Fig. 1). Again, few tracers are trapped near the walls. The same observations can be made for the case with St bulk = 2 and K 0 = 20 in Fig. 11c, which is in between favorable and unfavorable for mixing according to the results shown in section V A. Moreover, loops of the page 14 tracer trajectories begin to appear in the vicinity of the rotating cylinders.

Considering now the case with St bulk = 1 and K 0 = 20, which is clearly favorable for mixing according to the results shown in section V A, it is apparent in Fig. 11d that the ageing distribution of the seeding tracers exhibits irregularities inside the mixing zone, contrary to what is observed in the other cases. The filamental patterns of different colours in Fig. 11d are the mark of the sensitivity to the initial conditions, and hence the signature of chaotic advection. Some filamental patterns of different colors are also visible in Fig. 11e for the case St bulk = 1 and K 0 = 40, which is better for mixing efficiency than the case St bulk = 1 and K 0 = 20 according to Fig. 13. The relationship between mixing efficiency and the level of chaotic advection has to be further investigated, but we can already note that the case favorable for mixing is associated with chaotic advection, which is absent in the cases not clearly favorable for mixing. This supports a well-established general correlation between good mixing capacity and presence of chaotic advection (see e.g. [START_REF] Sundararajan | Transport phenomena in chaotic laminar flows[END_REF]). Incidentally, the results shown in Figs. 11d and 11e also demonstrate that the proposed RAW mixer is able to exhibit chaotic mixing. page 15 The basic difference in the flow dynamics between cases clearly favorable or not for mixing is evidenced in Fig. 12. The trajectories of 50 seeding tracers, uniformly distributed at the initial time over the lower half of the inlet line at the left of each frame, are coloured by the tracer age. It is shown that only in the case favorable for mixing (St bulk = 1 and K 0 = 20, see Fig. 12b), large transverse displacements of the fluid tracers occur on the whole mixing zone. For the direct simulations, at fixed values of K 0 , of the bulk Reynolds number Re bulk , and of the bulk Péclet number Pe bulk , the bulk Strouhal number St bulk dependence of the mixing efficiency represented by σ exhibits a maximum. When the bulk Strouhal number exceeds a critical value, the mixing efficiency decreases for large Strouhal number, and increases with the Strouhal number below the critical value. Thus, the direct simulations conform with theoretical suggestions by Rom-Kedar and Poje [START_REF] Rom-Kedar | Universal properties of chaotic transport in the presence of diffusion[END_REF] for flows with hyperbolic saddles, and with numerical and experimental observations by Horner et al. [START_REF] Horner | Transport enhancement mechanisms in open cavities[END_REF] for flows with parabolic-type saddles.

The variations of the mixing level as a function of the control parameters of the mixer obtained from the direct simulations are in good agreement with the results obtained from the heuristic conditions. Even if the mixing level itself is known only approximately, the heuristic model thus provides useful information on the control parameters to be chosen for reaching mixing efficiency without performing a large number of costly direct simulations. We believe that this is remarkable in itself, as our heuristic model was minimalistic and did not account for all the details of the the complex flow revealed by the numerical simulations.

VI. CONCLUSION

An active mixing device for open flows was presented, along with some numerical results on the mixing of a passive scalar concentration in a Newtonian fluid at low Reynolds number and high Péclet number in the two-dimensional case. The flow modulation protocol is imposed by three rotating circular arc-walls in a straight channel. The selected scenario of a smooth alternation of the direction rotations of the cylinders allows for the stretching of concentration areas and the folding of concentration filaments, which are ordered processes known to be suitable for good mixing of a passive scalar concentration. Two simple heuristic conditions based on phenomenological arguments were presented, page 17 6)) and K 0 (see Eq.( 8)), with ζ = 2 and Pe bulk = 10 6 . Each of the 90 circles (o) corresponds to the result of a direct simulation. Linear interpolation is applied subsequently between the results. The black and white lines are the lines shown in Fig. 3, based on the folding condition (8) and the condition [START_REF] Omari | Thermal chaotic mixing of power law fluids in a mixer with alternately-rotating walls[END_REF] for mixing.

one for folding and the other for mixing with scalar diffusion embedded in the model. These involve two control parameters: the Strouhal number based on the bulk flow velocity, and the strength of the cross flow relative to the transport flow. The heuristic conditions lead to values of these parameters that are in good agreement for mixing efficiency with those obtained from the direct simulations. It was thus shown that the conditions could provide useful informations on the control parameters to be selected for mixing efficiency, for example as a first attempt before mixing optimization. In the cases of good mixing efficiency, chaotic advection was detected in the flow both through the study of the tracer trajectories and the residence time distribution profiles. The beneficial role of the elliptic regions of the flow, generated by the rotating walls, for stretching and folding the concentration areas, was evidenced. The active inline RAW mixer is thus promising for producing efficiently chaotic mixing, in particular by the adaptability of the mixing protocols that can be chosen.

page 4 FIG. 2 :

 42 FIG. 2: Typical steady-state streamline flow topologies for the three cylinders rotating at constant angular velocity. The eight possible stirring configurations are 1(+ + +), 2(---), 3(+ + -), 4(--+), 5(-+ +), 6(+ --), 7(+ -+), and 8(-+ -). The symbol (+) indicates counter-clockwise direction of the angular velocity of the rotating cylinder (blue arrows), and (-) clockwise direction (red arrows). The elliptic zones are colored in light blue.Here we took Re bulk = 1 and K 0 = 10; see section III A.

) page 7 FIG. 3 :

 73 FIG. 3: Region of possible mixing (colored) in the parameter space (St bulk , K 0 ), based on the folding condition (8) (blue line) and the mixing condition (22) (black line). Here, we take ζ = 2 and Pe bulk = 10 6 . The region below the black line corresponds to values for which mixing is not expected. The region below the blue line corresponds to values for which folding is not expected. Note that the black line, which represents the mixing condition, has the vertical asymptote St bulk = 0.

  (a) t = 0. (b) t = T /8. (c) t = T /2. (d) t = T . (e) t = 10T . (f) t = 50T .

FIG. 5 :

 5 FIG. 5: Spatio-temporal plots of the concentration c as a function of the time t (abscissa) and the channel height y (ordinate), for St bulk = 1 and K 0 = 20. The colour code is: blue for c = 0, yellow for c = 0.5 and red for c = 1. Two cross-sections located at different axial coordinates x are considered along the channel of the RAW mixer: (a) cross-section passing by the lower cylinder axis (see Fig. 1); (b) outlet cross-section. The concentration data are recorded from time t = 0 to time t = 13.75 T for both frames. The black segment in Fig. 5b represents the duration of one period T .

FIG. 7 :FIG. 8 :

 78 FIG.7: Evolution of the concentration along the RAW mixer for St bulk = 2 and K 0 = 20. This case is in between favorable and unfavorable for mixing if we refer to the parameter planes in Figs.3 and 13. Based on the bulk velocity v bulk , the duration needed for a tracer to cross the mixer is 10T . The color code is: blue for c = 0, yellow for c = 0.5 and red for c = 1.

FIG. 9 :

 9 FIG. 9: Log-log plots of the residence time distribution of 10 000 tracers in the mixing channel without rotation of the cylinder (K 0 = 0) over the duration 50 in time unit W/v bulk . Each bin width represents the duration 0.5 W/v bulk .

  K 0 = 0: no cylinder rotation. 9 494 tracers in 10 000 have left the mixing zone at t = 25. -0.2 (b) St bulk = 2 and K 0 = 5. 9 516 tracers in 10 000 have left the mixing zone at t = 25. -0.2 (c) St bulk = 2 and K 0 = 20. 9 321 tracers in 10 000 have left the mixing zone at t = 25.

-0. 2 (

 2 FIG. 10: Semi-log plots of the residence time distribution of 10 000 tracers in the mixing channel over the duration 25 in time unit W/v bulk . Each bin width represents the duration 0.25 W/v bulk .

FIG. 11 :

 11 FIG. 11: Trajectories of 100 tracers in the mixing channel, colored by their ages. The duration is 25 in time unit W/v bulk .

FIG. 12 :

 12 FIG. 12: Trajectories of 50 tracers in the mixing channel, colored by their ages. The tracers are injected only over the lower half of the inlet line at the left of the frames. The duration is 25 in time unit W/v bulk .

FIG. 13 :

 13 FIG.13: Standard deviation of the concentration c as a function of St bulk (see Eq. (6)) and K 0 (see Eq.(8)), with ζ = 2 and Pe bulk = 10 6 . Each of the 90 circles (o) corresponds to the result of a direct simulation. Linear interpolation is applied subsequently between the results. The black and white lines are the lines shown in Fig.3, based on the folding condition (8) and the condition[START_REF] Omari | Thermal chaotic mixing of power law fluids in a mixer with alternately-rotating walls[END_REF] for mixing.

TABLE II :

 II Mean (RT) and standard deviation (σ RT ) of the residence time of the tracers, among 10 000, that have left the mixing channel at t = 25 (in time unit W/v bulk ).

VII. ACKNOWLEDGMENTS

The authors acknowledge the financial support from the french research program ANR-16-CE06-0003 through the NaiMYS project (Novel active inline Mixer for Yield Stress fluids). The computations were carried out on page 18

the Avakas cluster of the Mésocentre de Calcul Intensif Aquitain (MCIA) and the Occigen supercomputer of the Centre Informatique National de l'Enseignement Supérieur (CINES). The authors are grateful to A. Achhoud for his participation to the preliminary numerical tests as a part of his MsC studies at University of Pau and Pays de l'Adour.