Poroelastic model parameter identification using artificial neural networks: on the effects of heterogeneous porosity and solid matrix Poisson ratio

Predictive analysis of poroelastic materials typically require expensive and time-consuming multiscale and multiphysics approaches, which demand either several simplifications or costly experimental tests for model parameter identification.This problem motivates us to develop a more efficient approach to address complex problems with an acceptable computational cost. In particular, we employ artificial neural network (ANN) for reliable and fast computation of poroelastic model parameters. Based on the strong-form governing equations for the poroelastic problem derived from asymptotic homogenisation, the weighted residuals formulation of the cell problem is obtained. Approximate solution of the resulting linear variational boundary value problem is achieved by means of the finite element method. The advantages and downsides of macroscale properties identification via asymptotic homogenisation and the application of ANN to overcome parameter characterisation challenges caused by the costly solution of cell problems are presented. Numerical examples, in this study, include spatially dependent porosity and solid matrix Poisson ratio for a generic model problem, application in tumour modelling, and utilisation in soil mechanics context which demonstrate the feasibility of the presented framework.


Introduction
A poroelastic medium consists of a solid porous material (solid matrix) interacting with fluid percolating its pores. Having a pore scale much smaller than the average size of the body, even in the most idealised cases, it is computationally expensive to work with a model including all geometrical details and solve the solid-fluid interaction problem via Direct Numerical Simulation (DNS). This issue imposes the need to reduce the complexity/order of the model, which can be carried out employing multiscale approaches. Phenomenological methods either consider only B Hamidreza Dehghani hamidreza.dehghani@uni.lu Andreas Zilian andreas.zilian@uni.lu 1 Institute of Computational Engineering and Sciences, Department of Engineering, Faculty of Science, Technology and Medicine, University of Luxembourg, 6 Avenue de la Fonte, 4364 Esch-sur-Alzette, Luxembourg the macroscale behaviour (e.g. Biot's theory of poroelasticity [1]) or include upscaling techniques such as volume averaging [2,3]. They provide the same macroscale governing equations that include some characteristic parameters. However, direct relationships between macro-and microscale properties are not considered.
A detailed comparison between different upscaling approaches is provided in [4]. The average-field theory is based on volume average of microscale strains and stresses which results in the effective properties that characterise macroscale response of the medium. This approach can be extended to a nonlinear constitutive behaviour of the subphases, however, the coefficients are typically not entirely related to the underlying microstructure [5][6][7].
In contrast, asymptotic homogenisation [8][9][10], which results in the same macroscale equations as the mentioned phenomenological methods, provides analytical relationships between their coefficients and underlying microscale properties, allowing for much simpler material characterisation tests. The calculation of the mentioned mathematical relationships has recently been carried out for a wide range of microscale properties in three-dimensions via DNS in [5], which demonstrates the viability of this method in a wide range of applications. This framework is then extended in [10,11] by solving the actual macroscale system of PDEs via Finite Element Method to visualise the role of Microscale Solid Matrix Compressibility (MSMC). However, it is difficult to apply the above-mentioned framework in the latter studies to more complex problems, e.g. spatially dependent properties, non-linear problems, considering growth and remodelling [9], etc. The reason is that the macroscale coefficients should be calculated several times (because the microscale properties can be varied at each coordinate and time increment), which is a time-consuming process. This motivates us to find an improved alternative to such upscaling techniques.
Data science, together with powerful modern computers, have provided us with the "fourth paradigm" of science, namely, data-driven modelling and computing that its integration into mechanics is becoming increasingly successful. There are several Machine Learning (ML) approaches that are applicable within the field of mechanics and can be found in several descriptive pieces of literature such as [12][13][14]. Here, we make use of the well-known method Artificial Neural Network (ANN) [15] as it is a powerful tool based on interpolation and optimisation techniques aiming at describing the input-output relationship of a provided dataset. This method is inspired by the architecture of the human brain and consists of input and output and some hidden layers of neurons between them. Each layer adopts an arbitrary number (to be determined by sensitivity analysis) of neurons that are connected to the ones in the adjacent layers via some weights (different from weighted residual method) and biases. An activation function is also applied to each neuron to control if a neuron should be active or not so that a highly accurate approximation of a complex non-linear relationship is possible using simple linear or non-linear piece-wise hierarchical interpolations.
In mechanics, Machine Learning (ML) approaches are exploited to replace either the whole or a part (e.g. constitutive equations) of DNS to make the analysis more efficient and accurate (or one of them) [16][17][18][19][20]. Recently, the capabilities of this paradigm were applied to the multiscale analysis of heterogeneous materials in 2D [21][22][23] and 3D [24] [25]. Nevertheless, in the field of poroelasticity, to the best of our knowledge, the potential of ML is not well exploited. This technique is used to estimate some of the poroelastic material parameters based on their response under certain loading conditions via ultrasound tomography such as [26]. In a recent study [27], a Convolutional Neural Network (CNN) is used to bypassing the numerical homogenisation procedure. However, due to the unbounded randomness of this method, and as CNN is based on interpolation, there is the risk of extrapolation. Moreover, Biot's coefficient and modulus, which are essential properties of Biot's theory of poroelasticity, and have shown to be highly dependent on porosity and MSMC [5], are assumed to be constant.
In contrast, in asymptotic homogenisation, the unit pressure and dynamic viscosity are non-dimensionalised. Thus, the inputs of the problem could be only porosity and Solid Matrix Poisson Ratio (SMPR) which are non-dimensional and bounded. Moreover, in this technique, the effective properties are determined by solving specified problems derived from a robust analytical upscaling process and are based on microscopic properties.
In this study, we aim at bypassing the process of macroscale properties identification derived from asymptotic homogenisation by providing an ANN that efficiently computes the coefficients of the macroscale system of PDEs without solving the cell problems. For this purpose, first, we need to create a sufficiently rich dataset mapping (microscale properties) as the input to the model coefficients at the macroscale (as the output data), which is constructed by solving a certain number of the cell problems derived from asymptotic homogenisation via DNS. This dataset should cover a wide range of Poisson ratios and porosities so that the machine-learnt model be applicable in a wide range of scenarios of interest thus reducing the risk of extrapolation. Next, we choose an appropriate ANN architecture, activation function, and learning rule from sensitivity analysis, considering that the feed-forward process must be fast enough to be called several times. However, as the training procedure itself is a one-time calculation and its efficiency here is not of great importance, we do not investigate more complex and advanced ML techniques such as Long short-term memory (LSTM) for exploding and vanishing gradient problems [28]. The trained network, here, is introduced as a further constitutive relationship to the macroscale system of PDEs and its Finite Element (FE) discretisation via the open-source FE package FEniCS [29] so that, the parameters required to characterise the macroscale behaviour are microscale properties such as porosity and MSMC instead of coefficients computed with a homogenisation technique.
The provided framework is a general one which can be employed in a wide range of applications from soil mechanics to the modelling of biological tissues by choosing appropriate characteristic properties specific to the scenario of interest. This paper focuses on the modelling of poroelastic bodies under different Boundary Conditions (BCs) with spatially dependant porosity and SMPR. We study three types of problems. One with heterogeneous porosity under mechanical pressure and other Neumann and Dirichlet BCs. A second problem class, with similar BCs but different geometry, includes an impermeable object inside a poroelastic one with spatially dependent pore fluid volume fraction. The third problem type with heterogeneous porosity, as well as SMPR, under no mechanical pressure but a con-stant cavity/pore pressure (resembling drug injection into a tumour mass) is subject to specific displacement and drainage BCs. The presented results demonstrate the feasibility of the ANN-supported framework for such complex problems. Furthermore, the importance of flexibly considering additional aspects in poroelastic problems is highlighted. In the present discussion an isotropic solid matrix and rotational invariance of the cell geometry with respect to permutation of the three orthogonal axes are assumed.
The remainder of the work is organised as follows. In the next section, we introduce the governing equations derived from asymptotic homogenisation together with the variational formulation of the cell problems. We discuss the advantages and downsides of the described upscaling process together with the application of ANN to bypass the process of solving the cell problems derived from homogenisation in Sect. 3. This is followed by providing a training dataset, mentioning the ANN features and architecture, and testing the trained network. In Sect. 4, we apply the latter into the macroscale system of PDEs as an additional relationship which is followed by solving poroelastic problems using the provided framework. The concluding remarks and future perspective are presented in Sect. 5.

Governing equations
In this section, for the reader's convenience and understanding, the basic equations for poroelastic media derived from asymptotic homogenisation that are originally developed in [8] and revisited in [9,10] are recalled. in addition to the mentioned equations, we exploit the cell problems derived from asymptotic homogenisation to provide a training dataset for ANN.

Multiscale technique
Let us assume that the poroelastic medium has the average pore size d much smaller than the medium size L (e.g. as shown in Fig. 1). Its solid matrix is linear elastic interacting with incompressible Newtonian fluid that percolates its pores with no-slip boundary condition on the solid-fluid interface.
At this stage, non-dimensionalisation with respect to appropriate characteristic quantities in terms of leng-th scales and velocity fields is carried out as in [8]. For example, where C L, r , and μ c are characteristic pressure, pore radius, and fluid dynamic viscosity, respectively. A very small length scale separation parameter = d/L 1 allows to make use of asymptotic homogenisa-tion technique. We assume that x and y ( y = x/ ) indicate two formally independent macroscale and microscale spatial variables, respectively. Then, the differential operators transform into two independent ones for macroscale and microscale as which separates the problem into two scales. Equating the coefficients of 0 and 1 , by replacing every field by its power of series representation, ψ (x, y) = ∞ l=0 ψ (l) (x, y) l , yields zero-th and first order equations which are the principal ones of the upscaling process. The first result of these equations is that the leading order hydrostatic pressure p (0) and solid displacement u (0) are microscale independent (locally constant) and could be referred to as macroscale variables.
As the relative fluid velocity and the solutions of the cell problems (will be defined later) are not constant at the microscale, we apply the integral average to exploit the acquired equations in the homogenisation process and macroscale. Following complex mathematical machinery, assuming y-periodicity, and no-growth limit [9], the closed differential problem governing the macroscale is obtained as where τ E and v r f indicate effective stress and relative fluid velocity. ε is the solid strain define usually as This system of PDEs is characterised by some coefficients, namely, effective elasticity tensor, Biot coefficient and modulus, and hydraulic conductivity, that are defined, respectively, as where C and φ indicate the elasticity tensor of solid matrix and porosity (i.e. volume fraction of fluid phase), respectively. From the asymptotic homogenisation technique developed in [8] and revisited in [9], the fourth rank tensor M and the second rank tensors Q and W are the solutions of the following systems of PDEs, respectively.
Where Ω s , Ω f , Γ , A, and a represent the solid and fluid domains, their interface, a third rank tensor and a vector, respectively. n is the inward unit vector normal to the solidfluid interface Γ , and with P is an auxiliary vector that encodes microscale information by relating the first order hydrostatic pressure to ∇ x p (0) ( p (1) = −P ·∇ x p (0) ), and I is the second rank identity tensor. In the following, we take the approach that was developed in [5,9] to solve the above-mentioned cell problems via the finite element method and determine the coefficients of the macroscale system PDEs.

Solid cell problems
The system of Eqs. (9)-(11) is rewritten component-wise as follows The solution of the problem (M lkνγ ) can be obtained, using Voigt notation by concatenating the average solution (strain components) of six linear elastic-type cell problems that are achieved by fixing the couple of indices (ν, γ ).
Assuming that the elastic matrix is isotropic at the pore scale, we have the following Neumann boundary conditions on Γ for each problem νγ = 11: f = λn + 2μn 1 e 1 (24) Where n 1 , n 2 , n 3 are the components of the inward unit vector normal to the interface Γ , λ and μ are the Lamé constants derived from C, and e 1 , e 2 and e 3 are the standard unit vectors in the Cartesian coordinate system. At this stage, we can describe the auxiliary problems presented in system of PDEs The weak form of such an elastic problem reads where τ is the Cauchy stress tensor and δu is the virtual displacement. Moreover, the cell problem (12)-(13) reads as a linear elastic problem equipped with inhomogeneous Neumann interface conditions on Γ (a unit inner pressure) and periodic conditions on ∂Ω s \ Γ .
In fact, by choosing a cubic cell with rotational invariance with respect to the three orthogonal axes and an isotropic solid matrix, the effective elasticity tensor with cubic symmetry is obtained. This allows to conduct the analysis of the elastic moduli in terms of the effective Young's modulus, Poisson's ratio and shear modulus E p , ν p , and μ p , respectively, which are related to the independent components of C, namely,C 11 ,C 12 ,C 44 via the following relationships [5,6,30]

Fluid cell problems
The system of Eqs. (15)-(17) can be rewritten componentwise as follows This problem corresponds to three Stokes' problems for i = 1, 2, 3. However, since a geometry with rotational invariance with respect to the three orthogonal axes is chosen, the solution of this problem can be obtained by solving, for example, by i = 1 only, i.e.
where p = P 1 . The above problem formally reads as a periodic Stokes' problem for an incompressible fluid driven by a unit body force directed along e 1 and the solution v is related to the components of the tensor W by the following convention Note, that for this specific setting, the other components of W are very small and assumed to be zero. Finally, in order to solve this problem approximately, the variational formulation of Eqs. (37) to (39) is written as where δv and δ p are test functions.
Having solved the cell problems and determined the coefficients of the macroscale system of PDEs, the macroscale problems can be tackled as is done in [10,11]. However, it is crucial to know the advantages and limits of this framework and, more importantly, how to overcome the latter, which is the core and novelty of this study.

Application of ANN
In this section, we discuss the advantages and downsides of the described upscaling process and introduce the application of ANN to eliminate the limits imposed by the time-consuming process of macroscale properties identification.
The advantages of asymptotic homogenisation can be summarised as -Existence of analytical relationships between microscale and macroscale properties opening the possibility to study the mutual effects of microscale properties and macroscale response during the analysis; -Avoids unnecessary experimental tests to determine the macroscale properties; -Solid matrix elasticity tensor and fluid dynamic viscosity are non-dimensionalised, and the most important inputs (porosity and MSMC) are bounded.
On the other hand, the approach suffers from the following limits -In the case of heterogeneity of the underlying material properties, we need to determine the macroscale properties at each spatial coordinate which is neither practical not efficient; -Restriction to linear cases which means that we are limited to small deformations; -Although considerable effort has been made (e.g. in [9]) the effects of the macroscale response on the physical scale (and consequently microscale) mass, property, and shape change [31] (growth, remodelling and morphogenesis) are not fully considered in the upscaling process.
At this point, it seems that we need to find an alternative approach for determination of macroscale properties to overcome the mentioned limits and to determine the coefficients of macroscale system of PDEs in a fast and reliable way so that we can consider more details of a problem. ANNs are a promising tool for the purpose of determining the relationship between micro and macro scale properties.
In general, the identification of the coefficients of the macroscale system of equations can be summarised by indicating the inputs and outputs. The former consists of solid matrix Young modulus and Poisson ratio, porosity (representing the cell geometry shown in Fig. 1), and fluid dynamic viscosity. The outputs are effective Young modulus, Poisson ratio, and shear modulus as well as Biot modulus and coefficient, and hydraulic conductivity. This transfer can be written as where due to invariant geometry with respect to permutation of the three orthogonal axes, such a transfer can be constructed by means of ANN.
Considering that the elasticity tensor and viscosity are nondimensionalised before the upscaling process (as is shown in Sect. 2.1), one can reduce the above mapping tō and enforce the effects of the characteristic values on the outputs of Eq. (43). This not only renders the ANN smaller and more efficient but also avoids extrapolation, which is crucial because ANN is based on interpolation. An ANN, typically, consists of one input layer, one output layer, and some layers in between which act as an intermediate that determines the output from input data. Each layer has several processors called neurons. The value of each neuron is determined by adding a scalar called bias to the summation of all the neurons of the previous layer, where each one is multiplied by a coefficient called weight, and passing it to the activation function. The latter is a function that determines if a neuron should be active (has a value) or not (zero value) which, in turn, provides a somewhat piece-wise interpolation. In order to reach an accurate prediction of the output, we need to train the network using a dataset of the known values. The ANN training aims at achieving the best weights, and biases, which minimises a distance function called Cost function, of the ANN output and known values. Having a priori training procedure, this is called a supervised ANN.

Training dataset
In order to train an ANN, first, we need a training dataset which consists of a sufficient number of input sets and their respective output sets. In this study, we create the latter by performing a parametric study in which we change the porosity and solid matrix Poisson ratio. We choose 50 different porosities ranging from 0.082 to 0.783, that correspond to the pore radii of 0.1d and 0.4d (with the steps equal to 0.008), respectively. The same number of Poisson ratios are also chosen to consider the role of solid matrix compressibility spanning from 0.02 to 0.498. In this study, we vary neither solid matrix Young modulus nor fluid dynamic viscosity as the results are non-dimensionalised with respect to their characteristic values in the upscaling process according to Eq. (1). The Young modulus and viscosity are chosen to be the same as the ones in a previous study in the literature [5] (E = 13.5[−] and μ = 1[−]) which are non-dimensionalised with arbitrary unit pressure and viscosity. Then, we can validate our training dataset by comparing it with published results. Figure 2 provides the reader with a quick comparison of the results with the ones in [5].

ANN features and architecture
There are several methods for ANN architecture, features, optimisation, learning etc. in the literature (see, e.g. [32,33]). The choice of the relevant elements and methods is application dependent. Here, we make use of Deep Neural Network (DNN) with ReLU (x) = max(0, x) as the activation function. The ANN training process starts from, in this case, a linear Feed Forward which can be expressed as where a  The non-linear profile of Biot modulus and coefficient derived from microscale properties via asymptotic homogenisation, see [5] Having computed the output, we can now calculate the cost function. From a sensitivity analysis Itakura-Saito Distance (ISD) [34] is chosen as the Cost function which reads as where C, n, a (L) m , y m are cost function, the number of training pieces of data, ANN output, and known value. (a (L) m ) 2 is an element-wise square.
Next, an optimisation algorithm to update the weig-hts and biases according to the cost function and complete training cycle is needed. We choose the Adam optimiser as it is based on adaptive estimates of lower-order (first and second) moments [35]. In this algorithm, the weights and biases are updated in each iteration via where n w and n b are the learning rates or step size,¯ is a term to improve the numerical stability.m w ,v w ,m b , andv b are defined and updated in each iteration as where t is the number of iteration and β 1 and β 2 are exponential decay rates for the moment estimates. Finally, the terms m w , v w , m b , and v b are defined and updated in each iteration as We choose the initial values of m w , v w , m b , v b , and t equal to zero. In each iteration, after computing the cost function, the parameters are updated from Eq. (55) to (48) until we reach the target accuracy. As w (i) k j is not directly related to C we need to use the chain rule, for example, the following In this study, we use the "gradient" feature of the open-source framework PyTorch [36]. Having a comparative look at the outputs, which are the coefficients of macroscale system of PDEs, (see, for instance, Figs. 2, 3, 4, 5) we note that the values are of very different orders of magnitude. For example, the non-dimensional maximum Biot modulus is around 5220 while the one of hydraulic conductivity is 0.012, although they are both significant in material characterisation. Consequently, we normalise the values of all outputs by dividing them by their maximum values (which means they all could range between 0 and 1) and re-scale, after training the network.
Finally, choosing an appropriate ANN architecture is an essential factor that affects the efficiency and accuracy. Here, the architecture is identified by a sensitivity analysis. The training time is not a major concern, as it is a one-time effort, but the size of the network is important because it affects the output computation time. Consequently, we choose the smallest network architecture that provides us with the desired ISD (1.0 −6 ), which was three hidden layers with 200 neurons per layer.

Macroscale response
This study, so far, was intended to provide a basis for solving several complex problems such as macroscopic heterogeneity of microscale properties, growth, and remodelling of poroelastic media. Here, we focus on the first one and solve a simple problem to demonstrate the application of the provided framework. Heterogeneous porosity is known to affect the macroscale response of porous material [37][38][39], however, the full effects of it are not investigated yet.

Variational formulation
Variational formulation of Eqs. (4) to (7) is necessary to obtain the weakly continuous form of them to visualise the macroscale response via FE method. The weak form is obtained by integrating mentioned equations with respect to the corresponding volume and multiplying them with arbitrary test functions. The latter, respectively, consists of the variation of the macroscale displacement vector field δu (0) and the variation of the macroscale pressure δ p (0) applied to Eqs. (4) and (6), which then reads Then, to be able to prescribe or measure surface traction and relative fluid velocity (in the direction of the surface normal vector) we apply the divergence theorem which results in the following form Furthermore, free drainage can be introduced on a surface by enforcing where Δx is the distance between the surface of free drainage (S f d ) and the environment which is practically zero but here is a very small number. This could be introduced as a part of the fourth term on the right hand side of Eq. (59). Noteworthy is that if we do not introduce the latter on a boundary surface, an impermeable boundary condition is imposed. The weak formulation is the starting point for discretising and solving the macroscale problem using FEniCS, more detail regarding FE implementation of the macroscale problem can be found in [10].

Verification
The verification of the FE implementation in FEniCS is done by solving the same problem and acquisition of the same results as in the benchmark part of references [10,11].

Example: 1D problem with heterogeneous porosity
Let us assume a column of a poroelastic material with spatially dependent porosity (e.g. a column of consolidated soil under self-weight effects after a long time [39]). Note that, such properties could also vary with time, stress level, etc. due to macroscale response, growth, etc. [31] which are the subject of future research. Solid matrix Young modulus and Poisson ratio are assumed to be E = 15 (MPa) and ν = 0.3, respectively, and porosity is assumed where z is the axial coordinate varying from 0 at the bottom to 15 m at the top. Figure 6 shows the mesh and geometry setup of the model. The conventions that are used throughout this example such as sides, bottom, degrees of freedom etc. are shown in Fig. 7. For the sake of comparison, we solve simi- (the average). Moreover, the fluid is assumed to be water with dynamic viscosity μ f = 8.9 × 10 −4 Pa.s and the average cell dimension to be d = 10 −4 m. Table 1 shows the effective poroelastic properties derived from the underlying microstructural properties of the mentioned cases. The counterplot for the heterogeneous case shows a good agreement with the literature and also with Figs. 3, 4, 5. For example the definite minimum reported in [5] and Fig. 4b can be seen in the last row of Table 1. Moreover, in all the cases, drainage is allowed only from the top surface of the column where a constant mechanical pressure P m = 10 4 Pa in the direction −z is applied. Zero displacements boundary condition is imposed on the bottom surface. We only allow displacements in z (axial) direction of the column so that it resembles a 1D problem.
The load is applied instantly at the beginning of the problem, causing an overpressure which decreases at the top of the column with a high spatial gradient (due to the free drainage condition) as shown in Fig. 9a. The latter drives fluid to percolate through the pores towards the top surface (Fig. 10) resulting in a decrease in the maximum value of pore pres-  Fig. 9d after 1000 s). Figures 8,9,10 show the solution of all the mentioned problems at different times along the axial direction. The general finding of this example is that introducing non-uniform porosity has considerable effects such as a non-linear settlement profile at t = 1000 s shown in Fig. 8. This, as expected, does not follow the overall response of any of the cases with constant porosity which highlights the importance of such factors. Furthermore, the profile of settlement (−u (0) z ) of the heterogeneous case is close to the one with φ = 0. 25[−] at the bottom tending to the one with the highest porosity as we move towards the top of the column which causes a higher polynomial degree of the characteristic relation.
The profile of pore pressure is one of the most important factors to be studied as it has considerable effects on the critical phenomena observed in several disciplines ranging from fracture net pressure [40] and fault reactivation [41] in soil and rock mechanics to abnormal therapeutic agents flow in tumour [42] and rearrangement of its properties [43]. Figure 9 shows the considerable difference that spatially dependent porosity distribution induces. Although similar maximum interstitial hydrostatic pressure of all the cases is shown in Fig. 9a this quantity, as is discussed later, reduces at higher rate in the heterogeneous case than the other ones as shown in Fig. 9d. Figure 10 shows relative fluid velocity induced by the pore pressure gradient. A specific fluid flow profile can be the final goal of several applications such as controlled drug delivery in biology [44] which imposes the need of gaining a better understanding of the effects of different parameters. As shown in Fig. 10a, after a short time (t = 1 (s)), the heterogeneous case demonstrates similar profile to the case of constant high porosity (φ = 0.7[−]) as the considerable pore pressure gradient is located near the top of the column where the porosity is high (see Fig. 9a). However, for longer simulation times, e.g. t = 100 s in Fig. 10c, when pore pressure propagates to the height of the model, the profile of relative fluid velocity diverges from the one of φ = 0.7[−].
The observed profile of relative fluid velocity and pore pressure motivates to monitor the volume change of fluid. This is possible by computing the volumetric integral of ∇ · v r f which, using divergence theorem and for this case could be simplified as where Δ vr f is the relative fluid volume change in the body of the model in one time increment, A f d is the area of the drainage surface, v d r f = v r f ·n indicates relative fluid velocity at the mentioned surface in its normal direction which is, in this case, the axial direction of the column (the only nonzero component) and Δt is the time increment. Figure 11a shows the fact that although the fluid volume change rate of heterogeneous case at early times is very high and near the case with φ = 0.7[−], after passing more or less 500 s it becomes much smaller than for other cases (the fluid volume change is almost constant). The latter, considering Eq. 6, implies that the problem is approaching the steady-state condition which meansṗ (0) → 0. Having no external fluid source together with a free drainage surface (which results in no excess pore pressure), the mentioned condition means that the excess pore pressure in the poroelastic body is approaching to zero which is shown in Fig. 12.
Finally, it is important to note that although the stored volume of the fluid in both, heterogeneous case and the one with φ = 0.475[−], are equal, this does not necessarily mean that the fluid volume change and other results should be equal. The reason is that the spatial profile of the coefficients that characterise the macroscale response change non-linearly at different porosities. For example, from Fig. 4b, one can see that not only the profile of Biot modulus is non-linear, but also it can have a definite minimum.

Example: 3D model of a rock within soil
Let us assume a unit cube of fully saturated soil, with similar material properties as the one in Sect. 4.2, inside which a crystalline spherical rock with radius r = 0.2 m is located. We assume that the latter is impermeable (as its porosity is very low) and considerably stiff compared to the background soil, so that only rigid body movement is probable. Consequently, the rock is modelled with material properties α = 0, K = 0, M → ∞, E p = 60, 000 MPa, and ν p = 0.25 [45]). The soil has a heterogeneous porosity distribution φ = 0.25z + 0.  Figure 13 shows the geometry and the employed conventions in this example.
Early in the process, the applied load causes a high excess pore pressure in the body that decreases with a high gradient at the free drainage surface as shown in Fig. 14a. Note, that the latter is a Neumann BC which results in a pore pressure that is almost equal to the environment pressure (which is assumed zero) where it is applied. The maximum pore pressure decreases sharply over time in all models creating a profile as shown in Fig. 15. We also note that the heterogeneous case stands below other cases, thus leaving room for targeted optimisation. Figure 14a shows that, almost instantly after the load application, the pore pressure increases from z = 0 to z = 0.3 Fig. 9 The pore pressure profile of four cases, with constant and spatially dependent porosities. The heterogeneous case (solid line) shows higher polynomial degree of the characteristic relation. At quasi steady-state (t = 1000s) the gradient of the pore pressure in the heterogeneous case is very small compared with the other cases as discussed in the Sect. 4.2 which is due to the pore pressure concentration around the rock. This behaviour is not observed at later times, such as in Fig. 14b, which will be studied in detail. As time passes, the maximum gradient of the pore pressure decreases and its distribution in the depth of the model is more flat. The pore pressure in the homogeneous cases has inverse relationship with porosity (i.e. the higher the porosity the lower the pore pressure at a specific time after t = 1). This is due to the direct relation of the porosity with hydraulic conductivity (m 2 /Pa.s) which is one of the parameters controlling the reduction rate of the excess pore pressure. On the other hand, from Fig. 14d and Fig. 15, it is concluded that the heterogeneous porosity distribution imposes higher reduction rate of pore pressure (as in the 1D example) which is one of its essential effects.
In order to study the pore pressure concentration near the rock we provide Fig. 16 which shows the pore pressure along the circumference of the inclusion at the plane x = 0.5 m. Figure 16a shows the instant pore pressure profile before considerable consolidation occurs, showing that this parameter decreases moving from z = 0.3 m and 0.7 m to z = 0.5 m (from the top and bottom towards the middle of the rock) and from y = 0.5 m to y = 0.7 m and 0.3 m. This profile at t = 1 s is affected by the drainage from the top of the cube creating a profile, as shown in Fig. 16b, that is the highest at the bottom ((θ = 1.5 rad) and lowest at the top (θ = 0.5 rad) showing that, near the rock, fluid flows from the bottom to the sides and the top. The pore pressure concentration results in nominal stress alternation that can, in practical scenarios, induce seismic activities [46] as an example.
The profiles of the settlement (displacements in −z direction in cartesian coordinate system) shown in Fig. 17 have considerable differences with the one in 1D example without inclusion (see Fig. 8). The instant settlement, which reflects the undrained elastic response of the medium, has direct relationship with the porosity which agrees the equation of the Fig. 10 The relative pore fluid velocity profile of the heterogeneous case shortly after the application of load (t = 1s) is similar to the case with high porosity. However, time passing, it decreases at higher rates compared to other cases Fig. 11 The profile of maximum settlement, which represents the solid volume change, is similar to the one of relative fluid volume change. In the heterogeneous case, both of them become almost constant (indicating the steady-state condition) faster than other cases Fig. 12 The pore pressure in the body with spatially dependent porosity reduces at higher rates than in the cases with constant porosities  + λ p ) where L is the depth of the model, P m is the mechanical pressure, and μ p and λ p are the effective lame constants [11]. Note, that the latter have an inverse relationship with porosity in the cases in this section. At this stage the heterogeneous case remains near the case with φ = 0.375[−]. However, as a consequence of consolidation process, after 1 s the mentioned relationship becomes inverse such that the higher the porosity the smaller the settlement (see Fig. 17b). At later times, such as t = 10 s in Fig. 17c, when the settlement profile is similar to its profile in the steady-state, the case with spatially depen-dent fluid volume ratio at z = 0 m is similar to the case with low porosity and tends to the case of φ = 0.5[−] approaching to z = 1 m. We again highlight that the fluid volume ratio of the mentioned case (heterogeneous) is equal to φ = 0.25 [−] at z = 0 m and to φ = 0.5[−] at z = 1 m. Figure 18a shows that moving from y = 0.5 m to y = 0.7 m in the plane x = 0.5 m the effects of the rock exists but to a smaller extent. It is shown in Fig. 18b that at the edge of the mentioned plane (y = 1 m) this effect is not considerable and the profile is similar to a model without inclusion. This motivates us to study the displacements along the line with (x, z) = (0.5, 0.5) in y direction. Figure 19 shows that the displacement decreases monotonically from the circumference of the rock (y = 0.7 m and y = 0.3 m), with high gradient, to the edge of the plane (y = 1 m and y = 0 m), with small gradient, while it is constant inside the rock. It is noteworthy that the mentioned effect is more notable in the heterogeneous case.
Last but not least, as shown in Fig. 20, the maximum settlement of the model with spatially dependent porosity, initially follows the profile of the case with high porosity while approaching to the steady-state it stands close to the case with average porosity (φ = 0.375[−]).

Example: 2D model of a tumour tissue
In many studies in the literature, the effective properties of poroelastic tissues such as Poisson ratio (PR), Young modulus (YM), and permeability/hydraulic conductivity are either averaged (in macroscale) to be homogeneous and isotropic, or consist of some layers including different effective properties with sharp interfaces [47]. However, using poroelastography and atomic force microscopy methods together with powerful machine learning approaches it is seen that, for example, in tumours, these properties vary gradually from one point in the space to another [48,49]. There are different and sometimes contradictory reports of the material properties of healthy and malignant tissues. For instance, in some studies, in which they make use of the nanomechanical indentation and topography tests via atomic force microscopy, it is concluded that at the tumour core, where there exist more cancer cells, the tissue shows softer behaviour than the outer part (see, for example, [49] and references therein). It is also reported that the cancer cells are more compressible than the normal ones [50,51]. However, in several studies, the PR and YM of healthy tissues are considerably less than the cancerous ones [52][53][54]. We embrace the reports that the tissue at the core of a tumour mass have the properties of cancerous tissues while moving towards the outer boundaries its properties tend to the healthy ones [49,55,56].
In this study, we make use of the values measured via poroelastography in [48], which provides some contour plots  Fig. 14a shows the coordinates indicated on the horizontal axis in the model. Inside the rock there is no pore pressure, so that the lines, here, are only connected linearly from one side to the other side of the rock showing no physical information Fig. 15 The maximum pore pressure profile in time demonstrating non-uniform YM and PR of a tumour mass of a mouse animal. Based on the latter, we assume that PR ranges from 0.45 at the tumour core to 0.3 at the adjacent peripheral region and, for the sake of simplicity, solid matrix YM of 50 kPa is constant in the whole space and time domain. The reported values of the porosities of cancerous tissues are higher compared with the healthy tissues (see, for example, [58,59] and the references therein). The porosity of the highly cancerous tumour core is assumed 50%, which decreases to 20% at the peripheral region. Having assumed the microscale properties at the inner and outer boundaries of the homogenised model, we choose a linear relationship (for the sake of simplicity) describing how these properties vary depending on the spatial location. The identified microscale properties at every specific coordinate are the inputs of our ANN. The outputs of the latter are the coefficients of the macroscale system of PDEs at the associated coordinate. Furthermore, the adjacent peripheral region of a tumour, where     we impose zero displacement BC, is stiffened due to the collagen alignment [49]. In order to minimise the effects of the geometry and focus on the material properties, we choose a spherical geometry for the tumour mass with the radius R t = 1 cm from which a qualitatively small cavity with the radius R c = 0.05R t is cored out. A continuous injection with a constant cavity pressure P c = 1 kPa is applied on the inner boundary as in [11,60]. Figure 21 shows the geometry and labelling used throughout this example. Having a cubic microscale cell with rotational invariance with respect to the three orthogonal axes and with an isotropic solid matrix material, as highlighted in [5], the obtained effective elasticity tensor has a cubic symmetry. Considering symmetric geometry and BCs, we reduce the 3D model to a 2D model assuming that the model response in the azimuthal direction represents the polar one. Zero displacements and the cavity pressure applied at r = R t and r = R c are the Dirichlet BCs and the free drainage condition at r = R t is the Neumann BC.
We study six cases among which three cases, for the sake of comparison, have homogeneous material properties with the assumptions of tumour core, outer area and their average, respectively, while the fourth/fifth case has non-uniform porosity/Poisson ratio in the radial direction. The sixth case has both heterogeneous porosity and Poisson ratio in the radial direction. Table 2 shows the details of the chosen microscale properties of each model.
The resultant coefficients of macroscale system of PDEs are provided in Table 3 including colourmaps for the heterogeneous cases which agree with the results provided in [5]. Figure 22 shows the profile of pore pressure p (0) along the radius highlighting that at short times (e.g. 1 s) this pressure is distributed only very near to the cavity interface with high  21 A representative sketch of the geometry where R c = 0.05 cm and R t = 1 (cm). We impose zero displacements and free drainage on Γ outer and the cavity pressure P c = 1 kPa on Γ cavity spatial gradient where the differences between the models are not considerable. However, when the fluid percolates inside the poroelastic domain, the pore pressure at larger distances from the cavity (e.g. r = 0.2 cm) and the effects of different material properties increase. For example, if we consider p (0) at r = 0.2 cm in Fig. 22a the differences between the models are not considerable, however, at longer times such as in Fig. 22b they are easily noticeable. The latter highlights that the spatial profile of pore pressure depends on the porosity of the medium such that the models with the heterogeneous porosity (Case4 and Case6) have the highest p (0) followed by Case1 with greatest, Case5 and Case3 with average, and Case2 with the lowest homogeneous porosity. Case4 and Case6 are the first models that reach the steady-state (as in the 1D example in the last section) which is again followed by the same order of homogeneous models. For example, by comparing Fig. 22c and d, one would notice that Case2 with the lowest porosity is the last model that (almost) reaches the steady-state. We do not observe considerable effects of the various solid matrix Poisson ratio, which agrees with Eq. (7). According to the latter hydraulic conductivity is the only parameter relating relative fluid velocity to spatial gradient of pore pressure. This parameter depends on porosity, geometry of the pore (e.g. straight or tortuous), characteristic cell dimension, dynamic fluid viscosity, etc. which, except porosity (the lower porosity, the smaller Hydraulic conductivity [5]), are all the same in our models. Note, that in the previous example, we did not have an external source of cavity pressure so, unlike this example, the final state of the problem was zero pore pressure and relative fluid velocity. The considerable pore pressure difference between the cases with homogeneous and heterogeneous porosity distribution, shown in Fig. 22d, highlight the importance of considering more realistic porosity distribution for real-world problems as this factor has a significant impact on tumour growth and drug delivery inside the tumour mass [42,43,57]. Figure 23 shows that the norm of maximum radial relative fluid velocity v r r f decreases as time passes which agrees with the spatial gradient of the p (0) . In other words, at the beginning of the problem, when the pore pressure decreases sharply along the radius, we observe high relative fluid velocity (which is also in agreement with Eq. (7)) and, as time passes and we have smoother pore pressure spatial distribution, the maximum v r r f decreases and its minimum increases (smoother relative fluid velocity distribution). Moreover, at the beginning of the analyses, when we have ∇ p (0) only near the inner boundary, the profile of relative fluid velocity of Case6 and Case4 is similar to the Case1 which is the upper bound of porosity (as shown in Fig. 22a. This is due to the close values of Hydraulic conductivity of the mentioned cases in that region. However, close to the end of the process, such as in Fig. 22c and d the profile of Case6 and Case4 approach the profile of the average porosity (Case5 and Case3).
As solid matrix Poisson's ratio does not affect the Hydraulic conductivity, we do not observe a considerable effect of this parameter on Figs. 22 and 23. This parameter, instead, has a considerable effect on Biot coefficient, Biot modulus, and effective Poisson ratio and small impact on effective Young and shear moduli (see, for example, Table 3) which, in turn, affect the displacement response of the media as shown in Fig. 24. This plot shows that the displacements at short times such as t = 1 s in Fig. 24a are considerably smaller than the ones at t = 1000 s in Fig. 24d. Moreover, from the radial displacement distribution of the problem at θ = 0 rad (where θ is the azimuthal coordinate) along the radius at short times one can observe that the effect of hetero- Fig. 22 Pore pressure along the radius showing a high spatial gradient near the inner boundary at t = 1 s. The high difference between the cases with homogeneous porosity (Case1, Case2, Case3, and Case5) and heterogeneous one is an important parameter with special importance in the tumour growth and therapy [57] geneous solid matrix Poisson ratio is not considerable while at longer times it becomes considerable. In other words, the profiles of Case6 and Case4 (Case5 and Case3) are well separated in Fig. 24c and d.
One important point in Fig. 24 is that the radial displacement at the inner boundary has a negative value that, considering zero displacement BC on the outer boundary, is an indicator of the swelling of the model due to the injection of the fluid which is shown in Fig. 25. The higher pore pressures of Case4 and Case6 shown in Fig. 22 also agree with this observation as it is directly added to the volumetric part of the stress tensor in Eq. (5) (note thatα is a diagonal second rank tensor). In the models with homogeneous porosity (Case1, Case2, Case3, and Case5) which have the similar pore pressure profiles after 1000 s, as shown in Fig. 22d, the final value of the model volume change directly obeys the Biot coefficient value as it is the pore pressure-displacement coupling term (i.e. the more the Biot coefficient, the more the change in the volume due to the fluid injection).
Another effect of a compressible solid matrix is more "deviation from isotropy" (μ p /μ iso ) introduced in [11] (and in different notation in [6]) which is the ratio of the effective shear modulus derived from effective elasticity tensor (which is anisotropic with cubic symmetry) to an isotropic shear modulus calculated from effective Young modulus and Poisson ratio with isotropy assumption (μ iso = E p /2(1+ν p )). It is shown in [11] that the more compressible the solid matrix the more the deviation from isotropy which results in a more deviation from homogeneous displacement distribution in azimuthal direction. Note, that the porosity can also affect this parameter. Figure 26 shows the difference between the radial displacements at θ = 45 deg and θ = 0 deg at the inner boundary (r = 0.05 cm) along the time. From this plot and the ones provided in Fig. 24 and by comparing, on the one hand, Case6 and Case4 and, on the other hand, Case5 Fig. 23 Radial relative fluid velocity (cm/s) along the radius. The maximum value is at the inner boundary where pore pressure spatial gradient is the highest and Case3, it could be concluded that the effect of spatially dependent solid matrix Poisson ratio is larger when we have heterogeneous porosity (i.e. the difference between Case6 and Case4 is more than the one of Case5 and Case3).

Concluding remarks and future work
In this study, we have presented a computational framework using Artificial Intelligence (AI) to identify the macroscale parameters of poroelastic media in an efficient approach. In particular, we have provided an ANN which gives the coefficients of the macroscale system of PDEs from microscale properties in almost real-time (0.003 s). We have trained the network by a dataset of known values which is acquired by solving numerous cell problems of asymptotic homogenisation in 3D. This technique uses the non-dimensional form of parameters such as solid Young modulus and fluid dynamic viscosity so that, one can include their effects directly on the outputs of ANN which reduces the inputs to SMPR and porosity for this specific case. As these properties are bounded, we avoid extrapolation, which is a major concern of using an interpolation-based technique such as ANN. Although here, we assume that the solid phase is isotropic, the cell geometry is invariant with respect to permutation of three orthogonal axes, the pores are straight, etc. this framework is readily applicable to more complex cases by providing an appropriate dataset. We introduce the trained ANN as an additional relationship, established at the beginning of the macro-scale simulation, which determines the effective properties of macroscale governing equations of the model at every spatial coordinate.
Heterogeneous porosity is known to affect the macroscale response of porous material [37][38][39]. However, the full effects of it were not investigated. Hence, we study the macroscale response of two distinct problems with this type of poroelastic material and one with spatially dependent porosity and SMPR. The presented numerical results highlight significant effects of modelling heterogeneous porosity that is not possible to be captured by choosing a constant  have the most volume change as they have more pore pressure as shown in Fig. 22d one for the whole model. The contour plot of the coefficients of the macroscale system of equations show non-linear dependence of them on porosity distribution. We have also provided the profiles of displacements, pore pressure, and relative fluid velocity in space and time as they have significant effects on several phenomena in different disciplines from soil and rock to tumour mechanics [40][41][42][43][44]. In case of the first and second example a poroelastic material with and without inclusion, respectively, is modelled. In both cases, higher polynomial degree of the characteristic relation of the profile of displacements, pore pressure, and relative fluid velocity together with higher drainage temporal rates are the most important consequences of heterogeneous porosity among many others. Moreover, the role of fluid drainage, and consequently heterogeneous porosity, which alters the profile of pore pressure concentration around the spherical impermeable inclusion has been studied. Finally, we have modelled a spherical tumour mass with spatially dependent pore volume fraction and SMPR under constant cavity pressure. Our find- Fig. 26 The difference between the radial displacement at θ = 45 deg and θ = 0 deg and both at r = 0.05 cm is and indicator of heterogeneous displacement distribution along azimuthal direction. The maximum value is before reaching the steady-state while the most difference between the models with similar porosity distribution can be measured after reaching steady-state ings highlight that the value of pore pressure along the radius is significantly higher in the cases with non-uniform porosity (except the zones with prescribed conditions) which results in a greater volume change of the model. The heterogeneous SMPR affects the displacement profile, which is more notable when studying no-uniform radial displacement distribution along the azimuthal direction.
The present work provides a basis to develop a computational framework for multiscale appositional growth [9] and rearrangement of the microstructural features due to the macroscale mechanical and hydraulic response [10,31]. Extension of the presented data-driven approach to multiscale problems in multi-physics contexts, in particular fluid-solid interaction [61,62] and neural network supported rheological characterisation of interfaces [63], is the subject of ongoing research.