Libo Feng
email: fenglibo2012@126.com

Ian Turner
email: i.turner@qut.edu.au

Patrick Perré
email: patrick.perre@centralesupelec.fr

Kevin Burrage
email: kevin.burrage@qut.edu.au

The use of a time-fractional transport model for performing computational homogenisation of 2D heterogeneous media exhibiting memory effects

Keywords: control volume method, homogenisation theory, two-layered problems, time-fractional derivative, periodic boundary, heterogeneous medium Methodology, Conceptualization, Validation, Formal analysis, Software, Writing -original draft Conceptualization, Methodology, Writing -original draft, Writing -review & editing, Supervision

In this work, a two-dimensional time-fractional subdiffusion model is developed to investigate the underlying transport phenomena evolving in a binary medium comprised of two sub-domains occupied by homogeneous material. We utilise an unstructured mesh control volume method to validate the model against a derived semi-analytical solution for a class of two-layered problems. This generalised transport model is then used to perform computational homogenisation on various two-dimensional heterogenous porous media. A key contribution of our work is to extend the classical homogenisation theory to accommodate the new framework and show that the effective diffusivity tensor can be computed once the cell problems reach steady state at the microscopic scale. We verify the theory for binary media via a series of well-known test problems and then investigate media having inclusions that exhibit a molecular relaxation (memory) effect. Finally, we apply the generalised transport model to estimate the bound water diffusivity tensor on cellular structures obtained from environmental scanning electron microscope (ESEM) images for Spruce wood and Australian hardwood. A highlight of our work is that the computed diffusivity for the heterogeneous media with molecular relaxation is quite different from the classical diffusion cases, being dominated at steady-state by the material with memory effects.

Introduction

Originating from Fourier's and Darcy's law, numerous real-world problems involve transport processes in heterogeneous media, especially multi-scale transport processes, in which the medium properties differ spatially. One important multi-scale problem is the distributed microstructure model or double/dual-porosity model [START_REF] Carr | Dual-scale computational modelling of water flow in unsaturated soils containing irregular-shaped inclusions[END_REF][START_REF] Lewandowska | Modeling of unsaturated water flow in doubleporosity soils by the homogenization approach[END_REF][START_REF] Szymkiewicz | Micromechanical approach to unsaturated water flow in structured geomaterials by dual-scale computations[END_REF], where a large scale domain has a number of small isolated or disconnected inclusions embedded within it.

To model this type of problem, macroscopic and microscopic equations are used to describe the global transport in the connected domain and the local transport in the isolated domains, respectively. Some common approaches for the dual-scale problems include the heterogeneous multi-scale method [START_REF] Carr | Dual-scale computational modelling of water flow in unsaturated soils containing irregular-shaped inclusions[END_REF], distributed microstructure model [START_REF] Lewandowska | Modeling of unsaturated water flow in doubleporosity soils by the homogenization approach[END_REF] and "equation-free" approach [START_REF] Samaey | The gap-tooth scheme for homogenization problems[END_REF]. In [START_REF] Carr | The extended distributed microstructure model for gradient-driven transport: a two-scale model for bypassing effective parameters[END_REF], Carr et al. proposed an extended distributed microstructure model, which can effectively describe the non-equilibrium field in the inclusions for the water flow in unsaturated soils. However, a number of recent studies highlight that anomalous transport phenomena is evident in some complex heterogeneous media [START_REF] Turner | Modelling non-Fickian behavior in the cell walls of wood using a fractionalin-space diffusion equation[END_REF][START_REF] Bueno-Orovio | Fractional diffusion models of cardiac electrical propagation: role of structural heterogeneity in dispersion of repolarization[END_REF][START_REF] Perré | Coupled heat and mass transfer in biosourced porous media without local equilibrium: A macroscopic formulation tailored to computational simulation[END_REF], which deviates from the classical Fickian diffusion. Generally, asymptotic long-time behaviour of the mean square displacement is utilised to characterise anomalous diffusion, which has the form

x 2 (t) ∼ 2K γ Γ(1 + γ) t γ , t → ∞,
where K γ is the generalised diffusion coefficient and γ (0 < γ < 1) is the anomalous diffusion exponent. It turns out that mathematical models involving fractional derivatives are useful tools to treat the anomalous diffusive transport in heterogeneous media [START_REF] Turner | Modelling non-Fickian behavior in the cell walls of wood using a fractionalin-space diffusion equation[END_REF][START_REF] Fomin | Non-Fickian mass transport in fractured porous media[END_REF][START_REF] Metzler | The random walk's guide to anomalous diffusion: a fractional dynamics approach[END_REF], which motivates us to explore these generalised transport equations in this research.

In the seminal paper by Whitaker on the method of volume averaging, macroscopic equations describing the drying process were derived with nonlinear effective parameters within their definitions [START_REF] Whitaker | Coupled transport in multiphase systems: a theory of drying[END_REF]. This theory also included a set of closure equations that relates these parameters to the pore morphology of the medium, which facilitates their prediction. Homogenization theory can also be used to predict these parameters from the microstructure of the porous medium [START_REF] Allaire | Homogenization and two-scale convergence[END_REF][START_REF] Hornung | Homogenization and porous media[END_REF]. The process starts with a full description of the model at the scale of the heterogeneities and upscales the relevant information to derive the macroscopic conservation law by replacing the actual geometry with a periodic idealisation. The macroscopic equation for this homogenised medium is then derived in the limit as the heterogeneities tend to zero. In the classical diffusion model, the effective parameters are defined in terms of the solution of an elliptic equation on the period, subject to periodic boundary conditions (see for example, [START_REF] Carr | Dual-scale computational modelling of water flow in unsaturated soils containing irregular-shaped inclusions[END_REF][START_REF] Lewandowska | Modeling of unsaturated water flow in doubleporosity soils by the homogenization approach[END_REF][START_REF] Szymkiewicz | Micromechanical approach to unsaturated water flow in structured geomaterials by dual-scale computations[END_REF][START_REF] Carr | The extended distributed microstructure model for gradient-driven transport: a two-scale model for bypassing effective parameters[END_REF][START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF]). A key contribution of our work is the derivation of a time-evolutionary homogenisation theory that is applied to predict the diffusivity of porous media that exhibit memory effects. A good example of such a material is wood, in which the cellular structure undergoes a molecular relaxation phenomenon during drying [START_REF] Olek | Implementation of a relaxation equilibrium term in the convective boundary condition for a better representation of the transient bound water diffusion in wood[END_REF][START_REF] Olek | Non-Fickian moisture diffusion in thermally modified beech wood analyzed by the inverse method[END_REF].

To model this phenomenon we propose a generalised transport equation involving a time-fractional operator. A challenge is then the way to deal with the fluxes at the interfaces between two different media when one exhibits memory and the other does not. We show how to correctly model this scenario, paying careful attention to the treatment of the interfacial boundary between the different media.

The types of fractional operators utilised in these models mainly are time-fractional derivatives, which is generally used to describe the algebraic decaying or power-law waiting time in a Lévy process. For the time-fractional derivative, it has a wide application in many problems [START_REF] Yuste | Subdiffusion-limited A + A reactions[END_REF][START_REF] Schumer | Fractal mobile/immobile solute transport[END_REF][START_REF] Langlands | Solution of a modified fractional diffusion equation[END_REF][START_REF] Henry | Fractional cable models for spiny neuronal dendrites[END_REF][START_REF] Chang | Time-space fractional derivative models for CO 2 transport in heterogeneous media[END_REF][START_REF] Feng | Numerical methods and analysis for simulating the flow of a generalized Oldroyd-B fluid between two infinite parallel rigid plates[END_REF][START_REF] Feng | Novel numerical analysis of multi-term time fractional viscoelastic non-Newtonian fluid models for simulating unsteady MHD Couette flow of a generalized Oldroyd-B fluid[END_REF][START_REF] Feng | Finite difference/finite element method for a novel 2D multi-term time fractional mixed sub-diffusion and diffusion-wave equation on convex domains[END_REF]. We start with the following two-dimensional time-fractional subdiffusion equation with variable coefficients ∂u(x, y, t) ∂t

= R 0 D 1-γ t [∇ • (Q∇u(x, y, t))] , (x, y, t) ∈ Ω × (0, T ], (1)
subject to the initial and Neumann boundary conditions given, respectively, as u(x, y, 0) = φ 1 (x, y), (x, y) ∈ Ω, R 0 D 1-γ t [Q∇u(x, y, t)] • n = ψ(x, y, t), (x, y, t) ∈ ∂Ω × (0, T ],

where Ω = (0, L 1 ) × (0, L 2 ), ∇ = [ ∂ ∂x , ∂ ∂y ] T , Q = diag[q 1 (x, y), q 2 (x, y)], ψ(x, y, t) = ® ψ 1 (x, y), n = ± i, ψ 2 (x, y), n = ± j, is the unit vector normal to ∂Ω outward to Ω and 0 D 1-γ t (0 < γ < 1) is the Riemann-Liouville fractional derivative defined by

R 0 D 1-γ t u(x, y, t) = 1 Γ(γ) ∂ ∂t t 0 u(x, y, ξ) (t -ξ) 1-γ dξ.
As an application, we will extend the equation to a time-fractional model for simulating transport in heterogeneous systems, which is motivated by the models in [START_REF] Carr | The extended distributed microstructure model for gradient-driven transport: a two-scale model for bypassing effective parameters[END_REF] and [START_REF] Zeng | A discrete least squares collocation method for two-dimensional nonlinear time-dependent partial differential equations[END_REF]. In [START_REF] Carr | The extended distributed microstructure model for gradient-driven transport: a two-scale model for bypassing effective parameters[END_REF], Carr et al. considered an extended distributed microstructure model for gradient-driven transport on a binary medium without memory. In [START_REF] Zeng | A discrete least squares collocation method for two-dimensional nonlinear time-dependent partial differential equations[END_REF], Zeng et al. investigated a coupled system of Caputo time-fractional partial diffusion equations. Different from the models in [START_REF] Carr | The extended distributed microstructure model for gradient-driven transport: a two-scale model for bypassing effective parameters[END_REF] and [START_REF] Zeng | A discrete least squares collocation method for two-dimensional nonlinear time-dependent partial differential equations[END_REF], this paper deals with a time-fractional system based on the Riemann-Liouville fractional derivative on domains including irregularly shaped inclusions with memory. We derive the semi-analytical solution for a class of two-layered problems to validate the proposed numerical method and present an homogenisation theory and numerical calculation, for two-dimensional periodic structures with memory effects. The main contributions of this paper are summarised as follows:

• A two-dimensional time-fractional subdiffusion equation with variable coefficients is considered, in which an unstructured mesh control volume method (CVM) is applied to solve the problem. Given the non-smoothness of the solution, the modified weighted shifted Grünwald-Letnikov (WSGL) formula with starting weights is utilised to deal with the non-smooth problem, which is further extended to solve the time-fractional transport models in heterogeneous media.

• A novel time-fractional transport model on a binary medium (a porous medium comprised of two distinct phases) with regular or irregular-shaped inclusions is simulated using the CVM. To verify the accuracy of the computational model we derive a semi-analytical solution for a class of two-layered problems with quasiperiodic boundary conditions, in which the finite Fourier and Laplace transforms together with a numerical inverse Laplace transform technique are used. In addition, the mass balance equation for the layered medium is also developed.

• Homogenisation theory is extended from a classical diffusion equation to a time-fractional diffusion equation for carrying out computational homogenisation on some well-known test problems and then on media having inclusions that exhibit a molecular relaxation effect. The generalised transport model is also applied to estimate the bound water diffusivity tensor on cellular structures obtained from ESEM images for Spruce wood and Australian hardwood (Eucalyptus pilularis). An important finding is that, unlike the classical case, the fractional-order indices have a significant effect on the mass transfer and the equivalent diffusivity is dominated by the material with memory effects.

The structure of this paper is as follows. In Section 2, the modified WSGL formulae with starting weights for the time-fractional operator is proposed and the unstructured mesh control volume method for the two-dimensional subdiffusion problem (1) is introduced. We show that the proposed computational scheme offers second order convergence in time and space. In Section 3, a time-fractional transport model on a binary medium is proposed and the semi-analytical solution for a class of two-layered problems is derived. In Section 4, homogenisation theory for the time-fractional diffusion equation is developed. A series of numerical examples are considered in Sections 5 and 6 to verify the homogenisation theory, in which different heterogeneous media morphologies and wood cellular periodic structures are simulated. The numerical results show that the fractional-order indices have a significant effect on the mass transfer for the time-fractional transport equation, which is very different to the classical case. Finally, some conclusions are drawn in Section 7.

The control volume method for the subdiffusion equation

Discretisation for the time-fractional operator

Firstly, we present the grid partition used to discretise equation (1) in the temporal direction. Define t n = nτ , n = 0, 1, 2, . . . , N , where τ = T N is the uniform temporal step. The Riemann-Liouville time-fractional derivative can be rewritten in a more general form as

R 0 D α t v(t) = 1 Γ(1-α) d dt t 0 v(ξ) (t-ξ) α dξ, 0 < α < 1.
To approximate the Riemann-Liouville fractional derivative, the shifted Grünwald-Letnikov difference operator was proposed [START_REF] Meerschaert | Finite difference approximations for fractional advection-dispersion flow equations[END_REF][START_REF] Wang | Compact difference schemes for the modified anomalous fractional sub-diffusion equation and the fractional diffusion-wave equation[END_REF]:

A α τ,p v(t) = τ -α ∞ k=0 g (α) k v(t -(k -p)τ ), where g (α) k = (-1) k α
k and p is an integer. To improve the accuracy of the approximation formula, some high order schemes based on the WSGL formula have been developed [START_REF] Wang | Compact difference schemes for the modified anomalous fractional sub-diffusion equation and the fractional diffusion-wave equation[END_REF][START_REF] Tian | A class of second order difference approximations for solving space fractional diffusion equations[END_REF].

Here we focus on the second order WSGL formula, which has the form

B α τ,p,q v(t) = α -2q 2(p -q) A α τ,p v(t) + 2p -α 2(p -q) A α τ,q v(t),
where p and q are two non equal integers. According to [START_REF] Tian | A class of second order difference approximations for solving space fractional diffusion equations[END_REF], p and q should be chosen satisfying |p| ≤ 1 and |q| ≤ 1.

When (p, q) = (0, -1), we have α-2q

2(p-q) = 2+α 2 , 2p-α 2(p-q) = -α 2 . Then at t = t n , we obtain R 0 D α t v(t n ) = B α,n τ,p,q v + O(τ 2 ) = τ -α n k=0 ω (α) n-k v(t k ) + O(τ 2 ), where B α,n τ,p,q v = 2+α 2 A α τ,p v(t n ) -α 2 A α τ,q v(t n ) and the convolution weights ω (α) k are defined as ω (α) 0 = 2+α 2 g (α) 0 , ω (α) k = 2+α 2 g (α) k -α 2 g (α) k-1 , k ≥ 1. Define D α,n τ,p,q v = 1 2 (B α,n τ,p,q v + B α,n-1 τ,p,q v), then at t = t n-1 2 , we have R 0 D α t v(t n-1 2 ) = D α,n τ,p,q v + O(τ 2 ) = τ -α n k=0 D (α) n-k v(t k ) + O(τ 2 ), (3)
where

D (α) 0 = ω (α) 0 2 , D (α) k = ω (α) k +ω (α) k-1 2 , k ≥ 1.
Remark 2.1. To guarantee the second-order convergence of (3), the required conditions are v(t) ∈ L 1 (R) and R 0 D α+2 t v(t) and its Fourier transform belongs to L 1 (R) [START_REF] Tian | A class of second order difference approximations for solving space fractional diffusion equations[END_REF].

The implementation of the control volume method

To implement the control volume method, we partition the solution domain Ω into a mesh comprised of nonoverlapping triangles. Denote T h as the triangulation, N e the number of triangles and h the maximum diameter of the triangular elements. For a group of triangles with the same vertex, a control volume can be formed by joining the midpoints and barycenters of each triangle. Integrating (1) over each control volume V i (i = 1, 2, . . . , N p ) and applying Gauss's Divergence Theorem, yields

Vi ∂u ∂t dV i = R 0 D 1-γ t ∂Vi (Q∇u) • n dΓ i ,
where dΓ i is an infinitesimal small line segment on ∂V i . Utilising a lumped mass approach to approximate the time derivative, gives

V i ∂u i ∂t = R 0 D 1-γ t ∂Vi (Q∇u) • n dΓ i , i = 1, . . . , N p , (4)
where V i is the volume of the control volume. We assume that the integration path is anticlockwise and the triangle vertices are numbered counter-clockwise. We choose piecewise linear polynomials on the triangles. Within an element triangle e p , p = 1, 2, . . . , N e , the shape function can be defined in the form

N i (x, y) = 1 2∆e p (a i x + b i y + c i ), i = 1, 2, 3
, where a i , b i and c i are some constants and ∆ ep is the area of triangle element p. Combining these shape functions, we can construct the basis function l k (x, y) and approximate u(x, y, t) as u(x, y, t) ≈ u h (x, y, t) = N P k=1 u k l k (x, y). The line integral in (4) can be approximated by the midpoint formula at each control surface:

∂Vi (Q∇u) • n dΓ i = Γi q 1 (x, y, t) ∂u ∂x dy - Γi q 2 (x, y, t) ∂u ∂y dx = mi j=1 2 r=1 Å q 1 (x,
∆x i j,r , (5)
where (x r , y r ) is the mid-point of the control face and m i is the number of sub-control volumes associated with the node i. For more details, the readers can refer to [START_REF] Feng | An unstructured mesh control volume method for two-dimensional space fractional diffusion equations with variable coefficients on convex domains[END_REF]. Then from (4), we can derive the following ODE system

M du dt = R 0 D 1-γ t Ku + F b , (6)
where

M = diag[ V 1 , V 2 , . . . , V Np ], u = [u 1 , u 2 , .
. . , u Np ] T , K is the stiffness matrix derived from (5) and F b is the contribution from the boundary conditions. Now we discuss the treatment of the boundary conditions. For example, if a point i is on the boundary with a control volume face and k and j are its adjacent points on the boundary (see Figure 1); k 0 is the midpoint of k and i and j 0 is the midpoint of j and i; i 1 is the midpoint of j 0 and i and i 2 is the midpoint of k 0 and i. Then according to (2) and ( 5), we obtain

[F b ] i = R 0 D 1-γ t Å q 1 (x, y, t) ∂u ∂x ã (xi 1 ,yi 1 ) (y j0 -y i ) - Å q 2 (x, y, t) ∂u ∂y ã (xi 1 ,yi 1 ) (x j0 -x i ) + R 0 D 1-γ t Å q 1 (x, y, t) ∂u ∂x ã (xi 2 ,yi 2 ) (y i -y k0 ) - Å q 2 (x, y, t) ∂u ∂y ã (xi 2 ,yi 2 ) (x i -x k0 ) = ψ 1 (x i1 , y i1 , t)(y j0 -y i ) -ψ 2 (x i1 , y i1 , t)(x j0 -x i ) + ψ 1 (x i2 , y i2 , t)(y i -y k0 ) -ψ 2 (x i2 , y i2 , t)(x i -x k0 ).
At t = t n-1 2 , applying the Crank-Nicolson scheme to the integer time derivative and the WSGL formula to the time-fractional derivative, we obtain

M u n -u n-1 τ = τ γ-1 n k=1 D (1-γ) n-k K(u k -u 0 ) + F n-1 2 b + t γ-1 n-1 2 Ku 0 Γ(γ) . (7)

Modified numerical scheme with correction terms

Since the WSGL formulae require some restrictions on the solution ( R 0 D α+2 t v(t) ∈ L 1 (R)), they may not exhibit the expected high order accuracy when the solution is not smooth enough, especially if it involves terms of the form t σr . To overcome this issue, the modified WSGL formulae with correction terms were proposed to deal with problems exhibiting a non-smooth solution [START_REF] Zeng | Second-order numerical methods for multi-term fractional differential equations: smooth and non-smooth solutions[END_REF][START_REF] Yin | A class of shifted high-order numerical methods for the fractional mobile/immobile transport equations[END_REF]. At t = t n-1 2 , we have

R 0 D α t v(t n-1 2 ) = τ -α n k=0 D (α) n-k v(t k ) + τ -α m k=1 E (n,α) k v(t k ) + O(τ 2 ), (8)
where the starting weights E

(n,α) k are chosen such that ( 8) is exact for v(t) = t σr (1 ≤ r ≤ m), which leads to the following linear system:

m k=1 E (n,α) k k σr = Γ(σ r + 1) Γ(σ r + 1 -α) Å n - 1 2 ã σr-α - n k=0 D (α) n-k k σr , r = 1, . . . , m.
In addition, we also need to define the starting weights for the first order time derivative at t = t n-1 2 , i.e.,

dv(t n-1 2 ) dt = v(t n ) -v(t n-1 ) τ + τ -1 m k=1 P (n,α) k v(t k ) + O(τ 2 ), (9)
where the starting weights P (n,α) k are chosen such that (9) is exact for v(t) = t σr (1 ≤ r ≤ m), which can be obtained from the following equations:

m k=1 P (n,α) k k σr = σ r Å n - 1 2 ã σr-1 -n σr -(n -1) σr , r = 1, . . . , m.
In this paper, we choose the correction terms σ r = rα, r = 1, 2, . . . , m, such that (m + 1)α ≥ 2. For more details, the readers can refer to [START_REF] Zeng | Second-order numerical methods for multi-term fractional differential equations: smooth and non-smooth solutions[END_REF]. Then the modified numerical scheme of (7) with correction terms is

M u n -u n-1 τ + 1 τ m k=1 P (n,α) k M (u k -u 0 ) =τ γ-1 n k=1 D (1-γ) n-k K(u k -u 0 ) + m k=1 E (n,α) k K(u k -u 0 ) + F n-1 2 b + t γ-1 n-1 2 Ku 0 Γ(γ) , (10)
which also can be extended to solve the time-fractional transport model [START_REF] Whitaker | Coupled transport in multiphase systems: a theory of drying[END_REF].

Investigation of the accuracy and convergence order of CVM with the modified WSGL formulae

We complete this section by assessing the accuracy and convergence order of the CVM method with the modified WSGL formulae for the two-dimensional time-fractional subdiffusion model (1) on a square domain Ω = [0, 1] × [0, 1]. We treat the problem as isotropic Q = 1 2 I with initial condition u(x, y, 0) = sin x sin y and boundary conditions ψ 1 (x, y, t) = 1 2 t γ-1 E γ,γ (-t γ ) cos x sin y, ψ 2 (x, y, t) = 1 2 t γ-1 E γ,γ (-t γ ) sin x cos y, (x, y, t) ∈ ∂Ω × (0, T ],

where the Mittag-Leffler function is

E α,β (z) = ∞ k=0 z k Γ(kα+β)
. The exact solution of this problem is u(x, y, t) = E γ,1 (-t γ ) sin x sin y.

Table 1: The error and convergence order of numerical schemes [START_REF] Bueno-Orovio | Fractional diffusion models of cardiac electrical propagation: role of structural heterogeneity in dispersion of repolarization[END_REF] and [START_REF] Metzler | The random walk's guide to anomalous diffusion: a fractional dynamics approach[END_REF] for problem (1) for different h with γ = 0.5, τ = 1 × 10 -3 at t = 1, in which we choose correction terms m = 3. We apply the numerical scheme [START_REF] Bueno-Orovio | Fractional diffusion models of cardiac electrical propagation: role of structural heterogeneity in dispersion of repolarization[END_REF] and the modified numerical scheme [START_REF] Metzler | The random walk's guide to anomalous diffusion: a fractional dynamics approach[END_REF] to solve (1) and present the error and convergence order for different h and τ with γ = 0.5, τ = 1 × 10 -3 at t = 1 in Tables 1 and2. It is straightforward to observe that the regularity of the exact solution is low for this choice of γ, which can cause difficulty for the numerical method to attain the theoretical temporal convergence order. This finding is evident in the results presented in Tables 1 and2, where we can see that the numerical scheme [START_REF] Bueno-Orovio | Fractional diffusion models of cardiac electrical propagation: role of structural heterogeneity in dispersion of repolarization[END_REF] without correction terms fails to provide good convergence behaviour. However, when correction terms are added, the singularity of the solution can be captured effectively and global second-order convergence is achieved. We conclude that the modified numerical scheme [START_REF] Metzler | The random walk's guide to anomalous diffusion: a fractional dynamics approach[END_REF] with correction terms is a very effective computational method, thus justifying its choice for solving the generalised transport models considered throughout the following sections.

No correction

Application to a heterogeneous medium

A time-fractional transport model

In this section, we extend the control volume method to solve a time-fractional transport model on a binary medium consisting of a connected phase Ω 2 and an inclusion (disconnected phase) Ω 1 . The time-fractional indices are different for the two media (see Figure 2). Here we consider the following system:

® ∂u(x,y,t) ∂t = R 0 D 1-γ1 t [∇ • (Q 1 ∇u(x, y, t))] , (x, y, t) ∈ Ω 1 × (0, T ], ∂v(x,y,t) ∂t = R 0 D 1-γ2 t [∇ • (Q 2 ∇v(x, y, t))] , (x, y, t) ∈ Ω 2 × (0, T ], (11)
subject to

u(x, y, 0) = φ 1 (x, y), (x, y) ∈ Ω 1 ∪ Γ 1 , v(x, y, 0) = φ 2 (x, y), (x, y) ∈ Ω 2 ∪ Γ 1 ∪ Γ 2 , R 0 D 1-γ2 t [Q 2 ∇v(x, y, t)] • n = ψ(x, y, t), (x, y, t) ∈ Γ 2 × (0, T ].
In addition, we also need the following boundary conditions at the interface Γ 1 between the two media

u(x, y, t) = v(x, y, t), R 0 D 1-γ1 t [Q 1 ∇u(x, y, t)] • n = R 0 D 1-γ2 t [Q 2 ∇v(x, y, t)] • n, (x, y, t) ∈ Γ 1 × (0, T ]. ( 12
)
Integrating ( 11) over V i ∈ Ω 1 and V j ∈ Ω 2 , respectively, yields Furthermore, the matrix form can be derived as

® Vi ∂u ∂t dV i = R 0 D 1-γ1 t ∂Vi (Q 1 ∇u) • n dΓ i , Vj ∂v ∂t dV j = R 0 D 1-γ2 t ∂Vj (Q 2 ∇v) • n dΓ j .
® M 1 du dt = R 0 D 1-γ1 t K 1 u + F b1 , M 2 dv dt = R 0 D 1-γ2 t K 2 v + F b2 -F b1 , (13)
where

u = [u 1 , u 12 ] T , v = [v 1 , v 12 ] T , u 1 ∈ Ω 1 , u 12 = v 12 , for (x i , y i ) ∈ Γ 1 (
which is the shared part at the interface Γ 1 ) and v 1 ∈ Ω 2 , F b1 is the contribution from the boundary conditions at the internal interface Γ 1 and F b2 is the contribution from the boundary conditions at the external interface Γ 2 . Denote U = [u 1 , u 12 , v 1 ] T , then the system ( 13) can be recast as

M dU dt = R 0 D 1-γ t KU + F b2 ,
where K is written in terms of K 1 and K 2 by applying the continuity of the field over Γ 1 . The system can be solved using the same technique discussed in the previous section. In this section, we consider the semi-analytical solution for the following 2D time-fractional two-layered problem (see Figure 3).

The semi-analytical solution for a class of two-layered problems

® ∂u(x,y,t) ∂t = R 0 D 1-α1 t [D 1 ∆u(x, y, t)] , (x, y, t) ∈ Ω 1 × (0, T ], ∂v(x,y,t) ∂t = R 0 D 1-α2 t [D 2 ∆v(x, y, t)] , (x, y, t) ∈ Ω 2 × (0, T ], ( 14
)
where

Ω 1 = (l 1 , l 2 ) × (0, L 2 ) and Ω 2 = (l 0 , l 3 ) × (0, L 2 )\Ω 1 .
We assume the medium is homogeneous within each layer, impose zero flux boundary conditions at y = 0 and y = L 2 , and constrain the initial condition in each layer to vary only as a function of x to ensure we reduce the problem to that of studying the following one-dimensional diffusion model in a two-layered problem (see the Appendix for further details).

   ∂u(x,t) ∂t = R 0 D 1-α1 t D 1 ∂ 2 u(x,t) ∂x 2 , (x, t) ∈ (l 1 , l 2 ) × (0, T ], ∂v(x,t) ∂t = R 0 D 1-α2 t D 2 ∂ 2 v(x,t) ∂x 2 , (x, t) ∈ (l 0 , l 1 ) ∪ (l 2 , l 3 ) × (0, T ]. (15)
In the following, we will consider the equivalent form of ( 15), which can be written as

∂X i (x, t) ∂t = R 0 D 1-γi t Å D bi ∂ 2 X i (x, t) ∂x 2 ã , x ∈ (l i-1 , l i ), i = 1, 2, 3, (16)
with the initial conditions X i (x, 0) = X i0 (x) and quasi-periodic boundary conditions of the following form

R 0 D 1-γ1 t Å D b1 ∂X 1 ∂x (l 0 , t) ã = R 0 D 1-γ3 t Å D b3 ∂X 3 ∂x (l 3 , t) ã , X 3 (l 3 , t) = X 1 (l 0 , t) + q 0 , R 0 D 1-γi t Å D b,i ∂X i ∂x (l i , t) ã = R 0 D 1-γi+1 t Å D b,i+1 ∂X i+1 ∂x (l i , t) ã , X i (l i , t) = X i+1 (l i , t), i = 1, 2, where u = X 2 , v = X 1 ∪ X 3 , q 0 is a constant, D b2 = D 1 , D b1 = D b3 = D 2 , γ 2 = α 1 and γ 1 = γ 3 = α 2 .
These boundary conditions have been chosen to validate the time evolutionary homogenisation theory proposed in the proceeding section. The problem ( 16) will be solved using a combination of finite Fourier and Laplace transforms, which is an extension of the semi-analytical approach discussed in [START_REF] Carr | A semi-analytical solution for multilayer diffusion in a composite medium consisting of a large number of layers[END_REF]. First define d i = l i -l i-1 , i = 1, 2, 3 and then apply the finite Fourier transform in each layer by noting that the transform associated with the following Sturm-Liouville system in ϕ i , i = 1, 2, 3:

- d 2 dx 2 ϕ i = λ 2 i ϕ i , dϕ i dx (l i-1 ) = 0 and dϕ i dx (l i ) = 0, is given by ‹ X i (λ i,m , t) := X i , ϕ i,m = li li-1 X i (x, t) ϕ i,m (x) dx, where the eigenvalues λ 2 i,m = m 2 π 2 d 2 i , m ≥ 0. The corresponding eigenfunctions are ϕ i,m (x) = 1 √ di for m = 0 and ϕ i,m (x) = » 2 di cos [λ m (x -l i-1 )] for m > 0.
We introduce the time fractional potential for each layer as

X i = R 0 D 1-γi t X i , i = 1, 2, 3
and then apply the finite Fourier transform to each system within the i th layer to obtain

≠ ∂X i ∂t , ϕ i,m ∑ = D bi ≠ ∂ 2 X i ∂x 2 , ϕ i,m ∑ .
The term on the righthand side is integrated by parts to obtain

D bi ≠ ∂ 2 X i ∂x 2 , ϕ i,m ∑ = -D bi ß≠ - d 2 dx 2 ϕ i,m , X i ∑ + ï - ∂X i ∂x (l i , t)ϕ i,m (l i ) + ∂X i ∂x (l i-1 , t)ϕ i,m (l i-1 ) ò™ . ( 17
)
Next, define the interfacial flux terms as

D b1 ∂X1 ∂x (l 1 , t) = D b2 ∂X2 ∂x (l 1 , t):=v 12 (t), D b1 ∂X1 ∂x (l 0 , t) = D b3 ∂X3 ∂x (l 3 , t):=v 13 (t), D b2 ∂X2 ∂x (l 2 , t) = D b3 ∂X3 ∂x (l 2 , t):=v 23 (t).
Substituting the relevant boundary condition information into [START_REF] Yuste | Subdiffusion-limited A + A reactions[END_REF], the transformed layer equations are then given by

∂ ‹ X 1 ∂t = -D b1 λ 2 1,m ‹ X 1 + [v 12 (t)ϕ 1,m (l 1 ) -v 13 (t)ϕ 1,m (l 0 )] , ∂ ‹ X 2 ∂t = -D b2 λ 2 2,m ‹ X 2 + [v 23 (t)ϕ 2,m (l 2 ) -v 12 (t)ϕ 2,m (l 1 )] , ∂ ‹ X 3 ∂t = -D b3 λ 2 3,m ‹ X 3 + [v 13 (t)ϕ 3,m (l 3 ) -v 23 (t)ϕ 3,m (l 2 )] ,
together with the transformed initial conditions ‹ X i (λ i,m , 0) = X i,0 (x), ϕ i,m , i = 1, 2, 3. We now apply the Laplace transform in time and denote

‹ X i (λ i,m , s) = L ¶ ‹ X i (λ i,m , t) © , i = 1, 2, 3. Since L ¶ R 0 D 1-γi t ‹ X i © = s 1-γi ‹ X i (λ i,m ,
s), the Laplace transformation of the three layer equations given above can be rearranged as

‹ X 1 (λ 1,m , s) = ‹ X 1 (λ 1,m , 0) η 1,m (s) + v12 (s)ϕ 1,m (l 1 ) -v13 (s)ϕ 1,m (l 0 ) η 1,m (s) , (18)
‹ X 2 (λ 2,m , s) = ‹ X 2 (λ 2,m , 0) η 2,m (s) + v23 (s)ϕ 2,m (l 2 ) -v12 (s)ϕ 2,m (l 1 ) η 2,m (s) , (19)
‹ X 3 (λ 3,m , s) = ‹ X 3 (λ 3,m , 0) η 3,m (s) + v13 (s)ϕ 3,m (l 3 ) -v23 (s)ϕ 3,m (l 2 ) η 3,m (s) , (20)
where η i,m (s) = s + D bi s 1-γi λ 2 i,m . In order to determine the three unknown interfacial flux values v12 (s), v13 (s) and v23 (s), we need the boundary conditions at the interfaces:

X 1 (l 1 , s) = X 2 (l 1 , s), X 3 (l 3 , s) = X 1 (l 0 , s) + q0
s and X 2 (l 2 , s) = X 3 (l 2 , s). Noting that the inverse finite Fourier transform gives

X i (x, s) = ∞ m=0 ‹ X i (λ i,m , s)ϕ i,m (x), i = 1,
2, 3 and substituting the expressions ( 18)-( 20) we obtain the following 3 × 3 linear system:

∞ m=0 ‹ X 1 (λ 1,m , s)ϕ 1,m (l 1 ) = ∞ m=0 ‹ X 2 (λ 2,m , s)ϕ 2,m (l 1 ), ∞ m=0 ‹ X 3 (λ 3,m , s)ϕ 1,m (l 3 ) = ∞ m=0 ‹ X 1 (λ 1,m , s)ϕ 2,m (l 0 ) + q 0 s , ∞ m=0 ‹ X 2 (λ 2,m , s)ϕ 2,m (l 2 ) = ∞ m=0 ‹ X 3 (λ 3,m , s)ϕ 3,m (l 2 ), (21)
that can be solved for v12 (s), v13 (s) and v23 (s) at a given value of s using Cramer's rule. Finally, the solutions within each layer can be determined by applying the inverse Laplace transform resolved numerically:

‹ X i (λ i,m , t) = L -1 ‹ X i (λ i,m , s) = 1 2πi Γ e st ‹ X i (λ i,m , s) ds = 1 2πi Γ e z t ‹ X i (λ i,m , z/t) dz ≈ -2 Ñ K/2 k=1 c 2k-1 ‹ X i (λ i,m , z 2k-1 /t) t é ,
where z = st, c 2k-1 and z 2k-1 are the residues and poles of the best (K, K) rational approximation of e z on the negative real line as computed by the Carathéodory-Fejér method. Full details of this numerical approach can be found in Trefethen et al. [START_REF] Trefethen | Talbot quadratures and rational approximations[END_REF]. Then the solution in each layer can be computed using the inverse finite Fourier transform

X i (x, t) = ∞ m=0 ‹ X(λ i,m , t)ϕ i,m (x), i = 1, 2, 3.
We will use this semi-analytical solution to validate the numerical solutions we compute using the finite volume method summarised in Section 2.

Layered medium mass balance equations

We now determine mass balance equations to confirm the accuracy of our simulation results. Define L = 3 j=1 d j and the average value of the variable X as follows:

X = 1 L l3 l0 X(x, t) dx = 1 L 3 i=1 li li-1 X i (x, t) dx = 3 i=1 Å d i L ã X i ,
where the average value of X in the i th layer is given by

X i = 1 di li li-1 X i (x, t) dx.
Using the continuity of the fluxes at the interfaces, we now integrate [START_REF] Olek | Non-Fickian moisture diffusion in thermally modified beech wood analyzed by the inverse method[END_REF] across each layer to obtain:

d X dt = 0, (22)
subject to

X (0) = 1 L l3 l0 X 0 (x) dx = X 0 = 3 i=1 di L X i0
, where the i th layer average initial condition is

X i0 = 1 di li li-1 X i0 (x, 0) dx.
Then we can obtain that the average value of X must remain constant at its initial average for the duration of the simulation, namely

X = X 0 . ( 23
)

Homogenisation theory for time-fractional differential equation

In this section, we will extend the homogenisation theory from classical differential equations to fractional differential equations in heterogeneous media. Let Ω be a bounded domain in R 2 with (sufficiently smooth) boundary ∂Ω. At first, we note that the following two equations are equivalent.

∂u ∂t = R 0 D 1-γ t [∇(D • ∇u)] , (24)
C 0 D γ t u = [∇(D • ∇u)] . (25)
This equivalence can be verified by imposing the operator R 0 D γ-1 t on both sides of [START_REF] Feng | Finite difference/finite element method for a novel 2D multi-term time fractional mixed sub-diffusion and diffusion-wave equation on convex domains[END_REF] and using the relationship between the Riemann-Liouville fractional operator and the Caputo fractional operator [START_REF] Chen | Finite difference approximations for the fractional Fokker-Planck equation[END_REF] to obtain [START_REF] Zeng | A discrete least squares collocation method for two-dimensional nonlinear time-dependent partial differential equations[END_REF]. Due to this equivalence, for ease of mathematical exposition we focus on the homogenisation theory for the fractional differential equation [START_REF] Zeng | A discrete least squares collocation method for two-dimensional nonlinear time-dependent partial differential equations[END_REF] involving the Caputo fractional derivative. We consider the following generalised fractional transport model:

C 0 D γ t u ε (x, t) = ∇ • (D ε (x, t)∇u ε (x, t)), in Ω, t ∈ (0, T ], ( 26
) u ε (x, 0) = f (x), in Ω, u ε (x, t) = 0, on ∂Ω, t ∈ (0, T ],
where Ω = (0, L 1 ) × (0, L 2 ), the Caputo derivative is defined as

C 0 D γ t u(x, t) = 1 Γ(1-γ) t 0 (t -η) -γ ∂u(x,η)
∂η dη, and 0 < γ < 1 is the fixed time fractional derivative index. We assume the diffusivity tensor D ε (x, t) is symmetric with

D ε ij ∈ L ∞ (Ω) and D ε ij (x, t)ξ i ξ j ≥ Cξ 2 i a.e. in Ω, for some constant C > 0, ξ i ∈ R. We define D ε (x, t) := D x ε , t ε p , p > 0.
The choice of p will be discussed below and we assume there is no periodicity in the time variables. We follow the discussion from [START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF] by assuming further conditions on the structure of the functions D ε ij , namely:

1. D ε ij (x, t) := D ij x ε , t ε p , p > 0. 2. D ij (y, τ
) is Y -periodic as a function of y, where the fast variables in space and time are y = x ε and τ = t ε p , respectively.

D

ij (y, τ ) ∈ L ∞ (R 2 y × R τ ).
Assuming the fast variable y and slow variable x are independent as ε → 0, we seek an asymptotic expansion of the solution u ε (x, t) in the form u ε (x, t) = ∞ k=0 ε k u k (x, y, t, τ ), with u k being Y -periodic in y and we will add further constraints on u k regarding the time variable τ below. Before proceeding with the two-scale asymptotic expansion of the solution u ε , we first consider a general function f (t, τ ), with τ = t ε p , p > 0 being the local scale in time. Actioning the Caputo derivative gives

C 0 D γ t f (t, τ ) = 1 Γ(1 -γ) t 0 (t -η) -γ f (η, τ )dη. (27)
Next, the chain rule gives

d dt = ∂ ∂t + 1 ε p ∂ ∂τ , i.e., f (η, τ ) = ∂f ∂t + 1 ε p ∂f ∂τ .
Substituting this expression into [START_REF] Wang | Compact difference schemes for the modified anomalous fractional sub-diffusion equation and the fractional diffusion-wave equation[END_REF], we obtain

C 0 D γ t f (t, τ ) = 1 Γ(1 -γ) t 0 (t -η) -γ Å ∂f (η, ζ) ∂t + 1 ε p ∂f (η, ζ) ∂ζ ã dη = C 0 D γ t f (t, τ ) + 1 Γ(1 -γ)ε p t 0 (t -η) -γ ∂f (η, ζ) ∂ζ dη.
The second term on the RHS can be written as

1 Γ(1 -γ)ε p t ε p 0 (t -ε p ζ) -γ ∂f (η, ζ) ∂ζ ε p dζ = 1 Γ(1 -γ) τ 0 ∂f (η, ζ) ∂ζ dζ (τ -ζ) γ ε pγ ,
where we have set ζ = η ε p and t = ε p τ . Our investigation of different choices for the parameter p identified that

p = 2 γ gives C 0 D γ t f (t, τ ) = C 0 D γ t f (t, τ ) + 1 ε 2 C 0 D γ τ f (t, τ
). We will use this relation below to derive the time-fractional unit cell model.

Next, denote A ε = -∂ ∂xi î D ij (y, τ ) ∂ ∂xj ó
and then expand A ε as:

A ε = ε -2 A 1 + ε -1 A 2 + ε 0 A 3 ,
where

A 1 = -∂ ∂yi î D ij (y, τ ) ∂ ∂yj ó , A 2 = -∂ ∂yi î D ij (y, τ ) ∂ ∂xj ó -∂ ∂xi î D ij (y, τ ) ∂ ∂yj ó , A 3 = -∂ ∂xi î D ij (y, τ ) ∂ ∂xj ó
. Substituting the expressions for u ε and A ε into system [START_REF] Feng | Finite difference/finite element method for a novel 2D multi-term time fractional mixed sub-diffusion and diffusion-wave equation on convex domains[END_REF] gives: LHS:

Å C 0 D γ t + 1 ε 2 C 0 D γ τ ã (u 0 + εu 1 + ε 2 u 2 + . . .) = ε -2C 0 D γ τ u 0 + ε -1C 0 D γ τ u 1 + ε 0 C 0 D γ t u 0 + C 0 D γ t u 2 + . . . RHS: -A ε (u 0 + εu 1 + ε 2 u 2 + . . .) = -ε -2 A 1 u 0 -ε -1 (A 1 u 1 + A 2 u 0 ) -ε 0 (A 1 u 2 + A 2 u 1 + A 3 u 0 ) + . . .
Then by equating powers of ε, we obtain:

O ε -2 : C 0 D γ τ u 0 = -A 1 u 0 , (28)
O ε -1 : C 0 D γ τ u 1 = -(A 1 u 1 + A 2 u 0 ), (29)
O ε 0 : C 0 D γ t u 0 + C 0 D γ τ u 2 = -(A 1 u 2 + A 2 u 1 + A 3 u 0 ). ( 30
)
Due to the complexity introduced by the time-fractional operator and the time dependent periodic boundary conditions at the external and interfacial boundaries, we choose to solve these problems to a sufficiently large enough time in τ whereby the homogenised effective diffusion tensor has converged to a desired level of accuracy.

The simulation results for the case studies carried out in Section 5 indicate that τ ∼ O(10 4 )s is required to achieve a pseudo steady-state in the unit cell, at which point the effective diffusivity converges. This finding suggests that the time evolutionary behaviour of this effective diffusivity must be continuously updated during the simulation of the homogenised macroscopic model, and this is a research direction we will explore in the future. It is at this pseudo steady-state that we assume u 0 = u(x, t) is independent of y and τ and ( 29) can be expressed as

C 0 D γ τ u 1 + A 1 u 1 = ∂D ij ∂y i ∂u ∂x j .
Next, we introduce θ j as the steady state Y -periodic solution of the following time-fractional problem

C 0 D γ τ θ j = -A 1 (θ j + y j ), j = 1, 2.
The solution u 1 can be defined up to an additive constant as

u 1 (x, y, t, τ ) = θ j ∂u ∂x j + u 1 (x, t).
Hence, the unit cell problem becomes to solve: which gives the homogenised equation as

C 0 D γ τ θ j = ∂ ∂y i ï D ij (y, τ ) ∂(θ j + y j ) ∂y j ò , j = 1,
C 0 D γ t u = ∂ ∂x i Å D ij (y, τ ) ∂u ∂x j ã , where D = 1 |Y | Y D(I + J T θ )dy, D = ï D 11 D 12 D 21 D 22 ò , J θ = ñ ∂θ1 ∂y1 ∂θ1 ∂y2 ∂θ2 ∂y1 ∂θ2 ∂y2 ô .
We can express the unit cell problem in terms of the new variable ϕ j = θ j + y j , j = 1, 2, as

C 0 D γ τ ϕ j = ∂ ∂y i ï D ij (y, τ ) ∂ϕ j ∂y j ò . ( 31
)
This variable change affects the boundary conditions:

ϕ 1 (L 1 , y 2 , τ ) = ϕ 1 (0, y 2 , τ ) + L 1 , ϕ 1 (y 1 , L 2 , τ ) = ϕ 1 (y 1 , 0, τ ), ϕ 2 (L 1 , y 2 , τ ) = ϕ 2 (0, y 2 , τ ), ϕ 2 (y 1 , L 2 , τ ) = ϕ 2 (y 1 , 0, τ ) + L 2 .
Also we need the boundary conditions ( 12) at the interface. The initial conditions become ϕ j (y 1 , y 2 , 0) = θ j (y 1 , y 2 , 0)+ y j , j = 1, 2. In this case we obtain the effective diffusivity tensor as

D = 1 |Y | Y DJ T ϕ dy.

Numerical examples

In this section, we present the results obtained from the computational homogenisation simulations for binary media having the three different morphologies shown in Figure 7: Morphology 1 is a two-layered problem with a rectangular inclusion; Morphology 2 is a binary medium with a circular inclusion; and Morphology 3 is a binary medium with an L-shaped inclusion. These particular geometries have been analysed previously for estimating the equivalent permeability obtained from the numerical solution of the unit cell problem by Szymkiewicz [START_REF] Szymkiewicz | Modelling water flow in unsaturated porous media: accounting for nonlinear permeability and material heterogeneity[END_REF]. All of the computations were carried out using MATLAB R2018a on a DELL desktop with the configuration: Intel(R) Core(TM) i7-6700 CPU3.40GHz and RAM 16.0 GB.

= [ 3 8 , 5 8 ] × [0, 1], Ω 2 = [0, 1] × [0, 1]\Ω 1 . (b) Morphology 2 Ω 1 = {(x, y)|(x -0.5) 2 + (y -0.5) 2 < 1 4π }, Ω 2 = [0, 1] × [0, 1]\Ω 1 . (c) Morphology 3 Ω 1 = [ 1 2 , 3 4 ] × [0, 1 2 ] ∪ [ 1 4 , 3 4 ] × [ 1 2 , 3 4 ] , Ω 2 = [0, 1] × [0, 1]\Ω 1 .
To perform the homogenisation, we solve [START_REF] Trefethen | Talbot quadratures and rational approximations[END_REF] with periodic boundary conditions and nonzero initial conditions. The effective diffusivity (32) is computed once steady-state is achieved. For the homogenisation problems considered we take D 12 = D 21 = 0 and denote

D = ï D bx D bxy D byx D by ò . (32)
In order to ensure the correct interfacial boundary conditions are invoked, we revert to the following time-fractional model with the Riemann-Liouville fractional derivative:

® ∂u(x,y,t) ∂t = R 0 D 1-γ1 t [∇ • (D b1 ∇u(x, y, t))] , (x, y, t) ∈ Ω 1 × (0, T ], ∂v(x,y,t) ∂t = R 0 D 1-γ2 t [∇ • (D b2 ∇v(x, y, t))] , (x, y, t) ∈ Ω 2 × (0, T ],
with quasi-periodic boundary conditions

R 0 D 1-γ1 t [-D b1 ∇u(x, 0, t)] • n = R 0 D 1-γ1 t [D b1 ∇u(x, 1, t)] • n, R 0 D 1-γ2 t [-D b2 ∇v(0, y, t)] • n = R 0 D 1-γ2 t [D b2 ∇v(1, y, t)] • n, R 0 D 1-γ2 t [-D b2 ∇v(x, 0, t)] • n = R 0 D 1-γ2 t [D b2 ∇v(x, 1, t)] • n,
and v(1, y, t) = v(0, y, t) + 1, v(x, 1, t) = v(x, 0, t) + 1, and with additional boundary conditions at the interface Γ 1 are u(x, y, t

) = v(x, y, t), R 0 D 1-γ1 t [D b1 ∇u(x, y, t)] • n = R 0 D 1-γ2 t [D b2 ∇v(x, y, t)] • n.
We use the time evolutionary homogenisation theory from Section 4 to estimate the effective parameters for these three different morphologies as comparisons with the values of the diffusivity tensors can be made with the test problems given in [START_REF] Szymkiewicz | Modelling water flow in unsaturated porous media: accounting for nonlinear permeability and material heterogeneity[END_REF]. To calculate the equivalent diffusivity, we advance the numerical scheme in time until the steady-state is reached. We then calculate the equivalent diffusivity D bx in the x-direction and D by in the y-direction using the following formulae obtained by approximating [START_REF] Chen | Finite difference approximations for the fractional Fokker-Planck equation[END_REF] in Section 4 using a midpoint quadrature rule applied to the integral over the unit cell

D bx = i χD b1 ∂u ∂x + (1 -χ)D b2 ∂v ∂x i S i i S i , D by = i Ä χD b1 ∂u ∂y + (1 -χ)D b2 ∂v ∂y ä i S i i S i ,
where S i is the area of the i th triangular element and the indicator function is defined as

χ = ß 1, (x, y) ∈ Ω 1 , 0, (x, y) ∈ Ω 2 .
Throughout this section we will frequently make reference to the harmonic average (K 1 ) and arithmetic average (K 2 ), which can be calculated respectively by

K 1 = 1 D b 1 + 2 D b 2 -1 , K 2 = 1 D b1 + 2 D b2
, where 1 is the volumetric fraction of Ω 1 and 2 is the volumetric fraction of Ω 2 .

• Morphology 1: Rectangular inclusion (layered composite material)

For Morphology 1, we are able to use the derived semi-analytic solution from Section 3.2 to validate the CVM method when simulating time-fractional transport in a layered medium with quasi-periodic boundary conditions. This problem is chosen to simulate the CVM method when applied to the homogenisation theory presented in Section 4. We exhibit the numerical solution at the central line y = 0.5 to compare directly with the semi-analytical solution derived in Section 3.2. The initial condition is taken as constant in each layer u(x, y, 0) = v(x, y, 0) = u 0 . The maximum error between the numerical and semi-analytical solutions is presented in Table 3 for different γ 1 and γ 2 with h = 1.848 × 10 -1 , τ = 1 × 10 -3 , D b1 = 10, D b2 = 1, u 0 = q 0 = 1 at t = 1. We conclude from the results in Table 3 and Figure 5 that the numerical and semi-analytical solutions agree very well, with the maximum error in all cases ≈ 3.2 × 10 -4 .

Table 3: The maximum error between the numerical solution and semi-analytical solution at the central line y = 0.5 for different γ 1 and γ 2 with h = 1.848 Next, we investigate the impact of varying the fractional indices γ 1 and γ 2 on the solution profile at steady-state. At first, we fix γ 2 = 1 to reduce the diffusion on Ω 2 to normal diffusion and then change γ 1 to observe the solution behaviour at t = 10 4 s (see Figure 5a). At this time, the diffusion process within layer 1 (with γ 1 = 1) has reached its steady-state, which corresponds with the solution of the classical two-layered diffusion problem. The diffusion process for the case γ 1 = 0.1 also reaches its steady-state, however the solution profile is totally different to the classical case, being much steeper at the centre point x = 0.5. We conclude that the fractional index γ 1 has a significant impact on the diffusion process for the two-layered media with memory. Another interesting finding is that a larger value of the fractional index γ 1 will delay the diffusion process from reaching steady-state. Furthermore, the final steady-state profiles for all the fractional cases are identical when γ 2 = 1. To investigate this phenomenon further, we now fix γ 2 = 1 and choose γ 1 = 0.5 to observe the solution behaviour from t = 1 s to t = 10 7 s (see Figure 5b). We can see that the diffusion for γ 1 = 0.5 needs t = 10 7 s to reach its steady-state, which is much longer than required for the case γ 1 = 0.1 (t = 10 4 s). We again find that the steady-state solution profile for γ 1 = 0.5 is the same as exhibited in Figure 5a for γ 1 = 0.1. Moreover, when we fix γ 1 = 0.1 to see the impact of γ 2 on the diffusion, the reverse phenomenon to that observed above is apparent (see Figure 5c). It appears that the larger the choice of γ 2 , the quicker it reaches the steady-state solution. In Figure 6, the effect of the other model parameters D b1 , D b2 , q 0 and u 0 on the diffusion process is exhibited. Figure 6a shows the impact of the diffusivity coefficient for the layer exhibiting memory (D b1 ) with fixed D b2 = 1, γ 1 = 0.1, γ 2 = 1, from which we see that the smaller the value of D b1 , the easier it becomes for the process to attain steady state. Figures 6b and6c show the solution profiles for different values of q 0 and u 0 for both the numerical and semi-analytical solutions, which validates the accuracy of the computational model (CVM) and the fact that the overall mass balance remains constant over the simulation (as derived in equation ( 23)). In particular, it is evident in Figure 6b that as the value of q 0 is decreased, the steady-state solution profile becomes flatter. Furthermore, for all choices of u 0 in Figure 6c, the average solution value remains constant at u 0 during the simulation.

× 10 -1 , τ = 1 × 10 -3 , D b 1 = 10, D b 2 = 1, u 0 = q 0 = 1 at t = 1. γ 1 =0.2 γ 1 =0.5 γ 1 =0.8
We now give a comparison between the Riemann-Liouville (RL) and Caputo variants of the model discussed in Section 4:

ß C 0 D γ1 t u(x, y, t) = D b1 ∆u(x, y, t), (x, y, t) ∈ Ω 1 × (0, T ], C 0 D γ2 t v(x, y, t) = D b2 ∆v(x, y, t), (x, y, t) ∈ Ω 2 × (0, T ], (33)
with the initial and periodic boundary conditions for the Caputo model now expressed in terms of the classical flux u(x, y, 0) = v(x, y, 0) = u 0 , D b2 ∂v(0, y, t) ∂x = D b2 ∂v(1, y, t) ∂x , v(1, y, t) = v(0, y, t) + q 0 , and the boundary conditions at the interface Γ 1 between the two media u(x, y, t) = v(x, y, t), D b1 ∂u(x, y, t) ∂x As stated in Section 4, the fractional partial differential equations ( 24) and ( 25) are equivalent, however a most important observation is that the interfacial boundary conditions ( 12) and ( 34) are not. Following a similar approach as outlined in Section 3.2, we have derived the semi-analytical solution for the Caputo model [START_REF] Zeng | Second-order numerical methods for multi-term fractional differential equations: smooth and non-smooth solutions[END_REF] and also exhibited its behaviour along the central line y = 0.5. Figure 7 compares the steady-state solution profiles between the classical and fractional two-layered problems for the two different time-fractional transport models. A key finding is that for the Caputo model [START_REF] Zeng | Second-order numerical methods for multi-term fractional differential equations: smooth and non-smooth solutions[END_REF], there is no observable difference in the steady-state profiles computed for the classical and fractional cases whatever the choice of γ 1 . This means that the impact of the internal layer memory effect is diminished, which may not be a true accord of the physics associated with the anomalous transport phenomena. The main reason for this outcome is due to the incorrect treatment of the interfacial boundary conditions [START_REF] Yin | A class of shifted high-order numerical methods for the fractional mobile/immobile transport equations[END_REF], which also explains why we need to transform the Caputo fractional derivative to the Riemann-Liouville fractional derivative to perform the computational homogenisation. Otherwise, the interfacial boundary conditions ( 12) would need to be imposed on the Caputo model [START_REF] Zeng | Second-order numerical methods for multi-term fractional differential equations: smooth and non-smooth solutions[END_REF], which is more complicated to calculate.

= D b2 ∂v(x, y, t) ∂x , (x, y, t) ∈ Γ 1 × (0, T ]. ( 34
)
We complete this section with a summary of our findings for the homogenised diffusivities computed using the time-fractional homogenisation theory for the layered binary medium. To ease the mathematical exposition, we denote α 1 = 1 -γ 1 and α 2 = 1 -γ 2 . At first we consider the case α 1 = α 2 = 0, which reduces the problem to the classical two-layered problem. Figures 8a and8b show the equivalent diffusivity, the harmonic average and arithmetic average curves of the two-layered medium with different diffusivity ratios

D b 1 D b 2
. We can see that D bx and D by coincide with the harmonic average and arithmetic average, respectively, which is a well known result. When 4 shows the relative error of our numerical results with the theoretical equivalent diffusivity in [START_REF] Szymkiewicz | Modelling water flow in unsaturated porous media: accounting for nonlinear permeability and material heterogeneity[END_REF], which illustrates the effectiveness of the proposed numerical scheme. . An interesting observation from these figures is that for this two-layered composite material, the equivalent diffusivity tends to the diffusivity of the material with memory (i.e. D bx → D b1 = 10 when D b1 = 10 and D bx → D b1 = 0.1 when D b1 = 0.1), which is different to the classical case (D bx = 1.2903 or 0.3077) (see Figures 8a and8b). Another finding is that the fractional-order α 1 has a significant impact on the convergence of the effective diffusivity tensor entries and it delays the diffusion to reach its steady-state (at least t > 10 4 s) significantly compared to the time for the classical case (t < 1s). We also notice that the final equivalent diffusivity D bx tends to the same constant, which is independent of the fractional-order α 1 for a fixed α 2 = 0 (see Figure 9). , where D b 2 = 1 and α 2 = 0.

We believe these special properties for the composite material with memory are due to the impact of the fractional operators. In addition, it appears that the fractional-order α 1 has no influence on the equivalent diffusivity D by , which may be due to the homogeneity in the y direction (see Figures 8a and8b). Now, we fix D b2 = 1 and decrease the value D b1 to observe the evolution of the equivalent diffusivity for the classical and fractional cases (α 1 = 0.9) in Figure 10a using log-log axis scales. It can be seen that the equivalent diffusivity decreases with decreasing D b1 for both cases but they never meet (see Figure 10a). Figures 11a and11b depict the steady-state homogenised solutions for the classical two-layered problem and the time-fractional two-layered problem. An interesting outcome is that a very different diffusion profile is obtained for the case with memory. Again it can be concluded that the fractional-order has a significant effect on the diffusion process. • Morphology 2: Circular inclusion.

At first, we consider the case α 1 = α 2 = 0. Figures 8c and8d illustrate the equivalent diffusivity of the binary medium under different diffusivity ratios

D b 1 D b 2
. We determined that when The equivalent diffusivity for this model is between the harmonic average and arithmetic average, which is different to the two-layered problem considered on binary medium 1. In Table 5, the numerical equivalent diffusivity is compared with the theoretical value given in [START_REF] Szymkiewicz | Modelling water flow in unsaturated porous media: accounting for nonlinear permeability and material heterogeneity[END_REF] using the initial condition u(x, y, 0) = v(x, y, 0) = 1. We can see the numerical equivalent diffusivity agrees very well with the theoretical value.

Next, we observe the case when Ω 1 has memory, i.e., α 1 = 0, α 2 = 0 and the evolution of the equivalent diffusivity D b for different α 1 with increasing time under different diffusivity ratios

D b 1 D b 2
is depicted in Figures 9c and9d, respectively. When = 10 (see Figures 8c and8d). Also, the equivalent diffusivity D b is independent of α 1 with fixed α 2 = 0. , where we fix the parameters D b 2 = 1 and α 2 = 0.

The effects of the fraction-order α 1 on the diffusion are similar to that reported for Morphology 1. Furthermore, we study the effects of the diffusivity ratio

D b 1 D b 2
for fixed D b2 = 1 on the equivalent diffusivity and the related results are displayed in Figure 10b. It can be seen that, with a decreasing diffusivity for the memory part D b1 , the memory effect is weakened and the equivalent diffusivity D b coincides with the value for the classical case when D b1 is small (D b1 ≤ 10 -6 ), which tends to a constant D b → 0.5999. Finally, we give a comparison of the steady-state solution profile for Morphology 2 with and without memory in Figures 11c and11d. It can be observed that the memory of material Ω 1 again has a significant impact on the steady-state solution of the model.

• Morphology 3: L-shaped inclusion.

We first consider the classical α 1 = α 2 = 0 case. Figures 8e and8f illustrate the equivalent diffusivity for the problem on binary medium 3 under different medium diffusivity ratios

D b 1 D b 2
with initial condition is u(x, y, 0) = v(x, y, 0) = 1. The relative error between the numerical equivalent diffusivity and the theoretical diffusivity reported in [START_REF] Szymkiewicz | Modelling water flow in unsaturated porous media: accounting for nonlinear permeability and material heterogeneity[END_REF] is presented in Table 6.

When Ω 1 has memory (α 1 = 0, α 2 = 0), the evolution of the equivalent diffusivity D b for different α 1 with increasing time under different medium diffusivity ratios

D b 1 D b 2
is displayed in Figure 12. We see that the equivalent = 10, the values for the fractional case do not fall within this bound (see Figures 8e and8f). An interesting finding is that the sign of the cross-diffusion component D bxy has changed in the whole iterative process when Fig 12). Finally, the comparison of the steady-state in the x-direction and y-direction between the classical case and fractional case is shown in Figures 11d to 11h, from which some clear differences can be observed in the steady-state solution.

D b 1 D b 2 = 10 (see
Remark 5.1. In this section, we note that the computational homogenisation process involves the calculation of the time-fractional derivative over a long time. In the above examples, we have chosen the space step h ≈ 1 64 and time step τ = T 1000 for different T , for which the computational cost for each simulation ranges from around 10 ∼ 30 minutes on the hardware mentioned at the start of Section 5 we used to perform the simulations. When the time step is small or T is extremely large, the fast fractional linear multistep methods [START_REF] Schadle | Fast and oblivious convolution quadrature[END_REF][START_REF] Guo | Efficient multistep methods for tempered fractional calculus: Algorithms and simulations[END_REF] can be applied, which will be a focus of our future research.

In summary, the equivalent diffusivity for the classical two-layered medium can be bounded by the harmonic average and arithmetic average. However, the equivalent diffusivity for the binary medium with memory effect appears not to be bounded by these limits. For the fractional homogenisation problem, the equivalent diffusivity value tends to a fixed value after a long time iteration. In addition, the equivalent diffusivity deceases and recovers the classical case as the diffusivity coefficient of the memory part is decreased.

Application of the model to estimating the bound water diffusivity of wood

We now use the time-fractional transport model to perform computational homogenisation using different microscale images of Spruce earlywood, Spruce latewood and an Australian hardwood (blackbutt) to estimate the bound water diffusivity when the cellular structure exhibits anomalous diffusion or memory effects. This type of behaviour has been observed experimentally in [START_REF] Turner | Modelling non-Fickian behavior in the cell walls of wood using a fractionalin-space diffusion equation[END_REF][START_REF] Olek | Implementation of a relaxation equilibrium term in the convective boundary condition for a better representation of the transient bound water diffusion in wood[END_REF][START_REF] Olek | Non-Fickian moisture diffusion in thermally modified beech wood analyzed by the inverse method[END_REF][START_REF] Perré | Mass diffusivity of beech determined in unsteady-state using a magnetic suspension balance[END_REF], whereby the lignocellulosic structure of the cell wall induces a memory effect during transient hydric changes as a result of molecular relaxation. The ESEM imaging of these wood species were performed at the Laboratoire Génie des Procédés et Matériaux, Centrale-Supélec in France and are shown in Figure 13. Recent experimental work highlighted that traditional diffusion equations cannot adequately describe the absorption of water in a cell wall [START_REF] Olek | Implementation of a relaxation equilibrium term in the convective boundary condition for a better representation of the transient bound water diffusion in wood[END_REF][START_REF] Perré | Mass diffusivity of beech determined in unsteady-state using a magnetic suspension balance[END_REF]. These effects are exacerbated even further for a modified cell wall occurring due to heat treatment [START_REF] Olek | Non-Fickian moisture diffusion in thermally modified beech wood analyzed by the inverse method[END_REF]. The reduction in diffusivity is thought to be due to the degradation of hemicelluloses. The thermal treatment leads to a change in the mechanism of attaining hygroscopic equilibrium by the modified wood. For these applications, traditional conservation laws prove inadequate for describing the underlying physical processes, while models based on fractional operators succeed [START_REF] Turner | Modelling non-Fickian behavior in the cell walls of wood using a fractionalin-space diffusion equation[END_REF]. Fractional-order derivatives provide excellent alternatives to their classical counterparts for such applications by interpolating between the integer orders of differential equations to capture nonlocal relations in time using powerlaw memory kernels. A major attraction is that the dissipative and dispersive properties observed in experimental data are representative of a variety of different forms of generalised constitutive laws. Diffusion involving molecular relaxation in the complex porous microstructure induces a second time constant in the macroscopic physical phenomena that gives rise to 'sub-diffusive' transport behaviour (see chapter 17 in [START_REF] Klages | Anomalous transport: foundations and applications[END_REF]), which motivates why we have used a fractional-in-time derivative for the generalised conservation law. For the computational homogenisations performed over the given unit cell pore structures, the following modified transient diffusion equation can be used to model the transport of water within the wood cell, which is an extension of the wood drying problem discussed in [START_REF] Carr | A dual-scale modeling approach for drying hygroscopic porous media[END_REF][START_REF] Foy | The Meshfree Finite Volume Method with application to multi-phase porous media models[END_REF]

∂Ψ w ∂t + R 0 D α t [∇ • Q w ] = 0.
Introducing the indicator variable

χ = ß 1, (x, y) ∈ C (g) , 0, (x, y) ∈ C (s) ,
then we can define the fractional indexes and the mass fluxes as

α = (1 -χ)α 1 + χα 2 , Q w = (1 -χ)Q (s) w + χQ (g) w , where Q (s) w = -ρ s D b ∇X m ,Q (g) w = - ρgDv 1-ωv ∇ω v ,
where X m is the cell wall bound water moisture content and ρ v is the vapour density. The conserved quantity is defined as Ψ w = (1 -χ)ρ s X m + χρ v . The quasi-periodic boundary conditions are

X m (L 1 , y) = X m (0, y) + 1 L 1 ∂X B ∂x , 0 < y < L 2 , X m (x, L 2 ) = X m (x, 0) + 1 L 2 ∂X B ∂y , 0 < x < L 1 , R 0 D α t [Q w (0, y)] • n ∂C = R 0 D α t [Q w (L 1 , y)] • n ∂C , 0 < y < L 2 , R 0 D α t [Q w (x, 0)] • n ∂C = R 0 D α t [Q w (x, L 2 )] • n ∂C , 0 < x < L 1 ,
where n ∂C is the unit vector normal to ∂C outward to C, ∂X B ∂x and ∂X B ∂y are the gradients imposed over the cell, respectively. The closure conditions for the problem are P a [START_REF] Olek | Implementation of a relaxation equilibrium term in the convective boundary condition for a better representation of the transient bound water diffusion in wood[END_REF] , where the parameters are calculated using the same formulae in [START_REF] Carr | A dual-scale modeling approach for drying hygroscopic porous media[END_REF] at T = 20℃.

+ P v = P atm , ρ v = PvMv R(T +273.15) , ρ a = PaMa R(T +273.
Since wood has a periodic structure [START_REF] Perré | New modelling approaches to predict wood properties from its cellular structure: image-based representation and meshless methods[END_REF], a representative elementary volume (unit cell), needs to be chosen. We used ESEM images of the real wood cellular structure by extracting the boundary of the pores and then a digital representation of the cell was formed. Furthermore, we use the mesh generator Gmsh [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in pre-and postprocessing facilities[END_REF] to perform the triangulation, which can divide the cell domain [0, L 1 ] × [0, L 2 ] into two sub-sets: the cell lumens (Ω 1 or C (g) ) and the solid phase (Ω 2 or C (s) ). Here, we refer to the x-direction as the radial direction and the y-direction as the tangential direction. We denote the harmonic average as the series value and the arithmetic average as the parallel value to be consistent with [START_REF] Perré | New modelling approaches to predict wood properties from its cellular structure: image-based representation and meshless methods[END_REF]. For the treatment of the periodic boundary conditions, one can refer to [START_REF] Carr | A dual-scale modeling approach for drying hygroscopic porous media[END_REF]. Here we will consider the problem on two different species of wood cells: softwood Spruce earlywood and latewood and Australian hardwood (see Figure 13). Figure 14a shows the dimensionless equivalent diffusivity of Spruce earlywood for two different cases: the solid parts without memory (classical case) and with memory (fractional case). Since the solid cell walls have a low diffusivity of water, the equivalent diffusivity value is very close to the series model. As the cell walls in earlywood are aligned with the radial direction, the value in the radial direction is larger than in the tangential direction. For the classical case, the radial value (D bx ≈ 1.0628×10 -2 ) and tangential value (D by ≈ 1.1517×10 -2 ) are bounded by the series and parallel values, which agrees with the results presented in [START_REF] Foy | The Meshfree Finite Volume Method with application to multi-phase porous media models[END_REF][START_REF] Perré | New modelling approaches to predict wood properties from its cellular structure: image-based representation and meshless methods[END_REF]. While for the fractional case (where the solid part has memory), the radial value (D bx ≈ 1.8437 × 10 -3 ) and tangential value (D by ≈ 1.8534 × 10 -3 ) are not bounded by the series and parallel values, being smaller than the classical case.

In the following discussion, we consider the equivalent diffusivity of the Spruce latewood cell problem. Figure 14b displays the dimensionless equivalent diffusivity of Spruce latewood for two different cases: the solid parts without memory (classical case) and with memory (fractional case). In contrast to Spruce earlywood, the tracheids in latewood are more flattened in the tangential direction, which blocks the mass flux in the radial direction. An inverted anisotropic ratio is observed. For the classical case, the radial value is D bx ≈ 2.7068 × 10 -3 and tangential value is D by ≈ 2.3991 × 10 -3 . For the fractional case (the solid part has memory), the radial value is D bx ≈ 1.9415 × 10 -3 and tangential value is D by ≈ 1.9009 × 10 -3 , both of which are smaller than the classical case. Next, we consider the Australian hardwood case. The dimensionless equivalent diffusivity of hardwood for the classical case and the fractional case is presented in Figure 14c, from which a similar result can be seen. For the classical case, the radial value (D bx ≈ 7.0923 × 10 -3 ) and tangential value (D by ≈ 6.3437 × 10 -3 ) are bounded by the series and parallel values. While for the fractional case (the solid part has memory), the radial value (D bx ≈ 2.4440 × 10 -3 ) and tangential value (D by ≈ 2.3370 × 10 -3 ) are not bounded by the series and parallel values. Finally, Figure 15 compares the evolution of the equivalent diffusivity for the three different species of wood cells, namely Spruce earlywood, Spruce latewood and Australian hardwood. It is can be seen that, in accordance with the findings reporting in Section 5, the effective bound water diffusivities slowly approach the values of the solid phase (that undergoes molecular relaxation) diffusivities as the homogenised cell problem reaches steady state. Again, this process can take considerable time for the fractional model and tends to be much more rapid for classical diffusion.

In summary, for Spruce earlywood without memory, the tangential value is larger than the radial value while an inverted anisotropic ratio is observed for Spruce latewood. For Australian hardwood, the tangential and radial values are between those of Spruce earlywood and Spruce latewood, which is due to the cell structure. For both species of wood, the tangential and radial values are bounded by the series and parallel values. However, the tangential and radial values do not fall within these bounds when the wood is assumed to have a memory effect. In addition, the tangential and radial values for wood with memory are smaller than those for wood without memory, which we believe is due to the anomalous transport behaviour. To conclude, for the homogenisation of the timefractional transport problem, the fractional-order indices have a significant impact on the mass transfer and the equivalent diffusivity appears to be dominated by the lignocellulosic material with memory effects.

Conclusions

In this paper, we considered a two-dimensional time-fractional subdiffusion equation, in which the unstructured mesh control volume method is applied. Our chosen numerical examples showed that the method was stable and effective. As an application, we successfully simulated a time-fractional transport model in a binary medium consisting of regular and irregular inclusions and derived a semi-analytical solution for a class of two-layered problems subjected to quasi-periodic boundary conditions. An important contribution was the extension of the classical homogenisation theory to accommodate the new framework and to show that the effective diffusivity tensor can be computed once the cell problems reach steady state. We found for all morphologies considered that the effective diffusivity slowly converged to the diffusivity of the material in the unit cell with memory, which was a very different finding to that observed for homogenised diffusivity parameters for standard materials.

In order to harness the main findings of our work, we postulate that the time evolutionary behaviour of the effective parameters must be continuously updated during the simulation of the homogenised macroscopic model until steady-state in the unit cell is achieved, at which time the homogenised effective parameter is dominated by the diffusivity of the material exhibiting memory effects. Furthermore, it may be that the fractional index used in the macroscopic time-fractional equation must be averaged over the unit cell. In our future research, we plan to investigate this homogenisation strategy in more detail and verify it using experimental data. We will also use the time-fractional multi-scale problem to simulate more complicated heterogeneous systems and explore the use of the multistep methods to ensure accurate and efficient computation over long simulation times are achieved. 14) can be reduced to the one-dimensional variant [START_REF] Olek | Implementation of a relaxation equilibrium term in the convective boundary condition for a better representation of the transient bound water diffusion in wood[END_REF].

Figure 1 :

 1 Figure 1: An illustration of a point with control volume face on the boundary.

Figure 2 :

 2 Figure 2: An illustration of a binary medium comprised of two sub-domains occupied by homogeneous material. The square domain (medium 2) is connected and the circular inclusion (medium 1) is isolated or disconnected.

Figure 3 :

 3 Figure 3: An illustration of a two-layered medium.

2 ,

 2 subject to periodic boundary conditions in y and initially θ j (y, 0) is specified. Denote by I Y (•) = 1 |Y | Y •dy the average of a quantity in y. At the pseudo steady-state in τ (30) becomes C 0 D γ t u = -I y (A 2 u 1 + A 3 u),

Figure 4 :

 4 Figure 4: Schematics of the three binary media taken from [35], each with volumetric fraction for Ω 1 taken as 0.25. The triangulation and control volume partitions are also shown. (a) Morphology 1 Ω 1= [ 3 8 , 5 8 ] × [0, 1], Ω 2 = [0, 1] × [0, 1]\Ω 1 . (b) Morphology 2 Ω 1 = {(x, y)|(x -0.5) 2 + (y -0.5) 2 < 1 4π }, Ω 2 = [0, 1] × [0, 1]\Ω 1 . (c) Morphology 3 Ω 1 = [ 1 2 , 3 4 ] × [0, 1 2 ] ∪ [ 1 4 , 3 4 ] × [ 1 2 , 3 4 ] , Ω 2 = [0, 1] × [0, 1]\Ω 1 .

γ 2 =0.2 3 .

 23 2198E-04 3.2206E-04 3.2088E-04 γ 2 =0.5 3.2002E-04 3.1897E-04 3.1818E-04 γ 2 =0.8 3.2182E-04 3.1903E-04 3.1749E-04 γ 2 =1.0 3.2632E-04 3.2194E-04 3.1917E-04

  (a) Different γ 1 (b) Different t (c) Different γ 2

Figure 5 :

 5 Figure 5: A comparison between the numerical solution (symbol) and semi-analytical solution (line) for different γ 1 , γ 2 and t with D b 1 = 10, D b 2 = 1 and q 0 = u 0 = 1, in which the other parameters are: (a) γ 2 = 1 at t = 10 4 s; (b) γ 1 = 0.5 and γ 2 = 1; (c) γ 1 = 0.1 at t = 10 4 s.

Figure 6 :

 6 Figure 6: A comparison between the numerical solution (symbol) and semi-analytical solution (line) for different D b 1 , q 0 and u 0 with γ 1 = 0.1 and γ 2 = 1, in which the other parameters are: (a) D b 2 = 1, q 0 = u 0 = 1 at t = 1 s; (b) D b 1 = 10, D b 2 = 1, u 0 = 1 at t = 10 4 s; (c) D b 1 = 10, D b 2 = 1, q 0 = 1 and t from 0 s to 10 7 s.

Figure 7 :

 7 Figure 7: A comparison between the numerical solution (symbol) and semi-analytical solution (line) between the classical and fractional two-layered problems for the different time-fractional transport models, in which the the parameters are γ 1 = 1, γ 2 = 1 for the classical case and γ 1 = 0.1, γ 2 = 1 for the fractional case with D b 1 = 10, D b 2 = 1, q 0 = u 0 = 1 at t = 10 4 .

D b 1 D b 2 = 10 , 2 = 0 . 1 ,

 1210201 D bx = 1.2903, D by = 3.250; and for D b 1 D b D bx = 0.3077, D by = 0.775. Table

b 2 = 0. 1 Figure 8 : 1 D b 2 ,

 21812 Figure 8: The equivalent diffusivity for three morphologies with different diffusivity ratios D b 1 D b 2 , in which we fix D b 2 = 1.

Figures 9a and 9b

  Figures 9a and 9b display the evolution of the equivalent diffusivity D bx for different α 1 with increasing time and α 2 = 0 under different diffusivity ratios D b 1 D b 2

Table 4 :b 2 = 0. 1 Figure 9 :

 4219 Figure 9: The evolution of the equivalent diffusivity D bx using different α 1 for Morphologies 1 and 2 with increasing time under different diffusivity ratios D b 1 D b 2

Figure 10 :Figure 11 :

 1011 Figure 10: The effects of the diffusivity ratio D b 1 D b 2for fixed D b 2 = 1 on the equivalent diffusivity, where we use the parameter α 1 = 0.9 for the fractional case. The equivalent diffusivity is presented using a log-log axis (a) and semi-log axis (b), respectively.

D b 1 D b 2 = 10 , 2 = 0 . 1 ,

 210201 D b = 1.5148 and when D b 1 D b D b = 0.6602.

D b 1 D b 2 = 10 ,

 1210 D b ≈ 4.6006 and when D b 1 D b 2 = 0.1, D b ≈ 0.6399, in which the D b value is not bounded between the harmonic average and arithmetic average when D b 1 D b 2

Figure 12 :

 12 Figure 12: The evolution of the equivalent diffusivity D b for different α 1 for Morphology 3 with increasing time under different medium diffusivity ratios D b 1 D b 2

D b 1 D b 2 = 10 , 2 = 0 . 1 ,

 210201 D bx → 6.5561, D by → 4.5559, D bxy → 0.46969 and when D b 1 D b D bx → 0.4444, D by → 0.6444, D bxy → -0.047007. For this morphology, the equivalent diffusivity values D bx and D by for the classical case all lie between the harmonic average and arithmetic average, while when D b 1 D b 2

Figure 13 :

 13 Figure 13: The anatomical images of two different species of wood cells (100µm × 100µm) and their boundary contour and triangulation, in which the volumetric fraction of solid part (Ω 2 ) is approximately (a) 0.3427, (b) 0.8445, (c) 0.4566, respectively.

Figure 14 :

 14 Figure 14: The equivalent diffusivity for three different species of wood cells. Note that all the values have been scaled through division by Dv.

Figure 15 :

 15 Figure 15: The evolution of the equivalent diffusivity for three different species of wood cells. Note that all the values have been scaled through division by Dv.

2 L2 0 Y 2 n +µ 2 mC 2 n +µ 2 mf 2 L1 0 f 0 L1 0 f

 2022222000 subject to u(x, y, 0) = f (x), ∂u(x, 0, t) ∂y = ∂u(x, L 2 , t) ∂y = 0, with general boundary conditions applied in the x-direction. For the y-variable, the Sturm-Liouville problem for the given homogenous flux boundary conditions leads to the normalised eigenfunctions:Y (µ m , y) cos(µ m y), m > 0, and eigenvalues µ m = mπ L2 , m = 0, 1, 2, . . . Observe that L2 (µ m , y)dy = ß √ L 2 , m = 0, 0, m > 0.According to the initial condition, we suppose u(x, y, t) has the general solutionu(x, y, t) = ∞ m=0 ∞ n=1 C m,n X(λ n , x)Y (µ m , y)e -α(λ )t ,with X(λ n , x) and λ n the normalised eigenfunctions and eigenvalues, respectively, for the Sturm-Liouville problem in the x-variable. Rearranging the sum over m leads tou(x, y, t) = ∞ n=1 C 0,n X(λ n , x) m,n X(λ n , x) 2 L 2 cos(µ m y)e -α(λ (x)X(λ n , x)Y (µ 0 , y)dxdy = L (x)X(λ n , x)dx,and as observed above for m > 0C m,n = L2 (x)X(λ n , x)Y (µ m , y)dxdy = 0.Now we see the solution reduces tou(x, t) = ∞ n=1 ‹ C n X(λ n , x)e -αλ 2 n t , with ‹ C n = L1 0 f (x)X(λ n , x)dx, which recovers the solution of the one-dimensional problem. Using a similar argument, the two-dimensional problem (

Table 2 :

 2 The error and convergence order of numerical schemes (7) and[START_REF] Metzler | The random walk's guide to anomalous diffusion: a fractional dynamics approach[END_REF] for problem (1) for different τ with γ = 0.5, h = 2.3063 × 10 -2 at t = 1, in which we choose correction terms m = 3.

		No correction	Apply correction
	τ	error	Order	error	Order
	1/4 3.5427E-02	-	2.5124E-03	-
	1/8 2.5116E-02	0.50	8.6583E-04	1.54
	1/16 1.7903E-02	0.49	2.7189E-04	1.67
	1/32 1.2720E-02	0.49	7.7000E-05	1.82
	1/64 9.0164E-03	0.50	1.7889E-05	2.11

Table 5 :

 5 The relative error of the numerical results with the theoretical equivalent diffusivity in[START_REF] Szymkiewicz | Modelling water flow in unsaturated porous media: accounting for nonlinear permeability and material heterogeneity[END_REF] for Morphology 2.

Table 6 :

 6 The relative error of the numerical results with the theoretical equivalent diffusivity in[START_REF] Szymkiewicz | Modelling water flow in unsaturated porous media: accounting for nonlinear permeability and material heterogeneity[END_REF] for Morphology 3.

	D b1 D b2	D bx	Theoretical Relative error
	10	1	1.4841	1.4800	2.75E-03
	0.1	1	0.5328	0.5330	3.72E-04
	D b1 D b2	D by	Theoretical Relative error
	10	1	1.8828	1.8800	1.49E-03
	0.1	1	0.6758	0.6750	1.15E-03
	D b1 D b2	D bxy	Theoretical Relative error
	10	1	-7.9601E-02	-0.0796	1.80E-05
	0.1	1	-2.8599E-02	-0.0286	4.51E-05
	diffusivity tends to a fixed value at steady-state, which is consistent with the findings reported for Morphologies
	1 and 2. When

Acknowledgment

This work was supported financially by a visiting Professorial Fellowship that enabled Turner to work at the Université Paris-Saclay, France for a period of three months in 2017-18. We acknowledge the financial support for this research received through the Australian Research Council Discovery Grant DP150103675 and the financial support of collectivities to the Chair of Biotechnology of CentraleSupélec (Département de la Marne, Grand Reims, Région Grand Est and FEDER). All the authors thank the anonymous referees for their insightful comments that helped improve the final version of the paper.

Appendix

In this section, we show under the stated assumptions given in Subsection 3.1, that the two-dimensional problem [START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF] reduces to the one-dimensional problem [START_REF] Olek | Implementation of a relaxation equilibrium term in the convective boundary condition for a better representation of the transient bound water diffusion in wood[END_REF]. Here we take the following problem to illustrate:

Declaration of interests

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

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


Taking conviviality seriously (extended abstract)

Sophie Quinton and Jean-Bernard Stefani Univ. Grenoble Alpes, INRIA, CNRS, Grenoble INP, LIG, 38000 Grenoble, France

In his seminal book, "Tools for conviviality" [START_REF] Illich | Tools for conviviality[END_REF], Ivan Illich laid down a strikingly prescient analysis of the woes that can plague a technological society. He advocated the development of convivial technologies, i.e. technologies consistent with the development of a convivial society. A convivial society is one which would: support "the autonomous and creative intercourse among persons, and the intercourse of persons with their environment" ([13] p. 24); "guarantee for each member the most ample and free access to the tools of the community and limit this freedom only in favour of another member's equal freedom" ([13] p. 25); and "allow all its members the most autonomous action by means of tools least controlled by others" ( [START_REF] Illich | Tools for conviviality[END_REF] p. 33). The digital technologies that are so ubiquitous today bear the hallmarks of non-convivial tools:

1. Biological degradation: the disastrous ecological footprint of digital technologies is undeniable [START_REF] Flipo | La numérisation du monde -Un désastre écologique[END_REF], a prime example of the absence of decoupling between the growth of industrial production and ecosystemic degradation, despite exponential gains in energy efficiency (e.g. Koomey's law [START_REF] Koomey | Implications of historical trends in the electrical efficiency of computing[END_REF]).

2. Radical monopoly: the computerization of all sectors of economy and social life induces a radical monopoly on the handling of information. Digital technology has become a primary infrastructure on which all other critical infrastructures depend [START_REF] Tipper | Critical resilient interdependent infrastructure systems and processes[END_REF], and a mandatory recourse for any information handling problem, subject to ideological, technical, industrial and commercial injunctions [START_REF] Biagini | La tyrannie technologique. Critique de la société numérique[END_REF][START_REF] Morozov | Pour tout résoudre cliquez ici -L'aberration du solutionnisme technologique[END_REF].

3. Polarization: the advent of surveillance capitalism [START_REF] Zuboff | The Age of Surveillance Capitalism: The Fight for a Human Future at the New Frontier of Power[END_REF], the biases and negative feedback loops induced by algorithms [START_REF] Neil | Weapons of Math Destruction[END_REF] (further reinforced with the large scale deployment of machine learning artefacts), are just examples of the unbalances and abuses of power made possible by digital platforms and automated decision making systems.

4.

Over-programming and obsolescence: the ever-increasing sophistication of digital technologies, their incessant evolutions and accelerated obsolescence induce cultural, economic, generational and geographical divides that further exacerbate social inequalities [START_REF] Van Dijk | The Digital Divide[END_REF].

Getting out of the current non-convivial nightmare is a holistic, systemic, and radically democratic effort, as Illich himself recognized. For computer scientists, one way to contribute is to design computer systems and digital infrastructures for conviviality, recognizing that infrastructure design in the Anthropocene is "wickedly complex" [START_REF] Chester | Infrastructure as a wicked complex process[END_REF].

Designing for conviviality is not a new idea in computer science. Since the 1980s, Illich's conviviality has been invoked by HCI researchers in relation to questions of accessibility and inclusivity, in particular in work by Gerhard Fischer [START_REF] Fischer | Constrained design processes: Steps towards convivial computing[END_REF]. But it has only more recently been considered as an analytical framework for the design of technical artefacts [START_REF] Lizarralde | A framework for the integration of the conviviality concept in the design process[END_REF][START_REF] Carrey | Les low-techs comme objet de recherche scientifique. vers une société pérenne, équitable et conviviale[END_REF][START_REF] Voinea | Designing for conviviality[END_REF][START_REF] Couillet | Changer de métrique: vers un traitement du signal convivial ?[END_REF]. As such, designing for conviviality can be considered an instance of Value-sensitive design (VSD) [START_REF] Van Den Hoven | Moral methodology and information technology[END_REF], and related to designing for specific values such as human well-being and sustainability. VSD stems from the consideation that technologies and artefacts cannot be considered ethically neutral, but are inherently ambivalent -as Jacques Ellul made abundantly clear [START_REF] Ellul | Le système technicien[END_REF] and that they have moral and political bearings on humans and their environment. Ethical issues that arise when they get deployed are thus best avoided by making sure that relevant moral values are explictily taken into account early on in their design. VSD has developed widely in the past two decades [START_REF]Handbook of Ethics, Values, and Technological Design[END_REF], by looking at particular moral values such as transparency, justice, democracy, human well-being, privacy, sustainability, and at specific technology areas such as defense, healthcare, biotechnology, ICT. In computer science, it was introduced originally by HCI experts in the 1990s with a strong emphasis on practical methodologies 1 [START_REF] Friedman | Value sensitive design : shaping technology with moral imagination[END_REF]. VSD however has received several criticisms, notably by N. Manders-Huits [START_REF] Manders-Huits | What values in design? the challenge of incorporating moral values into design[END_REF], which highlight its lack of an overarching ethical framework, a form of relativism which seems somewhat contradictory given that VSD for example requires all stakeholders to be involved in the design process of a tool. More generally, along with other approaches such as Ethically-aligned design [START_REF]Ethically Aligned Design -Version II[END_REF] or Responsible Research and Innovation [START_REF] Macnaghten | The Making of Responsible Innovation[END_REF], VSD suffers from a lack of commitment to the "good society" technology should contribute to, and from an insufficient critical stance on technology. The convivial society provides this missing overarching commitment, and Illich's threats to conviviality constitute powerful diagnoses on an ailing technical system.

How then, do we go about designing convivial digital technology? Let us side-step the question for the time being and start by exploring two avenues of inquiry. The first one would question the goals: what is a convivial society? what is a convivial handling of information? The work of Illich questioned major technical systems of the industrialized world: medecine, education, transport. We are missing an Illichean analysis of information and communication technologies, to understand in particular the convivial balance in our field, the thresholds where information speed and automation threaten human creativity, autonomy and control (an example of what Illich calls counterfoil research [START_REF] Illich | Tools for conviviality[END_REF] p. 92). But we need also to question conviviality itself in at least three dimensions: principles, complexity and scale. First, conviviality summons multiple strands of reflection, around issues of justice ("In an age of scientific technology, the convivial structure of tools is a necessity for survival in full justice which is both distributive and participatory." [START_REF] Illich | Tools for conviviality[END_REF] p. 26), dignity, and autonomy. Each of these is the topic of ongoing debates and proposals, in an age of rediscovered planetary limits and humanity's embedding in its biosphere. Situating conviviality and the convivial society in this landscape, in particular in discussions about prospective scenarios, examining its possible role as an attractor for post-growth societies, holds for us a crucial importance. Second, complexity and scale both seem naively antithetic to conviviality. Yet Illich gives the example of the (1970s) telephone system as an eminently convivial tool. Does conviviality lie in unity or simplicity of purpose? But what to make of the industrial scale at which the production of ICT technologies apparently needs to operate? Is conviviality recovered at scale by alternative modes of productions [START_REF] Kostakis | Design global, maufacture local: Exploring the contours of an emerging productive model[END_REF][START_REF] Arvidsson | Changemakers -The industrious future of the digital economy[END_REF]? We need to disentangle this paradox to understand the multi-dimensional aspects of conviviality in digital technologies.

The second avenue focuses on a redesign of digital infrastructures as commons. Commons, as a form of collective action for the autonomous handling of shared interests, goods or resources [START_REF] Dardot | Commun -Essai sur la révolution au XXIe siècle[END_REF], are an integral part of Illich's vision of a convivial society, as discussed in [START_REF] De | Ivan illich's radical thought and the convivial solution to the ecological crisis[END_REF]. They are a necessary alternative to what Illich calls "managerial fascism" to get us out of our ecosystemic crisis ("The alternative to managerial fascism is a political process by which people decide how much of any scarce resource is the most any member of society can claim; a process in which they agree to keep limits relatively stationary over a long time, and by which they set a premium on the constant search for new ways to have an ever larger percentage of the population join in doing ever more with ever less." [START_REF] Illich | Tools for conviviality[END_REF] p. 116). They allow collective virtues of sobriety, modesty and temperance to exert their homeostatic influence. We believe that designing digital infrastructures as commons raises not only institutional concerns but also hardware and software architectural challenges. C. Voinea in discussing designing for conviviality highlights several design criteria, organized around the "axiological pillars identified by Ivan Illich, namely personal autonomy and social cohesion" [START_REF] Voinea | Designing for conviviality[END_REF]: flexibility, transparency, simplifiability and usability along the first one, sharedness, creativity and sociality along the second one. These constitute a useful initial set but other criteria need to be considered such as dependability, maintanability, compositionality, and accountability, which should take priority other performance, the overarching criterion in current computer system design. Architecturally, the commons imperative suggests designing digital infrastructures as federated systems, operating at different scales with communities at each scale in full administrative authority, in stark contrast to the administratively centralized character of current cloud systems -an effort at "disintermediation" as advocated by [START_REF] B Raghavan | Means and ends in human-computer interaction: Sustainability through disintermediation[END_REF].

Illich's ideas around convivial tools seem nowadays to resonate well within the computer science and engineering community. Is it because the nonconvivial nature of the tools we built is so obvious? Is it because performance increase, the overarching design criterion in our community, is so well-aligned with a productivist worldview which seems at odds [START_REF] Parrique | Decoupling debunked: Evidence and arguments against green growth as a sole strategy for sustainability[END_REF] with any hope of respecting a safe operating space for humanity [START_REF] Rockström | A safe operating space for humanity[END_REF]? In any case, taking Illich's conviviality seriously could start useful discussions and collaborations between computing and science and technology scholars.