Heni Dallagi
email: heni.dallagi@univ-amu.fr

Nassim Ait-Mouheb
email: nassim.ait-mouheb@inrae.fr

Audrey Soric
email: audrey.soric@centrale-marseille.fr

Olivier Boiron
email: olivier.boiron@centrale-marseille.fr

Nassim Ait

Keywords: Millifluidic, CFD, Clogging, Labyrinth channel, oscillatory shear index

Water scarcity is a global concern, with irrigation of food crops contributing significantly to freshwater depletion. Drip irrigation technology reduces water consumption but faces issues like clogging in narrow discharge sections, diminishing efficiency, and increasing costs. Accurate prediction of flow characteristics and understanding parameters affecting biofilm growth and particle deposition is crucial for effective anti-clogging strategies. Computational fluid dynamics (CFD) using turbulence models can be a valuable tool. This study evaluated the accuracy and efficiency of different turbulence models (standard k-ε, Reynolds Stress Model, and Large Eddy Simulation) in predicting the flow characteristics of a commercial emitter in a drip irrigation system. Results showed the standard k-ε model as a preferred choice for simulating mean flow characteristics and emitter discharge due to its balance between accuracy and computational efficiency. However, the Large Eddy Simulation model provided the most accurate results, considering the emitter discharge, unsteady flow behavior, wall shear stress distribution, and oscillatory index, despite requiring more computational resources. This model is valuable for understanding hydrodynamic effects on emitter clogging. The study also investigated the impact of velocity fluctuations, wall shear stress, and oscillatory shear index on biofilm growth and deposition in the emitter. Low shear stress in inlet and return zones reduced self-cleaning ability, leading to particle and microorganism attachment. Maintaining appropriate wall shear stress values in other regions proved crucial for improving anti-clogging ability. High oscillatory shear index values enhanced mass transfer, nutrient mixing, diffusion within the biofilm, and self-cleaning capacity.

In summary, this study greatly enhances our understanding of how flow dynamics and biofilm management impact drip irrigation systems. It provides practical insights for engineers and practitioners, aiding in the creation of more efficient and clog-resistant systems. By optimizing these dynamics and strategies, this research promotes sustainable water use in agriculture, while also minimizing maintenance costs and maximizing crop yields

Introduction

To date, different irrigation techniques such as surface, sprinkler, and drip irrigation are employed around the world. However, compared to other technologies, drip irrigation has proven to be an efficient technology for improving irrigation efficiency by up to 95%, using 30 -60% less water than surface irrigation, and improving crop yield (20 -50%) by ensuring a uniform wetting of the soil, and also convenient in distributing mineral fertilisers to cultivated fields [START_REF] Abduraimova | Methods for determining water flow from rice field to open drainage system[END_REF][START_REF] Van Der Kooij | The efficiency of drip irrigation unpacked[END_REF]. In addition, drip irrigation is considered the safest irrigation technology for wastewater reuse [START_REF] Tripathi | Impact of municipal wastewater reuse through micro-irrigation system on the incidence of coliforms in selected vegetable crops[END_REF]. Nevertheless, drippers that generate low flow rates (between 0.5 and 8 l h -1 ) have only narrow discharge sections that can be easily clogged thus perturbing the irrigation uniformity, reducing irrigation efficiency, and inducing additional cost [START_REF] Do Amaral | Dripper clogging by soil particles entering lateral lines directly during irrigation network assembly in the field[END_REF][START_REF] Shi | A Review of the Category, Mechanism, and Controlling Methods of Chemical Clogging in Drip Irrigation System[END_REF].

In general, the structure of the emitters in drip irrigation systems is complex and characterised by a small size flow path ranging from 0.2 to 1.2 mm [START_REF] Al-Agele | Testing Novel New Drip Emitter with Variable Diameters for a Variable Rate Drip Irrigation[END_REF]. Thus, clogging can occur when the irrigation water used contains particles or microorganisms. It has been found in the literature that clogging occurs due to biomass growth, chemical precipitations, and/or physical deposits [START_REF] Petit | Methods for drip irrigation clogging detection, analysis and understanding: State of the art and perspectives[END_REF][START_REF] Shi | A Review of the Category, Mechanism, and Controlling Methods of Chemical Clogging in Drip Irrigation System[END_REF]. Several factors control this process including biological processes (biofilm attachment, growth, and detachment) [START_REF] Rizk | Treated wastewater Reuse in micro-irrigation: Effect of shear stress on biofilm development kinetics and chemical precipitation[END_REF][START_REF] Zhang | Effects of microorganisms on clogging process and clogging substances accumulation of drip irrigation emitters using the high-sediment water sources[END_REF], water quality (pH, conductivity, temperature, and size and particle concentration) [START_REF] Liu | Emitter clogging characteristics under reclaimed wastewater drip irrigation: A meta-analysis[END_REF][START_REF] Oliveira | Influence of Concentration and Type of Clay Particles on Dripper Clogging[END_REF] and hydrodynamic conditions [START_REF] Lequette | Hydrodynamic effect on biofouling of milli-labyrinth channel and bacterial communities in drip irrigation systems fed with reclaimed wastewater[END_REF][START_REF] Li | Research on Hydraulic Properties and Energy Dissipation Mechanism of the Novel Water-Retaining Labyrinth Channel Emitters[END_REF].

Several studies have investigated the impact of hydrodynamic conditions on limiting the growth of biofilms and the deposition of particles within complex geometries, with particular attention being paid to emitters. These studies have revealed that both velocity and wall shear stress play a crucial role in determining the growth rate, detachment rate, and biofilm structure. Specifically, higher velocity values have been found to limit biofilm growth and particle deposition, while lower values can promote their accumulation. While high levels of shear stress can be effective in removing biofilms from surfaces (self-cleaning), there exists a threshold level below which the shear stress does not have a significant effect. Moreover, excessive shear stress can result in the compression of the biofilm and have negative consequences as the transport and renewal of microorganisms and nutrients can be facilitated [START_REF] Gamri | Biofilm development in micro-irrigation emitters for wastewater reuse[END_REF][START_REF] Liu | Microfluidic study of effects of flow velocity and nutrient concentration on biofilm accumulation and adhesive strength in the flowing and no-flowing microchannels[END_REF][START_REF] Paul | Effect of shear stress and growth conditions on detachment and physical properties of biofilms[END_REF][START_REF] Wang | Accumulation mechanism of biofilm under different water shear forces along the networked pipelines in a drip irrigation system[END_REF]. It is important to note that the threshold of shear stress required for biofilm detachment from the walls can depend on various factors, such as the strain of the biofilm, the age of the biofilm, the composition of the biofilm, and the properties of the irrigation emitter material [START_REF] Tsvetanova | STUDY OF BIOFILM FORMATION ON DIFFERENT PIPE MATERIALS IN A MODEL OF DRINKING WATER DISTRIBUTION SYSTEM AND ITS IMPACT ON[END_REF][START_REF] Yan | Biofilm structure and its influence on clogging in drip irrigation emitters distributing reclaimed wastewater[END_REF]. In general, the threshold required for loosening biofilm, without achieving complete detachment, is believed to fall within the range of 0.1 -10 Pa, as demonstrated in previous work [START_REF] Nejadnik | Determination of the Shear Force at the Balance between Bacterial Attachment and Detachment in Weak-Adherence Systems, Using a Flow Displacement Chamber[END_REF][START_REF] Wang | Accumulation mechanism of biofilm under different water shear forces along the networked pipelines in a drip irrigation system[END_REF].

In order to improve anti-clogging ability and to ensure high irrigation uniformity at a low discharge rate, current commercial emitters usually adopt a milli-labyrinth channel approach with discharge areas in the order of 1 mm 2 . This induces the development of a turbulent flow regime inside the channel but also creates some stagnant zones and zones with vorticity [START_REF] Al-Muhammad | Experimental and numerical characterization of the vortex zones along a labyrinth milli-channel used in drip irrigation[END_REF].

However, analysing the flow characteristics in such complex and small channels using traditional experimental methods remains a challenge. Thus, in the field of drip irrigation, Computational Fluid Dynamics (CFD) approaches have been used as a crucial tool for studying hydrodynamics and anticlogging performance. Several experimental studies coupled with numerical simulations were performed to analyse the effect of the labyrinth path on the hydraulic and anticlogging performance of different emitters. For example, a CFD simulations of flow in labyrinth emitters, using laminar or steady k-ε models, concluded that the movements of swirls as well as the stagnant flowing zones existing in the emitter channels were the main reason for emitters clogging [START_REF] Ait-Mouheb | Impact of hydrodynamics on clay particle deposition and biofilm development in a labyrinth-channel dripper[END_REF][START_REF] Li | Research on Hydraulic Properties and Energy Dissipation Mechanism of the Novel Water-Retaining Labyrinth Channel Emitters[END_REF][START_REF] Yavuz | The mechanism of emitter clogging analyzed by CFD-DEM simulation and PTV experiment[END_REF]. Recently, [START_REF] Al-Muhammad | Experimental and numerical characterization of the vortex zones along a labyrinth milli-channel used in drip irrigation[END_REF] analysed the flow in a three baffles-fitted labyrinth channel using the steady Reynolds stress model (RSM) modelling and the micro-PIV technique. Their results showed that the flow in the labyrinth channel is turbulent. Furthermore, the RSM model results have shown good concordance with the experimental data, with a difference of about 5% in the prediction of the mean velocity in the mainstream flow. However, a difference that possibly reaches 75% was observed in the estimation of the mean velocity. Such unsteady occurrences are particularly difficult to predict with the usual steady-state modeling approaches used in previous studies. The transition to an unsteady simulation technique is crucial to improve the accuracy of numerical simulations and to be able to better understand the flow characteristics in a labyrinth channel. In this context, the Large Eddy Simulation (LES) could provide a promising alternative for the common Reynolds-Averaged Navier-Stokes (RANS) turbulence models. Previous studies have utilised LES modeling to investigate the hydraulic performance of drip irrigation emitters, including those equipped with milli-labyrinth channels [START_REF] Feng | Effect of optimization forms of flow path on emitter hydraulic and anti-clogging performance in drip irrigation system[END_REF][START_REF] Li | Analysis of Tracing Ability of Different Sized Particles in Drip Irrigation Emitters with Computational Fluid Dynamics[END_REF][START_REF] Wu | Simulation of the flow characteristics of a drip irrigation emitter with large eddy methods[END_REF]. The results of these studies suggest that the LES model is superior to the standard k-ε (SKE) model in predicting flow velocities within the labyrinth region. However, these analyses have been limited to characterising the discharge curve of the emitter and identifying the vortex regions and main flow pathways within the labyrinth. Furthermore, the studies have not provided a comprehensive investigation of hydraulic parameters, including the unsteady behavior and fluctuation of the shear stress at the wall.

Additionally, there has been no exploration of the potential correlation between these hydraulic parameters and the clogging effect of the emitter. Therefore, further research is needed to fully understand the hydraulic behavior of drip irrigation emitters and optimize their design.

This study aims to simulate the flow in a commercial drip irrigation emitter using three different turbulence models: SKE, RSM, and LES. The simulation was performed using a 3D CFD approach. The objective is to determine the unsteady turbulence model that could provide the most accurate prediction of the flow in the industrial emitter by comparing the results of these simulations with experimental data obtained from the discharge tests. In addition to the comparison of turbulence models, the study also aims to understand the effects of the flow characteristics on the clogging process that occurred in the drip irrigation emitter presented in a previous study [START_REF] Lequette | Hydrodynamic effect on biofouling of milli-labyrinth channel and bacterial communities in drip irrigation systems fed with reclaimed wastewater[END_REF]. The simulation results were analysed in order to identify the regions of the emitter that are most susceptible to clogging and the main hydrodynamic parameters that cause the clogging effect.

Materials and Methods

Tested emitter

The labyrinth-channel geometry selected in this work consisted of a commercial non-pressurecompensating (NPC) emitter (model D2000, Rivulis Irrigation SAS, France). The channel geometry was described in previous work [START_REF] Lequette | Hydrodynamic effect on biofouling of milli-labyrinth channel and bacterial communities in drip irrigation systems fed with reclaimed wastewater[END_REF]. Table 1 summarises its characteristics. As shown in Fig. 1, the fluid flows through a T-shaped inlet (three inlets) and exits through a straight cylinder channel to avoid backflow effects that could affect the upstream flow simulation. It should be noted that this industrial emitter has the particularity of an alternate dental angle between 28° and 37°.

Additionally, this type of emitter includes a return zone, which serves to increase the length of the channel.

The Reynolds number (Re) of the studied emitter was around 305 at the labyrinth inlet which indicates that the flow behavior corresponds to a turbulent regime in such a labyrinth channel (Zhang et al., 2016). The Reynolds number was calculated as follows:

𝑅𝑒 = 𝜌𝑣 𝐷 ℎ µ (1)
Where 𝜌, µ, and 𝑣 are the water density (kg m -3 ), water viscosity (Pa s), and water bulk velocity at the inlet (m s -1 ), respectively. 𝐷 ℎ being the hydraulic diameter (m).

Turbulence model

In this study, the SKE and RSM models were used to simulate the steady-state flow, whilst the LES model was selected for the unsteady state. Liquid water was imposed without taking into account buoyancy and gravity. The exclusion of considerations for gravity and buoyancy is justified by the confined nature of the simulated flow system and the relatively uniform density of the fluid within. In confined conduit flow scenarios, where the cross-section is completely filled with fluid without a free surface, buoyancy effects resulting from significant density variations across the fluid are negligible.

Moreover, whilst gravity contributes to the hydrostatic pressure component in the confined conduit, its influence on altering the actual flow dynamics within the system remains minimal

The RANS governing equations of the continuity and momentum conservation for a steady incompressible flow are expressed as follows:

𝜕𝑢 ̅ 𝑖 𝜕𝑥 𝑖 = 0 (2) 𝜕(𝑢 ̅ 𝑖 𝑢 ̅ 𝑗 ) 𝜕𝑥 𝑗 = - 1 𝜌 𝜕𝑝̅ 𝜕𝑥 𝑗 + 𝜕 𝜕𝑥 𝑗 [ µ 𝜌 ( 𝜕𝑢 ̅ 𝑖 𝜕𝑥 𝑗 + 𝜕𝑢 ̅ 𝑗 𝜕𝑥 𝑖 ) -𝑢 𝑖 ′ 𝑢 𝑗 ′ ] (3)
Where 𝑢 ̅ 𝑖 and 𝑢 ̅ 𝑗 denotes the average velocities in i and j directions (Cartesian coordinate indices), 𝑝̅ is the mean pressure, 𝑢 𝑖 ′ and 𝑢 𝑗 ′ are the fluctuating velocities, and -𝑢 𝑖 ′ 𝑢 𝑗 ′ is the turbulent Reynolds Stress Tensor that needs to be modeled by solving additional transport equations (turbulence kinetic energy k, and/or turbulence dissipation rate ε or specific dissipation rate ω).

In the SKE model, the turbulence kinetic energy (k) and its rate of dissipation (ε), are obtained by solving the following equations [START_REF] Launder | Lectures in mathematical models of turbulence[END_REF]:

𝜕 𝜕𝑥 𝑖 (𝜌𝑢 ̅ 𝑗 𝑘) = 𝜕 𝜕𝑥 𝑗 [(µ + µ 𝑡 𝜎 𝑘 ) 𝜕𝑘 𝜕𝑥 𝑗 ] + 𝐺 𝑘 -𝜎𝜀 (4) 𝜕 𝜕𝑥 𝑖 (𝜌𝑢 ̅ 𝑗 𝜀) = 𝜕 𝜕𝑥 𝑗 [(µ + µ 𝑡 𝜎 𝜀 ) 𝜕𝜀 𝜕𝑥 𝑗 ] + 𝐶 𝜀1 𝜀 𝑘 𝐺 𝑘 -𝐶 𝜀2 𝜌 𝜀 2 𝑘 (5)
Where:

𝐺 𝑘 = µ 𝑡 𝑆 2 ; 𝑆 = √2𝑆 ̅ 𝑖𝑗 𝑆 ̅ 𝑖𝑗 ; µ 𝑡 = 𝜌𝐶 µ 𝑘 2 𝜀
The following default values are used for the model constants (σk, σε, Cε1, Cε2, Cµ): 1.0, 1.3, 1.44, 1.92, and 0.09, respectively [START_REF] Al-Muhammad | Experimental and numerical characterization of the vortex zones along a labyrinth milli-channel used in drip irrigation[END_REF].

In RSM, the eddy viscosity approach is avoided and the individual components of the Reynolds stress tensor are directly computed [START_REF] Gibson | Ground effects on pressure fluctuations in the atmospheric boundary layer[END_REF][START_REF] Launder | Second-moment closure: Present… and future?[END_REF]. Note that the RSM model does not assume that the flow is turbulent isotropic as the Reynolds stress tensor 𝑢 𝑖 ′ 𝑢 𝑗 ′ is usually anisotropic.

This anisotropy results from the properties of turbulent production, dissipation, transport, pressurestain-rate, and the viscous diffusive tensors. In this study, the Linear Pressure-Strain Model derived by [START_REF] Gibson | Ground effects on pressure fluctuations in the atmospheric boundary layer[END_REF] was also chosen to simulate the flow. Thus, the transport equations of the turbulent energy (k) and turbulent dissipation rate (ε) are as follows:

𝜕 𝜕𝑥 𝑖 (𝜌𝑢 ̅ 𝑗 𝑘) = 𝜕 𝜕𝑥 𝑗 [(µ + µ 𝑡 𝜎 𝑘 ) 𝜕𝑘 𝜕𝑥 𝑗 ] + 1 2 𝐺 𝑘 -𝜎𝜀 (6) 𝜕 𝜕𝑥 𝑖 (𝜌𝑢 ̅ 𝑗 𝜀) = 𝜕 𝜕𝑥 𝑗 [(µ + µ 𝑡 𝜎 𝜀 ) 𝜕𝜀 𝜕𝑥 𝑗 ] + 1 2 𝐶 𝜀1 𝜀 𝑘 𝐺 𝑘 -𝐶 𝜀2 𝜌 𝜀 2 𝑘 (7)
Where:

𝜀 𝑖𝑗 = 2 3 𝜌𝜀 𝛿 𝑖𝑗 ; 𝐺 𝑘 = µ 𝑡 𝑆 2 ; µ 𝑡 = 𝜌𝐶 µ 𝑘 2 𝜀
The following default values are used for the model constants (σk, σε, Cε1, Cε2, Cµ): 0.82, 1.0, 1.44, 1.92, and 0.09, respectively [START_REF] Lien | Assessment of turbulence-transport models including nonlinear rng eddy-viscosity formulation and second-moment closure for flow over a backwardfacing step[END_REF].

LES is a numerical simulation method for turbulence that lies between DNS (Direct Numerical Simulation) and RANS. Its basic concept is the direct solution of the transient Navier-Stokes equation of turbulent motion that is larger than the grid scale. For smaller-scale eddies, that tend to be isotropic, the subgrid-scale is introduced into the Navier-Stokes equation to reflect its effect on the larger-scale eddy. In the LES, the continuity and momentum equations for incompressible flow are calculated as follows:

𝜕𝜌 𝜕𝑡 + 𝜕 𝜕𝑥 𝑖 (𝜌𝑢 ̅ 𝑖 ) = 0 (8) 𝜕 𝜕𝑡 (𝑢 ̅ 𝑖 ) + 𝜕(𝑢 ̅ 𝑖 𝑢 ̅ 𝑗 ) 𝜕𝑥 𝑖 = - 1 𝜌 𝜕𝑝̅ 𝜕𝑥 𝑗 + 𝜈 𝜕 2 𝑢 ̅ 𝑖 𝜕𝑥 𝑖 𝜕𝑥 𝑗 - 𝜕𝜏 𝑖𝑗 𝜕𝑥 𝑗 (9)
Where 𝜏 𝑖𝑗 is the subgrid-scale stress tensor. It is computed from the following equation:

𝜏 𝑖𝑗 = -2µ 𝑡 𝑆 ̅ 𝑖𝑗 + 1 3 𝜏 𝑖𝑗 𝛿 𝑖𝑗 (10) 𝑆 ̅ 𝑖𝑗 = 1 2 ( 𝜕𝑢 ̅ 𝑖 𝜕𝑥 𝑗 + 𝜕𝑢 ̅ 𝑗 𝜕𝑥 𝑖 ) (11)
In this study, the eddy-viscosity is modeled using the Dynamic Smagorinsky-Lilly Model [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF][START_REF] Lilly | A proposed modification of the Germano subgrid-scale closure method[END_REF]:

µ 𝑡 = 𝜌𝐿 𝑠 2 |𝑆 ̅ | ; |𝑆 ̅ | = √2𝑆 𝑖𝑗 𝑆 𝑖𝑗
𝐿 𝑠 is the mixing length for subgrid scales and is computed, in ANSYS FLUENT, using: 𝐿 𝑠 = min (𝑘𝑑, 𝐶 𝑠 ∆),

where k is the von Kármán constant, d is the distance to the closest wall, Cs is the Smagorinsky coefficient, and ∆ is the local grid scale.

In commercial code ANSYS Fluent, the local grid scale was computed according to the volume of the computational cell (V) as ∆ = 𝑉 1/3 . Cs was dynamically computed based on the information provided by the resolved scales of motion; however, it was averaged over homogenous flow directions to avoid any numerical instability. The idea behind the dynamic process was to apply a second filter to the motion equations. The width of the new filter ∆ ̂ was equal to double the grid filter width ∆. More details can be found in [START_REF] Kim | Application of the localized dynamic subgrid-scale model to turbulent wall-bounded flows[END_REF].

Numerical approach and boundary conditions

software ANSYS/Fluent V2020 R2. The simulation of fluid flow in the flow path was performed in threedimensional modeling geometry. Two types of mesh were adopted. In the case of RANS models, 5,111,375 hexahedral cells were used. However, in order to improve the accuracy of the solution and CPU time (≈ 40 s time step -1 ), only 1,842,372 polygonal cells were used in the case of LES simulation.

In all simulations, the mesh quality and the wall functions were satisfied, and the independence of the results to the meshing was verified.

The boundary conditions used were represented by the velocity inlet, pressure outlet, and wall conditions. A velocity of 0.092 m s -1 (which corresponds to a flow rate of around 1 l h -1 and Re of 305)

and atmospheric pressure conditions were imposed at the inlet and the outlet respectively. The noslip boundary condition was employed at the walls where the velocity was set as zero, while the enhanced wall treatment approach was used to model the near-wall region [START_REF] Chen | Near-wall turbulence models for complex flows including separation[END_REF]. To compute the turbulence effect for the SKE, RSM, and LES models, a low turbulence intensity of 5% and hydraulic diameter of 1.02 mm were selected and applied in a normal to boundary direction. It should be noted that this parameter has not shown a significant effect on the results (data not shown). Indeed, the use of lower and higher turbulence intensity (1% and 15%, respectively) resulted in 0.1 and 0.5% deviation of the velocity profile inside a flow path.

The turbulent kinetic energy (k) and the turbulent intensity are selected as additional conditions for the LES, and RSM models. The vortex method was performed for LES simulation to enable the generation of synthetic turbulence at the inlets. The number of vortices was set to N/4, where N is the number of cell faces at the inlets, as was suggested [START_REF] Gerasimov | Quick Guide To Setting Up LES-type Simulation[END_REF].

In this study, the momentum and pressure are resolved as the primary variables using the pressurebased segregated solver. The SIMPLE and STANDARD methods were used for pressure-velocity coupling and pressure interpolation, respectively. Third-order MUSCL discretisation is selected to solve the momentum, Rij, and dissipation rate equations and the Least Squares Cell-Based method was selected to calculate gradients. The scheme of gradient calculation was presented in detail by [START_REF] Al-Muhammad | Experimental and numerical characterization of the vortex zones along a labyrinth milli-channel used in drip irrigation[END_REF]. Further, the bounded second-order implicit was used for the temporal discretization in the case of LES.

For the unsteady LES model, the simulation was performed for a time flow of 1 s with a time step of 10 -5 s. However, this model has a higher computational cost and required a long CPU time. For this purpose, the time step was changed with the simulation time, as was suggested in [START_REF] Gerasimov | Quick Guide To Setting Up LES-type Simulation[END_REF][START_REF] Murthy | Assessment of standard k-ε, RSM and LES turbulence models in a baffled stirred vessel agitated by various impeller designs[END_REF]. In fact, simulations were initially performed with the time step size of 10 -5 s where the corresponding CFL number was 1. As the solution progressed, the time step size gradually increased up to 0.0001 s leading to an important gain in CPU time. The convergence accuracy of the solutions was fixed at least 10 -5 , and the residual can reach 10 -12 in some cases. To ensure smooth and better convergence in the LES case, initially, the SKE simulation was performed until reaching the complete steady-state flow field, then the synthetic turbulence was superimposed on the mean flow in the computational domain (using "init-turb-vel-fluctuations" command) and then switch the solutions to the LES model [START_REF] Gerasimov | Quick Guide To Setting Up LES-type Simulation[END_REF]. The solution was saved at intervals of 100-time steps, whereas the reports for selected parameters such as velocity fluctuation and oscillatory shear index were recorded at every time step, with a frequency of 10,000 Hz. In the context of comparing steadystate (SKE and RSM) and unsteady (LES) flow regimes, it's crucial to recognize the significance of these differences. Even within the realm of unsteady-state flow simulations carried out using the LES model, consistency was aimed to be maintained by computing all data as average values across multiple time steps. This particular approach was adopted to harmonise the interpretation of unsteady-state data with the methodology employed for the RSM and k-ε models. These latter models operate on the principles of RANS, wherein the pertinent variables are derived as time-averaged values. By employing this strategy, the intent was to facilitate a fair and meaningful comparison of prediction accuracy amongst the various flow states.

Validation approach and selected results

Emitter discharge

In order to validate the simulations, emitter discharge was determined numerically and experimentally depending on the inlet pressure. For this purpose, a series of pressure flow rate tests on the D2000 drip irrigation emitter was piloted under working range pressures. Therefore, the discharge-pressure curves were presented and then the relation between flow rate and pressure drop was identified. This relationship is of particular importance as it defines the emitter's performance. It was described by

Karmeli's power law [START_REF] Karmeli | Classification and flow regime analysis of drippers[END_REF]:

𝑞 = 𝑘 𝑑 ∆𝑃 𝑥 ( 12
)
where q is the flow rate (l h -1 ), kd is the proportionality coefficient, ∆P is the working pressure head at the emitter (Pa), and x is the emitter discharge exponent. Note that for long-path emitters, the values of the emitter discharge exponent (x) are between 1 and 0.5 which are for laminar and fully turbulent flows. Noting that at 100 kPa pressure, the average flow rate measured for the emitter was 1.2 l h -1 .

The flow rate (q) was determined at the outlet channel and ∆P represents the average change in inlet minus outlet pressure drop. Experimentally, the discharge test was performed using Irstea's protocol according to the International Organization for Standardization (ISO 9260, 1991;ISO 9261, 2004). 25 emitters were sampled from a set of 500 emitters and divided into 5 lines as shown in Fig. 2. Different pressure values were set at the inlet of the pipe, then the flow rate was measured at the outlet of each emitter. The pressure drops (ΔP) were identified using a pressure sensor (CERABAR PMC 731; Endress+Hauser Huningue, France; accuracy of 0.2%), while the flow rate was determined by measuring the volume of water (dV) in a test tube as a function of time (q = dV/dt). The volume was calculated by multiplying the evolution of the water column length (dH) by a calibration coefficient (ai) corresponding to the section of the test tube (dV = ai.dH).

The calibration coefficient (ai) was determined by plotting a curve that represents the variation of volume with respect to the pressure at the bottom of the test tube (using a pressure sensor (CERABAR PMC 134; Endress+Hauser Huningue, France; accuracy of 0.2%)). During the calibration process, the volume increment for each measurement point was determined by weighing. It should be noted that the flow rate uniformity of drip emitters was checked. The coefficient of variation (CV), which refers to the degree of variation in the flow rate between the different emitters in the system (tested over 500 emitters), was considered an excellent coefficient of 3.11% [START_REF] Yavuz | The mechanism of emitter clogging analyzed by CFD-DEM simulation and PTV experiment[END_REF].

Fig. 2. Schematic representation of a pilot rig illustrates emitter discharge was evaluated by measuring flow rates and pressure losses in 25 emitters across 5 parallel lines.

Velocity fluctuation

Turbulence in fluid flows is often caused by the presence of large-scale eddies, which transfer energy to smaller eddies through a process known as energy cascade. These smaller eddies eventually dissipate the energy through frictional forces and generate a highly fluctuating flow (Stephen, 2001).

For this purpose, velocity data (x = 12 mm, y = 3.41 mm, z = 0.4 mm at the top and x = 12 mm, y = 1 mm, z = 0.4 mm at the bottom; the two positions are presented by the + character in Figure 1) were collected at each time step during a 1.6 s of time flow and then analyzed using the statistical analysis "OPEN_FPE_IFT", an open-source MATLAB package (Fuchs et al., 2022). The stationarity of the data was verified and plotted (Figure 8.a). For both positions, the turbulence intensity (Ti) was calculated as Ti = 100 σ/〈U〉; where σ is the standard deviation of the data and 〈U〉 is the mean value of a velocity time series. Furthermore, a comprehensive view of the data through the probability density function (PDF) distribution, the range (i.e., the difference between the maximum and minimum values), skewness, and flatness were presented according to (Fuchs et al., 2022).

Oscillatory Shear Index calculation

In this article, the Oscillatory Shear Index (OSI), is a fundamental element of our research and is introduced as representing an innovative tool for deepening our understanding of fluid flow dynamics.

OSI is a metric used to quantify the intensity of shear forces experienced by fluid or particles in a fluid flow system. Our study focuses on harnessing the OSI to advance our insights into biofilm growth and clogging within drip irrigation systems. Our primary goal is to evaluate the OSI's effectiveness in elucidating the processes involved in biofilm development within emitters. The calculation of this parameter is detailed as follows:

𝑂𝑆𝐼 = 1 2 [1 - |∫ 𝜏 𝑤 𝑑𝑡 𝑇 0 | ∫ |𝜏 𝑤 |𝑑𝑡 𝑇 0 ] (12)
Where

| 1 𝑇 ∫ 𝜏 𝑤 𝑑𝑡 𝑇 0
| is the average intensity values of the wall shear stress vector, and

1 𝑇 ∫ |𝜏 𝑤 |𝑑𝑡 𝑇 0
is the time-averaged of the wall shear stress vector. The range of values for the OSI is 0 to 0.5, where a value of 0 indicates the presence of unidirectional shear stress. A high OSI value indicates a significant degree of oscillatory flow behavior.

Results and discussion

Discharge of the emitter

The study aimed to evaluate the precision of different turbulence models in predicting the flow characteristics of an industrial emitter. To achieve this, the discharge-pressure loss curves (q = f(∆P))

for each turbulence method were compared with the experimental data (Fig. 3). Karmeli's power law was used to fit the obtained curves, but they showed some deviation for SKE and RSM models compared to experimental data. In contrast, the LES model coincided with the experimental results,

showing no significant difference (p-value <0.001). The relative error between the simulated and experimental flow rates ranged between 7.2 -9.1% and 25.5 -29.6% for SKE and RSM models, respectively. The differences in results obtained using various turbulence models could be attributed to the differences in the turbulent energy spectrum resolved by each model. These findings are consistent with previous studies that have also reported the performance of the LES model in simulating flow characteristics [START_REF] Li | Analysis of Tracing Ability of Different Sized Particles in Drip Irrigation Emitters with Computational Fluid Dynamics[END_REF][START_REF] Wu | Simulation of the flow characteristics of a drip irrigation emitter with large eddy methods[END_REF]Zhangzhong et al., 2015). For example, [START_REF] Li | Analysis of Tracing Ability of Different Sized Particles in Drip Irrigation Emitters with Computational Fluid Dynamics[END_REF] noted that compared to experimental results, the LES model showed better discharge prediction accuracy than the SKE model with deviations of 8.8% and 10.6% respectively, when tracing particles in the flow path of a drip irrigation emitter (3.19 l h -1 ). In addition, [START_REF] Wu | Simulation of the flow characteristics of a drip irrigation emitter with large eddy methods[END_REF] reported relative errors of 10.3% and 4.7% for SKE and LES, respectively, when simulating the flow characteristic in the flow path of an NPC cylindrical emitter. They also noticed that the LES model demanded more computational resources. It is important to note that accuracy is not the only factor to consider when selecting a turbulence model for drip irrigation system design and optimization. Computational resources and time requirements are also significant considerations. Indeed, in this study, the LES model provided better accuracy than the RANS models. However, it required more computational resources, taking 140 hours longer to solve the numerical computation than RANS models (using a powerful machine: Intel(R)

Xeon(R) Gold 6140 CPU with 36 cores and 192 GB RAM), which limits its practical application for largescale system optimization. Whilst the SKE turbulence model was more suitable for calculating the mean flow field in the emitter's flow path, the LES model may not be ideal for this purpose. However, for obtaining precise information on the flow characteristics inside the emitters, such as the unsteady behavior of the flow, oscillatory shear index, and clogging process understanding, LES could be a preferred choice.

Flow characterisation

Flow characterization: mean velocity and streamline

The mean velocity modulus fields along the labyrinth channel (Figs. 4 and S1 in the supplementary material) and the streamline (Fig. 5) details obtained at the middle of the channel (z=0. The streamlines (Fig. 5) show that the form, size, and centre of the vortex zone slightly vary from one model to another. For the RSM model, the vortex zone centre is predicted closer to the wall with a weak velocity value in comparison to SKE and LES where the vortex zone centre is in the middle part of the baffle. This difference between SKE and RSM is in agreement with the results obtained by (Al-Muhammad, 2016). In the same way, the streamlines obtained from LES revealed that the vortices present in the return zone exhibit larger sizes (more than double the size) compared to the ones obtained through RANS simulations. This observation can be attributed to the fact that LES captures unsteady flow phenomena and provides higher resolution of the turbulent eddies, whereas RANS models are based on time-averaged equations assuming that the turbulence is steady. In the literature, Figure S2 in the supplementary material illustrates the turbulent eddy dissipation obtained by three different models. It was observed that all models show an increase in energy dissipation in the mainstream flow after the second baffle. This can be attributed to the fact that the high mean flow velocity in the mainstream flow results in insufficient development of turbulence at first baffles. This phenomenon is consistent with previous studies on turbulent flow in confined geometries [START_REF] Clercx | Dissipation of coherent structures in confined twodimensional turbulence[END_REF][START_REF] Wei | Integral properties of turbulent-kinetic-energy production and dissipation in turbulent wall-bounded flows[END_REF]. When comparing the results of the three models, the LES model resolves the large-scale eddies and was capable of capturing the development of turbulence in the flow, resulting in less energy dissipation, especially in the first baffles at the inlet and the first baffles located after the return zone, but more energy was dissipated in the return zone where energetic eddies are difficult to dissipate. In contrast, the RSM model relies on a turbulence closure model that may not accurately capture the effects of large-scale eddies, resulting in an over-prediction of energy dissipation in certain regions of the flow. Furthermore, the observed negligible dissipation in the swirl region for the RSM model could be connected to the discharge curve of the emitter used in the simulation. These results are consistent with the previous study [START_REF] Al-Muhammad | Flow in millimetric-channel: Numerical and experimental study[END_REF].

In order to track the development of flow inside the labyrinth channel, the velocity profiles inside the first four baffles at the inlet and after the return zones (which had the same angle) were plotted in Fig. 6 (red lines in Fig. 1). It is worth noting that the results of LES were reported as the average of the velocity at all the time steps. The obtained results revealed that the velocity profile remained almost unchanged from the third baffle (l2) for both SKE and RSM models, with deviations of 3.2% and 3.4% between the velocity at baffles 3 and 6 (l2 and l4) for SKE and RSM models, respectively. However, the velocity profile continued to develop until the fourth baffle in the LES model, where the deviation reached 1% between baffles 4 and 6, while it was 13.4% between baffles 3 and 4. Interestingly, after the return zone, the flow was fully developed from the fourth baffle for both SKE and RSM models, with deviations of less than 2%. However, the LES model showed a deviation of 7% between the velocity at baffles 23 and 25 (l7 and l8). These results show that the LES model was more effective in capturing the turbulent behavior of the flow inside the labyrinth, which was not adequately captured by the other two models.

Therefore, in order to investigate with more accuracy this comparison and analyze the differences between the models, the mean velocity profiles along lines 4 and 8 (where the flow is fully developed at the inlet and return zones) were plotted in Fig. 7. The comparison revealed notable discrepancies between the models. For instance, compared to RSM, the SKE model predicts a higher velocity in the mainstream flow, with a maximum deviation of 11% at the inlet zone (increasing to 30% after the return zone). However, a significant difference can be observed in the vortex zone, where the flow velocity is lower and the deviation is more pronounced, reaching 62% at the inlet zone (68% after the return zone). On the other hand, the LES showed a deviation of 14% (13%) in the mainstream flow but a larger deviation of 32% (34%) in the vortex zone when compared to the SKE model. In addition, the LES model accurately predicted the highest velocity values in the mainstream zone, with maximum velocity deviations of 15.4 and 30.8% when compared to the SKE and RSM models, respectively. This finding is consistent with the existing literature. For instance, [START_REF] Feng | Effect of optimization forms of flow path on emitter hydraulic and anti-clogging performance in drip irrigation system[END_REF] recently conducted a numerical evaluation of the hydraulic performance of a flat irrigation emitter using LES, SKE, and RNG models. Their study demonstrated that the LES model yielded higher velocity predictions in the mainstream region, with maximum velocity errors of 7.12, 1.95, and 1.30%, respectively, compared to experimental characterization using an improved digital particle image velocimetry (DPIV) system. When interpreting the obtained results regarding emitter clogging, it is crucial to take into account the low-velocity regions close to the channel walls, where the flow slows down. These zones can have a considerable impact on the occurrence of clogging. As water flows through these stagnant zones, sediment, microorganisms, and other particles can settle out of the flow and accumulate, eventually leading to biofilm developments and blockages in the emitter. In the literature clogging of irrigation emitters due to the presence of low-velocity zones near the walls has been extensively studied [START_REF] Feng | Effect of optimization forms of flow path on emitter hydraulic and anti-clogging performance in drip irrigation system[END_REF][START_REF] Li | CFD and digital particle tracking to assess flow characteristics in the labyrinth flow path of a drip irrigation emitter[END_REF]Li et al., , 2019)). The researchers concluded that low-velocity zones near the walls of the emitter were more prone to sediment accumulation and clogging than zones with higher flow velocities. They also observed that increasing the flow velocity as well as increasing the labyrinth channel dimensions could reduce sediment accumulation and improve the performance of the irrigation system. The high turbulence present in the labyrinths creates vortices that serve as essential mechanisms for controlling fluid discharge and dissipating energy [START_REF] Wang | Influence of Emitter Structure on Its Hydraulic Performance Based on the Vortex[END_REF][START_REF] Wei | Study on hydraulic performance of drip emitters by computational fluid dynamics[END_REF].

These findings support the observations made in this study that, in the first and nearby return zone baffles, the mainstream flow is disrupted by stagnation zones (as shown in Fig. S1 in the supplementary material). These zones may promote biofilm development, as reported by [START_REF] Lequette | Hydrodynamic effect on biofouling of milli-labyrinth channel and bacterial communities in drip irrigation systems fed with reclaimed wastewater[END_REF].

Indeed, based on optical coherence tomography (OCT) analysis, the researchers have reported that the inlet (mostly in the first baffle) and return zones of this drip emitter were the most sensible zones for clogging as the biofouling volume tended to be higher (3.32 mm 3 ) at the inlet zones after four months, and it gradually decreased towards the return zone (2 mm 3 ). However, the presence of the low-velocity zone alone cannot fully account for the growth of biofilm or clogging of the emitter, as demonstrated in [START_REF] Lequette | Hydrodynamic effect on biofouling of milli-labyrinth channel and bacterial communities in drip irrigation systems fed with reclaimed wastewater[END_REF], where the bacterial attachment and the thickness of biofilm at the walls varied across different baffles sections of the emitters.

Velocity fluctuations

The instantaneous velocity profile obtained from the LES model exhibited a deviation from the timeaveraged velocity profiles as shown in Fig. 7. This deviation is indicative of an oscillating flow, particularly after the return zone, where eddy structures are overlapping. Results show that the turbulence intensity after the return zone was around 6 times higher than that before it which indicates that there is a higher degree of flow unsteadiness and fluctuation in the flow after the return zone. Additionally, the plot illustrates the Gaussian distribution with identical standard deviation and mean values as those of the data. In the first part of the labyrinth, the dataset has a narrow range of values, a slight positive skewness, and also a positive kurtosis (Fig. 8a). These properties suggest that the velocity fluctuations, in this labyrinth part, are more tightly clustered around the mean, potentially indicating a more stable flow with less variability. In contrast, the velocity fluctuation after the return zone has a wider range of values, a slight negative skewness, and a lower kurtosis than the normal distribution indicating a flatter distribution (Fig. 8b). These suggest that the velocity fluctuations in this region are more variable and potentially influenced by more complex and dynamic flow phenomena.

Such oscillatory behavior is thus observed, probably due to the instantaneous transfer of energy from smaller to larger scales [START_REF] Yasuda | A vortex interaction mechanism for generating energy and enstrophy fluctuations in high-symmetric turbulence[END_REF]. indicating a significant change in the flow characteristics. Indeed, the velocity data after the return zone exhibited a high peak corresponding to the value of the dominant frequency (747 Hz), which was 7 times higher than the frequency observed before the return zone indicating a significant change in the flow characteristics at this point. The high frequency of this peak suggests that these structures are small and short-lived, possibly due to the high turbulence intensity. In contrast, the first dominant frequency (101 Hz) observed before the return zone is likely associated with the flow behavior in the preceding flow path. This frequency is expected to represent a larger-scale motion. The ESD data was further analysed by applying the Kolmogorov f -5/3 law, which is commonly used to describe the inertial subrange of fully developed turbulence. This law indicates that the turbulent kinetic energy is dissipated through the action of small-scale turbulent eddies, giving rise to a characteristic length scale, known as the Kolmogorov length scale [START_REF] Wei | Integral properties of turbulent-kinetic-energy production and dissipation in turbulent wall-bounded flows[END_REF]. The analysis revealed that the Kolmogorov length scale was in the range of 0.2 -5.3mm and 0.1 -0.4mm before and after the return zone, respectively.

The ESD data before the return zone showed good agreement with the predictions of Kolmogorov's law, suggesting that the small-scale turbulence was almost isotropic and universal. However, after the return zone, the ESD data exhibited a deviation from the predictions of Kolmogorov's law, indicating that the dissipation became anisotropic [START_REF] Kolmogorov | Dissipation of Energy in the Locally Isotropic Turbulence[END_REF][START_REF] Wei | Integral properties of turbulent-kinetic-energy production and dissipation in turbulent wall-bounded flows[END_REF]. Table 3 reported the turbulence intensity (Ti) and dominant frequency at three specific positions for each of the following locations (presented by the x character in Fig. 1): the first baffle near the inlet, the return zone, and the outlet just before the last baffle. The chosen positions correspond to y coordinates of d5/4, d5/2, and 3d5/4. The results suggest that the inlet zone has the lowest turbulence intensity and a similar dominant frequency in all baffle sections. This is expected as the inlet zone is typically the most laminar region in a fluid flow system, and the flow develops and becomes more turbulent as it progresses through the system. However, in the remaining sections of the geometries, the fluctuation in mainstream flow was characterised by the highest frequency, suggesting that the flow is becoming more turbulent and that the energy in the flow is concentrated at a particular frequency. Notably, the dominant frequency observed in this region was found to be close to the fundamental frequency calculated for the system (625 Hz for a mean velocity of 1 m s -1 ).

The velocity fluctuations can cause variations in the wall shear stress exerted on the surface of the drip irrigation emitter, which can affect the deposition of biofilm and particles, therefore, a specific analysis such as the oscillatory shear index is required. As our study's findings indicate that the LES offers crucial information that other models, such as SKE and RSM, cannot capture, this paper will solely focus on the results obtained from LES in the following sections.

Shear stress at the walls and biofilm behavior

In this study, the aim was to investigate the impact of wall shear stress on biofilm growth and deposition in the studied drip irrigation emitter. To achieve this, the wall shear stress distribution was plotted. The results of the LES simulation showed that the distribution of wall shear stress was nonuniform, with values ranging from very low levels around 0.003 Pa in zones with low velocity to significantly higher levels in the mainstream flow. To facilitate a better understanding of the impact of wall shear stress on biofilm growth and deposition, the colour code was rescaled and the results are presented in Fig. 10. Our analysis focused on comparing regions of the emitter where the wall shear stress was below a general threshold with those that were more sensitive to biofilm deposition reported in [START_REF] Lequette | Hydrodynamic effect on biofouling of milli-labyrinth channel and bacterial communities in drip irrigation systems fed with reclaimed wastewater[END_REF]. As seen in Fig. 10 the inlet and return zone (including the first 3 baffles) are characterised by very low shear stress (<3 Pa) reducing their self-cleaning ability. This led to the attachment of solid particles and microorganisms to the surfaces which may explain the biofilm accumulation in these zones. Whilst most of the other regions have important wall shear stress values promoting selfcleaning and causing a dynamic growing-detaching process of the substances. As a result, biofilms formed in these zones are more likely to detach partially resulting in a lower thickness [START_REF] Lequette | Hydrodynamic effect on biofouling of milli-labyrinth channel and bacterial communities in drip irrigation systems fed with reclaimed wastewater[END_REF]. This suggests that maintaining an appropriate value of wall shear stress to enhance the selfcleaning ability of the emitter flow path is the key to improving its anti-clogging ability.

Apart from flow conditions, the availability of nutrients is also a significant factor in determining the thickness and structure of biofilms. According to [START_REF] Stoodley | Oscillation characteristics of biofilm streamers in turbulent flowing water as related to drag and pressure drop[END_REF][START_REF] Stoodley | Oscillation characteristics of biofilm streamers in turbulent flowing water as related to drag and pressure drop[END_REF], the presence and thickness of biofilm were found to be linked to nutrient availability under high shear rates. A high nutrient concentration can promote biofilm growth, but it may result in weak adhesive strength of the biofilm, as observed in the studies by [START_REF] Liu | Microfluidic study of effects of flow velocity and nutrient concentration on biofilm accumulation and adhesive strength in the flowing and no-flowing microchannels[END_REF] and Zhou et al., (2021). On the other hand, in zones where nutrients are lacking, the attached microorganisms have limited access to nutrients, and therefore, the production of extracellular polymeric substances (EPS) becomes insufficient. As a result, the clogging substances spread over a larger surface area with a smaller average thickness and roughness, forming a porous structure [START_REF] Li | Effects of average velocity on the growth and surface topography of biofilms attached to the reclaimed wastewater drip irrigation system laterals[END_REF]Zhou et al., 2019). This allows microorganisms to obtain nutrients and secreted EPS, thereby maintaining a stable structure.

In this case, the effect of nutrient supplementation is stronger than the self-cleaning effect, meaning that the reduction in biofilm thickness occurs when the detachment rate induced by the shear stress exceeds the growth rate. By combining advanced techniques such as biofilm thickness measurement, EPS characterisation, and CFD simulation methods, our understanding of biofouling in drippers can be enhanced. Such a multidisciplinary approach would enable us to delve deeper into the mechanisms involved in biofilm formation and growth, as well as the impact of the EPS on the process. Moreover, this could help develop better strategies for preventing and controlling biofouling in drip irrigation systems.

Oscillatory Shear Index (OSI)

Along the wall shear stress (𝜏 𝑤 ), the OSI is another parameter that can affect biofilm growth in drip irrigation emitters. The OSI represents the ratio of the time-averaged shear stress to the root mean square of the velocity gradient in a fluid flow. In other words, it describes the relative importance of the time-averaged shear stress and the turbulent fluctuations in the flow. This parameter is extensively applied in cardiovascular flow studies to quantify the spatial characteristics of oscillatory flow behavior [START_REF] Huo | Flow patterns in three-dimensional porcine epicardial coronary arterial tree[END_REF]. This index is particularly useful for analysing the impact of flow disturbances on the endothelial cells that line blood vessels. It is calculated by analyzing directional changes in shear stress that occur during the cardiac cycle. High degrees of flow disturbance, such as those near a bifurcation or stenosis, result in higher OSI values, whereas regions with a more laminar flow have lower OSI values. The OSI provides valuable insight into the hemodynamic forces acting on the endothelium, which can affect endothelial function and contribute to the development of cardiovascular diseases such as atherosclerosis.

In the context of biofilm cleaning processes, recent research has delved into the impact of shear stress in scenarios involving intermittent flows (Li et al., 2019). However, a notable gap exists in the lack of standardised index parameters to precisely quantify this phenomenon. Specifically, there remains a dearth of information concerning the influence of the OSI on biofilm control, particularly within the context of drip irrigation systems. Thus, this article aims to address this gap by investigating the potential of OSI as a tool for comprehending biofilm growth in emitters. Our study endeavors to provide plausible explanations, supported by simulations outlined in this article, alongside experimental measurements of biofilm growth in similar [START_REF] Lequette | Effects of the chlorination and pressure flushing of drippers fed by reclaimed wastewater on biofouling[END_REF] or diverse types of drippers documented in the literature [START_REF] Stoodley | Oscillation characteristics of biofilm streamers in turbulent flowing water as related to drag and pressure drop[END_REF]. Moreover, we aim to establish connections between these hypotheses and pertinent examples from the food industry to bolster our assertions [START_REF] Dallagi | Removal of biocontamination in the food industry using physical methods; an overview[END_REF]. This study marks the initial effort in exploring the potential correlation between OSI and biofilm control specifically within the domain of drip irrigation systems.

Figure 11 shows that in the first half of the emitter, the flow is quasi-steady. However, in the return zone, flow disturbances are observed, which affect the behavior of the wall shear stress vector. This results in a significant value of the OSI occurring not only in this zone but also in the four preceding and all following zones. These disturbances can lead to increased turbulence and recirculation, resulting in a higher degree of oscillatory flow behavior as quantified by the OSI. These observations confirm the results discussed earlier regarding velocity fluctuations. The dominant frequency after the return zone was found to be 7.5 times higher than that before this zone. It is crucial to note that in zones where the OSI is high, both the frequency and amplitude of oscillation are concurrently intensified, resulting in a pronounced enhancement of mass transfer, which promotes the mixing of nutrients and enhances the diffusion of nutrients within the biofilm. This can improve nutrient availability and distribution within the biofilm biomass. The impact of this factor can extend beyond the biofilm biomass, as it significantly alters the physical structure of biofilms as demonstrated by [START_REF] Tsagkari | The role of shear dynamics in biofilm formation[END_REF]. This concept is not in contradiction with our results since this heightened mass transfer is accompanied by an equivalent intensification in drag force, consequently improving the detachment process [START_REF] Pechaud | Influence of shear stress, organic loading rate and hydraulic retention time on the biofilm structure and on the competition between different biological aggregate morphotypes[END_REF]. Several researchers have demonstrated the role of the fluctuation of the wall shear stress in improving cleaning efficiency using intermittent flow, bubbles, or foam flow in the case of cleaning of standard pipes in food industries [START_REF] Dallagi | Wet foam flow: A suitable method for improving surface hygiene in the food industry[END_REF] as well as the use of intermittent fluctuated water pressure for the case of irrigation emitters (Li et al., 2019).

Therefore, the OSI could be a promising tool for characterising the cleaning processes in drip irrigation systems. By measuring the intensity of the shear forces generated by fluid flow, OSI provides a quantitative metric for assessing the cleaning efficiency of irrigation systems.

Conclusion

This study aims to analyze and evaluate the accuracy and computational efficiency of different turbulence models (LES, RSM, and SKE), in simulating the flow characteristics of an industrial emitter in drip irrigation systems. Our results indicate that the LES model offers the most accurate simulation results for understanding the impact of hydrodynamics on emitter clogging by solid particles and biofilm. However, it demands significant computational resources. However, the SKE model provides a good balance between accuracy and computational efficiency and is recommended for simulating the mean flow characteristics, such as the emitter discharge.

The results conclude that the accumulation of biofilm in a flowing microchannel is highly dependent on the flow velocity, which affects two crucial factors: mass transfer and wall shear stress:

-Low-velocity zones near the walls of the emitter are more prone to sediment accumulation and biofilm growth whilst increasing the flow velocity and labyrinth channel dimensions can reduce sediment accumulation and improve the emitter's anti-clogging performance.

-Maintaining an appropriate level of wall shear stress is crucial for enhancing the self-cleaning ability of the emitter by detaching the biofilm and improving its anti-clogging ability.

-High OSI values can improve nutrient availability and distribution within the biofilm by enhancing mass transfer and promoting the mixing and diffusion of nutrients. Nerveless the fluctuation of the wall shear stress could improve the self-cleaning ability. The OSI parameter can also be a valuable tool for evaluating the effectiveness of disinfection protocols in removing biofilms and other contaminants from the system.

Future research in this area should continue to explore the potential of OSI as a tool for characterising cleaning processes in drip irrigation systems and should investigate the impact of nutrient distribution on biofilm growth and control.

Declaration of Competing Interest

There is no conflict of interest for this work.

Fig. 1 :

 1 Fig. 1: Structure of the drip irrigation emitter. (a) 3D geometry of the industrial emitter. (b) 2D

Fig 3 .

 3 Fig 3. Discharge-pressure loss curves for the different turbulence models (CFD) and experimental data.

  4 mm) were plotted to understand the hydrodynamic in this complex geometry and to compare the turbulent models used in this study. Note that the reported results of LES consisted of the velocity average across all time steps. The results show that the flow inside this kind of emitter is characterised by two zones, the vorticity zones in channel corners and the mainstream flow characterized by a high velocity. This description of the velocity regions is in agreement with the description cited by several experimental as well as numerical works such as[START_REF] Ait-Mouheb | Impact of hydrodynamics on clay particle deposition and biofilm development in a labyrinth-channel dripper[END_REF][START_REF] Al-Muhammad | Experimental and numerical characterization of the vortex zones along a labyrinth milli-channel used in drip irrigation[END_REF][START_REF] Lequette | Hydrodynamic effect on biofouling of milli-labyrinth channel and bacterial communities in drip irrigation systems fed with reclaimed wastewater[END_REF].

Fig. 4 .

 4 Fig. 4. Mean velocity fields obtained at z= 0.4 mm using SKE (a), RSM (b), and LES (c) models. The

Fig. 5 .

 5 Fig. 5. Velocity streamlines at the inlet zone of the labyrinth (a, b, and c for the SKE, RSM, and LES

Feng

  et al. (2018) reported that the LES model demonstrated exceptional accuracy in predicting the position of the vortex distribution in an emitter flow path as compared to the RANS models, as evidenced by experimental results. The researchers explained that the vortex could be generated near the wall due to two reasons. Firstly, the partition motion of viscous fluid changes the flow regime when fluids with different velocities and flow directions come in contact with each other. This creates relative motion between the flows, causing the fluid in the non-mainstream zone to lag behind and change direction, resulting in a vortex. Secondly, the wall can block the motion of the fluid, causing inelastic collisions between fluid particles and the wall, leading to a reduction in kinetic energy and a change in direction.

Fig. 6 .

 6 Fig. 6. Velocity profiles of different turbulence models at key sections of the emitter: SKE (a), RSM (b),

Fig. 7 .

 7 Fig. 7. (a) comparison of the velocity profile at baffle number 6 (l4), (b) comparison of the velocity

Fig. 8 .

 8 Fig. 8. (a) velocity fluctuation before (red colour) and after (black colour) the return zone. Turbulence

Fig. 9 .

 9 Fig. 9. Energy spectral density (ESD) with respect to the frequency domain f before (a) and after the

Fig. 10 .

 10 Fig. 10. Wall shear stress distribution along the flow path. (a) Three-dimensional representation of the

Fig. 11 .

 11 Fig. 11. Oscillatory shear index (OSI) distribution along the flow path. (a) Three-dimensional

Fig. 2 .

 2 Fig. 2. Schematic representation of a pilot rig illustrates emitter discharge was evaluated by measuring

Figure 2

 2 Figure 2 Click here to access/download;Figure;Figure 2.docx

Figure 3

 3 Figure 3 Click here to access/download;Figure;Figure 3.docx

Fig. 4 .

 4 Fig. 4. Mean velocity fields obtained at z= 0.4 mm using SKE (a), RSM (b), and LES (c) models. The

Figure 4

 4 Figure 4 Click here to access/download;Figure;Figure 4.docx

Fig. 5 .=

 5 Fig. 5. Velocity streamlines at the inlet zone of the labyrinth (a, b, and c for the SKE, RSM, and LES

Fig. 6 .

 6 Fig. 6. Velocity profiles of different turbulence models at key sections of the emitter: SKE (a), RSM (b),

Figure 6

 6 Figure 6 Click here to access/download;Figure;Figure 6.docx

Fig. 7 .

 7 Fig. 7. (a) comparison of the velocity profile at baffle number 6 (l4), (b) comparison of the velocity

Figure 7

 7 Figure 7 Click here to access/download;Figure;Figure 7.docx

Fig. 8 .Figure 8

 88 Fig. 8. (a) velocity fluctuation before (red colour) and after (black colour) the return zone. Turbulence

Fig. 9 .

 9 Fig. 9. Energy spectral density (ESD) with respect to the frequency domain f before (a) and after the

Fig. 10 .

 10 Fig. 10. Wall shear stress distribution along the flow path. (a) Three-dimensional representation of the

Figure 10

 10 Figure 10 Click here to access/download;Figure;Figure 10.docx

Fig. 11 .

 11 Fig. 11. Oscillatory shear index (OSI) distribution along the flow path. (a) Three-dimensional

Figure 11

 11 Figure 11 Click here to access/download;Figure;Figure 11.docx

Table 1 :

 1 Geometry parameters of the studied emitter

	Parameters	Nomenclature	Value	Unit
	Inlet diameter	dinlet	1	mm
	Outlet diameter	doutlet	1	mm
	Depth	e	0.8	mm
	Hydraulic diameter	Dh	1.02	mm
	Unit length	d6	1.6	mm
		d1	8.2	mm
		d2	34.6	mm
		d3	4.3	mm
		d4	2.03	mm
		d5	0.95	mm
		d7	4.5	mm
	Angles	α1	28	°
		α2	37	°

Table 2 .

 2 Proportionality Coefficient (kd) and Emitter Discharge Exponent (x) Values for Discharge-Pressure Loss Curves in Laminar (x = 1) and Fully Turbulent (x = 0.5) Flows, based on Karmali Equation (q = kd ∆P x ).

		Experimental	SKE	RSM	LES
	kd	1.062	1.1514	1.3731	1.0714
	x	0.4777	0.4924	0.4501	0.4544

On the other hand, regardless of the turbulence model used, the emitter discharge exponent was close to 0.5, indicating that the flow was well-turbulent (Table

2

). The results demonstrate that all turbulence models used provided acceptable results, with the deviation not exceeding 5.78% in the RSM model. This is in agreement with the previous research

[START_REF] Al-Muhammad | Flow in millimetric-channel: Numerical and experimental study[END_REF]

, where the emitter discharge exponent obtained from SKE and RSM simulation were in good agreement with the experimental result, but the SKE model exhibited less deviation (4%) compared to the RSM model (22.9%).

Table 3 .

 3 The turbulent intensity and dominant frequency at the inlet, return, and outlet zones.

	Position		Inlet		Return		Outlet
		Intensity	Frequency	Intensity	Frequency	Intensity	Frequency
	Y=d5/4	0.2%	102 Hz	9.4%	714 Hz	10.6%	702 Hz
	Y=d5/2	1.1%	102 Hz	26.6%	95 Hz	26.1%	171 Hz
	Y=3d5/4	0.3%	102 Hz	12.9%	169 Hz	11.1%	179 Hz

Funding

This work was supported by A*MIDEX -Mechanical and Engineering Institute, Aix-Marseille University.

Abbreviations

with a labyrinth channel. Journal of Hydrodynamics,Ser. B,28(4),[610][611][612][613][614][615][616] https://doi.org/10.1016/S1001-6058(16)60665-0 Zhangzhong, L., Yang, P., Ren, S., Liu, Y., & Li, Y. (2015). Flow Characteristics and Pressure-Compensating Mechanism of Non-Pressure-Compensating Drip Irrigation Emitters. Irrigation andDrainage, 64(5), 637-646

Declaration of interests

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐ The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: