Homogenization of biomechanical models for plant tissues

In this paper homogenization of a mathematical model for plant tissue biomechanics is presented. The microscopic model constitutes a strongly coupled system of reaction-diffusion-convection equations for chemical processes in plant cells, the equations of poroelasticity for elastic deformations of plant cell walls and middle lamella, and Stokes equations for fluid flow inside the cells. The chemical process in cells and the elastic properties of cell walls and middle lamella are coupled because elastic moduli depend on densities involved in chemical reactions, whereas chemical reactions depend on mechanical stresses. Using homogenization techniques we derive rigorously a macroscopic model for plant biomechanics. To pass to the limit in the nonlinear reaction terms, which depend on elastic strain, we prove the strong two-scale convergence of the displacement gradient and velocity field.


Introduction
Analysis of interactions between mechanical properties and chemical processes, which influence the elasticity and extensibility of plant cell tissues, is important for better understanding of plant growth and development, as well as their response to environmental changes. Plant tissues are composed of cells surrounded by cell walls and connected by a cross-linked pectin network of middle lamella. Plant cell walls must be very strong to resist high internal hydrostatic pressure and at the same time flexible to permit growth. It is supposed that calcium-pectin cross-linking chemistry is one of the main regulators of plant cell wall elasticity and extension [51]. Pectin is deposited to cell walls in a methylesterified form. In cell walls and middle lamella pectin can be modified by the enzyme pectin methylesterase (PME), which removes methyl groups by breaking ester bonds. The de-esterified pectin is able to form calciumpectin cross-links, and thus stiffen the cell wall and reduce its expansion, see e.g. [50]. On the other hand, mechanical stresses can break calcium-pectin cross-links and hence increase the extensibility of plant cell walls and middle lamella. It has been shown that chemical properties of pectin and the control of the density of calcium-pectin cross-links greatly influence the mechanical deformations of plant cell walls [34], and the interference with PME activity causes dramatic changes in growth behavior of plant tissues [50].
To analyze the interactions between calcium-pectin dynamics and deformations of a plant tissue, we derive a mathematical model for plant tissue biomechanics at the length scale of plant cells. In the microscopic model we consider a system of reaction-diffusion-convection equations describing the dynamics of the methylesterified pectin, demethylesterified pectin, calcium ions, and calcium-pectin cross-links. Elastic deformations and water flow are modelled by the equations of poroelasticity for cell walls and middle lamella coupled with the Stokes system for the flow velocity inside cells. The interplay between the mechanics and the chemistry comes in by assuming that the elastic properties of cell walls and middle lamella depend on the density of the calcium-pectin cross-links and that the stress within cell walls and middle lamella can break the cross-links. Thus the microscopic problem is a strongly coupled system of the Stokes equations, reaction-diffusion-convection equations, with reaction terms depending on the displacement gradient, and equations of poroelasticity, with elastic moduli depending on the density of cross-links. To address the situations when a plant tissue is completely and not completely saturated by water, we consider both evolutional and quasi-stationary equations of poroelasticity.
To show the existence of a weak solution of the microscopic equations we use a classical approach and apply the Banach fixed-point theorem. However, due to quadratic non-linearities of reactionterms, the proof of the contraction inequality is not standard and relies on delicate a priori estimates for the L ∞ -norm of a solution of the reaction-diffusion-convection system in terms of the L 2 -norm of displacement gradient and flow velocity. The Alikakos [2] iteration technique is applied to show the uniform boundedness of some of components of solutions of the microscopic equations.
To analyze effective mechanical properties of plant tissues, we derive rigorously a macroscopic model for plant biomechanics using homogenization techniques. The two-scale convergence, e.g. [3,31], and the periodic unfolding method, e.g. [15], are applied to obtain the macroscopic equations. The main mathematical difficulty in the derivation of the macroscopic problem arises from the strong coupling between the equations of poroelasticity and the system of reaction-diffusion-convection equations. In order to pass to the limit in the nonlinear reaction terms we prove the strong two-scale convergence for the displacement gradient and fluid flow velocity, essential for the homogenization of the coupled problem considered here. Due to the dependence of the elasticity tensor on the time variable, in the proof of the strong two-scale convergence a specific form of the energy functional is considered.
Similar to the microscopic problem, to prove uniqueness of a solution of the macroscopic equations we derive a contraction inequality involving the L ∞ -norm of the difference of two solutions of the reactiondiffusion-convection equations. This contraction inequality also ensures also the well-posedness of the limit system.
The poroelasticity equations, modelling interactions between fluid flow and elastic stresses in porous media, were first obtained by Biot using a phenomenological approach [8,9,10] and subsequently derived by applying techniques of homogenization theory. Formal asymptotic expansion was undertaken by [5,13,23,42] to derive Biot equations from microscopic description of elastic deformations of a solid matrix and fluid flow in porous space. The rigorous homogenization of the coupled system of equations of linear elasticity for a solid matrix combined with the Stokes or Navier-Stokes equations for the fluid part was conducted in [17,19,24,32] by using the two-scale convergence method. Depending on the ratios between the physical parameters, different macroscopic equations were obtained, e.g. Biot's equations of poroelasticity, the system consisting of the anisotropic Lamé equations for the solid component, and the acoustic equations for the fluid component, the equations of viscoelasticity. The homogenization of a mathematical model for elastic deformations, fluid flow and chemical processes in a cell tissue was considered in [20]. In contrast to the problem considered in the present paper, in [20] the coupling between the equations of linear elasticity and reaction-diffusion-convection equations for a concentration was given only through the dependence of the elasticity tensor on the chemical concentration. The existence and uniqueness of a solution for equations of poroelasticity were studied in [45,53].
Compared to the many results for the equations of poroelasticity, there exist only few studies of interactions between a free fluid and a deformable porous medium. In [46] nonlinear semigroup method was used for mathematical analysis of a system of poroelastic equations coupled with the Stokes equations for free fluid flow. A rigorous derivation of interface conditions between a poroelastic medium and an elastic body was considered in [26]. Numerical methods for coupled Biot poroelastic system and Navier-Stokes equations were derived in [6]. The Nitsche method for enforcing interface conditions was applied in [12] for numerical simulation of the Stokes-Biot coupled system.
Several results on homogenization of equations of linear elasticity can be found in [7,21,33,42] (and the references therein). Homogenization of the microscopic model for plant cell wall biomechanics, composed of equations of linear elasticity and reaction-diffusion equations for chemical processes, has been studied in [39]. This paper is organized as follows. In Section 2 we derive the microscopic model for plant tissue biomechanics. A priori estimates as well as the existence and uniqueness of a weak solution of the microscopic problem are obtained in Section 3. In Section 4 we prove the convergence results for solutions of the microscopic problem. The multiscale analysis of the coupled poroelastic and Stokes problem is conducted in Sections 5. In Section 6 we show strong two-scale convergence of displacement gradient and flow velocity. The macroscopic equations for the system of reaction-diffusion-convection equations are derived in Section 7. The well-posedness and uniqueness of a solution of the macroscopic problem are proven in Section 8. In Section 9 we consider the incompressible and quasi-stationary cases for the equations of poroelasticity.

Microscopic model
In the mathematical model for plant tissue biomechanics we consider interactions between the mechanical properties of a plant tissue and the chemical processes in plant cells. A plant tissue is composed of cell interior (intracellular space), the plasma membrane, plant cell walls, and the cross-linked pectin network of the middle lamella joining individual cells together. Primary plant cell walls consist mainly of oriented cellulose microfibrils (that strongly influence the cell wall stiffness), pectin, hemicellulose, proteins, and water. It is supposed that calcium-pectin chemistry, given by the de-esterification of pectin and creation and breakage of calcium-pectin cross-links, is one of the main regulators of cell wall elasticity, see e.g. [51]. Hence in our mathematical model we consider the interactions and two-way coupling between calcium-pectin chemistry and elastic deformations of a plant tissue. To describe the coupling between the mechanics and chemistry we consider the dynamics of pectins, calcium, and calcium-pectin cross-links, water flow in a plant tissue and the poroelastic nature of cell walls and middle lamella.
To derive a mathematical model for plant tissue biomechanics, we denote a domain occupied by a plant tissue by Ω ⊂ R 3 , where Ω is a bounded domain with C 1,α boundary for some α > 0. Notice that all results also hold for a two-dimensional domain. Then the time-independent domains Ω f ⊂ Ω and Ω e ⊂ Ω, with Ω = Ω e ∪ Ω f and Ω e ∩ Ω f = ∅, represent the reference (Lagrangian) configurations of the intracellular (cell inside) and intercellular (cell walls and middle lamella) spaces, respectively, and Γ denotes the boundaries between the cell inside and cell walls and corresponds to the plasma membrane. Since Γ represents the interface between elastic material and fluid in the Lagrangian configuration, it is also independent of time.
Pectin is deposited into the cell wall in a highly methylesterified state and is modified by the wall enzyme pectin-methylesterase (PME), which removes methyl groups [50]. It was observed experimentally that pectins can diffuse in a plant cell wall matrix, see e.g. [18,38,48]. Thus in the balance equation for the density of the methylesterified pectin b e,1 and demethylesterified pectin b e,2 ∂ t b e,j + divJ b,j = g b,j in Ω e , j = 1, 2, we assume the flux to be determined by Fick's law, i.e. J b,j = −D be,j ∇b e,j , for j = 1, 2, and D be,j > 0. The term g b,j models chemical reactions, that correspond to the demethylesterification processes and creation and breakage of calcium-pectin cross-links. In general, diffusion coefficients for pectins and calcium depend on the microscopic structure of cell wall given by the cell wall microfibrils and hemicellulose network, which is assumed to be given and not to change in time, as well as on the density of pectins and calcium-pectin cross-links. For presentation simplicity we assume here that the diffusion coefficient does not depend on the dynamics of pectin and calcium-pectin cross-link densities. However the analysis can be conducted in the same way for the generalised model in which the diffusion of pectin, calcium, and cross-links depends on pectin and cross-link densities, assuming that diffusion coefficients are uniformly bounded from below and above, which is biologically sensible. The modification of methylesterified pectin by PME is modelled by the reaction term g b,1 = −µ 1 b e,1 with some µ 1 > 0. For simplicity we assume that there is a constant concentration of PME enzyme in the cell wall. By simple modifications of the analysis considered here, the same results can be obtained for a generalized model including the dynamics of PME and chemical reactions between PME and pectin, see [39] for the derivation of the corresponding system of equations.
The deposition of the methylesterified pectin is described by the inflow boundary condition on the cell plasma membrane. We also assume that the demetylesterified pectin cannot move back into the cell interior: D be,1 ∇b e,1 · n = P 1 (b e,1 , b e,2 , b e,3 ), D be,2 ∇b e,2 · n = 0 on Γ.
To account for mechanisms controlling the amount of pectin in the cell wall, we assume that the inflow of new methylesterified pectin depends on the density of methylesterified and demethylesterified pectin, i.e. b e,1 and b e,2 , and calcium-pectin cross-links b e,3 . We consider the diffusion and transport by water flow of calcium molecules in the symplast (in the cell inside) and diffusion of calcium in apoplast (cell walls and middle lamella), see e.g. [49]. Then the balance equations for calcium densities c f and c e in Ω f and Ω e , respectively, are given by where the chemical reaction term g f = g f (c f ) in Ω f describes the decay and/or buffering of calcium inside the plant cells, see e.g. [52], g e models the interactions between calcium and demethylesterified pectin in cell walls and middle lamella and the creation and breakage of calcium-pectin cross-links, and G is a bounded function of the intracellular flow velocity ∂ t u f . The condition that G is bounded is natural from the biological and physical point of view, because the flow velocity in plant tissues is bounded. This condition is also essential for a rigorous mathematical analysis of the model. We assume that as the result of the breakage of a calcium-pectin cross-link by mechanical stresses we obtain one calcium molecule and two galacturonic acid monomers of demethylesterified pectin. A detailed derivation of the chemical reaction term g e is given in [39]. See also Remark 2.3 for the detailed form of the reaction terms. We assume a passive flow of calcium between cell walls and cell inside and assume that the exchange of calcium between cell insides and cell walls is facilitated only on parts of the cell membrane The regulation of calcium flow by mechanical properties of the cell wall will be considered in future studies. Calcium-pectin cross-links b e,3 are created by electrostatic and ionic binding between two galacturonic acid monomers of pectin chains and calcium molecules. It is also known that these cross-links are very stable and can be disturbed mainly by thermal treatments and/or mechanical forces, see e.g. [35,36]. Thus assuming a constant temperature, the calcium-pectin chemistry can be described as an reaction between calcium molecules and pectins, where the breakage of cross-links depends on the deformation gradient of the cell walls. Hence we assume that the cross-links are disturbed by the mechanical stresses in cell walls and middle lamella, see [39] for a detailed description of the modelling of the calcium-pectin chemistry. A similar approach was used in [41] to model a chemically mediated mechanical expansion of the cell wall of a pollen tube. There are no experimental observations of diffusion of calcium-pectin cross-links b e,3 , however since most calcium-pectin cross-links are not attached to cell wall microfibrils [18] it is possible that cross-links can move inside the cell wall matrix by a very slow diffusion 3 in Ω e , where D be,3 > 0 and the reaction term g b,3 models the creation and breakage by mechanical stresses of calcium-pectin cross-links, see Remark 2.3 for a detailed form of g b, 3 . For the analysis presented here the diffusion term in the equations for calcium-pectin cross-link density is important. However the same results can be obtained if one assumes that calcium-pectin cross-links do not diffuse and that the reaction terms in equations for pectin, calcium and calcium-pectin cross-links depend on a local average of the deformation gradient, reflecting the fact that in a dense pectin network mechanical forces have a non-local effect on the calcium-pectin chemistry, see [39].
To describe elastic deformations of plant cell walls and middle lamella we consider the equations of poroelasticity reflecting the microscopic structure of cell walls and middle lamella permeable to fluid flow: ρ e ∂ 2 t u e − div(E(b e,3 )e(u e )) + α∇p e = 0 in Ω e , Here u e denotes the displacement from the equilibrium position, e(u e ) stands for the symmetrized gradient of u e , and ρ e denotes the poroelastic wall density. Since we consider the equations of poroelasticity, one more unknown function that should be determined is the pressure, denoted by p e . The mass storativity coefficient is denoted by ρ p and K p denotes the hydraulic conductivity of cell walls and middle lamella. In what follows, we assume that the Biot-Willis constant is α = 1.
It is observed experimentally that the load bearing calcium-pectin cross-links reduce cell wall expansion, see e.g. [51]. Hence elastic properties of cell walls and middle lamella depend on the chemical configuration of pectin and density of calcium-pectin cross-links, see e.g. [55]. This is reflected in the dependence of the elasticity tensor E of the cell wall and middle lamella on the density of calcium-pectin cross-links b e,3 . The differences in the elastic properties of cell walls and middle lamella result in the dependence of the elasticity tensor E on the spatial variables. Since the properties of calcium-pectin cross-links are changing during the deformation and the stretched cross-links have different impact (stress drive hardening) on the elastic properties of the cell wall matrix than newly-created cross-links, see e.g. [11,37,43], we consider a non-local in time dependence of the Young's modulus of the cell wall matrix on the density of calcium-pectin cross-links, see Assumption A1. A similar approach was used in [20] to model the dependence of cell deformations on the concentration of a chemical substance. We assume that the hydraulic conductivity tensor varies between cell wall and middle lamella and, hence, K p depends on the spacial variables.
In the cell inside, that is in Ω f , the water flow is modelled by the Stokes system where ∂ t u f denotes the fluid velocity, p f the fluid pressure, µ the fluid viscosity, and ρ f the fluid density. As transmission conditions between free fluid and poroelastic domains we consider the continuity of normal flux, which corresponds to mass conservation, and the continuity of the normal component of total stress on the interface Γ, i.e. the total stress of the poroelastic medium is balanced by the total stress of the fluid, representing the conservation of momentum, Also taking into account the channel structure of a cell membrane separating cell inside and cell wall, given by the presence of aquaporins, see e.g. [14], we assume that the water flow between the poroelastic cell wall and cell inside is in the direction normal to the interface between the free fluid and the poroelastic medium. Hence we assume the no-slip interface condition, which is appropriate for problems where at the interface the fluid flow in the tangential direction is not allowed, see e.g. [12], By Π τ w we define the tangential projection of a vector w, i.e. Π τ w = w − (w · n)n, where n is a normal vector and τ indicates the tangential subspace to the boundary. The balance of the normal components of the stress in the fluid phase across the interphase is given by Notice that the transmission conditions (1) and (2) imply E(b e,3 ) e(u e ) n · n = 0 on Γ. The transmission conditions are motivated by the models describing coupling between Biot and Navier-Stokes or Stokes equations considered in e.g. [6,12,27,28,46]. The coupling between elastic deformations and fluid flow is described in the Lagrangian configuration and hence Γ is a fixed interface between fluid domain and elastic material. Since in our model we consider only the linear elastic nature of the solid skeleton of the cell walls, the transmission conditions (1) and (2) are the corresponding linearizations of the fluid-solid interface conditions, i.e. | det(I + ∇u e )|(µ e(∂ t u f (t, x + u e )) − p f (t, x + u e )I)(I + ∇u e ) −T n is approximated by (µ e(∂ t u f (t, x)) − p f (t, x)I)n on Γ and the first Piola-Kirchhoff stress tensor is equal to the Cauchy stress tensor in the first order approximation. Then the model for the densities of calcium, pectins and calcium-pectin cross-links reads as in Ω e , t > 0 ∂ t c e = div(D e ∇c e ) + g e (c e , b e , e(u e )) in Ω e , t > 0, stands for the diagonal matrix of diffusion coefficients for b e,1 , b e,2 , and b e,3 .
For elastic deformations of cell walls and middle lamella and fluid flow inside the cells we have a coupled system of Stokes equations and poroelastic (Biot) equations: in Ω e , t > 0, in Ω e , t > 0, For multiscale analysis of the mathematical model (3)-(4) we derive the nondimensionalized equations from the dimensional model by considering t =tt * , . The dimensionless small parameter ε = l/L represents the ratio between the representative size of a plant cell l and considered size of a plant tissue L and reflects the size of the microstructure. For a plant root cell we can consider l = 10µm and L = 1m and, hence, ε is of order 10 −5 . We considerx = L, p = Λε, with Λ = 1MPa, andû = l. For the time scale we taket =μ/(Λε 2 ), which together witĥ µ = 10 −2 Pa·s corresponds approximately to 1.7min. We also considerÊ = Λ,K =x 2 ε/(pt) = l 2 /μ, ρ = (Λt 2 )/x 2 =μ 2 /(Λε 4 L 2 ),ρ p = 1/Λ,D =x 2 /t = l 2 Λ/μ, andR =xε/t = ε 3 LΛ/μ. Hydraulic conductivity K p is of order 10 −9 − 10 −8 m 2 · s −1 · Pa −1 and the minimal value of the elasticity tensor is of order 10MPa [55]. Hence the minimal value of the nondimensionalized elasticity tensor E * is approximately 10 and K * p ∼ 0.01 − 0.1. The parameters in the inflow boundary condition, i.e. in P (b e ), are of order 10 −7 m/s, and withR = 10 −7 m/s we obtain the nondimensional parameters in the boundary condition for b e to be of order 1. Here we assume that ρ j > 0, with j = e, p, f , are fixed. The case when the density ρ e and/or ρ p is of order ε 2 can be analyzed in the same way as the case when ρ e = 0 and ρ p = 0, considered in Section 9.
To describe the microscopic structure of a plant tissue, we assume that cells in the tissue are distributed periodically and have a diameter of the order ε. The stochastic case will be analyzed in a Figure 1: Schematic diagram of the geometry of a plant tissue and unit cell future paper. We consider a unit cell where Y e represents the cell wall and a part of middle lamella, and Y f , with Y f ⊂ Y , defines the inner part of a cell. We denote ∂Y f = Γ and let Γ be an open on Γ regular subset of Γ.
Then the time-independent domains Ω ε f and Ω ε e , representing the reference (Lagrangian) configuration of the intracellular (cell inside) and intercellular (cell walls and middle lamella) spaces, are defined by and defines the boundaries between cell inside and cell walls, Γ ε = ξ∈Ξ ε ε( Γ + ξ), see Figure 1. We shall use the following notation for time-space domains: . Neglecting * we obtain the nondimensonalised microscopic model for plant tissue biomechanics in Ω ε e,T , ∂ t c ε e = div(D e ∇c ε e ) + g e (c ε e , b ε e , e(u ε e )) in Ω ε e,T , and in Ω ε f .
On the external boundaries we prescribe the following force and flux conditions: The elasticity and permeability tensors are defined by Y -periodic functions where E(·, ξ) and K p (x, ·) are Y -periodic for a.a. ξ ∈ R and x ∈ Ω.
We emphasize that the diffusion coefficients D b , D e , and D f in equations (6) are supposed to be constant just for presentation simplicity. The method developed in this paper also applies to the case of non-constant uniformly elliptic diffusion coefficients.
We suppose the following conditions hold.
A3. G is a Lipschitz continuous function on R 3 such that |G(r)| ≤ R, for some R > 0 and all r ∈ R 3 .
Here and in what follows we identify the space of symmetric 3 × 3 matrices with R 6 .
Remark 2.2. Our approach also applies to the case when the initial velocity u 1 f 0 has the form u 1,ε f 0 (x) = U 1 f 0 (x, x/ε), where the vector function U 1 f 0 (x, y) is periodic in y, sufficiently regular, and such that div 3 , and c ε e can be considered in the following form: where µ 1 , µ 2 , r dc , r d , κ > 0, and R b (b ε e,3 ) is a Lipschitz continuous function of calcium-pectin crosslinks density, e.g. R b (b ε e,3 ) = r b b ε e,3 with some constant r b > 0. We assume that the concentration of the enzyme PME is constant and hence methylesterified pectin is de-esterified at a constant rate. The demethylesterified pectin is produced through the de-esterification of methylesterified pectin by PME, demethylesterified pectin can decay and through the interaction between two galacturonic acid groups of pectin chains and a calcium molecule a calcium-pectin cross-link is produced. If a cross-link breaks due to mechanical forces we regain two acid groups of demethylesterified pectin and one calcium molecule. We consider the decay of calcium inside the cells. The positive part of the trace of the elastic stress reflects the fact that extension rather than compression causes the breakage of calcium-pectin cross-links. See [39] for more details on the derivation of a microscopic model for the biomechanics of a plant cell wall.

A priori estimates, existence and uniqueness of a solution of the microscopic problem
We begin by proving the existence of a weak solution of the microscopic model (6)- (8) and uniform in ε a priori estimates. In order to obtain uniform in ε estimates we shall extend H 1 -functions from a perforated domain into the whole domain.
• There exists an extension c ε of c ε from L 2 (0, T ; Here the constant C is independent of ε and Proof Sketch. The assumptions on the geometry of Ω ε e and Ω ε ef and a standard extension operator, see e.g. [1,16], ensure the existence of extensions of b ε e , c ε e , and c ε satisfying estimates (12), (13), and (14), respectively.
Remark. Notice that we have a jump in c ε across Γ. Thus in order to construct an extension of c ε in H 1 (Ω) we have to consider c ε outside a δ-neighborhood of Γ. Also since we would like to have an extension of c ε f from Ω ε f to Ω, we have to consider Γ δ ∩ Y e , see Figure 1.
for δ > 0 sufficiently small Γ δ will satisfy the assumption of lemma.
for the densities we have and forT × Ω ε j , with j = e, f , and the constant C is independent of ε.
Proof. The non-negativity of c ε e , c ε f , and b ε e is justified in the proof of Theorem 3.3 on the existence and uniqueness of a weak solution of the microscopic problem (6)- (8).
To derive the estimates in (15), we first take (∂ t u ε e , p ε e , ∂ t u ε f ) as test functions in (9) and obtain ρ e ∂ t u ε e (s) 2 ) Ω ε e for s ∈ (0, T ]. As was defined just after formula (5), Ω ε j,s := (0, s) × Ω ε j for j = e, f . Using assumptions A1., A2., and A5. yields for s ∈ (0, T ]. Under our standing assumptions A1. on E we have Applying the trace and Korn inequalities [33] and using extension properties of u ε e we obtain Our assumptions A5. on the initial conditions ensure for s ∈ (0, T ]. Then applying the trace and Gronwall inequalities in (18) yields the following estimate where the constant C is independent of ε. Using the Korn inequality [33] for deformation and velocity, together with a scaling argument, we obtain Differentiating all equations in (7) with respect to time t and taking (∂ 2 t u ε e , ∂ t p ε e , ∂ 2 t u ε f ) as test functions in the resulting equations, we obtain for s ∈ (0, T ]. Here we used the following equality Assumptions A5. on the initial conditions together with the microscopic equations in (7) ensure that where the constant C is independent of ε. To justify (24), first we consider the Galerkin approximations of u ε e and ∂ t u ε f and a function φ k in the corresponding finite dimensional subspace, with φ k = 0 on ∂Ω and div φ k = 0 in Ω ε f , Taking t → 0 and using the regularity of u ε,k e , ∂ t u ε,k f , and b ε e,3 with respect to the time variable we obtain Then the integration by parts in the last two terms and the assumptions on the initial values ensure and hence In a similar way we also obtain the boundedness of ∂ t p ε,k e (0) L 2 (Ω ε e ) uniformly in k. Then the estimates similar to (23) for the Galerkin approximations of u ε e , p ε e , and . Then from the equations for u ε e and p ε e and the continuity of e(u ε e ), ∂ t u ε e , and ∇p ε e with respect to the time variable we obtain the continuity of ∂ 2 t u ε e and ∂ t p ε e with respect to the time variable. Then the assumptions on u ε e0 , u 1 e0 , and p ε e0 ensure the boundedness of Considering the continuity of e(u ε e ), ∂ 2 t u ε e , ∇p ε e and e(∂ t u ε f ) with respect to the time variable and taking t → 0 we obtain the continuity of ∂ 2 t u ε f and The integration by parts, the boundary conditions for u ε e0 , and the assumptions on φ imply From the assumptions on u ε e0 and p ε e0 we have that div Then considering assumptions A1.-A2. and applying the Hölder and Gronwall inequalities in (23), we obtain the estimates for ∂ 2 t u ε e , ∂ t p ε e and ∂ 2 t u ε f stated in (15). Here we used the fact that assumptions A1. on E imply the following upper bound: Testing the first and third equations in (7) with φ ∈ L 2 (0, T ; H 1 (Ω)) 3 and using the a priori estimates for u ε e , p ε e , and ∂ t u ε f , we obtain Here we used the properties of an extension of p ε e from Ω ε e to Ω, see Lemma 3.1, and the trace estimate where the constant C is independent of ε. This implies, by the definition of the L 2 -norm and the estimates for p ε e , that p ε To justify estimates (16) we take b ε e and c ε as test functions in (10) and (11), respectively. Using and c ε e (s) 2 . The Gagliardo-Nirenberg and trace inequalities, together with the extension properties of b ε e and c ε , see for an arbitrary δ > 0, and C δ depending on δ and independent of ε. Notice that since the Gagliardo-Nirenberg inequality is applied to the extension of b ε e and c ε defined in Ω, the constant in the Gagliardo-Nirenberg inequality is independent of ε. Then applying the Gronwall inequality and using the assump- The uniform boundedness of b ε e , i.e.
Since the derivation of estimate (28) is rather involved, we present the detailed proof of this estimate in Appendix, see Lemma 10.1. In the same Lemma in Appendix we also prove the estimate To justify the last estimate (17), we integrate the equation for b ε e in (6) over (t, t + h) and consider θ h b ε e − b ε e as a test function: . Then using the a priori estimates for u ε e , b ε e and c ε e in (15) and (16) together with the Hölder inequality implies the estimate for b ε Proof. We shall use a contraction argument to show the existence of a solution of the coupled system. We consider an operator K over L ∞ (0, s; ) in place of (u ε e , ∂ t u ε f ) and with external boundary conditions in (8), and then (u ε,j e , p ε,j e , ∂ t u ε,j f , p ε,j f ) are solutions of (7) with b ε,j e in place of b ε e . For each j = 2, 3, . . ., the proof of existence and uniqueness of (b ε,j e , c ε,j e , c ε,j f ) for given (u ε,j−1 e , ∂ t u ε,j−1 f ) follows the same arguments (with a number of simplifications) as the proof that K is a contraction for (u ε,j e , ∂ t u ε,j f ), i.e. using the Galerkin method and fixed point arguments. Notice that the fixed point argument for the system for b ε,j e and c ε,j allows us to consider the equations for b ε,j e and c ε,j recursively. Thus using the non-negativity of initial data b e0 , c e0 , and c f 0 and assumptions A4. on the reaction and boundary terms and applying iteratively the Theorem on positively invariant regions [40,47] we obtain the non-negativity of all components of b ε,j e and c ε,j . We choose the first iteration (u ε,1 e , p ε,1 e , ∂ t u ε,1 f , p ε,1 f ) to satisfy the initial and boundary conditions in (7) and (8). Then applying the Galerkin method (using the basis functions for H 1 (Ω ε e ) × H 1 (Ω \ Γ ε )) and fixed-point argument we obtain the existence of solutions (b ε,2 e , c ε,2 e , c ε,2 f ) of system (6) with external boundary conditions in (8) where the constant C depends only on e(u ε,1 e ) L ∞ (0,T ;L 2 (Ω ε e )) and the constants in assumptions A1.-A4. The estimates (29) can be justified in the same way as those in (16). Next we consider system (7) with b ε,2 e in place of b ε e . To show the existence result we use the Galerkin method with the basis functions {φ j , ψ j , η j } j∈N for the space and consider the approximate solutions in the form The linearity of equations for u ε,2 e , p ε,2 e , ∂ t u ε,2 f ensures the existence of unique solutions q k j (t) of the corresponding linear system of second order ordinary differential equations with initial conditions q k j (0) = α k j and where α k j and β k j are derived from the initial conditions in (7), and hence, the existence of a unique solution u ε,2 e,k , p ε,2 e,k , ∂ t u ε,2 f,k for k ∈ N. Then using the a priori estimates derived in the same way as in Lemma 3.2 (by considering assumptions A1., A2., and A5.), and taking the limit as k → ∞ we obtain the existence of u ε,2 in Ω ε f }, as test functions in the weak formulation we obtain the equations for u ε,2 e and p ε,2 e in (7) and f . Using first ψ = 0, φ = 0, and η ∈ L 2 (0, T ; H 1 (Ω ε f )) 3 , with Π τ η = 0 on (0, T ) × Γ ε , as a test function in the weak formulation of the equations for u ε,2 e , p ε,2 e , ∂ t u ε,2 f we obtain the transmission condition −n · ε 2 µ e(∂ t u ε,2 f ) n + p ε,2 f = p ε,2 e on (0, T ) × Γ ε , satisfied in the distribution sense. Choosing ψ = 0, φ ∈ L 2 (0, T ; H 1 (Ω ε e )) 3 and η ∈ L 2 (0, T ; H 1 (Ω ε f )) 3 , with φ = η on (0, T ) × Γ ε , as test functions and using the equations for u ε,2 e and ∂ t u ε,2 f ensure (ε 2 µ e(∂ t u ε,2 f ) − p ε,2 f I)n = (E ε (b ε,2 e,3 )e(u ε,2 e ) − p ε,2 e I)n on (0, T ) × Γ ε . Then, using the equations for u ε,2 e , p ε,2 e and ∂ t u ε,2 f and considering ψ ∈ L 2 (0, T ; H 1 (Ω ε e )), φ ∈ L 2 (0, T ; H 1 (Ω ε e )) 3 and η ∈ L 2 (0, T ; H 1 (Ω ε f )) 3 , with Π τ φ = Π τ η on (0, T ) × Γ ε and ψ = 0, φ = 0 on (0, T ) × ∂Ω, as test functions we obtain the transmission condition (−K ε p ∇p ε,2 e + ∂ t u ε,2 e ) · n = ∂ t u ε,2 f · n on (0, T ) × Γ ε in the distribution sense. Taking ψ ∈ L 2 (0, T ; H 1 (Ω ε e )), φ ∈ L 2 (0, T ; H 1 (Ω ε e )) 3 and η ∈ L 2 (0, T ; H 1 (Ω ε f )) 3 , with Π τ φ = Π τ η on (0, T )× Γ ε , as test functions we obtain the boundary conditions on (0, T )× ∂Ω. Hence we obtain that u ε,2 e , p ε,2 e , ∂ t u ε,2 f , p ε,2 f is a weak solution of (7), with b ε,2 e,3 in place of b ε e,3 , together with the corresponding external boundary conditions in (8). Standard arguments pertaining to the consideration of two solutions of (7) imply the uniqueness of a weak solution of (7), (8). The transmission condition −n · ε 2 µ e(∂ t u ε,2 f ) n + p ε,2 f = p ε,2 e on (0, T ) × Γ ε ensures that p ε,2 f is defined uniquely. Also, we obtain that the estimates similar to (15) are valid for the functions u ε,2 e , p ε,2 e , ∂ t u ε,2 f , p ε,2 f uniformly with respect to solutions of equations (6) with boundary conditions in (8): Iterating this step, we conclude the existence of a solution b ε,j e , c ε,j e , c ε,j f of (6) with (u ε,j−1 e , ∂ t u ε,j−1 f ) instead of (u ε e , ∂ t u ε f ) and a solution u ε,j e , p ε,j e , ∂ t u ε,j f , p ε,j f of system (7) with b ε,j e instead of b ε e , and that the estimates similar to (29) and (30)  Then the differences b ε,j e = b ε,j−1 e − b ε,j e , c ε,j e = c ε,j−1 e − c ε,j e , and c ε,j f = c ε,j−1 f − c ε,j f satisfy the following equations: together with the boundary conditions Using b ε,j e , c ε,j e and c ε,j f as test functions in the weak formulation of (31) and (32) we obtain, for any and for any δ > 0. Then combining (33) and (34) and applying the Gronwall inequality we obtain Notice that C 1 = C 2 e C 3 s ≤ C 2 e C 3 T and C δ = C 4 e C 5 s ≤ C 4 e C 5 T for s ∈ (0, T ] and we can consider C 1 and C δ to be independent of s. Considering | b ε,j e | p−1 , with p = 2 k , k = 2, 3, . . ., as a test function in the weak formulation of (31) and (32), applying the Gagliardo-Nirenberg inequality to | b ε,j e | p 2 , and using the iteration in p = 2 k with k ∈ N, see [2, Lemma 3.2], we derive the estimate for s ∈ (0, T ], an arbitrary 0 < δ < 1, and for any 0 < σ < 1/9. For more details see the proof of Lemma 10.2 in Appendix. Notice that C 1 and C δ depend on T and are independent of s. Now letting u ε,j e = u ε,j−1 e − u ε,j e , p ε,j e = p ε,j−1 e − p ε,j e , u ε,j f = u ε,j−1 f − u ε,j f , considering the equations for ( u ε,j e , p ε,j e , ∂ t u ε,j f ) and using (∂ t u ε,j e , p ε,j e , ∂ t u ε,j f ) as test functions in the integral formulation of these equations, we arrive at the following inequality Thus using a priori estimates for u ε,j e , ∂ t u ε,j e , ∂ t u ε,j f , b ε,j e , and b ε,j e we have that u ε,j e L ∞ (0,s;L 2 (Ω ε e )) + e( u ε,j e ) L ∞ (0,s;L 2 (Ω ε e )) + p ε,j e L ∞ (0,s;L 2 (Ω ε e )) for s ∈ (0, T ], 0 < δ < 1, and σ 1 > 10, with the constants C 1 , C 2 , and C δ depending on T and the model parameters but independent of the solutions, initial data, and of s ∈ (0, T ]. Considering δ < 1 and sufficiently small time intervals s in the inequality u ε,j e L ∞ (0,s;L 2 (Ω ε e )) + e( u ε,j e ) L ∞ (0,s;L 2 (Ω ε e )) + ∂ t u ε,j we obtain by the contraction arguments the existence of a fixed point of K and hence the existence of a unique weak solution of the microscopic problem (6)- (8) in (0, s). Since the constants C 2 and C δ depend only on T and the model parameters and do not depend on s, iterating over time intervals we obtain the existence and uniqueness result in the whole time interval (0, T ).

Convergence results
The a priori estimates proved in Lemma 3.2 imply convergence results for the components of solutions of the microscopic problem (6)-(8).
Lemma 4.1. There exist functions u e ∈ H 1 (0, T ; H 1 (Ω)) ∩ H 2 (0, T ; L 2 (Ω)), p e ∈ H 1 (0, T ; H 1 (Ω)), p ε e → p e strongly in L 2 (Ω T ), and for fluid velocity and pressure we have Additionally we have weak two-scale convergence ∂ t u ε e ⇀ ∂ t u e and ∂ t u ε f ⇀ ∂ t u f on Γ ε T . Proof. Applying standard extension arguments, see e.g. [1,16] or Lemma 3.1, and using the same notation for the original and extended sequences, from estimates (15) in Lemma 3.2 we obtain a priori estimates, uniform in ε, for u ε e , ∇u ε e , ∂ t u ε e , ∂ 2 t u ε e and ∇∂ t u ε e , as well as p ε e , ∇p ε e and ∂ t p ε e in L 2 (Ω T ). Then the convergence results for u ε e and p ε e follow directly from the compactness of the embedding of H 1 (0, T ; L 2 (Ω)) ∩ L 2 (0, T ; H 1 (Ω)) in L 2 (Ω T ), the a priori estimates (15), and the compactness theorems for the two-scale convergence, see e.g. [3,31]. The a priori estimates (15) and the compactness theorems for the two-scale convergence ensure the convergence results for ∂ t u ε f and p ε f . Using the trace inequality and a scaling argument together with a priori estimates (15) we obtain where the constant C is independent of ε. Then the compactness theorem for the two-scale convergence on oscillating surfaces [4,30] ensures the weak two-scale convergence of ∂ t u ε e and ∂ t u ε f on Γ ε T . In what follows we shall use the same notation for b ε e , c ε e and their extensions to Ω, whereas the extension of c ε from Ω ε ef to Ω will be denoted by c ε . Then for b ε e and c ε we have the following convergence results.
Proof. Using estimates (16) and the extensions of b ε e , c ε e , and c ε , defined in Lemma 3.1, we obtain where the constant C is independent of ε. The estimates (40), the compactness of the embedding of H 1 (Ω) in L 2 (Ω), along with the estimate (17) and the Kolmogorov compactness theorem [29] yield the strong convergence of b ε e → b e , c ε e → c e and c ε → c in L 2 (Ω T ).
Since Ω ε e,T ∩ Ω ε ef,T = ∅ and c ε e (t, x) = c ε (t, x) in Ω ε e,T ∩ Ω ε ef,T , along with the fact that c e and c are independent of the microscopic variables y, we obtain that c e (t, x) = c(t, x) in Ω T .
From the estimates for c ε , applying the compactness theorem for the two-scale convergence we obtain that there exists c 1 ∈ L 2 (Ω T ; H 1 per (Y \ Γ)/R) such that ∇c ε ⇀ ∇c + ∇ y c 1 weakly two-scale [54].

Derivation of macroscopic equations for the flow velocity and elastic deformations
This section focuses on homogenization of the microscopic problem (7)-(8). First we define the effective tensors E hom , K hom p , and K u . The macroscopic elasticity tensor E hom = (E hom ijkl ), permeability tensor K hom p = (K hom p,ij ), and K u = (K u,ij ) are defined by where w kl = w kl (b e,3 , ·) for k, l = 1, 2, 3, are Y -periodic solutions of the unit cell problems functions w k p = w k p (x, ·), for k = 1, 2, 3, are Y -periodic solutions of the unit cell problems div y K p (x, y)(∇ y w k p + e k ) = 0 in Y e , K p (x, y)(∇ y w k p + e k ) · n = 0 on Γ, Ye w k p dy = 0, and w k e = w k e (x, ·), for k = 1, 2, 3, are Y -periodic solutions of the unit cell problems div y (K p (x, y)∇ y w k e − e k ) = 0 in Y e , (K p (x, y)∇ y w k e − e k ) · n = 0 on Γ, Ye w k e dy = 0.
Here b kl = e k ⊗ e l and {e j } 3 j=1 is the canonical basis of R 3 . Proof Sketch. Assumptions A1. on E and the Korn inequality for periodic functions ensure the existence of a unique solution of the unit cell problems (42) for a given b e,3 ∈ L 2 (Ω T ), see e.g. [33]. Assumptions A2. on K p yield the existence of unique solutions of the unit cell problems (43) and (44). The positive definiteness of E and K p , the definition of E hom and K hom p , and the fact that w kl and w k p , for k, l = 1, 2, 3, are solutions of (42) and (43) ensure in the standard way (see [7]) that E hom and K hom p are positive definite. The definition of E hom implies that E hom satisfies the same symmetry assumptions in A1. as E.
Applying the method of the two-scale convergence and using the convergence results in Lemmas 4.1 and 4.2 we derive the homogenized equations for displacement gradient, pressure and flow velocity for a given {b ε e } such that b ε e → b e strongly in L 2 (Ω T ) 3 as ε → 0. It should be emphasized that we have not yet derived the equation for the limit function b e . We only use the strong convergence of {b ε e }. In the formations of the macroscopic problem for (u e , p e , ∂ t u f ) we shall use the function Q(x, ∂ t u f ) defined as where for (t, x) ∈ Ω T the function q is a Y -periodic solution of the problem Theorem 5.2. A sequence of solutions {u ε e , p ε e , ∂ t u ε f , p ε f } of microscopic problem (7) and (8) converges, as ε → 0, to a solution (u e , p e , ∂ t u f , π f ) of the macroscopic equations with boundary and initial conditions and the two-scale problem for the fluid flow velocity and pressure on Ω T × Γ, where ϑ e = |Y e |/|Y |, ϑ f = |Y f |/|Y |, and with w k p , w k e and q being solutions of (43), (44), and (46), respectively. We have u e ∈ H 2 (0, T ; L 2 (Ω)) ∩ H 1 (0, T ; H 1 (Ω)), p e ∈ H 1 (0, T ; , and π f ∈ L 2 (Ω T × Y f ) and the convergence in the following sense u ε e → u e in H 1 (0, T ; L 2 (Ω)), p ε e → p e in L 2 (Ω T ), ∇u ε e ⇀ ∇u e + ∇ y u 1 e , ∇p ε e ⇀ ∇p e + ∇ y p 1 e weakly two-scale, Remark. In the original microscopic problem the equations of poroelasticity and the Stokes system are coupled through the transmission conditions. The limit system shows the strong coupling in the whole domain Ω T . Namely, the equations for macroscopic displacement and pressure defined in the whole domain Ω T are coupled with the two-scale equations for the fluid flow defined on Ω T × Y f . This coupling in the limit problem can be observed through both the lower order terms in the equations and the boundary conditions. 3 , as test functions in the weak formulation of (7), with the corresponding boundary conditions in (8), we obtain
Thus we have p f = p e in Ω T . Taking (φ(t, x), ψ(t, x), η(t, x, x/ε)), where • φ ∈ C ∞ (Ω T ) 3 and ψ ∈ C ∞ (Ω T ), • η ∈ C ∞ (Ω T ; C ∞ per (Y f )) 3 with Π τ η = Π τ φ on Ω T × Γ and div y η(t, x, y) = 0 in Ω T × Y f , as test functions in the weak formulation of (7), with external boundary conditions in (8), yields Letting ε → 0 and using the two-scale convergence of u ε e , p ε e , and ∂ t u ε f we obtain Here we used the relation p e = p f a.e. in Ω T , as well as the fact that due to the relation div ∂ t u ε f = 0 and the two-scale convergence of ∂ t u ε f , we have To show the convergence of p ε e , η · n Γ ε T we use div y η = 0 and the fact that p 1 e is well-defined on Γ: Notice that n is the internal for Y f normal at the boundary Γ. Also, for an arbitrary test function η 1 ∈ C ∞ 0 (Ω T ; C ∞ 0 (Y f )), from the two-scale convergence of ∂ t u ε f and the fact that ∂ t u ε f is divergence-free it follows that Considering φ ≡ 0 and ψ ≡ 0, and taking first η ∈ C ∞ 0 (Ω T ; C ∞ 0 (Y f )) 3 , with div y η = 0 and then η ∈ C ∞ 0 (Ω T ; C ∞ per (Y f )) 3 with Π τ η = 0 on Ω T × Γ we obtain the two-scale problem (49) for ∂ t u f . From the boundary conditions Π τ ∂ t u ε e = Π τ ∂ t u ε f on Γ ε T and the two-scale convergence of ∂ t u ε e and ∂ t u ε f on Γ ε T , see Lemma 4.1, we obtain x, y) ψ(t, x, y) dγ y dxdt for all ψ ∈ C 0 (Ω T ; C per (Y )). Thus Π τ ∂ t u e = Π τ ∂ t u f on Ω T × Γ. 3 and Π τ η = Π τ φ on Γ, and using the equations (49) for ∂ t u f , we obtain the limit equations for u e and p e : where ϑ e = |Y e |/|Y | and the effective elasticity tensor E hom is defined by (41), and with p 1 e defined by the two-scale problem (55). Considering the weak formulation of (49) with the test function η = 1 yields Using the Y -periodicity of p 1 e we obtain Thus we can rewrite the equation for u e as where ϑ f = |Y f |/|Y |. Considering the structure of problem (55) we represent p 1 e in the form where w k p and w k e are solutions of unit cell problems (43) and (44), and q is a solution of the twoscale problem (46). Incorporating the expression (63) for p 1 e into equations (61) and considering (62) we obtain that p e and u e satisfy the macroscopic problem (47)- (48), where E hom , K hom p , and K u are defined by (41). The coupling with the flow velocity ∂ t u f is reflected in the interaction function Q, defined by (45). Notice that since div ∂ t u f = 0 in Y f , we have that Γ ∂ t u f · n dγ = 0 and the problem (46) is well posed, i.e. the compatibility condition is satisfied.
Proof. To show the strong two-scale convergence we prove the convergence of the energy functional related to equations (7) for u ε e , p ε e , and ∂ t u ε f . Because of the dependence of E on the temporal variable we have to consider a modified form of the energy functional. We consider a monotone decreasing function ζ : R + → R + , e.g. ζ(t) = e −γt for t ∈ R + , and define the energy functional for the microscopic problem (7)-(8) as Ω ε e,s + K ε p ∇p ε e ζ, ∇p ε e ζ Ω ε e,s as a test function in (9) we obtain the equality Due to the assumptions on E and ∂ t E there exists a positive constant γ such that all continuous bounded functions ξ, and y ∈ Y.

Derivation of macroscopic equations for reaction-diffusion-convection problem
The homogenized coefficients in the reaction-diffusion-convection equations, which will be obtained in the derivation of the macroscopic problem, are defined by with ω b and ω being Y -periodic solutions of the unit cell problems and div y (D(y)(∇ y ω j + e j )) = 0 in Y \Γ, D e (∇ y ω j e + e j ) · n = 0, D f (∇ y ω j f + e j ) · n = 0 on Γ, where ω j e (y) = ω j (y) for y ∈ Y e and ω j f (y) = ω j (y) for y ∈ Y f , and z is a Y -periodic solution of Here in Ω T , where W (b e,3 , y) = W klij (b e,3 , y) 3 k,l,i,j=1 = b ij kl + (e y (w ij (b e,3 , y))) kl with w ij being solutions of the unit cell problems (42), and b kl = e k ⊗ e l , {e k } 3 k=1 is the canonical basis of R 3 .
Here ϑ e = |Y e |/|Y |, ϑ f = |Y f |/|Y |, and ϑ Γ = |Γ|/|Y |. We have the convergence in the following sense ∇b ε e ⇀ ∇b e + ∇ y b 1 e , ∇c ε ⇀ ∇c + ∇ y c 1 weakly two-scale. Proof. We can rewrite the microscopic equation for b ε e as , and χ Ω ε e is the characteristic function of Ω ε e . Taking into account the strong convergence of b ε e and c ε e and the two-scale convergence of ∇b ε e and ∇c ε e , see Lemma 4.2, together with the strong two-scale convergence of e(u ε e ) we obtain − |Y e |b e , ∂ t φ 1 Ω T + D b (∇b e + ∇ y b 1 e ), ∇φ 1 + ∇ y φ 2 Ω T ×Ye − |Y e |b e0 , φ 1 Ω T = g b (c, b e , e(u e ) + e y (u 1 e )), Here we used the fact that due to the strong two-scale convergence of e(u ε e ) we have lim ε→0 T * ε (e(u ε e )) − e(u e ) − e y (u 1 e ) L 2 (Ω T ×Ye) = 0, where T * ε is the periodic unfolding operator for the perforated domain Ω ε e , see e.g. [15]. Assumptions on g b in A4. and the a priori estimates for c ε , b ε e , and u ε e ensure where , b e L 2 (Ω T ) and the constants C 1 and C 2 are independent of ε. Combining the estimates in (79), the definition of Ω ε e , and the strong convergence of c ε e and b ε e in L 2 (Ω T ) and of T * ε (e(u ε e )) in L 2 (Ω T × Y e ), along with the properties of the unfolding operator, we obtain for all ψ ∈ C ∞ 0 (Ω T ; C per (Y )). Thus using the estimate for g b (c ε e , b ε e , e(u ε e )) L 2 (Ω ε e,T ) in (79) we conclude g b (c ε e , b ε e , e(u ε e )) → g b (c e , b e , e(u e ) + e y (u 1 e )) two-scale. To show the convergence of the boundary integral over Γ ε we used the Lipschitz continuity of P and the trace estimate Then due to the strong convergence of b ε e in L 2 (Ω T ), the regularity of b e , i.e. b e ∈ L 2 (0, T ; H 1 (Ω)), and the boundedness of ∇b ε e in L 2 (Ω ε e,T ), uniformly in ε, we obtain Taking in (78) first φ 1 ≡ 0 and then φ 2 ≡ 0 and considering such φ 1 that φ 1 (0) = 0 we obtain macroscopic equations for b e in (75). The standard arguments for parabolic equations imply that ∂ t b e ∈ L 2 (0, T ; H 1 (Ω) ′ ). Combining this with the fact that b e ∈ L 2 (0, T ; H 1 (Ω)), see Lemma 4.2, we conclude that b e ∈ C([0, T ]; L 2 (Ω)). Then from (78) we obtain that b e satisfies the initial condition.
The properties of Ω ε f and of the unfolding operator T * ε,f for the domain Ω ε f yield lim ε→0 for all ψ ∈ C ∞ 0 (Ω T ; C per (Y )). Using the Lipschitz continuity of G and the strong convergence of T * ε,f (∂ t u ε f ), ensured by the strong two-scale convergence of ∂ t u ε f , we obtain Thus taking into account the boundedness of G(∂ t u ε f ) we conclude In the same way as for g b , the assumptions in A4. ensure where C = C T * ε (e(u ε e )) L 2 (Ω T ×Ye) , e(u e ) + e y (u 1 . The a priori estimates for c ε , b ε e , and u ε e and assumptions on g in A4. imply g e (c ε e , b ε e , e(u ε e )) L 2 (Ω ε e,T ) ≤ C, with a constant C independent of ε. Then estimate (82) and the strong convergence of c ε e and b ε e in L 2 (Ω T ) and of T * ε (e(u ε e )) in L 2 (Ω T × Y e ), together with calculations similar to (80), yield g e (c ε e , b ε e , e(u ε e )) → g e (c, b e , e(u e ) + e y (u 1 e )) two-scale.
, as a test function in (11) we obtain − c ε e χ Ω ε e , ∂ t ϕ 2 Ω T + D c ∇c ε e , ∇ϕ 2 χ Ω ε e Ω T − g c (c ε e , b ε e , e(u ε e )), , ϕ 2 (∂Ω) T . The two-scale and the strong convergences of c ε e and c ε f together with strong two-scale convergence of e(u ε e ) and ∂ t u ε f ensure that − |Y e |c, ∂ t ψ 1 Ω T + D c (∇c + ∇ y c 1 ), where c 1 l (t, x, y) = c 1 (t, x, y) for y ∈ Y l and (t, x) ∈ Ω T , with l = e, f . Taking into account the structure of equation (83) we represent c 1 in the form x, y) = 8 Well-posedness of the limit problem. Uniqueness of a weak solution To ensure that the whole sequence of solutions of microscopic problem converges we shall prove the uniqueness of a solution of the limit problem (47)-(49), (75). In fact we are going to prove, using the contraction arguments, that the limit problem is well-posed and in particular has a unique solution.
The macroscopic equations for elastic deformation and pressure are coupled with the two-scale problem for fluid flow velocity. Thus the derivation of the estimates for u e and ∂ t u f is non-standard and is shown in the following lemma.
From the compactness of the embedding for any δ > 0. Adding (92) and (93), and applying the Hölder and Gronwall inequalities yield for s ∈ (0, T ) and any 0 < σ < 1/9, where C is independent of s and iterative solutions of the limit problem. Considering a time interval (0,T ), such that CT 9 Incompressible case. Quasi-stationary poroelastic equations in Ω ε e Problem (6)- (8) was derived under a number of assumptions on plant tissue. In some cases these assumptions should be changed and system (6)-(8) should be modified accordingly.
In this section we consider two possible modifications of problem (6)-(8): (i) the incompressible case, when the intercellular space is completely saturated with water and we have the elliptic equation for p ε e ; (ii) the quasi-stationary case for the displacement u ε e . In this case we can consider both compressible and incompressible fluid phase in the elastic part Ω ε e . In the first case the equation for p ε e in (7) is replaced with the following elliptic equation: In the second situation we consider in (7) the quasi-stationary equations for u ε e − div(E ε (b ε e,3 )e(u ε e )) + ∇p ε e = 0 in Ω ε e,T .
The analysis of the quasi-stationary problems considered in this section is very similar to the analysis of (6)- (8) presented in the previous sections. The only part that should be slightly modified is the derivation of a priori estimates.
For the incompressible case, in the same way as in the proof of Lemma 3.2, but now with equation (95) for p ε e , we obtain for s ∈ (0, T ] and arbitrary δ > 0. Then, as in the proof of Lemma 3.2, applying the trace and Korn inequalities [33] and using extension properties of u ε e and assumptions A5. on initial data u ε e0 , u 1 e0 , and u 1 f 0 , we obtain estimates (19), (20), and (22). The trace and Poincare inequalities together with the constraints in (97) and properties of an extension of p ε e from Ω ε e to Ω, see Lemma 3.1, ensure that for s ∈ (0, T ]. Then applying the Gronwall inequality we obtain from (98) the estimates for u ε e , ∂ t u ε e , p ε e and ∂ t u ε f in (21). Differentiating the equations in (7) and (95), with respect to time t and taking (∂ 2 t u ε e , ∂ t p ε e , ∂ 2 t u ε f ) as test functions in the weak formulation of the resulting equations, we obtain 3 )e(u ε e (s)), e(∂ t u ε e (s)) Ω ε e + E ε (b ε e,3 )e(∂ t u ε e (0)) − 2∂ t E ε (b ε e,3 (0))e(u ε e (0)), e(∂ t u ε e (0)) Ω ε e − 2∂ 2 t E ε (b ε e,3 )e(u ε e ) + ∂ t E ε (b ε e,3 )e(∂ t u ε e ), e(∂ t u ε e ) Ω ε e,s for s ∈ (0, T ]. As before, applying the Korn inequality and the Poincare inequality together with the constraint in (97), we obtain the estimates for ∂ 2 t u ε e , ∂ t p ε e , and ∂ 2 t u ε f stated in (15). The equations for ∂ t u ε f and u ε e and calculations similar to those in the proof of Lemma 3.2 ensure the estimate for p ε f . To derive the a priori estimates in the second case, when u ε e satisfies the quasi-stationary equations (96), we have to check that the Korn inequality holds for u ε e .
Lemma 9.1. For u ε e (s) ∈ H 1 (Ω ε e ), with s ∈ (0, T ], we have the following estimate u ε e (s) H 1 (Ω ε e ) ≤ C e(u ε e (s)) L 2 (Ω ε e ) + ε Proof. Consider first Y e and V = {v ∈ H 1 (Y e ) 3 : Considering scaling x = εy and summing over ξ ∈ Ξ ε we obtain v 2 whereΩ ε e = Int ξ∈Ξ ε ε(Y e + ξ) . Using the fact that Π τ ∂ t u ε e = Π τ ∂ t u ε f on Γ ε and estimating u ε e by ∂ t u ε e and the initial value u ε e (0) we obtain Hence applying (103) to u ε e and using the fact that ε u ε e (0) 2 Then considering the extension of u ε e from Ω ε e to Ω, see e.g. [33], and applying the Korn inequality in Ω yield the estimate stated in the Lemma.
Then, in the same way as in the proof of Lemma 3.2, applying the Korn inequalities proved in Lemma 9.1, using extension properties of u ε e and the regularity of the initial data u 1 f 0 ∈ H 2 (Ω) 3 we obtain the following a priori estimates for solutions of the quasi-stationary problem where the constant C is independent of ε. Notice that in the incompressible and quasi-stationary case, i.e. in the case of equations (95) and (96) for p ε e and u ε e , respectively, problem (7), (8), (95), and (96) is well-posed without the initial conditions for u ε e and p ε e . In this case u ε e (0, ·) and ∂ t u ε e (0, ·) are determined from the corresponding elliptic equations and the initial values for the fluid flow u 1 f 0 . In contrast with the limit equations given by (47), in the quasi-stationary and incompressible case the macroscopic equations for effective displacement and pressure do not contain time derivatives and take the form together with the two-scale equations (49) for u f and π f .
where the constant C is independent of ε.
Now we consider (b ε e ) p−1 as a test function in (10)  |b ε e | p 2 2 L 1 (Ω ε e ) dt, where the constants C 1 and C 2 are independent of ε. Then the Alikakos iteration lemma implies the boundedness of b ε e , uniformly in ε. We turn to c ε . Considering first (c ε e,N ) p−1 and (c ε f,N ) p−1 , where c ε j,N (t, x) = min{c ε j (t, x), N } for (t, x) ∈ Ω ε j,T with j = e, f and N > 0, as test functions in (11) and performing calculations similar to those in the derivation of (107), we obtain The boundary integral can be estimated in the same way as in (112) (Ω ε f,s ) ≤ C, where the constant C depends on p and ε and is independent of N . Letting N → ∞ we obtain that (c ε j ) p−1 ∈ L 2 (0, T ; H 1 (Ω ε j )) with j = e, f and p = 3, . . . , 6. Thus we can consider (c ε e ) p−1 and (c ε f ) p−1 , with p = 3, 4, as test functions in (11): In the same way as in (115), applying the Gagliardo-Nirenberg inequality to |c ε j | p 2 in L 2 (Ω ε j ) and L 4 (Ω ε j ) and using properties of the extension of c ε e from Ω ε e to Ω and of c ε from Ω ef to Ω we obtain for s ∈ (0, T ), any δ > 0 and 0 < σ < 1/9, where the constants C and C δ are independent of s and j.