An upscaled model for permeable biofilm in a thin channel and tube

In this paper, we derive upscaled equations for modelling biofilm growth in porous media. The resulting macro-scale mathematical models consider permeable multi-species biofilm including water flow, transport, detachment and reactions. The biofilm is composed of extracellular polymeric substances (EPS), water, active bacteria and dead bacteria. The free flow is described by the Stokes and continuity equations and the water flux inside the biofilm by the Brinkman and continuity equations. The nutrients are transported in the water phase by convection and diffusion. This pore-scale model includes variations of the biofilm composition and size due to reproduction of bacteria, production of EPS, death of bacteria and shear forces. The model includes a water-biofilm interface between the free flow and the biofilm. Homogenization techniques are applied to obtain upscaled models in a thin channel and a tube, by investigating the limit as the ratio of the aperture to the length $\varepsilon$ of both geometries approaches to zero. As $\varepsilon$ gets smaller, we obtain that the percentage of biofilm coverage area over time predicted by the pore-scale model approaches the one obtained using the effective equations, which shows a correspondence between both models. The two derived porosity-permeability relations are compared to two empirical relations from the literature. The resulting numerical computations are presented to compare the outcome of the effective (upscaled) models for the two mentioned geometries.


Introduction
Biofilms are sessile communities of bacteria housed in a self-produced adhesive matrix consisting of extracellular polymeric substances (EPS), including polysaccharides, proteins, lipids and DNA ( [1]). The proportion of EPS in biofilms is 50% to 90% of the total organic matter ( [9,33]). Water is by far the largest component of the matrix, giving biofilms the nickname 'stiff water' ( [11]). Biofilms provoke chronic bacterial infection, infection on medical devices, deterioration of water quality and the contamination of food ([15]). On the other hand, biofilms can be used for wastewater treatment and bioenergy production ([20]). In microbial enhanced oil recovery (MEOR), one of the strategies is selective plugging, where bacteria are used to form biofilm in highly permeable zones to diverge the water flow and extract the oil located in low-permeability zones ([24]).
In wastewater treatment, one of the strategies consists of using biofilms to break down compounds which are not desirable to discharge into the natural environment ( [5]). Two of the motivations to derive upscaled models are to accurately describe the average behaviour of the system with relatively low computational effort compared to fully detailed calculations starting at the microscale ( [31]) and to determine effective parameters ( [12]). The values of these effective parameters can be determined using known values from pore-scale experiments. Recent works have been carried out to derive upscaled models, e.g., [7] obtained a mathematical model describing macroscopic tumour growth, transport of drug and nutrient through homogenization and [14] upscaled a pore-scale model for primary fluid recovery and showed that the macroscopic equation for the water flux is fundamentally different from Darcys' law. We also refer to [23] and [25], where the authors upscaled various pore-scale models for biofilm formation in perforated and strip geometries.
The present work builds on [17], where a pore-scale model is discussed. The model includes permeable biofilm and evolution of different biofilm components: active bacteria, dead bacteria and EPS. The importance of including biofilm permeability is underlined by the fact that the dominated mechanism of nutrient transport within some biofilms is convection ([18]). This mathematical model is based on laboratory experiments performed by [19], where the biofilm was grown in micro-channels. Here we upscale this pore-scale model to derive effective equations, by investigating the limit as the ratio of the height to the length of the micro-channel approaches to zero.
In this general context, the objective of the research reported in the present article is to obtain corescale models (also known as Darcy-scale or macro-scale models) for permeable biofilm in two different pore geometries: a thin channel and a tube. The motivation for choosing these two geometries is because experiments are performed in the laboratory in micro-channels ([19]) and tubes ( [3]), they may represent a fracture in a core sample ( [4]) and some porous media can be modelled as a stack of micro-tubes or micro-channels ( [31]).
properties t = T Initially,

Water and nutrients
Water flux: Nutrients: Biomass: Water-biofilm interface: Inside and outside the difussion.
convection and Made of EPS, active and dead bacteria.
Changes in time due to EPS production, bacterial reproduction, death of bacteria and detachment.
Biomass, water and nutrients biofilm.  To summarize, the novel aspect in this work is the derivation of core-scale models from a pore-scale model for a biofilm which is permeable to the flow and has a variable (in time and space) height. The fluid flow in the biofilm is modelled by the Brinkman equation, whereas in the remaining pore space the Stokes model is adopted. This is done for two different geometries. We derive analytical expressions for the upscaled quantities and provide numerical simulations for the upscaled models in both cases.
The structure of this paper is as follows. In Sec. 2, we describe the pore-scale biofilm model. In Sec. 3, we present the dimensionless pore-scale biofilm model. In Sec. 4, we perform formal homogenization on the model equations and obtain upscaled equations. In Sec. 5, we compare the upscaled models with the upscaled model of [31] and with the well-known core-scale model of [6]. We compare the derived porosity-permeability relations to empirical porosity-permeability relations from the literature. Also, we perform numerical simulations in the upscaled models and we compare the results for the biofilm height and nutrient concentration for the two different effective models. Finally, in Sec. 6 we present the conclusions.

Pore-scale model
The pore-scale mathematical model considered here follows ideas from [2] , [31] and [8]. A detailed description of this model can be found in [17], where a comparison of laboratory measurements and numerical simulations is also presented.
The biofilm has four components: water, EPS, active and dead bacteria ( j = {w, e, a, d}). Let θ j and ρ j denote the volume fraction and the density of species j. The sum of volume fractions is constraint to 1 (θ w + θ e + θ a + θ d = 1), where the volume fraction of water θ w is assumed constant. The biomass phases and water are assumed to be incompressible (∂ t ρ j = 0) and the biofilm layer is attached to the pore walls. Fig. 1 shows schematically the phenomena considered for the biofilm formation.
We consider two different pore geometries: a tube in cylindrical coordinates r r r = (r, ϕ, z) and a thin channel in Cartesian coordinates x x x = (x, y, z). The z direction is taken along the length L of the tube and thin channel (see Figs. 2 and 6). In the first case, the pore has circular cross-section and in the second a rectangular one. In both cases, the length is much larger than the cross-sectional aperture. In both cases, we assume a certain symmetry. For the cylindrical pore we assume that the processes are radially symmetric, hence there is no angular dependence (see Fig. 2). For the thin channel, there are no changes in the x direction, i.e., the width of the channel, so it can be reduced to a two-dimensional strip (see Fig. 6). This assumption is based on experiments, showing that when the width of the channel is much smaller than its height, the growth of the biofilm occurs only at the upper and lower walls along the channel [19]. We present in detail the upscaling of the model equations on the tube geometry. The upscaling on the channel geometry is shown in Appendix A. Fig. 2 shows the different domains, boundaries and interface in the pore with tubular geometry. We consider a thin tube of radius and length L. We denote the biofilm height by d which only depends on z and time as a result of the symmetry assumption. The domain is occupied by the water Ω T,w (t) = {r r r | r ∈ [0, − d(z, t)), ϕ ∈ [0, 2π), z ∈ (0, L)} and biofilm Ω T,b (t) = {r r r | r ∈ ( − d(z, t), ), ϕ ∈ [0, 2π), z ∈ (0, L)} phases with the biofilm located along the tube wall Γ T,s = {r r r | r = , ϕ ∈ [0, 2π), z ∈ (0, L)}. Clearly r = − d(z, t) separates the water and biofilm regions. The water domain has three boundary parts: the inflow Γ T,iw (t) = {r r r | r ∈ [0, − d(z, t)), ϕ ∈ [0, 2π), z = 0}, the outflow Γ T,ow (t) = {r r r | r ∈ [0, − d(z, t)), ϕ ∈ [0, 2π), z = L} and the interface between the water and the biofilm Γ T,wb (t) = {r r r | r = − d(z, t), ϕ ∈ [0, 2π), z ∈ (0, L)}. Similarly, for the biofilm we have the inflow Γ T,ib (t) = {r r r | r ∈ ( − d(z, t), ), ϕ ∈ [0, 2π), z = 0} and outflow Γ T,ob (t) = {r r r | r ∈ ( − d(z, t), ), ϕ ∈ [0, 2π) z = L} boundary parts, the water-biofilm interface Γ T,wb and the solid tube wall Γ T,s . Although the tube is a three-dimensional domain, recalling the rotational symmetry, we only write the rand z-components of the vectors in order to reduce the length of the mathematical expressions.
The unit normal ν ν ν pointing into the biofilm and the normal velocity of the interface ν n can be written in terms of the biofilm height d as [31] The water flux outside the biofilm Ω T,w (t) is described by the Stokes and continuity equations µ∆q q q w = ∇p w , ∇ · · · q q q w = 0, while the water flux inside the biofilm Ω T,b (t) is described by the Brinkman and continuity equations Here p w and p b are the water pressures and q q q w and q q q b are the water velocities in the water domain and biofilm domain respectively; µ is the water viscosity (constant, not dependent on biofilm species) and k is the permeability of the biofilm (assumed isotropic). At the interface Γ T,wb (t) one has the continuity of the velocity and of the normal stress tensor where I is the identity matrix. At the wall Γ T,s we consider the no-slip boundary condition q q q b = 0 0 0. To model the nutrient transport and consumption, we let c α (α ∈ {w, b}) stand for the nutrient concentration in water or biofilm (mass per total volume of biofilm) respectively and D is the nutrient diffusion coefficient in water. Then, the nutrients in the water Ω T,w (t) and the biofilm Ω T,b (t) satisfy the convection-diffusion equation where J J J w and J J J b are given by Further at Γ T,wb (t) we impose the mass conservation and we assume continuity of nutrients At the solid wall Γ T,s the normal flux is υ υ υ · · · J J J b = 0, where υ υ υ is the normal vector at the pore wall. The reaction term R b for the consumption of nutrients is given by where µ n is the maximum rate of nutrient consumption and k n is the Monod-half nutrient velocity coefficient. The movement of the biomass components θ i in Ω T,b (t) due to reproduction, production of EPS and death of active bacteria can be modelled as a Darcy flow [2]. We denote by u u u the velocity of the biomass and Φ the growth velocity potential. Then, we consider the following equations [2,17] The growth velocity potential is set to zero Φ = 0 at the interface Γ T,wb (t) and homogeneous Neumann boundary condition υ υ υ · · · ∇Φ = 0 at the wall Γ T,s . For each of the biomass components θ i in Ω T,b (t), we assume mass conservation [2] and Neumann condition ν ν ν · · · ∇θ i = 0 at the interface Γ T,wb (t) and υ υ υ · · · ∇θ i = 0 at the wall Γ T,s . The reaction terms for the biomass components are given by where Y e and Y a are yield coefficients and k res is the bacterial decay rate. The water-biofilm interface changes in time due to changes inside the biofilm and the water flux provoking detachment of components, which is known as erosion. Thus, the normal velocity of the interface Γ T,wb (t) is given by [31] where we ensure that the biofilm-water interface does not overlap by taking the positive cut f Here k str is the stress coefficient and S is the tangential shear stress, given by [31] This pore-scale model can be extended to consider more complex systems. For example, one can add different kind of nutrients, different active bacteria species in the biofilm or bacterial attachment.

Non-dimensional model
Before seeking an effective model, we bring the mathematical equations to a non-dimensional form. To this aim, we introduce the reference time t re f , length L re f , radius re f , water velocity q re f = L re f /t re f , biomass velocity u r e f , pressure p r e f and concentration c re f . The thin tube is characterized by the ratio of its radius to the length ε = r e f /L r e f , which is called the dimensionless aspect ratio. We define dimensionless coordinates and time asr = r/ r e f ,z = z/L re f andt = t/t re f . The non-dimensional biofilm height is given byd = d/ r e f . The non-dimensional unit normal (1) is given byν ν ν(r,z) = (1, ε∂zd)/(1 + (ε∂zd) 2 ) 1/2 . We notice that a factor of ε appears in the second component of the non-dimensional unit normal, as a result of the transformation of the coordinates Note that we have omitted the dependence of the vector variables on ϕ (ν ν ν(r, ϕ,z) =ν ν ν(r,z)). This is justified by our assumption of the radial symmetry. The non-dimensional nutrient concentrations and densities are given The water velocities are given byq q q w (r,z) The biomass velocity is given byū ū u(r,z) = (ūr,ūz) = (u r /(εu re f ), u z /u re f ). We assume that the velocities in the radial direction are of the order ρ re f /t re f (see [31]). Hence, they scale by 1/ε when compared to the longitudinal velocities. The biomass volume fractions are dimensionless; therefore, in the non-dimensional model we simply defineθ i = θ i , i ∈ {w, e, a, d}. Finally, the pressures and growth velocity potential We observe that the growth velocity potential Φ is scaled by 1/ε 2 in order to have the biomass velocities in the radial direction of the order ρ r e f /t r e f (see [31]). We define the following dimensionless parameters Pe = q re f L re f /D,μ n = t re f µ n , k n = k n /c r e f ,k = k/ 2 r e f ,μ = µL r e f q r e f /( 2 re f p re f ),k str = p re f k str /u re f andk res = t re f k res . In this way, the dimensionless system of equations for the water flux (2-4) is given by where (14-16) are the dimensionless Stokes and continuity equations, (17-19) are the dimensionless Brinkman and continuity equations, (20-22) are the dimensionless interface conditions and (23) is the dimensionless condition on the wall. The dimensionless equations for the transport of nutrients (5)(6)(7) in the water and biofilm are given by where (24) (28) is the dimensionless condition on the wall. The dimensionless reaction rate (8) for the consumption of nutrients is given byR b = −μ nθaρacb /(k n +c b ). The equations for the growth velocity potential (9) become where (29)(30) are the dimensionless equations for the biomass growth velocity potential, (31) is the dimensionless reference potential at the interface and (32) is the dimensionless condition on the wall. We define the dimensionless sum of the biomass reaction terms asΣ = (Y eρa /ρ e +Y a )μ nθacb /(k n +c b )+(ρ a /ρ d −1)k resθa .
The equations for the biomass components (10) become where (33)(34)(35) are the dimensionless conservation of mass equations for the biomass components, (36) is the dimensionless condition at the interface and (37) is the dimensionless condition on the wall. The dimensionless biofilm height (12) is given by The dimensionless tangential shear stress (13) is given bȳ where the matrixM is given byM

Upscaling
The pore-scale mathematical model describes the biofilm formation in a three-dimensional domain. Under some model assumptions, when the length of the tube is much larger than its radius, it is possible to reduce the dimensionality of the problem from three to one dimension, letting the aspect ratio ε approach to zero. We perform a formal asymptotic expansion of the variables depending on ε, In [31], [16] and [4], the authors present upscaled models of pore-scale mathematical models for reactive flows. Following the same ideas, we can obtain the corresponding upscaled model in the tube pore geometry. We define the average water velocity q as the following integral Notice that we divide by the cross-sectional area of the tube. We consider the following space regions in the tube:Ξ These regions are a disk of radius 1 −d and a ring of thicknessd respectively; both of length δz. Integrating (14) and (17) over the previous regions and using the Gauss's theorem, we obtain 0 = ∫Ξ w∇ · · ·q q q w dV + ∫Ξ b∇ · · ·q q q b dV = 2π Recalling the no-slip condition for the water flux on the wall (23) and the continuity of fluxes at the interface (22), the previous equation becomes Dividing the previous equation by δz and in the limit where δz approach to zero, we obtain for the lowest-order terms in ε ∂z q = ∂z q w (z,t) + ∂z q b (z,t) = 0, where we have used the definition of the average water velocity q (41). The lowest order terms in the Stokes model (14-16) lead to ∂r (rq w,r,0 )/r + ∂zq w,z,0 = 0, ∂rp w,0 = 0,μ∂r (r∂rq w,z,0 )/r = ∂zp w,0 . (42a, b, c) From (42b), we conclude thatp w,0 does not depend on ther coordinate. Analogously, for the Brinkman model (17-19), the lower-order terms in ε give From (43b), we conclude thatp b,0 also does not depend on ther coordinate. Since neitherp w,0 norp b,0 depend on ther coordinate and from the lowest order terms in (20) we have thatp w,0 =p b,0 at the biofilm-water interface, we obtain thatp w,0 (z,t) =p b,0 (z,t) =p 0 (z,t). We turn our attention to equations (42c) and (43c). It is possible to find solutions for q w,z,0 and q b,z,0 integrating twice with respect tor both equations and in addition using the symmetry, interface and boundary conditions (21-23). After integration, we get where J ν (z) and Y ν (z) are the Bessel functions of order ν of first and second kind respectively (see [22]). The coefficients appearing in (44) are , where w = 1 −d 0 , ξ = i(θ w /k) 1/2 and i is the imaginary number. We remark that most of the mathematical commercial software include Bessel functions; therefore, it is easy to use the above expression. Although the Bessel functions are evaluated with complex numbers, both fluxes q w,z,0 and q b,z,0 are real numbers.
To obtain the water velocity defined in (41), we integrate (44) as follows

This gives the Darcy's law q
is the effective permeability given by which changes according to the biofilm heightd 0 .

Discussion and comparison with other biofilm models
The extension from a channel (tube) to a porous medium is done by considering a stack of channels (tubes) of void space and solid material [31], where we denote by φ the porosity of the porous medium. As [26] and [31], we assume that all tubes have the same diameter and are equally separated. Therefore, multiplying the upscaled model equations by φ, the corresponding core-scale mathematical models are obtained. Table  2 shows the core-scale equations of the van Noorden model [31], the porous medium formed by channels and the porous medium formed by tubes, where v = φq is the Darcy velocity. The van Noorden model accounts for water flux, transport of nutrients and bacteria, bacterial attachment, detachment of biomass due to erosion, growth of biomass due to nutrient consumption and death of bacteria. In our model we do not include transport of bacteria and bacterial attachment. The reason is because the pore-scale model was built based on laboratory experiments, where only nutrients were injected after inoculation of bacteria [17]. For the details of the upscaling on the thin channel domain, see Appendix A.
From Table 2, we observe that for the Darcy flow, the permeability is different for the three models. A discussion of these relations is given at the end of this Section. For the equations describing the transport of nutrients, the difference between the van Noorden and channel model is on the reaction term. For the van Noorden model the nutrient consumption depends on the biofilm height, while in the channel model depends also in the volume fraction of active bacteria. Comparing the tube and channel model, we observe a different function of the biofilm height, due to the cylindrical geometry. For the biofilm height, the difference between the van Noorden and channel model is on the Σ term, where for the van Noorden model the total sum of reactions accounts for bacterial reproduction and decay, while in the channel structure it also accounts for EPS. Finally, for the bacterial, EPS and dead bacterial volume fractions, the model equations are the same for the channel and tube models, while for the van Noorden model the active bacterial volume fraction is constant with value 1.
In [31], the authors compared their upscaled model with a well-known macro-scale model by [28], where a mathematical model for an impermeable single-species biofilm including flow, transport and reactions is built. We compare the two derived uspscaled models with a macro-scale model by [6], where a mathematical model for impermeable multi-species biofilm including water flow, transport of nutrients and reactions is built (Table 3). In this model the authors considered a biofilm formed by two bacterial species which consume nutrients, reproduce and die. The porosity of the core decreases as the biofilm grows which decreases the permeability of the core.
In Table 3, K is the permeability, θ B and θ K the volume fractions of the contaminant-degrading microbe and the strong biofilm-forming microbe, φ O and K O the clean surface porosity and initial permeability respectively and η an experimentally determined parameter. In [6], the porosity decreases as the component concentrations increases. In our models, the porosity in the porous medium decreases as the biofilm height increases. However, in [6] the authors do not include the detachment effects. The permeability in both tube Table 2: Core-scale equations for the a channel, b tube and c van Noorden models Name Upscaled equation Name Equation and channel models have different functions as a result of the different geometries and also because of the water flow inside the biofilm. Notice that there is a quartic function of the biofilm height in one of the permeability terms in the porous medium formed by tubes, as proposed in [27] and [21]. For the nutrient transport, those models, in addition to the [6] model, have the same form with different reaction rates. Porosity-permeability relations for evolving pore space is an active research field (see [13] for a review of these relations). [29] present the following relation which includes the biofilm permeability k where φ crit is the critical porosity at which the permeability becomes zero. [32] proposed the following relation of permeability and porosity for a plugging model where B is a relative porosity given by B = 1− φ/φ O and B c is the critical point where biofilm begins to detach and form plugs. Fig. 3 shows our derived porosity-permeability relations, the one derived by van Noorden and the two proposed relations by Thullner and Vandevivere for different values of biofilm permeability. The values of parameters are φ crit = 0, η = 1.76, B c = 0.1 [13] and θ w = 0.1. For a biofilm with high permeability k = 10 −1 , we observe a faster reduction of permeability for Vandevivere. For a biofilm with low permeability k = 10 −3 , the channel and tube relations approach the one for van Noorden in contrast with the relations of Thullner and Vandevivere. In general, we observe different behaviours of the relations as the porosity decreases. We perform numerical simulations considering both effective models (channel and tube) to compare the biofilm height over time. We consider two different porous media of length L = 0.1 m: the first one has pores formed by thin channels of height 2l = 0.2 mm and the second ones with tubes of diameter 2 = 0.2 mm. For the inlet boundary, we set p i = 4 Pa. The injected nutrient concentration is c i = 1 kg m −3 . The porosity Φ is set to 0.4. Recalling that biofilms are mostly composed by water, we set the water volume fraction in the biofilm equal to 90%. We set the initial EPS and active bacterial volume fraction equal to 5%; thus the initial dead bacterial volume fraction is 0. In Table 4, the values of parameters for the numerical simulations are presented.
We implement the model equations in the commercial software COMSOL Multiphysics (COMSOL 5.2a, Comsol Inc, Burlington, MA, www.comsol.com). A decoupled finite element algorithm is used to solve the mathematical model equations. Firstly, we solve for the pressure and concentration. Then, we compute the volume fractions and biofilm height. We iterate between both steps until the error (the difference between successive values of the solution) drops below a given tolerance ζ. We perform numerical simulations and we compare the results of the two upscaled mathematical models.
To check if there is a correspondence between the pore-scale and upscaled models as ε is close to zero, numerical simulations can be done for both models to compare the average solution of one of the variables. Fig. 4 compares the upscaled model with the pore-scale model in the channel for different values of ε, where the percentage of biofilm on the whole domain is plotted over time. We called this coverage area a and for the channel is given by For all numerical simulations, we fixed the value of the height of the channel 2l, where we set the initial biofilm height as d = l/5. Then, the length of the channel was changed accordingly to match the value of   ε. We observe that the coverage area in the pore-scale simulations approaches the one computed from the upscaled models as ε gets smaller. Fig. 5 shows the changes over time of biofilm height and nutrient concentration for both porous media. Initially, the left part (0 < z < L/2) has a biofilm height of d = l/2 (d = /2 for the tubular pores), while the right part (L > z > L/2) has a height of d = l/4 (d = /4 for the tubular pores). This initial condition is given to study the biofilm development after clogging. We observe that the biofilm height increases faster for the pore channels than in the pore tubes. The explanation of this result is that the right-hand side of the equation for the biofilm height in the channel is two times greater than the right-hand side of the equation for the biofilm height in the tube which is obtained after upscaling in the two different geometries. For the porous medium formed by channels, we observe that the biofilm keeps growing even though the left part of the pore is clogged and after 4 days it reaches a stationary state. This result cannot be observed using the van Noorden model because the water flux stops once the channel is clogged. For the nutrients, we observe that the biofilm consumes the nutrients in the porous medium formed by tubes faster, but the biofilm in the porous medium formed by channels consumes more nutrients.

Conclusions
In this work, we upscale a mathematical model for permeable biofilm considering a thin channel and tubular pore geometries. The upscaled models differ mainly in the effective permeability terms which are functions of the biofilm height. As ε gets smaller, we obtain that the percentage of biofilm coverage area over time predicted by the pore-scale model approaches the one obtained using the effective equations, which shows a correspondence between both models. After comparing with the model proposed by [31], it is possible to derive this model as a particular case of the channel model. The derived upscaled models and the [6] model are very similar. In this manner, the upscaling provides additional support for this model. The numerical simulations show that the biofilm height increases faster in the porous medium formed by channels than the one formed by tubes. These two upscaled models could be used to model porous media where the geometries of the fractures are similar to thin channels or tubes. To validate the core-scale upscaled models, designed laboratory experiments are necessary which is the subject of our future research.  Figure 6: Pore of length L, height 2l and width w in Cartesian coordinates.
The dimensionless system of equations for the water flux is given by The equations for the nutrients become whereR b = −μ nθaρacb /(k n +c b ). The dimensionless equations for the growth velocity potential are given by (ũỹ,ũz) = −(∂ỹΦ, 2 ∂zΦ), The equations for the biomass components become For the biofilm height we have q r e f (−ũỹ + ∂zdũz)),d = 1, We define the average water velocity q as the following integral q (z,t) = q w (z,t) + q b (z,t) = 1 2 We define the following space regions in the channel Integrating (57) and (60) over the previous regions and using the Gauss's theorem yield Recalling the no-slip condition for the water flux on the wall (66) and the continuity of fluxes at the interface (65), the previous equation becomes Dividing the previous equation by δz and letting δz approach zero, we obtain for the lowest-order terms in where we have used the definition of the water velocity q (84). The lowest order terms in the Stokes model (57-59) leads to ∂ỹq w,ỹ,0 + ∂zq w,z,0 = 0, ∂ỹp w,0 = 0,μ∂ 2 yq w,z,0 = ∂zp w,0 .
Letting δ go to zero in order to return to the nonregularized formulation, we get