Local sensitivity via the complex-step derivative approximation for 1D Poro-Elastic and Poro-Visco-Elastic models

Poro-elastic systems have been used extensively in modeling fluid flow in porous media in petroleum and earthquake engineering. Nowadays, they are frequently used to model fluid flow through biological tissues, cartilages, and bones. In these biological applications, the fluid-solid mixture problems, which may also incorporate structural viscosity, are considered on bounded domains with appropriate non-homogeneous boundary conditions. The recent work in [ 12 ] provided a theoretical and numerical analysis of nonlinear poro-elastic and poro-viscoelastic models on bounded domains with mixed boundary conditions, focusing on the role of visco-elasticity in the material. Their results show that higher time regularity of the sources is needed to guarantee bounded solution when visco-elasticity is not present. Inspired by their results, we have recently performed local sensitivity analysis on the solutions of these fluid-solid mixture problems with respect to the boundary source of traction associated with the elastic structure [ 3 ]. Our results show that the solution is more sensitive to boundary datum in the purely elastic case than when visco-elasticity is present in the solid matrix. In this article, we further extend this work in order to include local sensitivities of the solution of the coupled system to the boundary conditions imposed on the Darcy velocity. Sensitivity analysis is the first step in identifying important parameters to control or use as control terms in these poro-elastic and poro-visco-elastic models, which is our ultimate goal.

1. Introduction. The science of soil mechanics was initiated with the theory of Terzaghi in 1925 [47], who mathematically explained soil consolidation in one dimension with remarkable success for engineering applications. Then Biot generalized the mathematical theory of consolidation to three dimensions and time dependent loads around 1940s [11]. It is his pioneering work that gave a firm mathematical formulation for coupling deformations to flow in porous media, and started the changes in sources can induce pathological changes in the dynamics of the coupling if the visco-elasticity is not present in the structure.
The results of [12] prompted us to investigate numerically the sensitivity analysis of the (elastic displacement, pressure, discharge velocity)-solution with respect to the boundary source of traction (i.e., the Neumann boundary condition for the elastic displacement), and compare the results between the poro-elastic and the poro-visco-elastic cases. In [3], we worked in the framework of the numerical test cases considered in [12], considering both stationary examples, where the data and the permeability are constants, and dynamical ones, with permeability depending nonlinearly on the porosity. In order to compute local sensitivities, we employed the complex-step method, described below. All our numerical results in [3] show that the solution is more sensitive to the boundary traction in the elastic case than in the visco-elastic scenario. This is in line with the well-posedness analysis provided in [12], where higher regularity in time was necessary for the boundary source in the poro-elastic case. Moreover, our simulations show that the effects of the boundary source are most significant for the discharge velocity (and implicitly the fluid energy), especially in the purely elastic case. This might explain why the fluid energy became unbounded as the boundary source of traction lost sufficient smoothness in time, and visco-elasticity was no longer present in [12].
Interestingly, the analysis in [12] does not point out any differences in the regularity requirements for the boundary conditions associated with the fluid, which we will refer to as BC f . This motivates the following questions: "Is the solution of the system not/less sensitive to the Neumann boundary condition associated with the fluid in comparison to the boundary source of traction?" and "Is there a difference in sensitivities w.r.t BC f in the two cases, purely elastic vs. poro-visco-elastic models?" Our motivation for performing sensitivity analysis is due to the fact that sensitivity analysis provides valuable information on how the fluid-mechanical responses are affected by changes in boundary data, and also reveals which parameters are most important in the system and can effectively be used as controls. Sensitivity analysis is the precursor to optimization and optimal control problems, and our ultimate goal is to develop and address relevant control and optimization problems for the poro-visco-elastic models in order to inform novel strategies to improve experimental approaches in bioengineering and medicine. Therefore, in this paper we complete the work started in [3] by including sensitivity analysis with respect to BC f , and thus have a full picture of the influence of boundary data on the fluid-solid mixture solution.
The rest of the paper is organized as follows. In Section 2, we present the one dimensional poro-visco-elastic PDE system considered. In Section 3 we summarize the complex-step method used for the approximation of sensitivity derivatives. Section 4 outlines the computational procedure. The results of sensitivity computations are reported in Section 5. In Section 6, we conclude with our findings and recommendations.
2. The 1D poro-visco elastic model. Following the numerics in [12], we use a 1D model of poro-visco-elasticity for our analysis. Let Ω = (x 0 , x f ) be the domain occupied by the fluid-solid mixture with boundary ∂Ω = {x 0 , x f } and length L = x f − x 0 . We also define the time domain t ∈ (0, T ), T ∈ R + . Assuming negligible inertia, small deformations, and incompressible mixture components, the equations for the balance of mass of the fluid component and the balance of linear 626 BANKS, BEKELE-MAXWELL, BOCIU, NOORMAN AND GUIDOBONI momentum for the fluid-solid mixture which describe the motion of the poro-viscoelastic material are given respectively by ∂ξ ∂t where v is the discharge velocity, σ is the total stress of the mixture, F is the body force, S is the fluid production rate, and ξ( To complete the balance equations, we also include the following constitutive equations. The total stress includes the elastic and visco-elastic stresses (σ e and σ v respectively), which are given by where u is solid displacement, λ e , µ e and λ v , µ v are Lamé elastic and viscoelastic parameters, respectively, and ℘ is the elastic pressure parameter which is discussed in more detail below. The total stress then is given by where p is the fluid pressure and δ ≥ 0 is a parameter indicating the extent to which the model includes viscoelastic effects for the solid component. The discharge velocity is given by where is the permeability which is dependent on porosity. In our computations, we consider permeability k that is either constant (f k = 1) or nonlinearly dependent on porosity, in which case we use the Carman-Kozeny formula for Newtonian fluid flow through spherical particles [24]: Fluid content is given by where c 0 is the constrained specific storage coefficient and α is the Biot-Willis coefficient. However, due to incompressibility of the fluid and solid component, we have c 0 = 0 and α = 1, thus reducing the fluid content to ξ = ∂u ∂x . Since we also have ξ = φ − φ 0 ⇒ φ = φ 0 + ∂u ∂x , we can consider porosity as a function of ∂u ∂x . This means we can also consider permeability k(φ) as a function of dilation ∂u ∂x . The elastic pressure parameter ℘ is defined using the equation and is introduced to avoid displacement differentiation in the evaluation of the permeability. That is, in computation we evaluate k −℘ λe rather than k ∂u ∂x . This prevents degradation of computational accuracy (see [12] for more details).
We let ∂Ω = Γ D ∪ Γ N and Γ D = Γ D,p ∪ Γ D,v where the subscripts N and D correspond to conditions imposed on the stress and displacement, respectively, and the subscripts p and v indicate conditions imposed on the fluid pressure and discharge velocity, respectively. Note that Γ N ∪ Γ D,p ∪ Γ D,v = ∂Ω = {x 0 , x f }, and that Γ N , Γ D,p and Γ D,v can be empty (but not all of them simultaneously). With this notation, we complete the model with the following initial condition and boundary conditions where n is the outward unit normal vector. Further, we note that no boundary condition is needed to impose on ℘ since the total stress is already prescribed on Γ N in (9). In summary, we solve the balance equations (1a) and (1b) with constitutive equations (3), (4), (6), and (7) along with initial condition (8) and boundary conditions (9)- (11). The principle unknowns are the solid displacement u and fluid pressure p.
In the numerical results, we consider the sensitivity of solutions (u, p) w.r.t to (g, ψ) in both elastic and visco-elastic cases. We will also compute the sensitivity of the Darcy velocity v w.r.t. (g, ψ), as the numerical simulations in [12] suggest blow up in the fluid energy in the case of non-smooth boundary traction g when no viscosity (δ = 0) is present in the system. There are many different techniques of approximating sensitivity derivatives, finite difference approximation and use of sensitivity equations being the most common ones (see [7,8,9] and the references therein). Sensitivity equations are accurate and computationally inexpensive for reasonably small systems. However, in our case, we are dealing with a nonlinear fluid-solid mixture, and the sensitivity system will result in a coupled system of linear equations for the sensitivities, coupled with the original nonlinear system. Hence the numerical investigation of this new coupled system will be challenging and not very efficient. The finite difference method is relatively easy and efficient to implement but suffers from cancellation error for small step size.
More recently, a method referred to as the complex-step has gained popularity to calculate sensitivities, especially among aerodynamic engineers ( [1,2,31,32,44]). The complex-step estimate is easy to implement, very accurate and extremely robust while retaining a reasonable computational cost.
In recent efforts [5,6], we demonstrated the efficiency of the complex-step method for computing sensitivities for biological models. There the method is applied to examples of various complexity, ranging from time delayed differential equations (DDEs) to Lamé systems, and the results are compared with solutions of traditional sensitivity equations. Our results showed that the computational effort for the complex-step algorithm is much less than that of sensitivity equations when the size of parameters vectors is not too large. In addition we observed that the method provides a consistent approximation of the derivative starting from h as large as 10 −2 up to h crit = 10 −300 (approximately the machine zero at h = 10 −324 ) with a true second order accuracy. Moreover, we showed that, even though the complex-step formula is derived assuming analyticity of the solution function, the approximation provides accurate one-sided derivatives for functions with far less smoothness (essentially functions that are at least C 2 ).
In the following we summarize complex-step method following [5] and [31]. Let z = x + iy, x, y ∈ R be a complex number and f (z) = f (x, y) = r(x, y) + iw(x, y) be a function of a complex variable. If f is analytic, we have the following Cauchy-Riemann equations that establish the relationship between the real and imaginary parts of the function: For a given step size h we can derive a finite difference-like first derivative estimate for real functions using complex calculus. From the first equation in (12), and the definition of derivatives we have If f takes a real-valued function, then Therefore, for small h, we have the complex-step derivative approximation This formula can also be obtained by approximating an analytic function f with a complex variable using a Taylor series expansion: Taking the imaginary parts of both sides and dividing by h gives Terms of order h 2 and higher can be ignored because the step size h can be chosen up to machine precision. Thus the complex-step derivative is given by The method is accurate down to a specific step size we call h crit . Below h crit , underflow occurs and the approximation becomes meaningless.
We note that the derivative estimate (16) constitutes a big advantage over the finite-difference approach. This is because the finite-difference approximation is subject to cancellation error due to the differencing operation. On the other hand, the accuracy of the complex-step estimates is only limited by the numerical precision of the algorithm that evaluates the function f . For detailed explanations, and remarks on implementation procedures, we refer the reader to [5], [31], [32], and the references therein.
4. Computational procedure. The numerical procedure for computing the sensitivities for the PDE model described in Section 2 consists of three main steps: (1) Semi -discretization in time, (2) Discretization in space and (3) Implementation of the complex-step method for computation of sensitivity derivatives. The first two steps are following the work in [12].
We summarize each step briefly below.

Semi-discretization in time using the Backward Euler (BE) scheme.
For discretization in time we consider a uniform partition of the time interval [0, T ] into subintervals of length ∆t = T /N t . We use the backward Euler scheme for time discretization and obtain the following semi-discrete system: for x in Ω, with The backward Euler scheme is an implicit scheme, i.e. the new approximation of the solution at t = i + 1 appears in both sides of the equations. Therefore, we use the fixed-point iteration to solve the algebraic equation for the unknown solution at t = i + 1.

Discretization in space using the Dual Mixed Hybridized (DMH) finite elements.
To avoid degradation in accuracy due to numerical approximations of derivatives, a dual mixed method is used. In particular a hybridized lowest order Raviart-Thomas finite element technique is adopted. In this technique the dual variables, the total stress σ and the discharge velocity v, as well as the primal variables, the solid displacement u and the fluid pressure p, are treated as independent variables. For more details we refer to [12] and the references therein. Let {I h } be a partition of the domain Ω into subintervals We denote the mid point of I k byx k and its boundary by ∂I k . The associated unit normal vector is given by Let P q (I k ) be the set of polynomials of degree less than or equal to q in I k . For the spatial discretization of the weak formulation of (17)-(24) (see [12]), we use the finite elements and we define spaces: In addition, we set The variables σ h , u h , v h , p h are the approximations of σ, u, v, p in the interior of I k ∈ I h andû h andp h are approximations of u and p at each node of I h .
The finite element approximation of (17)- (24) is: where µ := 2(µ e + δµ v /∆t), m p := (1 + δλ v /(∆tλ e )), The equations (29)-(34) constitute a linear algebraic system for seven scalar variables U h , P h and ℘ h . Including the nodal variablesû andp, the large number of unknowns may result in expensive computation. However, almost all of the equations are local and, consequently, for each element I k ∈ I h the internal elements σ h , u h , v h and p h can be eliminated in favor of the nodal variablesû h andp h and the problem data. This elimination procedure is called static condensation and one of the steps that makes the hybridized method efficient. We refer the reader to [12] for the details.
4.3. Implementation of the complex-step method. For computing sensitivity functions, the complex-step method is built into the DMH finite element algorithm and implemented as follows: Step 1: Define all functions and operators that are not defined for complex arguments such as for example max, min and abs in the program used. For this work, we use Matlab programming language, hence we redefine the Matlab functions appropriately.
Step 2: Add a small complex step ih to the desired sensitivity parameter 'x', and run the DHFEM/BE solver to obtain U (·, x + ih).
Step 3: Compute ∂U/∂x using the formula

5.
Results. For our numerical study, we consider two of the test cases presented in [12]. These cases include two poro-visco-elastic models, one where the permeability of the tissue is constant and the other where the permeability depends nonlinearly on the dilation.

Constant permeability.
In this example, we consider the problem outlined in Section 2 on (0, L) × (0, T ) = (0, 1) × (0, .1), with boundary conditions given by u = p = 0 at x 0 = 0, and σn = g 1 (t) and vn = ψ 1 (t) at x f = 1. We prescribe the following volumetric terms where the spatial and temporal shape functions are given by We assume that the porosity and permeability are constant and given by φ = φ 0 = 0.5 and k = k ref = 1 cm 3 sg −1 , respectively.
The boundary sources (g 1 , ψ 1 ) are functions of time. Therefore, assuming that the data-to-state map is Gâteux differentiable, we compute the directional derivatives of the state (u, p) and v with respect to the functions (g 1 , ψ 1 ) in an arbitrary direction (ḡ,ψ). We approximate both the data and the directions using linear splines. We partition the time interval [0, .1] into ten subintervals of equal length ∆t = .01. We define φ i : [−.01, 0.11] → R for i = 1, · · · , 11 as follows: with the understanding that t 1 = 0 and t 11 = .1. We use these eleven splines to approximate g 1 and ψ 1 , with α i = g 1 (t i ) and β i = ψ 1 (t i ) (see Figure 1): By the chain rule, we have that D α u = (D g 1h u)(D α g). Note that (D α g 1h )[ᾱ 1ᾱ2 ...ᾱ 11 ] T =ḡ. Therefore which means that the action of the directional derivative of u w.r.t. g 1h (or ψ 1h ) is completely defined through the vector D α u (D β u, respectively), and the structure of the direction. Similar results can be obtained for the directional derivatives of p and v. Thus it is sufficient to numerically compute the derivatives of the states with respect to the parameters α i (β i ), which represent the values of g 1 (ψ 1 ) at specific times t i .
In the figures below, we display the graphs of the directional derivatives in the direction of a particular spline (φ 4 ). For all the sensitivities, we include both viscoelastic (δ = 1) and purely elastic (δ = 0) cases, and compare the results.
Comparing Figures 2 and 3 (note the different scales for δ = 0 and δ = 1), we see that u appears to be more sensitive to α 4 than it is to β 4 , regardless of the values of δ. We also note that both ∂u ∂α4 and ∂u ∂β4 are larger in magnitude when δ = 0 than when δ = 1. However, this difference in magnitude is more significant for ∂u ∂α4 where we see a full order of magnitude difference.
Further, we notice that the sensitivity of u with respect to α 4 is more prevalent throughout the spatial domain than the sensitivity of u with respect β 4 , regardless of the value of δ (though, when δ = 1, the sensitivity of u with respect to α 4 steadily decreases as we get farther away from the point where g 1h is applied, the sensitivity of u with respect to α 4 stays consistent throughout the spatial domain when δ = 0). In contrast, the sensitivity of u throughout the temporal domain is most prevalent  when δ = 1, p appears to be more sensitive to ψ than it is to g. However, when δ = 0, the opposite is true. Similar to the sensitivity of u with respect to g, the sensitivity of p with respect to g is about one order of magnitude larger when δ = 0 than when δ = 1. In contrast, the sensitivity of p with respect to ψ reacts in the opposite manner; that is, p appears to be more sensitive to ψ when δ = 1 than when δ = 0. Similar to the sensitivities for u, we see that p is more sensitive to g throughout the spatial domain than to ψ. However, we do not see much sensitivity throughout the temporal domain. This is most likely due to the lack of time derivatives on p. Similar to the sensitivities of p, we see from Figures 6 and 7 that v is more sensitive to ψ than it is to g when δ = 1, but it is more sensitive to g than to ψ when δ = 0. We again see that, when δ = 0, the sensitivity of v with respect to g is one order of magnitude larger than when δ = 1. On the other hand, the magnitude of the sensitivity of v with respect to ψ is the same for δ = 0 and δ = 1.
Moreover, these sensitivities are similar to the sensitives for p in that we do not see the sensitivity continuing throughout the temporal domain. Unlike the sensitivities for u and p, however, the sensitivity of v is more prevalent throughout the spatial domain when δ = 1.
Overall, we see that u is the least sensitive to both g and ψ and v is the most sensitive. We also note that the behavior of the sensitivities with respect to g when we vary δ between 1 and 0 are more consistent and more drastic than the sensitivities with respect to ψ.

Nonlinear permeability.
We consider the problem outlined in Section 2 in the space-time domain (−1, 1) × (0, T ), so that L = 2cm and T = 2s. The boundary conditions are given by Porosity and permeability are now nonlinear functions of the solution. The volumetric terms are given by: where: with ω t = 2π/T and all the other parameter values given as in the above example. Again, since the boundary sources g 2 and ψ 2 are time dependent, we approximate the function using linear splines (as described in the example above): where α i = g 2 (t i ) and β i = ψ 2 (t i ) (see Figure 8).
We computed the sensitivities with respect to all the parameters α i and β j , for i = 1, . . . , 11 and j = 1, . . . , 13. These represent the directional derivatives of the solution with respect to boundary source in the directions of the splines φ i . We display here one representative case, i.e., for α 4 and β 4 . All the results are displayed for both δ = 0 and δ = 1.
In Figures 9 and 10, we see that unlike the sensitivities in the case of constant permeability, here u appears to be more sensitive to ψ than g. While there does not appear to be a large difference in magnitude of the sensitivity of u with respect to g when δ = 0 and δ = 1, we note that when the sensitivity is taken in the direction φ 1 , there is a difference of three orders of magnitude (see [4]). We also see that the sensitivity of u with respect to ψ is one order of magnitude larger when δ = 0 than when δ = 1.   Further, when δ = 0, the sensitivity of u with respect to ψ is much more irregular than when δ = 1. This behavior was seen in the sensitivities of u with respect to ψ in the direction of each of the splines φ i , i = 1, . . . , N ψ . On the other hand, the sensitivities of u with respect to g are irregular in the direction of some of the earlier splines (i.e. j < 4), but smooth out for directions φ j , j = 4, . . . , N g .   In Figures and 11 and 12, we see that p appears to be more sensitive to ψ than to g and both ∂p ∂α4 and ∂p ∂β4 are one order of magnitude larger when δ = 0 than when δ = 1. When δ = 0, the sensitivities of p with respect to g is somewhat irregular when taken in the directions φ j , j ≤ 8, however they smooth out for j = 8, . . . , N g (not shown here). In contrast, the sensitivities of p with respect to ψ are irregular in the directions of each spline and for both δ = 0, 1 [4].
In Figures 13 and 14, we again see that v appears to be more sensitive to ψ than it is to g. Similarly to the sensitivity of u with respect to g, there does not appear to be a large difference in magnitude of the sensitivity of v with respect to g when δ = 0 and δ = 1. However, when the sensitivity is taken in the direction φ 1 , there is a difference of three orders of magnitude [4]. We also see a drastic difference in magnitude of the sensitivities of v with respect to ψ.
Remark. The significant difference in magnitude of sensitivities between the cases δ = 1 and δ = 0 shown in Figure 14 is not due to the particular direction we have chosen. Results in [4] also show similar differences in most "spline directions" as well as a linear combination of them, which we have included here as well (see Figure 15).      Overall, we note that, in this case the sensitivities are much more complex than those of the previous case where we have constant permeability. This is at least in part due to the interaction between the nonlinear permeability with the non-zero body force F and net volumetric fluid production rate S. In fact, we saw in one case when F and S were set to zero, the sensitivities were much more similar to those in our constant permeability and nonzero F , S case [4]. 6. Conclusions. In this paper, we perform numerical sensitivity analysis on the poro-visco-elastic system outlined in Section 2 to further investigate the effects of the boundary data on the solution of the poro-visco-elastic equations. While this work and the work in [3] both expand on the computational results obtained in [12], they also act as a precursor to a control problem by identifying appropriate control parameters. The results in [3] and [12] are concerned solely with the effect of the boundary traction g; in this work we also investigate the effects of the boundary condition ψ, associated with the fluid.
Our numerical results in [3] suggest that the solution (u, p, v) is more sensitive to the boundary traction g in the elastic case than in the visco-elastic scenario. In this work, however, we haven't observed consistent results for the sensitivities with respect to boundary datum ψ. In particular, the elastic displacement u, the pressure p, and the Darcy velocity v are all more sensitive to ψ in the purely elastic case, in the case of nonlinear permeability, dependent on dilation. For constant permeability, the elastic displacement is again more sensitive w.r.t. ψ for δ = 0. However, the pressure appears to be more sensitive to ψ in the poro-visco-elastic case (δ = 1) than in the poro-elastic one, and the magnitudes of the sensitivities of the Darcy velocity w.r.t ψ are about the same in both scenarios. These observations are made based on the additional examples found in the technical report [4].
For the most part, the effects of the boundary sources (both g and ψ) are most significant for the discharge velocity v, especially in the δ = 0 case. This is inline with the reports of the numerical investigation in [12], which hinted that the fluid energy (which is dependent on the discharge velocity) becomes unbounded as the boundary source of traction loses H 1 -smoothness in time, and viscosity is no longer present.
When comparing the sensitivities in terms of the presence/absence of viscoelasticity, the solutions consistently appear to be more sensitive to g in the absence of viscosity (δ = 0) than when viscosity is present (δ = 1), which is not always the case for ψ. This implies that g could be an important control parameter in future control and optimization problems governed by poro-visco-elastic models.
When we compare the sensitivities of each variable to the two boundary terms g and ψ, the difference in magnitude depends on the specific variable when permeability is constant [4]. However, in the case when permeability is nonlinearly dependent on dilation, all of the variables (u, p, v) are consistently more sensitive to ψ than they are to g, regardless of the presence or absence of visco-elasticity in the system. This implies that, according to this preliminary analysis, while the theoretical results in [12] focus on g and F , ψ should not be disregarded from the control point of view.
However, we would like to point out that apart from the presence/absence of viscosity and nonlinearity of permeability, the sensitivity results were also dependent on other characteristics of the problem, such as the volumetric body force F and net volumetric fluid production rate S. This suggests that further investigation may be necessary to make definitive conclusions about the effects of the boundary sources g and ψ.