Numerical Modeling of Mass Transfer in the Interaction between River Biofilm and a Turbulent Boundary Layer

Falilou Coundoul^1*#, Evrad M. D. Ngom^1*, Frédéric Y. Moulin^2*

¹Institut Polytechnique de Saint-Louis, Université Gaston Berger, Saint-Louis, Senegal.
²Institut de Mécanique des Fluides de Toulouse, UMR 5502-Université de Toulouse, Occitanie, France.

DOI: 10.4236/ojfd.2024.141001                

Abstract

In this article dedicated to the modeling of vertical mass transfers between the biofilm and the bulk flow, we have, in the first instance, presented the methodology used, followed by the presentation of various results obtained through analyses conducted on velocity fields, different fluxes, and overall transfer coefficients. Due to numerical constraints (resolution of relevant spatial scales), we have restricted the analysis to low Schmidt numbers (S_c=0.1, S_c=1, and S_c=10) and a single roughness Reynolds number (Re_*=150). The analysis of instantaneous concentration fields from various simulations revealed logarithmic concentration profiles above the canopy. In this zone, the concentration is relatively homogeneous for longer times. The analysis of results also showed that the contribution of molecular diffusion to the total flux depends on the Schmidt number. This contribution is negligible for Schmidt numbers S_c≥0.1, but nearly balances the turbulent flux for S_c=0.1. In the canopy, the local Sherwood number, given by the ratio of the total flux (within or above the canopy) to the molecular diffusion flux at the wall, also depends on the Schmidt number and varies significantly between the canopy and the region above. The exchange velocity, a purely hydrodynamic parameter, is independent of the Schmidt number and is on the order of 10% of in the present case. This study also reveals that nutrient absorption by organisms near the wall depends on the Schmidt number. Such absorption is facilitated by lower Schmidt numbers.

Keywords

Epilithic Biofilm, Passive Scalar Transport, Direct Numerical Simulation, Navier-Stokes

1. Introduction

To enhance the modeling of epilithic biofilm dynamics, it is necessary to properly parameterize exchanges between the canopy and the flow above it. The aim of this article is to analyze the evolution of a passive scalar diffusing from a hydraulically rough bottom subjected to turbulent boundary layer flow. The underlying idea is to study the spatio-temporal evolution of a mass boundary layer to quantify the vertical fluxes of chemical species encountered in turbulent river-type flows, particularly in the presence of an epilithic biofilm during colonization and initial growth. The study of interactions between river biofilms and turbulent boundary layers is an essential area of ecological research, significantly contributing to our understanding of aquatic ecosystems. These interactions are pivotal in nutrient cycling and maintaining the health of riverine environments. Recent studies, such as Font et al. [1] , have highlighted how environmental factors like low flow and heatwaves can alter ecosystem functioning in stream mesocosms. Koch et al. [2] explored the impact of chronic elevated nitrate concentrations on river biofilm structure and function, providing insights into nutrient loading effects. Additionally, Boulêtreau et al. [3] investigated how temperature influences denitrification processes in phototrophic river biofilms, adding to the knowledge of temperature-dependent biogeochemical cycles. Complementing these studies, Nepf [4] provided a comprehensive review of flow and transport in regions with aquatic vegetation, elucidating the hydrodynamic interactions within these habitats. In this context, our work aims to bridge knowledge gaps in the field by examining the influence of varying Schmidt numbers on mass transfers in turbulent boundary layers, offering new perspectives on these complex hydrodynamic-biological interactions. There has been relatively limited research on mass transfer in hydraulically rough turbulent boundary layer flows. The work Calmet [5] clearly quantified the importance of the Schmidt number in transfer processes, notably in turbulent flow over a smooth bottom with the presence of a free surface. Scale laws were obtained, particularly for the evolution

of the global transfer coefficient $K_L$ : $\frac{K_L}{u_{\mathrm{<mo>\; *\; </mo>}}}\propto S_c^{-2/3}$ and $\frac{K_L}{u_{\mathrm{<mo>\; *\; </mo>}}}\propto S_c^{-1/2}$ for a smooth

wall and a gas-liquid interface, respectively. For rough canopies, Nikora [6] describes the vertical exchange of nutrients or chemicals in the flow using the double mean. According to these last ones, these transfer or uptake processes are driven by the interaction between turbulence, viscous effects, and the complexity of the bottom geometry or the organisms living there. Nikora et al. [7] and Larned et al. [8] postulate that for mass transfer and absorption in a canopy colonized by periwinkles, there are at least three distinct absorption regimes defined by the interactions between the flow and the periphyton. The transfer regimes at the plant canopy-flow interface are classified based on the ratio of canopy

height ( $h_c$ ) to diffusive sublayer thickness ( ${\delta }_{Ds}=S_c^{-1/3}{\delta }_{\nu }$ ), where $S_c=\frac{\nu }{D}$ represents the Schmidt number, $\nu$ is the kinematic viscosity, and D is the mass diffusivity. Vertical transfers are essential for modeling biofilm dynamics, as they are responsible for nutrient and oxygen supply to the biofilm and higher trophic levels. It is necessary to quantify them properly for different situations, as shown by Coundoul et al. [9] . The work presented here involves modeling biofilm activity by the instantaneous consumption of chemical species with concentration denoted as “C” and, for different conditions, studying the hydrodynamic processes involved in matter transfer within this hydraulically rough turbulent flow. The manuscript is structured as follows: The “Material and Methods” section details our numerical modeling approach, emphasizing specific parameters and conditions. The “Results and Discussion” section presents and analyzes our findings, with a focus on the impact of varying Schmidt numbers on mass transfer and contextualizes our results within the existing literature. The “Conclusions” section summarizes our key contributions and suggests avenues for further study.

2. Material and Methods

2.1. General Solved Equations

In our case, which involves the unsteady flow of an incompressible and Newtonian fluid with gravity effects neglected, the equations for mass conservation, momentum, and passive scalar transport (mass concentration) can be written as follows, respectively, by Equations (1)-(3):

$\nabla \cdot u=0,$ (1)

$\frac{\partial u}{\partial t}+u\nabla \cdot u=-\frac{1}{\rho }\nabla p+\nu \nabla \left(\nabla \cdot u+{\left(\nabla \cdot u\right)}^{\text{T}}\right)\mathrm{,}$ (2)

$\frac{\partial C}{\partial t}+\nabla \cdot \left(Cu-D\nabla C\right)=\mathrm{0,}$ (3)

where $\nu \left(\nabla \cdot u+{\left(\nabla \cdot u\right)}^{\text{T}}\right)$ represents the viscous stress tensor, u is the fluid velocity vector, p is the pressure, C is the concentration of the passive scalar, $\rho$ is the fluid mass density, $\nu$ is the kinematic viscosity of the fluid, and D is the mass diffusivity of the passive scalar.

2.2. Presentation of the Digital Tool Used

The numerical code used throughout this study is a community code known as Jadim. This research tool was initially developed by Magnaudet [10] at the Institute of Fluid Mechanics in Toulouse (IMFT). Jadim is a finite volume method that solves the three-dimensional, unsteady, and incompressible Navier-Stokes equations (cf Tougma [11] , Badahmane [12] ). In recent years, new modules have been added to the core of the code, expanding its scope to study increasingly complex configurations. Jadim consists of a central core that enables the structured mesh resolution of the three-dimensional, unsteady, and incompressible Navier-Stokes equations. For over two decades, the research code Jadim has been used at IMFT in various fluid mechanics research activities, including biomechanics, environmental studies, energy, and processes. In addition to the immersed boundary method module (cf. Coundoul et al. [9] ) used in this study, the code also includes several other modules. These modules include one for simulating turbulent flows using the Large Eddy Simulation (LES) approach with a mixed dynamic subgrid-scale model. Another module deals with mass transfer (passive scalar) through the resolution of an additional advection-diffusion equation. There is also a Volume Of Fluid module without interface reconstruction developed by Calmet [5] and Legendre [13] . Efforts have been made in recent years to parallelize the Jadim code, leading to the parallelization of the central core and most of the modules. The code can now run on several tens to hundreds of processors.

2.2.1. Spatial and Temporal Discretization

The spatial discretization of Equations (1)-(3) is performed on a three-dimensional, structured, orthogonal, and staggered Cartesian grid in velocity-(pressure/passive scalar). The discretization follows a finite volume approach. Using this type of grid enhances the accuracy of flux calculations on the facets of integration volumes and provides good properties for the conservation of hydrodynamic quantities and the passive scalar. The three velocity components have different integration volumes. Pressure, concentration, and other scalars are integrated within the volume centered at p. A more detailed description of the spatial discretization in the code is provided in Calmet [5] . The numerical schemes used for time advancement of the solution to the system of Equations (1)-(3) are of the Runge-Kutta/Crank-Nicolson type. The advective terms and source terms are handled with a third-order Runge-Kutta scheme. The viscous and diffusive terms are treated using a semi-implicit Crank-Nicolson scheme. These scheme choices are particularly suitable for turbulent flows with very low diffusion. In Jadim, the numerical resolution of the conservation equations follows a projection method. It involves advancing the advective and diffusive terms in the first step and the pressure in the second step to satisfy the incompressibility condition (cf. Calmet [5] ).

2.2.2. The Immersed Boundary Method (I.B.M.) Used in Jadim

In the JADIM code is implemented an immersed boundary method (IBM) based on the method proposed by Yuki et al. [14] . This method employs a body force proportional to a solid volume fraction for coupling the solid and the fluid. It allows for the simulation of high-Reynolds-number flows in the presence of fixed or moving objects of arbitrary shape. Here, the immersed boundary is defined by the volume fraction of the solid denoted $\alpha$ which is equal to 1 if the cell is completely in the object, 0 otherwise and $0<\alpha <1$ if the cell is crossed by the immersed boundary. The expression of $\alpha$ depends on the shape of the object. Here, a hyperbolic-tangent function is used as a surface digitizer for computing the volume fraction. The detailed calculation $\alpha$ used in the present work can be found in Yuki et al. [14] . In the present case, we assume a no-permeability and no-slip condition at the motionless immersed boundaries (representing the cobbles), so the mass and momentum equations can be written in the new form

$\nabla \cdot u=0$ (4)

$\frac{\partial u}{\partial t}=-\nabla P+H+f$ (5)

$H=-u\nabla u+\nu \nabla \cdot \left[\nabla u+{\left(\nabla u\right)}^{\text{T}}\right]$ (6)

where $P=\frac{p}{\rho }$ , and $\nu =\mu /\rho$ . Note that (5)-(6) are similar to the Navier-Stokes

Equations (2) with the exception of the last (body force) term f. This new term allows us to modify the velocity field so as to take into account the presence of the immersed boundary. The simplified time-advancement scheme for $u$ between time $t=n\Delta t$ and $t=\left(n+1\right)\Delta t$ , n being the time increment and $\Delta t$ the time-step, is the following.

We first calculate a fluid predictor velocity $\tilde{u}$

$\tilde{u}=u^n+\Delta t\left(-\nabla P^n+H^n\right)$ (7)

We then compute the body force f as,

$f=-\alpha \tilde{u}/\Delta t$ (8)

and update the fluid predictor velocity by

${\tilde{u}}^{n+1}=\tilde{u}+\Delta tf$ (9)

Note that for cells inside the solid hemispheres $\alpha =1$ , thus according to (8)-(9) $f=-\tilde{u}/\Delta t$ and ${\tilde{u}}^{n+1}=0$ so as expected there is no motion inside the hemisphere, while for cells far from the hemispheres $f=0$ , the velocity field is unaffected by the immersed objects. Finally, the projection method is applied in order to get the new velocity field $u^{n+1}$ satisfying the incompressibility constraint, $u^{n+1}$ being obtained ${\tilde{u}}^{n+1}$ .

2.3. Results Analysis

As with the velocity field, according to Nikora [6] , it is possible to decompose the local instantaneous concentration field C into a doubly averaged part $\langle \bar{c}\rangle$ , a dispersive part $\tilde{c}=\bar{c}-\langle \bar{c}\rangle$ , and a fluctuating part $c^{\prime }=c-\bar{c}$ . The transport equation for c can then be expressed using Equation (10):

$\frac{\partial \langle \bar{c}\rangle }{\partial t}+\langle \bar{u}\rangle \frac{\partial \langle \bar{c}\rangle }{\partial x_j}=\frac{\partial }{\partial x_j}\langle \overline{D\frac{\partial c}{\partial x_j}}\rangle -\frac{\partial \langle \widetilde{u_j}\tilde{c}\rangle }{\partial x_j}-\frac{\partial \langle \overline{{u^{\prime }}_j{c}^{\prime }}\rangle }{\partial x_j}-\frac{1}{V_f}\stackrel{¯}{{\displaystyle {\iint }_{S_{int}}}D\frac{\partial c}{\partial x_j}{n}_j\text{d}S}$ (10)

where $J=-\frac{1}{V_f}\stackrel{¯}{{\displaystyle {\iint }_{S_{int}}}D\frac{\partial c}{\partial x_j}{n}_j\text{d}S}$ quantifies the absorption of nutrients by organisms near the bottom. D is the mass diffusivity of the passive scalar, $S_{int}$ represents the exchange surface with the rough wall, and $V_f$ represents the volume occupied by the fluid. The vertical transfer is driven by the vertical turbulent flux $\langle \overline{w^{\prime }{c}^{\prime }}\rangle$ and the dispersive flux $\langle \tilde{w}\tilde{c}\rangle$ .

The diffusive flux near a fixed wall $J_p$ is given by Equation (11):

$J_p=\frac{1}{V_0}\stackrel{¯}{{\displaystyle {\iint }_{S_{int}}}D\frac{\partial c}{\partial x_j}{n}_j\text{d}S}$ (11)

While the turbulent flux is given by Equation (12):

$J_{tur}=\langle \tilde{w}\tilde{c}\rangle +\langle \overline{w^{\prime }{c}^{\prime }}\rangle$ (12)

The viscous flux is given by Equation (13):

$J_v=\langle \overline{D\frac{\partial c}{\partial z}}\rangle$ (13)

And the total flux by Equation (14):

$J_t=\langle \tilde{w}\tilde{c}\rangle +\langle \overline{w^{\prime }{c}^{\prime }}\rangle +\langle \overline{D\frac{\partial c}{\partial z}}\rangle$ (14)

For mass transfer in a hydraulically smooth turbulent boundary layer, two quantities are generally defined to normalize turbulent quantities: a reference concentration $c_{\mathrm{<mo>\; *\; </mo>}}$ given by Equation (15):

$c_{\mathrm{<mo>\; *\; </mo>}}=\frac{D}{u_{\mathrm{<mo>\; *\; </mo>}}}\frac{\partial \langle \bar{c}\rangle }{\partial z}$ (15)

And $\Delta c=c_p-c_{ref}$ , the concentration difference between the wall and another reference concentration typically taken in the free stream. Note that these choices do not account for the dispersive contribution $\langle \widetilde{u_j}\tilde{c}\rangle$ in the calculation of the turbulent concentration flux.

2.4. Presentation of Different Simulations

To study the evolution of a passive scalar diffusing from a hydraulically rough bed subjected to turbulent flow, we conducted three numerical simulations for three different Schmidt numbers: $S_c=0.1$ , $S_c=1$ , and $S_c=10$ . In all cases, the initial instantaneous velocity field is obtained from the simulation of turbulent flow established under the conditions of Coundoul et al. [9] , with the geometry shown in Figure 1, and some parameters are listed in Table 1. This simulation was chosen for the quality of its mesh, allowing spatial scales to be resolved at approximately 3 times the Kolmogorov scale η in the vertical direction and 4.65η in the transverse and longitudinal directions. The concentration field c is initialized uniformly to zero: $c\left(\stackrel{\to }{x},t=0\right)=0$ .

In the study of biofilm evolution in rivers, we are interested in oxygen and carbon dioxide transfer. The mass diffusivities of these chemical compounds in water result in Schmidt numbers of 490 and 660, respectively. This leads to describing transfers with a diffusive sublayer ${\delta }_D$ approximately ten times smaller than the viscous sublayer ${\delta }_{\nu }~S_c^{-1/3}{\delta }_{\nu }$ . Therefore, it would be necessary to increase the spatial resolution by a factor of 10, making the computational times prohibitively long. For this reason, we restricted ourselves to cases where $S_c\le 10$ , for which the same spatial resolution could be adopted. In any case, the three

Table 1. Control parameters for different numerical simulations.

Figure 1. 3D view and top view of the computational domain. The flow is primarily directed along x.

simulations allow us to observe the influence of the Schmidt number $S_c$ in the development of the diffusive boundary layer, with a reasonable spatial resolution. In the subsequent work of Coundoul et al. [9] , these results will be used to extend the study of passive scalar transfer processes to cases of turbulent boundary layer flow over hydraulically rough beds. Periodic boundary conditions are imposed in the longitudinal and transverse directions, and a frictionless slip condition is imposed at the top of the domain, with a no-slip condition imposed on the obstacles and the bed. For the passive scalar, a constant concentration $c_p$ is imposed on the obstacles and the bed, with periodic boundary conditions in the longitudinal and transverse directions. A zero-flux condition is imposed at the free surface. In our study, the boundary condition of constant concentration on the bed and obstacles, leaving the wall flux free, can lead to significant variations in $c_{\mathrm{<mo>\; *\; </mo>}}$ . On the other hand, the imposed zero flux condition at the free surface leaves the concentration at the upper boundary-free, resulting in large variations in concentrations at the bed and free surface over time. As we will see later $c_{\mathrm{<mo>\; *\; </mo>}}$ and $\Delta c$ vary significantly between the different cases $S_c=0.1$ , $S_c=1$ , and $S_c=10$ , which poses challenges for normalizing statistical quantities. The choice of a concentration flux at the free surface has the consequence of not guaranteeing a steady state for the average concentration. However, it is possible to estimate a characteristic time at which all the scalar is consumed and to obtain the spatiotemporal evolution of the mass boundary layer. The non-stationarity of the average concentration fields makes it impossible to directly calculate the dispersive concentration fluctuations $\tilde{c}$ , which are necessary to calculate the dispersive flux following the double-averaging technique $\langle \tilde{w}\tilde{c}\rangle$ . Additionally, as Coundoul et al. [9] have shown, it is necessary to have a fully established flow and a large number of samples to obtain statistically converged results. To address this problem, we considered that each temporal sample $c\left(x\mathrm{,}y\mathrm{,}z\mathrm{,}t\right)$ of the concentration field is a fluctuation of the mean concentration field. This allows us to decompose the fluctuating field $c_f$ as:

$c_f\left(x\mathrm{,}y\mathrm{,}z\right)=c^{\prime }+\tilde{c}=c\left(x\mathrm{,}y\mathrm{,}z\mathrm{,}t\right)-\langle \bar{c}\rangle$ (16)

The dispersive flux is then given by:

$c_f\cdot w_f=\left(c-\langle c\rangle \right)\left(w-\langle w\rangle \right)\approx \left(c^{\prime }+\tilde{c}\right)\left(w^{\prime }+\tilde{w}\right)$ (17)

$\langle c_f\cdot w_f\rangle =\langle c^{\prime }{w}^{\prime }\rangle +\langle \tilde{c}\tilde{w}\rangle +\langle c^{\prime }\tilde{w}\rangle +\langle w^{\prime }\tilde{c}\rangle$ (18)

The last two terms are negligible compared to the rest of Equation (18), allowing us to directly deduce the turbulent flux at time t from the equation. The flux through the boundaries formed by the hemispheres is calculated by moving from a Cartesian coordinate system to a spherical coordinate system. The hemispherical shape of the obstacles that make up the bed makes the calculation of fluxes in Cartesian coordinates more complex.

3. Results and Discussions

3.1. Mass Transfer at the Surface of a Cylinder

To validate the mass transfer (or heat transfer) part of the immersed boundary module in Jadim, a 2D numerical simulation of flow around a cylinder with a diameter d in the presence of mass/heat transfer was performed. The results obtained with Jadim are compared with the numerical results of Kim and Choi [15] and the experimental results of Eckert and Soehngen [16] . The Reynolds number

$Re=\frac{u_{\infty }d}{\nu }$ , based on the cylinder diameter d and the upstream velocity $u_{\infty }$ , is

120. The computational domain size is $-20<x/d<50$ in the longitudinal direction and $-50<y/d<50$ in the vertical direction, with flow directed along the x direction and the cylinder’s center located at $x=y=0$ . The total number of cells is 530 × 334. The mesh is orthogonal and refined near the cylinder to accurately describe the thermal and dynamic boundary layers. The Prandtl number,

defined as $Pr=\frac{\nu }{k}$ (with k being the thermal diffusivity and $\nu$ the kinematic

viscosity of the fluid), is 0.7. Note that in the case of mass transfer, k is replaced by the mass diffusivity, and the Prandtl number becomes the Schmidt number. For the flow, a constant velocity is imposed at the domain’s inlet ( $u=u_{\infty }$ and $v=0$ ), an outlet boundary condition is imposed at the downstream edge in the longitudinal direction, and slip conditions are applied in the vertical direction on the boundaries located at the top and bottom of the domain. An adhesion condition is imposed on the cylinder’s surface. For temperature (or concentration or passive scalar), a Dirichlet condition is imposed at the domain’s inlet ( $T=T_{\infty }$ ), and on the cylinder’s surface ( $T=T_p$ ). A Neumann condition (zero flux) is applied on the upper and lower boundaries. In this case, the volume fraction $\alpha$ of the immersed cylinder is calculated according to Yuki et al. [14] :

$\alpha \left(x\mathrm{,}y\right)=\frac{1}{2}\left[1-\tanh \left(\frac{{\delta }_s}{\sigma \lambda \Delta }\right)\right]$ (19)

where: - ${\delta }_s$ is the distance from the center of the cylinder, - $\stackrel{\to }{n}=\left(n_x\mathrm{,}{n}_y\right)$ is the outward normal to the cylinder’s boundary, - $\lambda =\left|n_x\right|+\left|n_y\right|$ , - $\Delta =\sqrt{\Delta x^2+\Delta y^2}$ , and - $\sigma =0.05\left(1-{\lambda }^2\right)+0.3$ .

This choice allows a smooth transition of $\alpha$ from 1 to 0 over an approximately $\Delta$ thickness through the immersed boundary. At short times, the velocity field shows the separation of the fluid along the cylinder and a perfectly symmetrical recirculation zone. For longer times, the velocity field loses its symmetry in favor of the generation of Von Karman vortices in the wake of the cylinder (see Figure 2). The instantaneous temperature (concentration) fields shown in Figure 2 are in good agreement with the numerical simulations of Kim and Choi [15] . For comparison, we have included one of their results in Figure 3, corresponding to the same initial parameters ( $Re=120$ , $T=T_p$ on the cylinder). The distribution of the Nusselt number Nu along the heated cylinder for the simulation with Jadim, the numerical results of Kim and Choi [15] , and the experimental results of Eckert and Soehngen [16] are shown in Figure 4. The Nusselt number is defined by Equation (20):

$Nu=\frac{{\frac{\partial T}{\partial r}|}_{z=0}}{\frac{T_p-T_{\infty }}{d}}$ (20)

where $T_p$ is the temperature (concentration) at the cylinder’s wall, $T_{\infty }$ is the

temperature (concentration) at infinity, d is the cylinder’s diameter, and $\frac{\partial T}{\partial r}$ is

the temperature (concentration) gradient at the cylinder’s wall. The results are in good agreement with the numerical and experimental results. The instantaneous temperature fields effectively illustrate the organization of flow structures in the wake of the cylinder in both time and space.

Figure 2. Temperature field at different instants. For (a), $\left(t{u}_{\infty }\right)/d\approx 36$ and the time separating each successive frame is approximately $\left(\Delta t{u}_{\infty }\right)/d\approx 3.5$ .

Figure 3. Isocontours of temperature obtained by Kim and Choi [15] .

Figure 4. Nusselt number distribution along the heated cylinder. The results from Jadim are shown as triangles ( $\Delta$ ), those from Kim and Choi [15] as circles (o), and the experimental results of Eckert and Soehngen [16] as squares ( $\square$ ).

3.2. Instantaneous Concentration Fields

Figure 5 presents the instantaneous concentrations obtained for $S_c=0.1$ , $S_c=1$ , and $S_c=10$ at $t=50T$ . The influence of the Schmidt number $S_c$ on the spatial distribution of concentration in the flow is quite evident. In case of $S_c=1$ , the mean concentration field is relatively homogeneous, suggesting that mixing occurs primarily through diffusion. For $S_c=1$ and $S_c=10$ , the concentration distribution becomes more heterogeneous throughout the domain, unlike the case of $S_c=0.1$ . It can be concluded that mass transfer is driven by turbulent fluctuations, in line with the results obtained for a hydraulically smooth bottom by Calmet [5] . In the case of a hydraulically rough bottom, the viscous sublayer is submerged within the rough sublayer characterized by high shear rates and strong turbulence, as discussed in Florens [17] and Coceal et al. [18] . For Schmidt numbers greater than 1, mass transfer is no longer controlled by molecular diffusion but by turbulent structures in the rough sublayer. Therefore, it appears more relevant to compare the thickness of the diffusive sublayer with the thickness of the rough sublayer (rather than the viscous sublayer) to study the evolution of mean velocity profiles and mean concentration profiles in the case of a hydraulically rough bottom.

3.3. Mean Concentration Fields between Patterns

Figures 6-11 display the spatial mean concentration fields between patterns for $S_c=0.1$ , $S_c=1$ , and $S_c=10$ at different values of T. The first three figures provide a view at the $\left(x\mathrm{,}z\right)$ plane at $y/h=1$ , while the latter three offer a top-down view at the $\left(x\mathrm{,}y\right)$ plane at $z/h=0.5$ .

For all three Schmidt numbers, mixing occurs in the rough sublayer, and concentration maxima are located in the canopy region ( $z/h<1$ ). In the case of $S_c=0.1$ , concentration iso-contours follow the upper part of the hemisphere and elongate in the longitudinal direction, unlike the cases of $S_c=1$ and $S_c=10$ . Notably, there is significant development of a concentration boundary layer downstream of the hemisphere for all three cases. However, in the case of

Figure 5. Instantaneous concentration at $t=50T$ for $S_c=0.1$ , $S_c=1$ , and $S_c=10$ . The iso-surface represented is $c=0.97$ . The flow is in the x direction, and the rough Reynolds number is $R{e}_{\mathrm{<mo>\; *\; </mo>}}=156$ .

Figure 6. Contours of the spatial mean concentration between patterns for $S_c=0.1$ . View on the median plane of the basic pattern at $y/h=1$ .

Figure 7. Contours of the spatial mean concentration between patterns for $S_c=1$ . View on the median plane of the basic pattern at $y/h=1$ .

Figure 8. Contours of the spatial mean concentration between patterns for $S_c=10$ . View on the median plane of the basic pattern at $y/h=1$ .

Figure 9. Contours of the spatial mean concentration between patterns for $S_c=0.1$ . Top-down view at $z/h=0.5$ .

Figure 10. Contours of the spatial mean concentration between patterns for $S_c=1$ . Top-down view at $z/h=0.5$ .

$S_c=0.1$ , the concentration boundary layer contours are smoother and thicker compared to the other cases.

These observations suggest that mixing is more efficient in the case of $S_c=0.1$ than in the cases of $S_c=1$ and $S_c=10$ . For $S_c=10$ , the presence of concentration streaks, coexisting with lower concentration regions, is notable. These streaks are known as “concentration streaks” and have been highlighted by Calmet [5] .

3.4. Vertical Concentration Profiles

In Figure 12, spatially-averaged vertical concentration profiles at various characteristic times are presented for the three Schmidt numbers studied: $S_c=1$ , $S_c=0.1$ , and $S_c=10$ . As previously discussed, the mean concentration profiles do not converge towards a steady profile due to the imposition of a zero flux

Figure 11. Contours of the spatial mean concentration between patterns for $S_c=10$ . Top-down view at $z/h=0.5$ .

Figure 12. Temporal evolution of spatially-averaged vertical concentration profiles $\langle c\rangle$ at different times t for the three Schmidt numbers: $S_c=1$ , $S_c=0.1$ , and $S_c=10$ . Time t is normalized by the characteristic eddy turnover time T.

condition for concentration on the upper boundary. The imposition of $c=1$ at $z/h=0$ , zero flux at $z/h=1$ , and periodic boundary conditions in the lateral and longitudinal directions leads to an accumulation of material in the fluid column. This results in a progressive increase in the concentration field throughout the domain. It is noticeable in Figure 12 that, for all three Schmidt numbers, the concentration gradient is maximal at the height of the hemispheres ( $z/h=1$ ). This is clearly visible in Figure 14 and Figure 15, where the same quantities as in Figure 12 are presented, but normalized by $\Delta c$ and $c_{\mathrm{<mo>\; *\; </mo>}}$ , respectively, and in a semi-logarithmic scale. Furthermore, mixing in the water column evolves differently depending on the Schmidt number, as suggested by the observations in Figure 5. The most surprising aspect is that this evolution is not consistent as the Schmidt number increases. To illustrate this point, we plotted the temporal evolution of the mean instantaneous concentration in the canopy region ( $z/h<1$ ) and outside the canopy region ( $z/h>1$ ) in Figure 13(b)). Based on these results, mixing is favored when the Schmidt number $S_c\ll 1$ or $S_c\gg 1$ outside the canopy region. For $z/h<1$ , mixing is almost the same for $S_c=1$ and $S_c=10$ . In general, mixing is more pronounced for $S_c=0.1$ . As previously mentioned, there is a challenge in using $\Delta c$ and $c_{\mathrm{<mo>\; *\; </mo>}}$ to normalize the statistical quantities. The vertical concentration profiles demonstrate that the concentration is nearly stationary within the canopy, indicating that the values of $\Delta c$ and $c_{\mathrm{<mo>\; *\; </mo>}}$ vary significantly over time and with the Schmidt number. Despite our normalization choices, the concentration profiles do not overlap (see Figure 14 and Figure 15). Nevertheless, a logarithmic behavior is observed for the mean concentration profiles obtained for Schmidt numbers $S_c=1$ and $S_c=10$ . This can be explained by the fact that for Schmidt numbers greater than 1, transfer is governed by turbulent agitation, whereas for Schmidt numbers less than 1, purely diffusive effects dominate.

Figure 13. (a) Temporal evolution of the diffusive concentration flux near the wall $J_p$ (kg·m⁻²·s⁻¹) for the three Schmidt numbers: $S_c=0.1$ , $S_c=1$ , and $S_c=10$ . (b) Temporal evolution of instantaneous concentration mean in the canopy region ( $z/h<1$ ) and outside the canopy region ( $z/h>1$ ) for the three Schmidt numbers: $S_c=1$ , $S_c=0.1$ , and $S_c=10$ . Time t is normalized by the characteristic eddy turnover time T.

Figure 14. Temporal evolution of spatially-averaged vertical concentration profiles $\langle c\rangle$ at different times t for the three Schmidt numbers: $S_c=1$ , $S_c=0.1$ , and $S_c=10$ . Time t is normalized by the characteristic eddy turnover time T. The profiles are normalized by $\Delta c$ .

Figure 15. Temporal evolution of spatially-averaged vertical concentration profiles $\langle c\rangle$ at different times t for the three Schmidt numbers: $S_c=1$ , $S_c=0.1$ , and $S_c=10$ . Time t is normalized by the characteristic eddy turnover time T. The profiles are normalized by $c_{\mathrm{<mo>\; *\; </mo>}}$ .

3.5. Concentration Diffusive Flux

In Figure 13(a)), the temporal evolution of the diffusive concentration flux $J_p$ (kg·m⁻²·s⁻¹) within the canopy region ( $z/h<1$ ) is presented. It is evident that $J_p$ is higher at $t=0$ and gradually decreases over time to approach a constant flux. Moreover, it decreases as the Schmidt number increases. The fluxes for Schmidt numbers $S_c=1$ and $S_c=10$ appear to converge, unlike the case with $S_c=0.1$ . The initial peak in diffusive flux can be attributed to a significant initial concentration gradient, which decreases over time due to the imposition of a constant concentration on the domain’s bottom and hemispheres. The convergence to constant flux values is observed for cases with Schmidt numbers $S_c=1$ and $S_c=10$ is a result of diffusive transfer being limited by turbulent transfer. Therefore, the plateau is reached more quickly $S_c=10$ than for $S_c=1$ . In the case of $S_c=0.1$ , diffusive effects dominate over turbulent effects, leading to a continuous decrease in diffusive flux. It’s worth noting that for longer times, such as $t/T=50$ , the ratios between the diffusive fluxes obtained for $S_c=0.1$ , $S_c=1$ , and $S_c=10$ are approximately 4 and 2, respectively.

3.6. Turbulent Concentration Flux

For a turbulent boundary layer flow, it is expected that, far from the wall, the

viscous flux $J_v=\langle \overline{D\frac{\partial c}{\partial z}}\rangle$ is negligible compared to the turbulent flux $\langle \overline{w^{\prime }{c}^{\prime }}\rangle$

and the dispersive flux $\langle \tilde{w}\tilde{C}\rangle$ . In Figure 16, vertical profiles of the total flux $J_t$ (kg·m⁻²·s⁻¹), turbulent flux $J_{tur}$ (kg·m⁻²·s⁻¹), and viscous flux $J_v$ (kg·m⁻²·s⁻¹) are presented for the three Schmidt numbers $S_c=1$ , $S_c=0.1$ , and $S_c=10$ at $t/T=50$ . It can be observed that the contribution of the viscous flux to the total flux is very low for $S_c=1$ and $S_c=10$ . However, $S_c=0.1$ , this contribution is higher, and it is not negligible, even comparable to the turbulent flux within the canopy region ( $z/h<1$ ) and in the upper part of the domain. Figure 17 depicts the temporal evolution of the total flux $J_t$ (kg·m⁻²·s⁻¹) normalized by the diffusive flux near the wall, $J_p$ (kg·m⁻²·s⁻¹), for the three cases of $S_c=1$ , $S_c=0.1$ , and $S_c=10$ . The same quantities are presented at different times in Figure 18. For short times ( $t/T<1$ ), these results indicate that the total flux is equivalent to the diffusive flux, and the maxima of turbulent flux are localized in the roughness sublayer (between $1<z/h<2$ ). In the case where diffusive transfer dominates over turbulent transfer, i.e., $S_c=0.1$ , the ratio between the total flux and diffusive flux decreases over time, reaching 0.25 at $t/T=50$ . Although the diffusive flux $J_p$ decreases, this is not primarily responsible for the significant reduction in the total flux $J_t$ . For the cases $S_c=1$ and $S_c=10$ , an increase in the total flux $J_t$ is observed, as shown in Figure 13(a)). In these cases, the diffusive flux remains constant for longer times. The accumulation of concentration flux observed mainly in the roughness sublayer is more significant for $S_c=1$ than for $S_c=10$ , with $J_t/J_p=65$ and 30, respectively. This difference is also visible in the curve of the evolution of the instantaneous concentration

Figure 16. Vertical profiles of total flux $J_t$ (kg·m⁻²·s⁻¹), turbulent flux $J_{tur}$ (kg·m⁻²·s⁻¹), and viscous flux $J_v$ (kg·m⁻²·s⁻¹) for three Schmidt numbers $S_c=1$ , $S_c=0.1$ , and $S_c=10$ at $t/T=50$ .

Figure 17. Temporal evolution of total flux $J_t$ (kg·m⁻²·s⁻¹) normalized by the diffusive flux near the wall $J_p$ (kg·m⁻²·s⁻¹) for three Schmidt numbers $S_c=1$ , $S_c=0.1$ , and $S_c=10$ .

Figure 18. Temporal evolution of total flux $J_t$ (kg·m⁻²·s⁻¹) normalized by the diffusive flux near the wall $J_p$ (kg·m⁻²·s⁻¹) at different times scaled by T for three Schmidt numbers $S_c=1$ , $S_c=0.1$ , and $S_c=10$ .

averaged within the canopy and above it, where a faster evolution of concentration is observed for $S_c=10$ compared to $S_c=1$ . The appearance of streaks in the case of $S_c=10$ could explain this difference. According to Calmet [5] , the dynamics of the concentration field is related to streaks in the case of sheared walls, and structures with large scale parallel to the flow direction (but often small in the vertical direction) govern the transfer.

3.7. Global Transfer Coefficients and Exchange Velocity

For given hydrodynamic conditions, the diffusive flux of concentration at the wall, $J_p$ (kg/m²·s), is related to the transfer velocity, $K_L$ (m/s), and the concentration disequilibrium, $\Delta c$ (kg/m³), by Equation (21):

$K_L=\frac{J_p}{\Delta c}$ (21)

The temporal evolution of $K_L$ , normalized by $u_{\mathrm{<mo>\; *\; </mo>}}$ , is presented in Figure 19, along with that of the local Sherwood number, $S_h$ , defined by Equation (22):

$S_h=\frac{J_{p}h}{D\Delta c}$ (22)

As expected, the transfer coefficient decreases as $S_c$ increases. In the case of a wall, Calmet [5] found values of $K_L/u_{\mathrm{<mo>\; *\; </mo>}}\approx S_c^{-0.633}$ close to the theoretical value

Figure 19. Temporal evolution of the local Sherwood number $S_{h\mathrm{,}C}$ (a), $S_{h\mathrm{,}E}$ (b) and the transfer coefficient $K_L/u_{\mathrm{<mo>\; *\; </mo>}}$ (c) for the three Schmidt numbers $S_c=1$ , $S_c=0.1$ , and $S_c=10$ . $S_{h\mathrm{,}C}$ is the Sherwood number calculated using the concentration disequilibrium within the canopy $\Delta c=c_p-\langle c_c\rangle$ , and $S_{h\mathrm{,}E}$ is calculated using the concentration disequilibrium between the canopy bottom and the outside $\Delta c=c_p-c_{sl}$ , where $\langle c_c\rangle$ and $c_{sl}$ are the instantaneous average concentration within the canopy and the concentration at the free surface level.

of $K_L/u_{\mathrm{<mo>\; *\; </mo>}}\approx S_c^{-n}$ with $n=-2/3$ . When considering the temporal evolution of $K_L$ , the found values of n encompass −2/3. The temporal evolution of $S_h$ for the three Schmidt numbers shows an increase over time for $S_c=0.1$ , unlike the cases of $S_c=1$ and $S_c=10$ . The choice to calculate the concentration disequilibrium $\Delta c$ within the canopy and outside is justified by the fact that even though concentration profiles do not converge to a steady-state profile in the upper part of the domain, they seem to do so within the canopy. $S_{h\mathrm{,}C}$ and $S_{h\mathrm{,}E}$ can be interpreted as the ratios between total mass transfer and mass transfer by diffusion within the canopy and in the free flow. It is observed that within the canopy, the values of $S_{h\mathrm{,}C}$ are relatively close (ranging between 2 and 3) for the three Schmidt numbers. In the upper part of the domain, mass transfer by diffusion is lower for all Schmidt numbers. For $S_c=1$ and $S_c=10$ , $S_{h\mathrm{,}E}$ quickly evolves towards a constant value (around 20) for both, contrary to $S_c=0.1$ , where $S_{h\mathrm{,}E}$ increases indefinitely. In the latter case, $S_{h\mathrm{,}E}$ ‘s gradual increase over time can be explained by a less significant concentration disequilibrium due to greater mixing. In addition to the values of $K_L$ and $S_h$ , we have also calculated the exchange velocity between the canopy and the flow above. This parameter is purely hydrodynamic and does not depend on the Schmidt number. For the three cases studied, we found very close values (see Table 2).

Table 2. Exchange velocity for the three Schmidt numbers $S_c=1$ , $S_c=0.1$ , and $S_c=10$ .

4. Conclusion

This comprehensive study of mass transfer in the interaction between river biofilm and a turbulent boundary layer unveils complex and varied dynamics. Our numerical simulations, constrained to relatively low Schmidt numbers and a specific roughness Reynolds number, have highlighted several key aspects of these mass transfers. Firstly, we observed that the concentration profile above the canopy follows a logarithmic trend, with a relative homogenization of concentration over extended durations. This characteristic suggests a dynamic interaction between the biofilm and the flow above, influenced by turbulence and the morphology of the biofilm surface. Secondly, our analysis revealed that the contribution of molecular diffusion to the total flux varies significantly with the Schmidt number. While this contribution is almost negligible for Schmidt numbers equal to or greater than 1, it becomes comparable to the turbulent flux for a Schmidt number of 0.1. This observation underscores the importance of molecular diffusion in low Schmidt number regimes. Moreover, the local Sherwood number, representing the ratio between the total flux and the molecular diffusion flux at the wall, also depends on the Schmidt number. This dependency indicates that nutrient absorption by organisms near the wall is facilitated by lower Schmidt numbers. Interestingly, the exchange velocity, a purely hydrodynamic parameter, proved to be independent of the Schmidt number. In our case, it is on the order of 10% of the friction velocity ( $u_{\mathrm{<mo>\; *\; </mo>}}$ ). This finding suggests a fundamental interaction between the flow structure and the mass transfer process, regardless of specific diffusion properties. In conclusion, our study significantly contributes to the understanding of mass transfer mechanisms in river biofilm systems, highlighting the complex interaction between flow dynamics, molecular diffusion, and biofilm structure. These results pave the way for future research that could explore broader ranges of Schmidt and Reynolds numbers, to expand our understanding of these crucial ecological systems.

[1] Font, R.A., Khamis, K., Milner, A.M., Smith, G.H.S. and Ledger, M.E. (2021) Low Flow and Heatwaves Alter Ecosystem Functioning in a Stream Mesocosm Experiment. Science of the Total Environment, 777, Article 146067. https://doi.org/10.1016/j.scitotenv.2021.146067

[2] Koch, R., Kerling, D. and O’Brien, J.M. (2018) Effects of Chronic Elevated Nitrate Concentrations on the Structure and Function of River Biofilms. Freshwater Biology, 63, 1199-1210. https://doi.org/10.1111/fwb.13126

[3] Boulêtreau, S., Salvo, E., Lyautey, E., Mastrorillo, S. and Garabetian, F. (2012) Temperature Dependence of Denitrification in Phototrophic River Biofilms. Science of the Total Environment, 416, 323-328. https://doi.org/10.1016/j.scitotenv.2011.11.066

[4] Nepf, H.M. (2012) Flow and Transport in Regions with Aquatic Vegetation. Annual Review of Fluid Mechanics, 44, 123-142. https://doi.org/10.1146/annurev-fluid-120710-101048

[5] Calmet, I. (1995) Analyse Par Simulation Des Grandes Échelles Des Mouvements Turbulents et Du Transfert de Masse Sous Une Interface Plane. Ph.D. Thesis, Institut National Polytechnique de Toulouse. https://www.theses.fr/1995INPT110H

[6] Nikora, V. (2007) Hydrodynamics of Aquatic Ecosystems: Spatial-Averaging Perspective. Acta Geophysica, 55, 3-10. https://doi.org/10.2478/s11600-006-0043-6

[7] Nikora, V., Larned, S. and Biggs, B.J.F. (2003) Hydrodynamic Effects in Aquatic Ecosystems with a Focus on Periphyton. Recent Research Developments in Fluid Dynamics, 4, 41-70.

[8] Larned, S., Nikora, V. and Walter, R. (2004) Mass-Transfert-Controlled Nitrogen and Phosphorus Uptake by Stream by Stream Periphyton: A Conceptual Model and Experimental Evidence. Limnology and Oceanography, 49, 1992-2000. https://doi.org/10.4319/lo.2004.49.6.1992

[9] Coundoul, F., Bonometti, T., Graba, M., Sauvage, S., Sanchez Pérez, J.-M. and Moulin, F. (2015) Role of Local Flow Conditions in River Biofilm Colonization and Early Growth. River Research and Applications, 31, 350-367. https://doi.org/10.1002/rra.2746

[10] Bonometti, T. and Magnaudet, J. (2007) An Interface-Capturing Method for Incompressible Two-Phase Flows. Validation and Application to Bubble Dynamics. International Journal of Multiphase Flow, 33, 109-133. https://doi.org/10.1016/j.ijmultiphaseflow.2006.07.003

[11] Tougma, J.L.W. (2023) Method of Analytical Resolution of the Navier-Stokes Equations. Open Journal of Fluid Dynamics, 13, 226-231. https://doi.org/10.4236/ojfd.2023.135017

[12] Badahmane, A. (2023) Linear System Solutions of the Navier-Stokes Equations with Application to Flow over a Backward-Facing Step. Open Journal of Fluid Dynamics, 13, 133-143. https://doi.org/10.4236/ojfd.2023.133011

[13] Legendre, D. (1996) Quelques aspects des forces hydrodynamiques et des transferts de chaleur sur une bulle sphérique.

[14] Yuki, Y., Takeuchi, S. and Kajishima, T. (2007) Efficient Immersed Boundary Method for Strong Interaction Problem of Arbitrary Shape Object with the Self-Induced Flow. Journal of Fluid Science and Technology, 2, 1-11. https://doi.org/10.1299/jfst.2.1

[15] Kim, J. and Choi, H. (2004) An Immersed-Boundary Finite-Volume Method for Simulation of Heat Transfer in Complex Geometries. Journal of Mechanical Science and Technology, 18, 1026-1035. https://doi.org/10.1007/BF02990875

[16] Eckert, E.R.G. and Soehngen, E. (1952) Distribution of Heat Transfer Coefficients Around Ciculair Cylinders in Crossflow at Reynolds Numbers from 20 to 500. Transactions of the American Society of Mechanical Engineers, 75, 343-347. https://doi.org/10.1115/1.4015778

[17] Florens, E. (2010) Couche limite turbulente dans les écoulements à surface libre: Etude expérimentale d’effets de macro-rugosités. Ph.D. Thesis, Université Paul Sabatier, Toulouse III.

[18] Coceal, O., Thomas, T.G., Castro, I.P. and Belcher, S.E. (2006) Mean Flow and Turbulence Statistics over Groups of Urban-Like Cubical Obstacles. Boundary-Layer Meteorology, 121, 491-519. https://doi.org/10.1007/s10546-006-9076-2

NOTES

*These authors contributed equally to this work.

^#Corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Font, R.A., Khamis, K., Milner, A.M., Smith, G.H.S. and Ledger, M.E. (2021) Low Flow and Heatwaves Alter Ecosystem Functioning in a Stream Mesocosm Experiment. Science of the Total Environment, 777, Article 146067. https://doi.org/10.1016/j.scitotenv.2021.146067
[2]	Koch, R., Kerling, D. and O’Brien, J.M. (2018) Effects of Chronic Elevated Nitrate Concentrations on the Structure and Function of River Biofilms. Freshwater Biology, 63, 1199-1210. https://doi.org/10.1111/fwb.13126
[3]	Boulêtreau, S., Salvo, E., Lyautey, E., Mastrorillo, S. and Garabetian, F. (2012) Temperature Dependence of Denitrification in Phototrophic River Biofilms. Science of the Total Environment, 416, 323-328. https://doi.org/10.1016/j.scitotenv.2011.11.066
[4]	Nepf, H.M. (2012) Flow and Transport in Regions with Aquatic Vegetation. Annual Review of Fluid Mechanics, 44, 123-142. https://doi.org/10.1146/annurev-fluid-120710-101048
[5]	Calmet, I. (1995) Analyse Par Simulation Des Grandes échelles Des Mouvements Turbulents et Du Transfert de Masse Sous Une Interface Plane. Ph.D. Thesis, Institut National Polytechnique de Toulouse. https://www.theses.fr/1995INPT110H
[6]	Nikora, V. (2007) Hydrodynamics of Aquatic Ecosystems: Spatial-Averaging Perspective. Acta Geophysica, 55, 3-10. https://doi.org/10.2478/s11600-006-0043-6
[7]	Nikora, V., Larned, S. and Biggs, B.J.F. (2003) Hydrodynamic Effects in Aquatic Ecosystems with a Focus on Periphyton. Recent Research Developments in Fluid Dynamics, 4, 41-70.
[8]	Larned, S., Nikora, V. and Walter, R. (2004) Mass-Transfert-Controlled Nitrogen and Phosphorus Uptake by Stream by Stream Periphyton: A Conceptual Model and Experimental Evidence. Limnology and Oceanography, 49, 1992-2000. https://doi.org/10.4319/lo.2004.49.6.1992
[9]	Coundoul, F., Bonometti, T., Graba, M., Sauvage, S., Sanchez Pérez, J.-M. and Moulin, F. (2015) Role of Local Flow Conditions in River Biofilm Colonization and Early Growth. River Research and Applications, 31, 350-367. https://doi.org/10.1002/rra.2746
[10]	Bonometti, T. and Magnaudet, J. (2007) An Interface-Capturing Method for Incompressible Two-Phase Flows. Validation and Application to Bubble Dynamics. International Journal of Multiphase Flow, 33, 109-133. https://doi.org/10.1016/j.ijmultiphaseflow.2006.07.003
[11]	Tougma, J.L.W. (2023) Method of Analytical Resolution of the Navier-Stokes Equations. Open Journal of Fluid Dynamics, 13, 226-231. https://doi.org/10.4236/ojfd.2023.135017
[12]	Badahmane, A. (2023) Linear System Solutions of the Navier-Stokes Equations with Application to Flow over a Backward-Facing Step. Open Journal of Fluid Dynamics, 13, 133-143. https://doi.org/10.4236/ojfd.2023.133011
[13]	Legendre, D. (1996) Quelques aspects des forces hydrodynamiques et des transferts de chaleur sur une bulle sphérique.
[14]	Yuki, Y., Takeuchi, S. and Kajishima, T. (2007) Efficient Immersed Boundary Method for Strong Interaction Problem of Arbitrary Shape Object with the Self-Induced Flow. Journal of Fluid Science and Technology, 2, 1-11. https://doi.org/10.1299/jfst.2.1
[15]	Kim, J. and Choi, H. (2004) An Immersed-Boundary Finite-Volume Method for Simulation of Heat Transfer in Complex Geometries. Journal of Mechanical Science and Technology, 18, 1026-1035. https://doi.org/10.1007/BF02990875
[16]	Eckert, E.R.G. and Soehngen, E. (1952) Distribution of Heat Transfer Coefficients Around Ciculair Cylinders in Crossflow at Reynolds Numbers from 20 to 500. Transactions of the American Society of Mechanical Engineers, 75, 343-347. https://doi.org/10.1115/1.4015778
[17]	Florens, E. (2010) Couche limite turbulente dans les écoulements à surface libre: Etude expérimentale d’effets de macro-rugosités. Ph.D. Thesis, Université Paul Sabatier, Toulouse III.
[18]	Coceal, O., Thomas, T.G., Castro, I.P. and Belcher, S.E. (2006) Mean Flow and Turbulence Statistics over Groups of Urban-Like Cubical Obstacles. Boundary-Layer Meteorology, 121, 491-519. https://doi.org/10.1007/s10546-006-9076-2