Journal of Scientific Computing (2022) 92:92 https://doi.org/10.1007/s10915-022-01944-2 Finite Element Methods for Large-Strain Poroelasticity/Chemotaxis Models Simulating the Formation of Myocardial Oedema N. A. Barnafi1 · B. Gómez-Vargas2 ·W. J. Lourenço3 · R. F. Reis3 · B. M. Rocha3 ·M. Lobosco3 · R. Ruiz-Baier4,5,6 · R. Weber dos Santos3 Received: 7 November 2021 / Revised: 5 July 2022 / Accepted: 7 July 2022 / Published online: 27 July 2022 © The Author(s) 2022 Abstract In this paper we propose a novel coupled poroelasticity-diffusion model for the formation of extracellular oedema and infectiousmyocarditis valid in large deformations, manifested as an interaction between interstitial flowand the immune-driven dynamics between leukocytes and pathogens. The governing partial differential equations are formulated in terms of skeleton displacement, fluid pressure, Lagrangian porosity, and the concentrations of pathogens and leukocytes. A five-field finite element scheme is proposed for the numerical approximation of the problem, and we provide the stability analysis for a simplified system emanating from linearisation. We also discuss the construction of an adequate, Schur complement based, nested preconditioner. The produced computational tests exemplify the properties of the new model and of the finite element schemes. Keywords Poroelasticity · Reaction-diffusion · Finite-strain regime · Cardiac applications · Oedema formation · Finite element discretisation Mathematics Subject Classification 92C10 · 65M60 · 74L15 · 35K57 1 Introduction Poroelastic structures are found in many applications of industrial and scientific relevance. Examples include the interaction between soft permeable tissue and blood flow, and the study of biofilm growth and distribution near fluids (see [34, 71]). Recent applications of poroelastic consolidation theory to the poromechanical characterisation of soft living tissues relate with oxygen diffusivity in cartilage [54], swelling of hydrogels [83], feather and scale development [10], tumour localisation and biomass growth [70], cardiac perfusion [11, 24, 30], chemically-controlled cell motion [56], lung characterisation [14], traumatic brain injury B R. Ruiz-Baier ricardo.ruizbaier@monash.edu Extended author information available on the last page of the article 123 http://crossmark.crossref.org/dialog/?doi=10.1007/s10915-022-01944-2&domain=pdf http://orcid.org/0000-0001-9818-0342 http://orcid.org/0000-0001-6068-4510 http://orcid.org/0000-0002-0508-8959 http://orcid.org/0000-0003-3144-5822 http://orcid.org/0000-0002-0633-1391 92 Page 2 of 40 Journal of Scientific Computing (2022) 92 :92 [32], the formation of inflammatory oedema in the context of blood-brain barrier failure [50], and immune systems for small intestine [75, 82] as well as for myocarditis [39]. In this work we are concerned with the latter application involving the formation and evolution of oedema, a build up of excess of fluid content in the intercellular space, in myocardial tissue. Local infection of tissue by pathogens contribute to an inflammatory reaction driven by the immune system with the aim of protecting the compromised region against invading organisms. Driven by oncotic and hydrostatic pressure gradients, there is a transvascular water flux: capillaries transport fluid into and out from the interstitium (the intercellular region), and lymphatic nodes contribute to returning fluid out from that space. This work also considers that, before infection occurs, the body is under osmotic equilibrium since the capillary barrier is permeable to ions. In this case, ions have no significant effect on the interstitial flow under normal conditions, otherwise, a previous disease may exist. The presence of pathogens inside a cell as well as tissue damage caused by them may change this scenario, i.e., ion-concentration may change, impacting local pressure. However, in this work, we assume, for simplification purposes, that the main cause of oedema formation is the presence of the pathogen in the tissue and the immune response to it. In the process of inflammatory response, local vasodilation occurs, which not only allows leukocytes to leave the bloodstream and to access the infection site but also increases the accumulation of fluid in the intercellular space, see Fig. 1. However such abnormal accumulation of interstitial fluid content may lead to local stress generation, tissue deformation and swelling, eventually producing a number of complications that can impair the normal function of the tissue. Oedema may be found in several types of tissue. Heart oedemas are commonly observed in myocarditis and they constitute one of the main criteria used in medical practice to its diagnosis [41, 62]. Deriving refined continuum-based models for cardiac oedema formation valid in the large strain regime and that include 3D personalised geometries have the potential to describe mechanistic processes that are observable by clinical evidence and to provide deeper insight into the complex multiscale mechanisms that are inherent to the impairing of the healthy cardiovascular function. We extend the preliminary results from [39, 40] and develop a simplified phenomenolog- ical model for the dynamic interaction between poroelastic finite strain deformations and the chemotaxis of leukocytes towards pathogens. The present framework assumes that the porous skeleton admits heterogeneities at the macro scale, and that the pores are all interconnected. We also consider that inertial forces are negligible. The solid-fluid interaction describing the distribution of flow within the deformable skeleton is then cast in Lagrangian coordinates using the solid displacement, the fluid pressure, and the nominal porosity (ratio between the pore volume and the total volume), whereas the immune system dynamics are represented by the concentrations of leukocytes and a pathogen. The model for the immune system has also been modified to better represent the dynamics that occur after the pathogen has been elimi- nated. Furthermore, the mechanical feedback into pathogens and leukocyte populations are now taken into account, and we have explored these aspects further in our recent work [52], which uses the present formulation for the coupled poromechanical/chemotaxis equations. Even if many numerical methods for the approximation of linear poroelasticity and their convergence analysis can be found in the recent literature (see [12, 17–19, 55] and the ref- erences therein), much less attention has been given to formulations addressing large-strain poromechanics. In this regard, recent works include an enriched Galerkin framework [27], stabilised finite elements [15, 84] and hybrid finite elements as well [83]. Although most contributions address the problem in a Lagrangian frame of reference, some alternative for- mulations include descriptions in Eulerian [65] and ALE [22] coordinates. Besides the stress response, nonlinear models differ from the linear ones in that pressure is not a primary 123 Journal of Scientific Computing (2022) 92 :92 Page 3 of 40 92 variable, but instead it is given by constitutive modelling [25, 31]. Instead, we follow the framework of [30, 53] (and other similar models) where the issue is circumvented through solid phase incompressibility, which yields a Lagrange multiplier acting as hydrostatic pres- sure. The extension of such solution techniques for large scale simulations is highly non-trivial due to the requirement of efficient preconditioned iterative methods for the poroelasticity problem. Preconditioning strategies are mainly given by operator preconditioning [9, 46] and block factorisations [42, 79]. In this regard, block preconditioners are more easily gen- eralisable, but an efficient implementation requires an adequate approximation of the Schur complement, usually approximated through a fixed-stress formulation in poroelasticity [38, 42, 79]. Despite all of this, an efficient numerical solver for a linearisation of large-strain poromechanics with nominal porosity as a primary variable remains largely unaddressed. In the present case we formulate the set of nonlinear equations using as primary variables the solid displacement, the fluid pressure, the Lagrangian porosity, and the concentrations of leukocytes and pathogens. A finite element method is employed for the space discretisation [67], combined with a first-order backward Euler time advancing scheme. The method uses a MINI element for the displacement-porosity pair [6] (see also [8, 23] for the case of hyperelasticity), and linear Lagrangian elements for the remaining unknowns. As the stability of the fully nonlinear system is still a paramount task (and an open problem even for the segregated nonlinear poroelasticity without the reaction-diffusion coupling, see, e.g., [12]), we address the solvability of the linearised set of equations, focusing on the semidiscrete in- time formulation. For the corresponding analysis, we propose a fixed-point strategy, and then, Schauder fixed point theorem, together with the Fredholm alternative, the Babuška–Brezzi theory, the Lax–Milgram lemma, and classical theory of quasi-linear equations, allow us to assert unique solvability robustly with respect to the material parameters of the linearised poroelastic solid. Then, we address the numerical approximation of the linearised problem through a Krylov subspace method with a Schur complement preconditioner. The main advantages of the proposed mathematical model and the associated computa- tional methods are i) a consistent theoretical framework valid for finite strains that stems as the natural generalisation of three-field linear poroelasticity from [60] (which is written in terms of displacement, fluid pressure, and total pressure), ii) the versatility of the formulation to accommodate 2D or 3D geometries, iii) the accuracy of the numerical scheme and effi- ciency of the Krylov subspace iterative methods under consideration, and iv) the potential of the methodology in replacing current invasive methods for the detection of interstitial fibrosis and myocarditis (such as endo-myocardial biopsy) by techniques hinging only on MRI data. Wehave organised the contents of this paper in the followingmanner. Section 2 outlines the main details of the model problem, motivating each component in the balance equations and statingmain assumptions together with typical boundary conditions. In Sect. 3 we present the weak formulation of the nonlinear coupled problem, as well as the derivation of the linearised form of the weak formulation. Then, in Sect. 4, we address the solvability of the linearised set of equations. Section 5 contains the derivation of a finite element scheme, including the fully discrete problem in matrix form, and we address the construction of efficient solvers tailored for the multiphysics coupled system. In Sect. 6 we present a number of computational results, consisting in a simple numerical study concerning the sensitivity analysis of relevant model parameters, the verification of spatio-temporal convergence, and an investigation of different cases on simplified and physiologically accurate geometries. We close with some remarks and a discussion on model extensions collected in Sect. 7. 123 92 Page 4 of 40 Journal of Scientific Computing (2022) 92 :92 2 ContinuumModel and Proposed Set of Field Equations The different scales and structural components in the tissue suggest to employ a continuum model. We assume that the tissue is a poroelastic medium with interstitial fluid completely filling the void spaces. Adopting usual kinematic conventions and notation, let us consider a domain Ω ⊂ R d , d = 2, 3 representing the volume occupied by a deformable porous structure in its reference configuration. We will denote by n the outward unit normal vector on the boundary ∂Ω . We also disjointly partition the boundary into the sub-boundaries Γ and Σ (also in the reference undeformed state) where different types of boundary condi- tions, associated with the equations of motion, will be applied. A material particle in the undeformed domain Ω is identified by its position x (and, unless specified otherwise, all dif- ferential operators such as gradients and divergences will be evaluated with respect to these undeformed coordinates), whereas u : Ω → R d will denote the displacement field defining its new position x+ u in the deformed configuration. The tensor F := I+∇u is the gradient (applied with respect to the fixed material coordinates) of the deformation map and I denotes the identity tensor; its Jacobian determinant, denoted by J = det F = det(I + ∇u), mea- sures the dilation (solid volume change during the deformation); and C = FtF is the right Cauchy–Green deformation tensor, where the superscript ()t denotes the transpose operator. Thenon-viscousfiltrationflow through thenon-deformedporous skeleton canbedescribed by Darcy’s law [2, 3]. Adopting its parabolic formulation solely in terms of pressure head (which condenses momentum and mass conservation equations for the fluid) and combining it with the conservation of linear momentum in the solid, we obtain the extension of Biot consolidation equations to the regime of large strains. This can be carried out using the descriptions already available in, e.g., [15, 53, 84]. In this paperwe choose towrite themusing the first Piola–Kirchhoff stress tensor P = Peff − α pJF−t and the effective (hyperelastic) stress Peff = ∂Ψs ∂F . For incompressible neo–Hookean and Holzapfel–Ogden solids one has the following forms (see, e.g., [45, 63]) Ψ NH s (C) = μs 2 (I1(C) − d), Ψ HO s (C) = a 2b eb(I1(C)−d) + a f s 2b f s [ eb f s I8, f s (C)2 − 1 ]+ ∑ i∈{ f ,s} ai 2bi [ ebi (I4,i (C)−1)2+ − 1 ] , (2.1) respectively, where μs , a, b, ai , bi with i ∈ { f , s, f s} are material parameters and we have used the notation (ζ )+ := ζ if ζ > 0 or zero otherwise, and the invariants and pseudo- invariants of Cmeasuring isotropic and direction-specific stretches I1(C) = trC, I4, f (C) = f0 ·(C f0), I4,s(C) = s0 ·(Cs0), I8, f s(C) = f0 ·(Cs0). Here f0 and s0 denote local alignment of fibres and sheetlets which is inherent to the local microstructure of the myocardial tissue [58]. Histological studies based on diffusion tensor imaging indicate a fibre architecture that rotates 120◦ from the epicardial to the endocardial surface. A volumetric part of the energy is added for the case of nearly incompressible materials, for example, the term U (J ) = 1 2λs(lnJ )2 (with λs another material parameter) which is designed to enforce convexity or polyconvexity of the strain density [85]. Provided that the volumetric stiffness of the solid is substantially smaller than that of the constituent (which is reasonable for soft living tissues, see, e.g., [50, 54]), volume changes of solid constituent are negligible, implying that J = 1 + φ f − φ0, 123 Journal of Scientific Computing (2022) 92 :92 Page 5 of 40 92 where φ f − φ0 is the nominal porosity change in the reference configuration [53] (here we are assuming that the mixture is saturated, that is, φ f + φs = 1, where φs is the volume fraction of the solid). Recall that the Lagrangian porosity is the natural state variable being work-conjugate to the fluid pressure p (see, e.g., [59]). Should pre-stress be considered, we need to take instead J = J0+φ f −φ0, where J0 is the volume associated with the pre-stress. Note that the specific form of the strain energy density Ψs is considered irrespective of the fluid saturating the solid. It can also be derived from subtracting the volumetric free energy of the fluid from the total Helmholtz free energy of the mixture Ψs = Ψ − φ f Ψ f , where φ f is the nominal (Lagrangian) porosity measuring the current fluid volume per unit reference total volume [59, 84]. Note also that rearrangement of tissue components as a result of change in fluid content will imply an additional solid deformation as well as a stress modification [27, 50, 61]. The problem is then stated in a Lagrangian framework and in terms of the solid displace- ment u, the fluid pressure p, and the nominal porosity φ f , as −Div[Peff − α pJF−t] = ρsb in Ω × (0, tfinal], (2.2a) ρ f Dtφ f − 1 J Div ( φ f ρ f JF−1 κ(J , φ f ) μ f F−t∇ p ) = ρ f � in Ω × (0, tfinal], (2.2b) J − φ f = 1 − φ0 in Ω × (0, tfinal], (2.2c) where Div and Div are the divergence operators (with respect to the coordinates in the refer- ence configuration) for vector and tensor fields, respectively. Equation (2.2a) is the balance of linear momentum for the multiphase system expressed in terms of the true total stress; (2.2b) states the mass conservation (or continuity) of fluid content in a material framework, written in terms of the nominal flux that includes the movement of the fluid phase due to pressure gradient under Darcy’s law (which in particular discards shear flow effects and is valid for small Reynolds numbers); and (2.2c) represents the aforementioned kinematic constraint relating volume and porosity. It states the material incompressibility of the con- stituents, but it still allows for compression of the medium since the pore volume can change locally in connection with modifications in the macroscopic mass density of the poroelastic material. A more complete discussion can be found in [53]. This description is independent of the material constitutive model for the solid matrix, and the balance of angular momentum simply furnishes the condition FPt = PFt. Here Dt indicates the material time derivative (with respect to skeleton particles), b is a vector of external body loads, � is a distributed source/sink strength term (local source volume of fluid added per volume of tissue per unit time) that accounts for fluid being added or removed from the interstitium by the blood and lymph capillaries and its specific form will be made precise later on, ρs is the reference density of the porous matrix, ρ f is the reference density of the interstitial fluid, α is the Biot–Willis modulus assuming values between 0 and 1, and μ f is the dynamic viscosity of the interstitial liquid. The porosity-dependent heterogeneous tensor of intrinsic permeability (or hydraulic con- ductivity) κ can be defined by using normalised Kozeny–Carman (KC), exponential (E), or by power-law (PL) relations defined, respectively, as κKC iso (J , φ f ) = κ0I + κ0 φ3 0 (1 − φ0) 2 Jφ3 f (J − φ f ) 2I, κE iso(J , φ f ) = κ0 ( J − φ3 f 1 − φ3 f )3 exp(M0(J 2 − 1)/2)I, κPL iso(J , φ f ) = κ0 ( Jφ f φ0 )2/3 I, 123 92 Page 6 of 40 Journal of Scientific Computing (2022) 92 :92 where κ0 = d2 mφ3 0/180(1 − φ0) 2, dm is the typical diameter of pores or of solid grains, and M0 is the slope of the normalised permeability curve [7, 14, 84]. As porosity is a volumetric quantity, these forms of permeability are isotropic and therefore indifferent to rotation and shear effects. A relatively simple phenomenological extension that accommodates anisotropy due to tissue microstructure is given as follows (using as example the Kozeny–Carman law) κKC(J , φ f ) = κ0 ( 1+ 1 φ3 0 (1−φ0) 2 Jφ3 f (J −φ f ) 2 )( κ f0 I4, f F f0 ⊗ F f0+ κs0 I4,s Fs0 ⊗ Fs0 ) , (2.3) where κ f0 , κs0 are weights for the principal hydraulic conductivities such that the ratio κ f0/κs0 ≈ 2.5 coincides with the anisotropy ratio typically used for the conductivity ten- sor in the extracellular domain of bidomain models [29, 74] (and much milder than the ratio observed in the intracellular domain, or in other porous skeletons such as sand or cartilage [21, 37]). This allows to reproduce the preferential perfusion direction of the micro-vascular sys- tem following the orientation of cardiac fibres f0, s0. Other transformations that go beyond the kinematic nonlinearity induced by change of frame can be found in models exploiting strain-dependence [7] or stress-assistance mechanisms [26]. Unless otherwise specified, we will use (2.3). The infection process starts when an invader (pathogen) enters the poroelastic tissue. Once there, it starts to propagate. The exact mechanism adopted for increasing the pathogen population (replication, in the case of viruses, or reproduction, otherwise) is not essential to our model. In our simplification, γp represents the rate at which the pathogen propagates. Fig. 1 Sketch of local infection of tissue by pathogens, triggering inflammatory and immune responses. In the inflammatory response process, local increment vascular permeability occurs, which allows leukocytes to leave the bloodstream to access the infection site which also leads to fluid accumulation in the interstitial space. Driven by oncotic and hydrostatic pressure gradients, blood capillaries transport fluid into and out of the interstitium, and lymphatic capillaries contribute to the resolution of the oedema 123 Journal of Scientific Computing (2022) 92 :92 Page 7 of 40 92 The innate immune cells (leukocytes) use their receptors to sense the presence of non-self molecules characteristic of pathogens. After the presence of pathogens is identified, leuko- cytes start to produce pro-inflammatory cytokines, which increase vascular permeability and recruit more leukocytes to the site of infection. Themodel does not include cytokines because the concentration of pathogens and leukocytes can indirectly represent their effect on vascular permeability and recruitment of leukocytes. Additionally, we assumed a simplified reaction term for the leukocytes since among these cells some have a very long life span when com- pared to the timescale of the simulations considered in this work. For example, macrophages can live from months to years [76]. Evidently, leukocyte death from apoptosis and pathogen phagocytosis create a sink that was not considered in the current work. Leukocytes adhere to the endothelial cells that line the vessels, a phenomenon called margination (Fig. 1). The increase in vascular permeability facilitates leukocyte extravasation to the tissue. The model considers that leukocytes leave the bloodstream to the tissue with an extravasation rate equals to λpl . The leukocytes move towards the site of infection due to the presence of chemoat- tractants. The chemotactic rate is given by χcl . Finally, leukocytes can start the phagocytosis of pathogens, eliminating them at a rate equal to λlp (Fig. 1). We denote by cp and cl the concentrations of pathogens and leukocytes, respectively. Both populations can also diffuse in the interstitial fluid. Their dynamics, in the current configuration and with derivatives in terms of the deformed coordinates, are governed by ∂t (φ f cp) + ( ∂tu − κ μ f ∇ p ) · ∇cp − div(Dp∇cp) = r p(φ f , cp, cl ) in Ωt ×(0, tfinal], ∂t (φ f cl ) + ( ∂tu − κ μ f ∇ p ) · ∇cl − div(Dl∇cl − χcl∇cp) = rl (φ f , cp, cl ) in Ωt ×(0, tfinal], whereDl ,Dp are positive-definite diffusion tensors for the species within the interstitial fluid, and rl , rp are reaction terms adopting the following specific forms (see, e.g., [39]) rp(φ f , cp, cl) = φ f (γ̄p − λ̄lpcl)cp, rl(φ f , cp, cl) = λ̄plφ f cpcl , with γ̄p, λ̄lp , and λ̄pl being constant coefficients (pathogen growth rate, leukocyte phagocy- tosis rate, and leukocyte migration rate, respectively) [39]. Note that the advecting velocity corresponds to the filtration velocity (that between the solid and fluid velocities, as the species are supposed to be transported in the interstitial fluid). Pulling back these advection-reaction- diffusion equations to the reference frame we obtain a modified system with kinematic nonlinearity Dt (φ f cp) − 1 J Div(φ f JF−1DpF−t∇cp) = r p(φ f , cp, cl ) in Ω ×(0, tfinal], (2.4a) Dt (φ f cl ) − 1 J Div(φ f JF−1DlF −t∇cl − χφ f cl JF−1F−t∇cp) = rl (φ f , cp, cl ) in Ω ×(0, tfinal]. (2.4b) Note that, in order to make the advection completely absorbed by the total derivative we have considered the total fluid-solid velocity. In turn, this implies that we are assuming that pathogens and leukocytes are able to move in the solid also. Next we return to the continuity Eq. (2.2b) to specify the source (rate of fluid movement) �. This term will encode the immune response feedback into the hydromechanics by setting a Starling-Hill function of the type �(p, cp) = c f (cp)[pc − p − σ(cp)(πc − πi )] − �0 [ 1 + vmax(p − p0)n kn m + (p − p0)n ] , 123 92 Page 8 of 40 Journal of Scientific Computing (2022) 92 :92 specified by the nonlinear functions of the pathogen concentration (capillary conductivity and reflection coefficient, respectively) c f (cp) = (S/V )L p0(1 + Lbpcp), σ (cp) = σ0(1 + Lbr cp) −1. The remaining coefficients �0, pc, πc, πi , vmax, km, n, p0, S/V , L p0, Lbp, Lbr (normal lymph flow, capillary basal pressure, capillary oncotic pressure, interstitial oncotic pressure due to plasma protein, maximal lymph flow, half-live of pressure to lymph flow, Hill coeffi- cient, normal interstitial fluid pressure, vessel area per volume unit, healthy tissue hydraulic permeability of the capillary wall, intensity of the pathogen infection into permeability, and influence of pathogens in the reflection coefficient; respectively) are positive model constants [39]. To close the system composed by (2.2a)–(2.2c) and (2.4a)–(2.4b), we need to provide suitable initial data and boundary conditions. We suppose that the system is initially at rest and with an homogeneous distribution of the leukocytes and pathogens in the domain φ f (0) = φ0, p(0) = 0, u(0) = 0, cl(0) = cl,0, cp(0) = cp,0 in Ω. (2.5) In addition, on the boundaries Γ and Σ we prescribe zero traction and displacement, respec- tively; and zero flux everywhere on the boundary for the species interacting in the immune system. That is, we consider the following set of boundary conditions u = 0, and φ f JF−1 κ(J , φ f ) μ f F−t∇ p · n = 0 on Γ × (0, tfinal], (2.6a) Pn = JF−t tΣ, and p = p0 on Σ × (0, tfinal], (2.6b) φ f JF−1DpF−t∇cp · n = [φ f JF−1DlF−t∇cl − χφ f cl JF−1F−t∇cp] · n = 0 on ∂Ω × (0, tfinal], (2.6c) where tΣ is a given boundary traction applied in the current configuration. 3 Weak Formulation and Consistent Linearisation We derive a weak form for the coupled system, carry out a Newton–Raphson linearisation, and then discuss the solvability of a simplified system arising from the linearisation with respect to conditions at rest. That resulting problem coincides in structure with the coupled Biot/reaction-diffusion equations. 3.1 GeneralWeak form and in-time Discretisation In view of the essential boundary conditions for displacement and fluid pressure (2.6a) and (2.6b), we define the Hilbert spaces H1 Γ (Ω) = {v ∈ H1(Ω) : v|Γ = 0}, H1 Σ(Ω) = {q ∈ H1(Ω) : q|Σ = 0}, associated with the classical norms inH1(Ω), and H1(Ω), respectively. The nonlinear weak form of the equations ofmotion and of the transport for the immune system can be established following a standard approach, that is, multiplying each field equation by a suitable test function, integrating overΩ , and invoking the divergence theorem whenever appropriate. To lighten the presentation, we drop the measures dx and ds from the integrals, which are to be understood throughout the manuscript as the Lebesgue measure. This leads to the following 123 Journal of Scientific Computing (2022) 92 :92 Page 9 of 40 92 continuous, five-field weak formulation: For almost all t > 0, find (cp, cl) ∈ [H1(Ω)]2, u ∈ H1 Γ (Ω), p ∈ H1 Σ(Ω), and φ f ∈ L2(Ω), such that ∫ Ω J Dt (φ f cp)wp + ∫ Ω φ f JF−1DpF−t∇cp · ∇wp = ∫ Ω Jrp(φ f , cp, cl)wp ∀wp ∈ H1(Ω), ∫ Ω J Dt (φ f cl)wl + ∫ Ω φ f JF−1(DlF−t∇cl − χclF−t∇cp) · ∇wl = ∫ Ω Jrl(φ f , cp, cl)wl ∀wl ∈ H1(Ω), ∫ Ω (Peff − α pJF−t) : ∇v − ∫ Σ JF−t tΣ · v = ∫ Ω ρsb · v ∀v ∈ H1 Γ (Ω), ∫ Ω ρ f J Dtφ f q + ∫ Ω ρ f φ f JF−1 κ(J , φ f ) μ f F−t∇ p · ∇q = ∫ Ω ρ f J�(p, cp) q ∀q ∈ H1 Σ(Ω), ∫ Ω (J − φ f )ψ = ∫ Ω (1 − φ0)ψ ∀ψ ∈ L2(Ω). (3.1) We do not address time regularity in this work, and simply assume that the solution is weakly differentiable and continuous in time, the latter in order to obtain well-defined ini- tial conditions. Next we partition the interval [0, tfinal] into N evenly spaced non-overlapping sub-intervals of fixed lengthΔt and apply a time semi-discretisation of (3.1) using the uncon- ditionally stable, backward Euler’s method. Starting from initial data, at each time iteration n = 0, . . . , N − 1 we have a nonlinear problem with solution (cn+1 p , cn+1 l ) ∈ [H1(Ω)]2, un+1 ∈ H1 Γ (Ω), pn+1 ∈ H1 Σ(Ω), and φn+1 f ∈ L2(Ω). Then we perform a Newton–Raphson method linearisation. This means that at each time step n + 1 we start the Newton iter- ates using as initial guess the solution at the previous time step ck=0 p = cn p , ck=0 l = cn l , uk=0 = un , pk=0 = pn , φk=0 f = φn f , and for k = 0, 1, . . . we look for the solution incre- ments (δck+1 p , δck+1 l ) ∈ [H1(Ω)]2, δuk+1 ∈ H1 Γ (Ω), δ pk+1 ∈ H1 Σ(Ω), δφk+1 f ∈ L2(Ω) satisfying a1(δck+1 p , wp) + b∗ 1(δck+1 l , wp) + b∗ 2(δu k+1, wp) + b∗ 3(δφ k+1 f , wp) = F1(wp) ∀wp ∈ H1(Ω), b1(wl , δck+1 p ) + a2(δck+1 l , wl) + b∗ 4(δu k+1, wl) + b∗ 5(δφ k+1 f , wl) = F2(wl) ∀wl ∈ H1(Ω), a3(δuk+1, v) + c∗ 1(δ pk+1, v) = F3(v) ∀v ∈ H1 Γ (Ω), c2(q, δck+1 p ) + c1(q, δuk+1) + a4(δ pk+1, q) + c∗ 3(δφ k+1 f , q) = F4(q) ∀q ∈ H1 Σ(Ω), c4(ψ, δuk+1) − a5(δφ k+1 f , ψ) = F5(ψ) ∀ψ ∈ L2(Ω). (3.2) Then the solution is updated ck+1 p = ck p + δck+1 p , ck+1 l = ck l + δck+1 l , uk+1 = uk + δuk+1, 123 92 Page 10 of 40 Journal of Scientific Computing (2022) 92 :92 pk+1 = pk + δ pk+1, φk+1 f = φk f + δφk+1 f , until either the increments or the residuals drop below a given tolerance. In order to specify these variational forms, as usual we use the directional derivatives of all nonlinearities at a particular trial solution in the direction of the increments. After rearranging terms and omitting the superscript (·)n+1 whenever clear from the context, the bilinear forms and linear functionals in (3.2) are defined as a1(δcp, wp) := ∫ Ω J k (2φk f − φn f Δt − φk f ∂rp(ck p, ck l ) ∂ck p ) δcpwp + ∫ Ω J kFk,−1DpFk,−t∇δcp · ∇wp, b∗ 1(δcl , wp) := − ∫ Ω J k ∂rp(φ k f , ck p, ck l ) ∂ck l δcl wp, F1(wp) := ∫ Ω Rk,n 1 wp, F2(wl) := ∫ Ω Rk,n 2 wl , b∗ 2(δu, wp) := ∫ Ω J kFk,−1([Fk,−t :∇δu]Dp −∇δuFk,−1Dp − DpFk,−t(∇δu)t ) Fk,−t∇ck p · ∇wp + ∫ Ω J kFk,−t : ∇δu ( φk f − φn f Δt ck p + ck p − cn p Δt φk f − rp(φ k f , ck p, ck l ) ) wp, b∗ 3(δφ f , wp) := ∫ Ω J k ( 2ck p − cn p Δt − ∂rp(φ k f , ck p, ck l ) φk f ) δφ f wp, F3(wl) := ∫ Ω Rk 3,Ω · v + ∫ Σ Rk 3,Σ · v, a2(δcl , wl) := ∫ Ω J k (2φk f − φn f Δt − ∂rl(φ k f , ck p, ck l ) ∂ck l ) δclwl + ∫ Ω J kFk,−1DlFk,−t∇δcl · ∇wl − ∫ Ω χ J kδclFk,−1Fk,−t∇ck p · ∇wl , b1(wl , δcp) := − ∫ Ω χ J kck l F k,−1Fk,−t∇δcp · ∇wl − ∫ Ω J k ∂rl(φ k f , ck p, ck l ) ∂ck p δcpwl , b∗ 4(δu, wl) := ∫ Ω J kFk,−1([Fk,−t : ∇δu]Dl − ∇δuFk,−1Dl − DlFk,−t(∇δu)t ) Fk,−t∇ck l · ∇wl + ∫ Ω J kFk,−t : ∇δu ( φk f − φn f Δt ck l + ck l − cn l Δt φk f − rl(φ k f , ck p, ck l ) ) wl + ∫ Ω χck l J k ( Fk,−1∇δuFk,−1 + Fk,−1Fk,−t(∇δu)t − Fk,−t : ∇δuFk,−1 ) × Fk,−t∇ck p · ∇wl , 123 Journal of Scientific Computing (2022) 92 :92 Page 11 of 40 92 b∗ 5(δφ f , wl) := ∫ Ω J k ( 2ck l − cn l Δt − ∂rl(φ k f , ck p, ck l ) φk f ) δφ f wl , c∗ 1(δ p, v) = − ∫ Ω αδ pJ kFk,−t : ∇v, a3(δu, v) := ∫ Ω ( ∂Peff (Fk) ∂Fk ∇δu + α pk J k[Fk,−t : ∇δu I − Fk,−t(∇δu)t]Fk,−t ) : ∇v − ∫ Σ J k[Fk,−t : ∇δu I − Fk,−t(∇δu)t]Fk,−t tΣ · v, c1(q, δu) := ∫ Ω φk f J kFk,−1 ( Fk,−t : ∇δu κ(J k, φk f ) μ f − ∇δuFk,−1 κ(J k, φk f ) μ f − κ(J k, φk f ) μ f Fk,−t(∇δu)t ) × Fk,−t∇ pk · ∇q + ∫ Ω ρ f J kFk,−t : ∇δu ( φk f − φn f Δt − �(pk, ck p) ) q a4(δ p, q) := ∫ Ω φk f J kFk,−1 κ(J k, φk f ) μ f Fk,−t∇δ p · ∇q − ∫ Ω ρ f J k ∂�(pk, ck p) ∂ pk δ p q, a5(δφ f , ψ) := ∫ Ω δφ f ψ, c2(q, δcp) := − ∫ Ω ρ f J k ∂�(pk, ck p) ∂ck p δcp q, F4(q) := ∫ Ω Rk,n 4 q, F5(ψ) := ∫ Ω Rk 5ψ, c∗ 3(δφ f , q) := ∫ Ω J k μ f Fk,−1 ∂[κ(J k, φk f ) · φk f ] ∂φk f δφ f Fk,−t∇ pk · ∇q + ∫ Ω ρ f J k δφ f Δt q, c4(ψ, δu) := ∫ Ω J kFk,−t : ∇δuψ, (3.3) where Fk = I+∇uk , J k = det Fk , Fk,−1 = (Fk)−1, Fk,−t = (Fk)−t, andRk 3,Ω,Rk 3,Σ ,Rk 5 are vector and scalar residuals coming from the linearisation of (2.2a) and (2.2c) and con- taining terms associated with the previous Newton step k as well as body and traction loads. The terms Rk,n 1 ,Rk,n 2 ,Rk,4 4 are scalar residuals from the linearisation of (2.4a), (2.4b), and (2.2b), respectively; which in addition contain contributions from the previous time step n. 123 92 Page 12 of 40 Journal of Scientific Computing (2022) 92 :92 4 Stability and Solvability of the Linearised Problem 4.1 Definition and Preliminaries For this section we draw inspiration from the stability analysis of a system similar to (2.2), recently proved in [12]; and from the study of the linearised hyperelasticity equations, recently studied in [36]. We restrict the description to neo–Hookean poroelastic solids with material parameters μs, λs , for which the effective first Piola–Kirchhoff stress tensor is Peff = μs(F − F−t) + λs ln(J )F−t. Note that a simplified version of the time-discrete tangent system (3.2) is readily obtained by focusing on the first Newton–Raphson iterate when starting from the following initial guess uk = 0, pk = p0, φk f = φ0, ck p = ck l = 0, for displacement, fluid pressure, nominal porosity and species concentrations, which gives Dtφ k f = 0, Fk = I, J k = 1, ∇ pk = 0. In addition, the residuals can be conveniently rearranged and we rescale the system in an appropriate manner using the parameters ρ f , Δt , α, μs , φ0, and λs . Next we introduce the total pressure ζ := α p − λs div u, (4.1) and the system can be recast in terms of u, ζ, p, cp, cl (that is, using ζ instead of φ f ), as follows ã1(cp, wp) + a1(cp, wp) + β1b∗ 3(p, wp) + β2b∗ 3(ζ, wp) = F̃1(wp) ∀wp ∈ H1(Ω), ã2(cl , wl) + a2(cl , wl) + β1b∗ 5(p, wl) + β2b∗ 5(ζ, wl) = F̃2(wl) ∀wl ∈ H1(Ω), a3(u, v) + b1(ζ, v) = F3(v) ∀v ∈ H1 Γ (Ω), b4(cp, q) ã4(p, q) + a4(p, q) − 1 α b2(δtζ, q) = F̃4(q) ∀q ∈ H1 Σ(Ω), b1(ψ, u) + b2(ψ, p) − a5(ζ, ψ) = 0 ∀ψ ∈ L2(Ω), (4.2) where for a function scalar Xn+1, the expression δt Xn+1 stands for Xn+1−Xn Δt , and where we have renamed the bilinear forms, which are now defined as ã1(cp, wp) = ∫ Ω δt cpwp, a1(cp, wp) = ∫ Ω Dp∇cp · ∇wp, b1(ψ, v) = − ∫ Ω div v ψ, b∗ 3(ζ, wp) = − 1 Δt ∫ Ω ζwp, ã2(cl , wl) = ∫ Ω δt clwl , a2(cl , wl) = ∫ Ω Dl∇cl · ∇wl , a5(ζ, ψ) = 1 λs ∫ Ω ζψ, F̃4(q) = ∫ Ω ( 1 Δt R̃n 4 − R4 ) q, b∗ 5(ζ, wl) = − 1 Δt ∫ Ω ζwl , b2(ψ, q) = α λs ∫ Ω qψ, b4(wp, q) = −1 Δt ∫ Ω wpq, a3(u, v) = 2μs ∫ Ω ε(u) : ε(v), ã4(p, q) = ( 1 + α λs )∫ Ω δt pq, a4(p, q) = ∫ Ω κ μ f ∇ p · ∇q, 123 Journal of Scientific Computing (2022) 92 :92 Page 13 of 40 92 c∗ 3(ζ, q) = −ρ f Δt ∫ Ω ζq, F3(v) = ∫ Ω (ρsb + R3) · v, F̃1(wp) = ∫ Ω ( 1 Δt R̃n 1 + R1 ) wp, F̃2(wl) = ∫ Ω ( 1 Δt R̃n 2 + R2 ) wl . (4.3) Here we have also set tΣ = 0. Note that the linear variational problem (4.2) corresponds to the linear and coupled reaction-diffusion / three-field Biot equations δt cp − 1 Δt (β1 p + β2ζ ) − div(Dp∇cp) = 1 Δt R̃n 1 + R1, δt cl − 1 Δt (β1 p + β2ζ ) − div(Dl∇cl) = 1 Δt R̃n 2 + R2, −div(2μsε(u) − ζ I) = ρsb + R3, − 1 Δt cp + ( 1 + α λs ) δt p − 1 λs δtζ − div ( κ μ f ∇ p ) = 1 Δt R̃n 4 − R4, α λs p − 1 λs ζ − div u = 0, (4.4) where the Ri ’s and R̃i ’s are scalar and vector residuals, and theβi > 0 are additional constants arising from the linearisation of the reaction terms and from the change of variables (4.1). A variant of system (4.4) including advection and nonlinear reaction, has been studied in [77]. We follow that reference and in the remainder of this section we present the solvability, stability, and convergence analysis for the in-time problem associated with (4.2). We will assume that the anisotropic permeability κ and the diffusion matrices Dp, Dl are uniformly bounded and positive definite in Ω . The latter means that, there exist positive constants κ1, κ2, and Dmin i , Dmax i , i ∈ {p, l}, such that κ1|w|2 ≤ wtκw ≤ κ2|w|2, and Dmin i |w|2 ≤ wtDiw ≤ Dmax i |w|2 ∀w ∈ R d . We also recall the following version of Korn’s inequality, valid for all u ∈ H1 Γ (Ω) Ck,1‖u‖21,Ω ≤ ‖ε(u)‖20,Ω, where Ck,1 is a positive constant (see, e.g., [20]). This bound, together with the assump- tions above on the permeability and diffusivity, imply the following coercivity and positivity properties a1(wp, wp) ≥ Dmin p |wp|21,Ω, a2(wl , wl) ≥ Dmin l |wl |21,Ω, a3(v, v) ≥ 2μsCk,1‖v‖21,Ω, a4(q, q) ≥ κ1 μ f ‖q f ‖21,Ω, a5(ζ, ζ ) = λ−1 s ‖ζ‖20,Ω, (4.5) for all v ∈ H1 Γ (Ω), ζ ∈ L2(Ω), wp, wl ∈ H1(Ω), q ∈ H1 Σ(Ω). Moreover, thanks to the structure of the formulation, the bilinear form b1 (which coincides with the usual non- diagonal bilinear form in the Stokes equations) satisfies the following inf-sup condition (see, e.g., [44]): For every ζ ∈ L2(Ω), there exists β > 0 such that sup v∈H1 Γ (Ω) b1(v, ζ ) ‖v‖1,Ω ≥ β‖ζ‖0,Ω . (4.6) 123 92 Page 14 of 40 Journal of Scientific Computing (2022) 92 :92 Finally, we recall an important discrete-in-time identity and introduce the discrete-in-time norm ∫ Ω Xn+1δt Xn+1 = 1 2 δt‖Xn+1‖2 + 1 2 Δt‖δt Xn+1‖2, ‖X‖2 �2(V ) := Δt n∑ m=0 ‖Xm+1‖2V , respectively, which will be useful for the subsequent analysis. 4.2 Fixed-point Scheme and Unique Solvability of the Decoupled Problems With the aim to facilitate the comprehension and clarity of the forthcoming analysis, through this section and Sects. 4.3–4.5, we take up again the notation (·)n+1 given in Sect. 3.1 to denote the current time step of the involved variables. Moreover, in order to apply a fixed-point scheme, we modify the variational formulation (4.2) in an equivalent one: find (cn+1 p , cn+1 l , un+1, pn+1, ζ n+1) ∈ H1(Ω)×H1(Ω)×H1 Γ (Ω)×H1 Σ(Ω)×L2(Ω) such that ã1(c n+1 p , wp) + a1(c n+1 p , wp) = F1,pn+1,ζ n+1(wp), (4.7a) ã2(c n+1 l , wl) + a2(c n+1 l , wl) = F2,pn+1,ζ n+1(wl), (4.7b) a3(un+1, v) + b1(ζ n+1, v) = F3(v), (4.7c) ã4(pn+1, q) + a4(pn+1, q) − 1 α b2(δtζ n+1, q) = F4,cn+1 p (q), (4.7d) b1(ψ, un+1) + b2(ψ, pn+1) − a5(ζ n+1, ψ) = 0, (4.7e) for each (wp, wl , v, q, ψ) ∈ H1(Ω) × H1(Ω) × H1 Γ (Ω) × H1 Σ(Ω) × L2(Ω), where the functionals F1,pn+1,ζ n+1 , F2,pn+1,ζ n+1 , F4,cn+1 p are defined as F1,q,ψ (wp) := F̃1(wp) − β1b∗ 3(q, wp) − β2b∗ 3(ψ,wp), F2,q,ψ (wl) := F̃2(wl) − β1b∗ 3(q, wl) − β2b∗ 3(ψ,wl), F4,wp (q) := F̃4(q) − b4(wp, q), (4.8) respectively. We then start with the fixed-point strategy. Let us define the operator S : H1(Ω) → H1 Σ(Ω) × L2(Ω) as S(cn+1 p ) := (S1(cn+1 p ),S2(cn+1 p )) := (pn+1, ζ n+1) ∀ cn+1 p ∈ H1(Ω), (4.9) where pn+1 ∈ H1 Σ(Ω) and ζ n+1 ∈ L2(Ω) are the second and third components of the solution of problem: Find (un+1, pn+1, ζ n+1) ∈ H1 Γ (Ω) × H1 Σ(Ω) × L2(Ω) such that a3(un+1, v) + b1(ζ n+1, v) = F3(v) ∀v ∈ H1 Γ (Ω), (4.10a) ã4(pn+1, q) + a4(pn+1, q) − 1 α b2(δtζ n+1, q) = F4,cn+1 p (q) ∀q ∈ H1 Σ(Ω), (4.10b) b1(ψ, un+1) + b2(ψ, pn+1) − a5(ζ n+1, ψ) = 0 ∀ψ ∈ L2(Ω), (4.10c) for a given cn+1 p ∈ H1(Ω). In turn, let S̃ : H1 Σ(Ω) × L2(Ω) → H1(Ω) be the operator defined by S̃(pn+1, ζ n+1) := cn+1 p ∀ (pn+1, ζ n+1) ∈ H1 Σ(Ω) × L2(Ω), (4.11) 123 Journal of Scientific Computing (2022) 92 :92 Page 15 of 40 92 where cn+1 p ∈ H1(Ω) is the first component of the solution of the problem: Find (cn+1 p , cn+1 l ) ∈ [H1(Ω)]2 such that ã1(c n+1 p , wp) + a1(c n+1 p , wp) = F1,pn+1,ζ n+1(wp) ∀wp ∈ H1(Ω), (4.12a) ã2(c n+1 l , wl) + a2(c n+1 l , wl) = F2,pn+1,ζ n+1(wl) ∀wl ∈ H1(Ω), (4.12b) for a given pair (pn+1, ζ n+1) ∈ H1 Σ(Ω) × L2(Ω). Therefore, we define the map T : H1(Ω) → H1(Ω) as T(cn+1 p ) := S̃(S1(cn+1 p ),S2(cn+1 p )) ∀ cn+1 p ∈ H1(Ω), (4.13) and one readily realises that solving (4.7a)–(4.7e) is equivalent to seeking a fixed point of the solution operator T. Now, in order to prove that the operator T is well-defined, we first need to investigate whether the uncoupled problems defined by the operators S and S̃ are well-posed. The fact that the poroelastic problem is well-posed for a given cn+1 p is given next. Its proof can be deduced by employing the Fredholm alternative approach and classical tools frequently used for showing the well-posedness of elliptic/parabolic equations. We refer to [77, Lemmas 2.1, 2.2 and 2.3] for further details. Lemma 4.1 For a given cn+1 p ∈ H1(Ω), the problem (4.10a)–(4.10c) has a unique solution (un+1, pn+1, ζ n+1)∈ H1 Γ (Ω)×H1 Σ(Ω)×L2(Ω). Moreover, there exists C > 0, independent of λs , such that for each n, ‖un+1‖21,Ω + ‖pn+1‖20,Ω + ‖ζ n+1‖20,Ω + ‖p‖2 �2(H1(Ω)) ≤ exp(C) {‖u0‖21,Ω + ‖p0‖20,Ω + ‖ζ 0‖20,Ω + ‖R3‖20,Ω + ‖b‖20,Ω + Δt‖R4‖20,Ω + n∑ m=0 ‖R̃m 4 ‖20,Ω + n∑ m=0 ‖cm+1 p ‖20,Ω } . (4.14) In turn, the well-posedness of the uncoupled diffusion-reaction problem is given next. Lemma 4.2 For fixed pn+1 ∈ H1 Σ(Ω) and ζ n+1 ∈ L2(Ω), the system of nonlinear parabolic Eqs. (4.12a)–(4.12b) has a unique solution (cn+1 p , cn+1 l ) ∈ [H1(Ω)]2, that is continuously dependent on data, that is, there exists C > 0, such that for each n, ‖cn+1 p ‖20,Ω + ‖cn+1 l ‖20,Ω + ‖∇cp‖2�2(L2(Ω)) + ‖∇cl‖2�2(L2(Ω)) ≤ exp(C) { ‖c0p‖20,Ω + ‖c0l ‖20,Ω + n∑ m=0 (‖R̃m 1 ‖20,Ω + ‖R̃m 2 ‖20,Ω) + Δt‖R1‖20,Ω + Δt‖R2‖20,Ω + n∑ m=0 (‖ζm+1‖20,Ω + ‖pm+1‖20,Ω) } . (4.15) Proof For each n, the unique solvability of (4.12a)–(4.12b) can be obtained by using the coercivity of the diffusivities Di i ∈ {p, l}, and classical tools for time-discrete approxi- mation schemes of parabolic equations (see for instance [64]). Regarding stability, we focus first on the equation for cn+1 p . Consider wp = cn+1 p in (4.12a), to get ∫ Ω δt c n+1 p cn+1 p + ∫ Ω Dp∇cn+1 p · ∇cn+1 p = ∫ Ω ( 1 Δt R̃n 1 + R1 ) cn+1 p + β1 1 Δt ∫ Ω ζ n+1cn+1 p 123 92 Page 16 of 40 Journal of Scientific Computing (2022) 92 :92 +β2 1 Δt ∫ Ω pn+1cn+1 p , and then, applying the property (4.5), and the classical Cauchy–Schwarz inequality, we obtain 1 2 δt‖cn+1 p ‖20,Ω + 1 2 Δt‖δt c n+1 p ‖20,Ω + Dmin p ‖∇cn+1 p ‖20,Ω ≤ 1 Δt ‖R̃n 1‖0,Ω‖cn+1 p ‖0,Ω + ‖R1‖0,Ω‖cn+1 p ‖0,Ω + 1 Δt ‖β1ζ n+1 + β2 pn+1‖0,Ω‖cn+1 p ‖0,Ω . Applying Young’s inequality, it is possible to get 1 2 δt‖cn+1 p ‖20,Ω + 1 2 Δt‖δt c n+1 p ‖20,Ω + Dmin p ‖∇cn+1 p ‖20,Ω ≤ 1 2Δt ‖R̃n 1‖20,Ω + ( 1 Δt + 1 2 ) ‖cn+1 p ‖20,Ω + 1 2 ‖R1‖20,Ω + β1 2Δt ‖ζ n+1‖20,Ω + β2 2Δt ‖pn+1‖20,Ω, and summing over n − 1 and multiplying by Δt , we finally deduce that 1 2 ‖cn p‖20,Ω + 1 2 Δt2 n−1∑ m=0 ‖δt c m+1 p ‖20,Ω + Dmin p Δt n−1∑ m=0 ‖∇cm+1 p ‖20,Ω ≤ ‖c0p‖0,Ω + 1 2 n−1∑ m=0 ‖R̃m 1 ‖20,Ω + (1 + Δt 2 ) n−1∑ m=0 ‖cm+1 p ‖20,Ω + Δt 2 n−1∑ m=0 ‖R1‖20,Ω + β1 2 n−1∑ m=0 ‖ζm+1‖20,Ω + β2 2 n−1∑ m=0 ‖pm+1‖20,Ω . (4.16) The same result can be derived for cn l . Finally, the estimate (4.15) follows by adding the aforementioned result to (4.16) and then, applying Gronwall’s inequality, and using the resulting approach for n + 1. �� 4.3 Existence of aWeak Solution to the Coupled Linearised System Let us define a closed subset of the Banach space L2(Ω) as K := {cn+1 p ∈ L2(Ω) : ‖cn+1 p ‖1,Ω ≤ r}, (4.17) for a given r > 0. We then restrict the space in (4.13) to the ball (4.17), and therefore, we look for a cn+1 p ∈ K such that T(cn+1 p ) = cn+1 p . In what follows, we prove the existence of that cn+1 p by means of the Schauder fixed-point approach. We start with a preliminary result. Lemma 4.3 For a given cn+1 p ∈ H1(Ω), the problem defined by the operator S (cf. (4.9)) satisfies the following estimate ‖S(cn+1 p )‖0,Ω ≤ ‖pn+1‖1,Ω + ‖ζ n+1‖0,Ω ≤ C { ‖R3‖0,Ω + ‖bm+1‖0,Ω + ‖R4‖0,Ω + ‖R̃n 4‖0,Ω + ‖ζ n‖0,Ω } + C̃‖cn+1 p ‖0,Ω, (4.18) 123 Journal of Scientific Computing (2022) 92 :92 Page 17 of 40 92 where C and C̃ are positive constants depending on Δt, μ f , μs, κ1, β, α, with C̃ := (3αμ f ((2Δt2κ1)min{μs/Δt, ακ1/2μ f , β 2})−1)1/2. (4.19) Proof The proof starts by multiplying the Eqs. (4.10a), (4.10b) and (4.10c) by 1 Δt , α, and − 1 Δt , respectively. Then, the estimate (4.18) follows after adding these resulting equations with v := un+1, q := pn+1 andψ := ζ n+1, respectively, using Cauchy–Schwarz and Young inequalities. We omit further details. �� Now, we are in a position to establish the existence of a fixed point of the operator T. This result is abridged in the following two lemmas. Lemma 4.4 Assume that C̃Ĉ ≤ 1 4 , and Ĉ { ‖cn p‖0,Ω + ‖R̃n 1‖0,Ω + ‖R1‖0,Ω +C { ‖R3‖0,Ω + ‖bm+1‖0,Ω + ‖R4‖0,Ω + ‖R̃n 4‖0,Ω + ‖ζ n‖0,Ω }} ≤ r 2 , where C̃, Ĉ and C are the constants specified in (4.19), (4.22) and (4.18), respectively. Then, the fixed-point operator T maps from K into itself. Proof We begin by obtaining an energy estimate for the problem defined by S̃ (cf. (4.11)). In fact, by taking wp = cn+1 p in (4.12a), we readily see that 1 Δt ‖cn+1 p ‖20,Ω + Dmin p ‖∇cn+1 p ‖20,Ω ≤ 1 Δt ∫ Ω cn+1 p cn p + ∫ Ω ( 1 Δt R̃n 1 + R1)c n+1 p + β1 1 Δt ∫ Ω ζ n+1cn+1 p + β2 1 Δt ∫ Ω pn+1cn+1 p , for a given pair (pn+1, ζ n+1) ∈ H1 Σ(Ω) ×L2(Ω), from which, we straightforwardly obtain ‖̃S(pn+1, ζ n+1)‖1,Ω = ‖cn+1 p ‖1,Ω ≤ Ĉ { ‖cn p‖0,Ω + ‖R̃n 1‖0,Ω + ‖R1‖0,Ω + ‖pn+1‖0,Ω + ‖ζ n+1‖0,Ω } , (4.20) and therefore, by using the definition ofT (cf. (4.13)), and applying the result given by (4.20), we get ‖T(cn+1 p )‖1,Ω := ‖̃S(S1(cn+1 p ),S2(cn+1 p ))‖1,Ω ≤ Ĉ { ‖cn p‖0,Ω + ‖R̃n 1‖0,Ω + ‖R1‖0,Ω +‖S1(cn+1 p )‖0,Ω + ‖S2(cn+1 p )‖0,Ω } , (4.21) where Ĉ := (Δt)−1 max{1, (Δt)−1, β1, β2} · min{(Δt)−1, Dmin p }−1. (4.22) Finally, the desired result follows after substituting the estimate (4.18) into (4.21), and then, applying the assumption on data given at the statement of the present lemma. �� 123 92 Page 18 of 40 Journal of Scientific Computing (2022) 92 :92 We remark that there is no inconvenience with the second assumption on data given in Lemma 4.4 since it depends on a constant r that can be appropriately chosen (cf. (4.17)). Moreover, the first smallness assumption given inLemma4.4 ismerely a theoretical condition to guarantee solvability of the continuous problem. However, even though for the constants C̃ and Ĉ defined by (4.19) and (4.22), respectively, the condition might not always be satisfied (especially for extreme values of the model constants), the numerical experiments obtained in Sect. 6 show satisfactory results even even without verifying the theoretical constraints. The present analysis could be improved by using a different approach that circumvents the current restrictions on data, but we are not aware of such a strategy being applicable directly to this context. Lemma 4.5 The map T : K → K is continuous and T(K) is relatively compact in L2(Ω). Proof We notice from the previous lemma, that T(K) is bounded in H1(Ω). In this way, the compact embedding of H1(Ω) into L2(Ω) together with boundedness of T(K) conclude that T(K) is relatively compact in L2(Ω). On the other hand, for the continuity property, we let {cn+1 p,k }k ∈ K be a sequence such that cn+1 p,k → cn+1 p in L2(Ω) as k → ∞, and as is known from (4.13), cn+1 p,k = T(cn+1 p,k ). In this way, proceeding as in [77, Lemma 2.8], we obtain that cn+1 p,k → T(cn+1 p ) in L2(Ω) as k → ∞, concluding the proof. �� Finally, the following lemma concerning the existence of solution for problem (4.7a)– (4.7e), is merely an application of Lemmas 4.4 and 4.5, and the Schauder fixed-point theorem. Lemma 4.6 The semi-discrete in-time formulation (4.7a)–(4.7e) possesses at least one solu- tion. 4.4 Uniqueness ofWeak Solutions We establish the uniqueness of weak solutions through the following result. Lemma 4.7 The semi-discrete weak formulation (4.7a)–(4.7e) has a unique solution. Proof We begin by defining two solutions (un+1 1 , pn+1 1 , ζ n+1 1 , cn+1 p,1 , cn+1 1,l ) and (un+1 2 , pn+1 2 , ζ n+1 2 , cn+1 p,2 , cn+1 l,2 ) associatedwith initial data b1, u01, p01, ζ 0 1 , c 0 p,1, c0l,1, R1,1, R̃n 1,1, R2,1, R̃n 2,1, R3,1, R4,1, R̃n 4,1 and b2, u 0 2, p02, ζ 0 2 , c0p,2, c0l,2, R1,2, R̃n 1,2, R2,2, R̃n 2,2,R3,2,R4,2, R̃4,2, respec- tively, and then Un+1 = un+1 1 − un+1 2 , Pn+1 = pn+1 1 − pn+1 2 , χ n+1 = ζ n+1 1 − ζ n+1 2 , Cn+1 p = cn+1 p,1 − cn+1 p,2 , Cn+1 l = cn+1 l,1 − cn+1 l,2 . Therefore, taking advantage of the linearity of the involved forms, choosing v = δtUn+1 in (4.7c), and then applying Cauchy–Schwarz and Young inequalities, and multiplying the resulting inequality by Δt and summing over n − 1, we finally obtain 123 Journal of Scientific Computing (2022) 92 :92 Page 19 of 40 92 μsCk,1‖Un‖21,Ω + μsCk,1Δt2 4 n−1∑ m=0 ‖δtUm+1‖21,Ω ≤ C { ‖U0‖21,Ω + n−1∑ m=0 ‖χm+1‖20,Ω + n−1∑ m=0 ‖b1 − b2‖20,Ω + n−1∑ m=0 ‖R3,1 − R3,2‖20,Ω } , (4.23) where C is a constant independent of Δt . On the other hand, by taking q = Pn+1 and ψ = χ n+1 in (4.7d) and (4.7e), respectively, and then, applying again the linearity of the bilinear forms and proceeding as in [77, Lemma 2.2], we arrive at the following estimate 1 2λs ‖χ n‖20,Ω + 1 2 ( 1 + α λs )( ‖Pn‖20,Ω + Δt2 n−1∑ m=0 ‖δtPm+1‖20,Ω ) + κ1Δt 2μ f n−1∑ m=0 ‖Pm+1‖21,Ω ≤ 1 2λs ‖χ 0‖20,Ω + 1 2 ( 1 + α λs ) ‖P0‖20,Ω + ( (1 + α)2 2λ + 1 ) n−1∑ m=0 ‖Pm+1‖20,Ω + 1 2 n−1∑ m=0 ‖R̃m 4,1 − R̃m 4,2‖20,Ω + μ f Δt 2κ1 n−1∑ m=0 ‖R4,1 − R4,2‖20,Ω + 1 2 n−1∑ m=0 ‖Cm+1 p ‖20,Ω + 1 2μsCk,1 ‖χ n+1‖20,Ω + μsCk,1 2 ‖Un+1‖21,Ω + 1 μCk,1 n−1∑ m=0 ‖χm+1‖20,Ω + μsCk,1Δt2 4 n−1∑ m=0 ‖δtUm+1‖21,Ω . (4.24) Now, for the diffusion-reaction problem,we apply the linearity of the bilinear forms in (4.7a)– (4.7b), with test functions Cp and Cl , respectively, and then, proceeding similar to Lemma 4.1, we get ‖Cn p‖20,Ω + ‖Cn l ‖20,Ω + Δt2 n−1∑ m=0 (‖δtCm+1 p ‖20,Ω + ‖δtCm+1 l ‖20,Ω) + DminΔt n−1∑ m=0 (|Cm+1 p |21,Ω + |Cm+1 l |21,Ω) ≤ C ( ‖C0p‖20,Ω + ‖C0l ‖20,Ω + n−1∑ m=0 ( ‖R̃m 1,1 − R̃m 1,2‖20,Ω + ‖R̃m 2,1 − R̃m 2,2‖20,Ω ) + Δt 2 n−1∑ m=0 (‖R1,1 − R1,2‖2 +‖R2,1 − R2,2‖20,Ω )+ ( 1 + Δt 2 ) n−1∑ m=0 ( ‖Cm+1 p ‖20,Ω + ‖Cm+1 l ‖20,Ω ) + n−1∑ m=0 ‖χm+1‖20,Ω + n−1∑ m=0 ‖Pm+1‖20,Ω ) , (4.25) 123 92 Page 20 of 40 Journal of Scientific Computing (2022) 92 :92 where Dmin := min{Dmin p , Dmin l }. Finally, thanks to the inf-sup condition (4.6), it is possible to obtain a bound for ‖χ n+1‖0,Ω independent of λs , and then, combining (4.23)–(4.25), and then using Gronwall’s lemma, it can be deduced that ‖Un+1‖1,Ω + ‖Pn+1‖0,Ω + ‖χ n+1‖0,Ω + ‖Cn+1 p ‖0,Ω + ‖Cn+1 l ‖0,Ω + ‖P‖l2(H1(Ω)) + ‖∇Cp‖l2(L2(Ω)) + ‖∇Cl‖l2(L2(Ω)) ≤ C √ exp(C1) ( ‖U0‖1,Ω + ‖P0‖0,Ω + ‖χ 0‖0,Ω + ‖C0p‖0,Ω + ‖C0l ‖0,Ω + n∑ m=0 ‖b1 − b2‖0,Ω + n∑ m=0 ( ‖R̃m 1,1 − R̃m 1,2‖2 + ‖R̃m 2,1 − R̃m 2,2‖2 ) + Δt‖R1,1 − R1,2‖20,Ω + Δt‖R2,1 − R2,2‖20,Ω + ‖R3,1 − R3,2‖20,Ω + n∑ m=0 ‖R̃m 4,1 − R̃m 4,2‖20,Ω + Δt‖R4,1 − R4,2‖20,Ω ) , from which, we can ensure the existence of at most one weak solution to the system (4.7a)– (4.7e). �� 4.5 Stability of the Linearised Coupled Problem The following lemma establishes the continuous dependence on data for problem (4.7a)– (4.7e). Its proof follows similar arguments to those used in Lemmas 4.1, 4.2 and 4.7, and therefore we omit it here. Lemma 4.8 The solution (cn+1 p , cn+1 l , un+1, pn+1, ζ n+1) ∈ H1(Ω) × H1(Ω) × H1 Γ (Ω) × H1 Σ(Ω) × L2(Ω) of problem (4.7a)–(4.7e) satisfies ‖un+1‖1,Ω + ‖pn+1‖0,Ω + ‖ζ n+1‖0,Ω + ‖p‖�2(H1(Ω)) + ‖cn+1 p ‖0,Ω + ‖cn+1 l ‖0,Ω + ‖∇cp‖�2(L2(Ω)) + ‖∇cl‖�2(L2(Ω)) ≤ √exp(C) { ‖u0‖1,Ω + ‖p0‖0,Ω + ‖ζ 0‖0,Ω + ‖c0p‖0,Ω + ‖c0l ‖0,Ω + ‖b‖0,Ω + n∑ m=0 ( ‖R̃m 1 ‖0,Ω + ‖R̃m 2 ‖0,Ω + ‖R̃m 4 ‖0,Ω ) + ‖R3‖0,Ω + Δt‖R1‖0,Ω + Δt‖R2‖0,Ω + Δt‖R4‖0,Ω } , where C > 0 is a constant independent of λs . Remark 4.1 – Although carrying out the solvability analysis of the linearised fully- continuous problem goes beyond the scope of this work, it is possible to extend the ideas presented here to establish that well-posedness using appropriate choices of Sobolev spaces, and performing a passage to the limit adapting to our problem the results from, e.g., [5] (which focus on reaction–diffusion–Brinkman problems). 123 Journal of Scientific Computing (2022) 92 :92 Page 21 of 40 92 – The solvability analysis of the fully-discrete problem associated with (4.7a)–(4.7e) can be established similarly to the semidiscrete-in-time case. More precisely, with the finite elements spaces specified in (5.1), we can define an appropriate fixed-point operator and prove that it is well-defined thanks to the well-posedness of each uncoupled problem (see a general form in, e.g., [64]). Then, by employing the continuous dependence in the fully-discrete case (which follows exactly as in Lemma 4.8) in combination with the analysis from [5, Sect. 5.3], we can prove the continuity of the fixed-point operator described above. Finally, the desired result follows from an application of Brouwer’s fixed–point theorem. – Given the approximation properties of the finite element spaces specified in (5.1) (and recalled in, e.g., [60]), and following the steps given by [77, Sect. 4], it is possible to derive an asymptotic O(h) convergence for the proposed method. This is also verified for the nonlinear case, as shown numerically in Sect. 6. – The convergence properties are also robust with respect to incompressibility of the solid phase. – We stress that the arguments presented in this section are unfortunately not readily generalised to arbitrary hyperelastic materials, nor a generic initial guess for the Newton– Raphson loop. In all of the numerical tests performed below, this has not shown to be a problem in practice. Still, the lack of a global existence result makes it is possible for the model to be increasingly difficult in the simulation of other scenarios of interest, such as pathological ones. 5 A Finite Element Formulation 5.1 Discretisation of the Nonlinear Problem In order to define a Galerkin finite element method we denote by {Th}h>0 a shape-regular family of partitions of Ω̄ , conformed by tetrahedra (or triangles in 2D) K of diameter hK , with mesh size h := max{hK : K ∈ Th}. Given an integer k ≥ 1 and a subset S of R d , d = 2, 3, by Pk(S) we will denote the space of polynomial functions defined locally in S and being of total degree up to k. Let us also denote by bK := ϕ1ϕ2ϕ3 a P3 bubble function in K , where ϕ1, ϕ2 , ϕ3 are the barycentric coordinates of the triangle K (in 3D the bubble is a quartic polynomial). Then the finite-dimensional subspaces for the pathogen and leukocyte concentrations Wh ⊆ H1(Ω), displacement Vh ⊆ H1 Γ (Ω), porous fluid pressure Qh ⊆ H1 Σ(Ω), and nominal porosity Φh ⊆ L2(Ω) are defined, respectively, as follows Wh := {wh ∈ C(Ω) : wh |K ∈ P1(K )d , ∀K ∈ Th}, Vh := {vh ∈ C(Ω) : vh |K ∈ [P1(K ) ⊕ span{bK }]d ∀K ∈ Th, vh |Γ = 0}, Qh := {qh ∈ C(Ω) : qh |K ∈ P1(K ), ∀K ∈ Th, qh |Σ = 0}, Φh := {ψh ∈ C(Ω) : ψh |K ∈ P1(K ), ∀K ∈ Th}. (5.1) The pair (Vh, Φh) is the well-known MINI-element, which is inf-sup stable in the context of saddle-point Stokes equations in their velocity-pressure formulation [28]. Then the fully discrete problem arises from (3.2) and for each time step n, we perform inner Newton–Raphson iterations from k = 0, . . . seeking (δck+1 p,h , δck+1 l,h , δuk+1 h , δ pk+1 h , δφk+1 f ,h ) ∈ 123 92 Page 22 of 40 Journal of Scientific Computing (2022) 92 :92 Wh × Wh × Vh × Qh × Φh =: Hh solutions to the unsymmetric matrix system ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ A1 B′ 1 B′ 2 0 B′ 3 B1 A2 B′ 4 0 B′ 5 0 0 A3 C′ 1 0 C2 0 C1 A4 C′ 3 0 0 C4 0 −A5 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ δCk+1 p δCk+1 l δU k+1 δPk+1 δΦk+1 f ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ Fk,n 1,h Fk,n 2,h Fk 3,h Fk,n 4,h Fk 5,h ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (5.2) where the entries δCk+1 p , δCk+1 l , δU k+1, δPk+1 and δΦk+1 f in the independent vector vari- able are the vectors containing all internal degrees of freedom associated with the discrete incremental solutions for all fields, and the operators in calligraphic letters appearing in the coefficient matrix and load vector from (5.2) are induced by the corresponding bilinear forms and linear functionals in (3.3). The functionalswith superscript n on the right-hand side vector indicate that they also receive contributions from the backward Euler time-discretisation. 5.2 Schur Complement Based Robust Preconditioner The numerical solution of problem (5.2) through direct methods is not feasible for large systems, which is a natural consequence of considering finer meshes for 3D geometries, needed to better capture anatomical details. This is the standard scenario in which Krylov space methods are the most useful, more specifically a GMRES method due to the indefinite and non-symmetric nature of the problem [69]. Still, scalable solvers require the construction of a robust preconditioner that allows for the Krylov iterations to remain roughly constant whenever (i) the amount of processors is increased and (ii) the mesh is refined. In this section we propose a preconditioner based on a two-level nested Schur complement inspired by the multiphysics nature of the model. If we consider a block matrix M given by the general structure M = [ A B1 B2 C ] , with A invertible, a Gauss elimination procedure yields M = [ I 0 B2A−1 I ] [ A 0 0 C − B2A−1B1 ] [ I A−1B1 0 I ] . (5.3) We note that ifC is invertible, the same procedure can be appliedwith respect to it. Schur com- plement based preconditioners enjoy excellent theoretical properties, as the preconditioned system possesses at most 3 distinct eigenvalues [57], implying that it converges in at most 3 iterations of aKrylov subspacemethod, but this is seldom true in practice as the Schur comple- ment blockS = C−B2A−1B1 is computationally expensive to compute. Two complementary strategies to circumvent this problem are (i) to consider the full/lower triangular/upper trian- gular/diagonal part of (5.3) and (ii) to approximate S. In turn, two standard approaches for approximating S are the SIMPLE preconditioner given by S ≈ C−B2diag (A)−1B1 and the block diagonal approximation S ≈ C (for other possibilities and their comparison see, e.g., the comprehensive review [35]). On the first level of the proposed preconditioner, we split the variables into poroelastic and chemotaxis, denoted by poro = (u, p, φ f ) and chem = (cp, cl) respectively, and use 123 Journal of Scientific Computing (2022) 92 :92 Page 23 of 40 92 this to write the linearised problem as J = [ Jporo Jporo, chem Jchem, poro Jchem ] , (5.4) where each block is given by Jporo = ⎡ ⎣ A3 C′ 1 0 C1 A4 C′ 3 C4 0 −A5 ⎤ ⎦ , Jchem = [ A1 B′ 1 B1 A2 ] , Jporo,chem = ⎡ ⎣ 0 0 C2 0 0 0 ⎤ ⎦ , Jchem,poro = [ B′ 2 0 B′ 3 B′ 4 0 B′ 5 ] . Using (5.3) we can write the inverse of (5.4) as J−1 = [ I −J−1 poroJporo,chem 0 I ] [ J−1 poro 0 0 (Jchem − Jchem,poroJ−1 poroJporo,chem)−1 ] [ I 0 −Jchem,poroJ−1 poro 0 ] . Still, the application of this preconditioner requires applying J−1 poro, for which we consider an additional splitting between displacement and the pair (pressure, nominal porosity) denoted as u − Π , with Π = (p, φ f ). We can apply the same argument to the block Jporo, now written as Jporo = [ Ju Ju,Π JΠ,u JΠ ] , by making use of the invertibility of the displacement block Ju, where the sub-blocks are given by Ju = A3, JΠ = [ A4 C′ 3 0 −A5 ] , Ju,Π = [ C′ 1 0 ] , JΠ,u = [ C1 C4 ] . This idea has been applied in [33] for an FSI problem, and in [48] for the Boussinesq equations. The equations governing (cp, cl) are of parabolic type and so their approximation does not pose a challenge, which indicates that the main difficulty lies in the poroelastic block approximation. Motivated by this, the proposed preconditioner is given by the following components: – On the first level we consider a lower-triangular Schur decomposition and on the second level we consider a full one. – The “chem” block is approximated by the action of the HYPRE-BoomerAMG [81] preconditioner. – Both blocks in the second level are approximated by an inexact GMRES solver precon- ditioned by an additive Schwarz method with a direct solver in each subdomain (using the MUMPS library [4]). The inexactness is given by a relative tolerance of 0.1, and we highlight that an inexact solver in the sub-blocks gives rise to a preconditioner that changes from one iteration to the next one. This requires the use of a flexible GMRES algorithm (fGMRES) [68]. The use of an additive Schwartz preconditioner, or in general a domain decomposition one, alleviates 123 92 Page 24 of 40 Journal of Scientific Computing (2022) 92 :92 the deteriorating scalability of AMG preconditioners for higher-order elements, in this case required by the displacement to satisfy an appropriate discrete inf-sup condition. We close this section by noting that the choice of the type of Schur decomposition is not arbitrary. For this, we consider exact, lower triangular Schur complement preconditioners at both levels. In this way, the preconditioner adopts the following form P = [ Pporo 0 Jchem,poro Schem ] , where Pporo = [ Ju 0 JΠ,u SΠ ] , SΠ = JΠ − JΠ,uJ−1 u Ju,Π , Schem = Jchem − Jchem,poroJ−1 poroJporo,chem. After algebraic manipulations, it can be seen that the preconditioned system is given by P−1J = [ P−1 poroJporo P−1 poroJporo,chem S−1 chemJchem,poro(I − P−1 poroJporo) S −1 chem(Jchem − Jchem,poroP−1 poroJporo,chem) ] ,(5.5) where P−1 poroJporo = [ I J−1 u Ju,Π 0 I ] . The preconditioned system does not exhibit the expected property of having a small number of eigenvalues. This is mainly due to the use of a preconditioner instead of an exact inverse, where Pporo �= Jporo. If such blocks coincide, then we would have P−1J = [ I J−1 poroJporo,chem 0 I ] , (5.6) which justified the use of an exact Schur preconditioner at the second level. This results in a preconditioned problem with only one eigenvalue. 6 Computational Experiments We now turn to the presentation of numerical examples serving to illustrate the performance of the finite element scheme and to examine further the main features of the model. All routines have been implemented using the open source finite element library FEniCS [1]. A fixed tolerance of 10−6 is used on the residuals for the convergence criterion of the Newton– Raphson iterative algorithm. We highlight that the error control of multiphysics problems is not trivial, for example in [16] as an error they use the norm of the increment, whereas in [80] they use the norm of the residual. We prefer the norm of the residual, as having an increment going to zero might mean stagnation and not convergence. Still, this can be improved by checking the error in each physics separatedly. 6.1 Sensitivity Analysis of Model Parameters This section presents the influence of some parameters on the dynamics of the proposed model to highlight the coupling mechanisms between the poroelastic component (2.2) and 123 Journal of Scientific Computing (2022) 92 :92 Page 25 of 40 92 the advection-reaction-diffusion component (2.4) describing pathogens and immune system dynamics, respectively. The analyses were performed in a one dimensional version of the fully coupled model given by Eqs. (2.2)–(2.4) considering a domain Ω ∈ [0, 8] cm. An initial pathogen concentration of cp,0 = 0.001 was considered in the middle of the domain x ∈ [3.95, 4.05] to start the infection, whereas in the remaining domain cp,0 = 0 was used. Throughout the domain we also consider the following initial conditions: cl,0 = 0.003, φ0 = 0.2, and p0 = 0. The following boundary conditions were considered: u = 0 on the left at x = 0, and p = 0 was prescribed at the right end at x = 8cm of the domain. A simplified sensitivity analysis varying one parameter at a time with respect to the refer- ence parameters reported in Table 1 was performed.We limit the presentation of these results to only a few parameters and variables divided into three cases that highlight the couplings of the model. Young’s modulus E , associated with the mechanical part, was investigated for the first case. In the second and third cases, associated with the immune response of the model, both phagocytosis rate, λlp , and pathogen’s reproduction rate, γP , were studied. These parameters were chosen because they better highlighted the couplings in the model. They are also related to critical biological responses depending on their values. Changes in the value of γP , for example, can be associated with pathogens with distinct proliferation abilities; changes in the value of λlp can describe the ability of leukocytes in dealing with the invading pathogen; whereas changes in E represent tissues with different stiffness. In all the three cases, the investigated parameter (E , λlp , or γP ) had its value doubled and halved with respect to the reference value reported in Table 1. To study these three scenarios, we assume that pathogens are responsible for triggering the inflammatory response. For this reason, the instant in which pathogens’ concentration reaches its peak is used as a reference to collect data from the other variables of interest. Figure 2 shows the results of the sensitivity analysis. The time in which pathogens con- centration reaches its maximum was distinct for each scenario only for the analysis of the influence of pathogen’s reproduction rate (cf. Fig. 2, third line), occurring after 12 days, six days, and 24 days, respectively, since the start of infection. For the analysis of λlp and λpl , the instant of time at which the pathogens reach the maximum was 12 days since the start of infection. The plots in the first row of Fig. 2 presents the results for varying tissue stiffness. A stiffer tissue will imply a lower capacity to deform and accumulate fluid. With a smaller amount of fluid in the infected region, the concentration of pathogens and leukocytes in the tissue increases (Fig. 2, top left), i.e., the same amount of pathogens in the region is divided by a smaller fluid volume. The second row of panels in Fig. 2 presents the results for varying the phagocytosis rate, λlp . A reduction in λlp value can be related to a scenario in which the immune system is inef- fective in dealing with the invading pathogen. This could be explained by immunodeficiency, for example. In this context, the concentration of pathogens will reach a higher and faster peak of infection for smaller values of λlp . The lower capacity of phagocytosis by leukocytes leads to a considerable increase in the concentration of pathogens (Fig. 2, centre left). The increase in pathogens’ concentration induces an increase in the capillary permeability, which allows a more significant accumulation of fluid in the infected region (Fig. 2 centre middle). The accumulation of fluid in the region causes the pressure to increase. As a consequence of these events, the deformation suffered by the tissue is greater than the one observed for higher values of λlp (Fig. 2 centre right). Note that an increase in γP can be interpreted as a pathogen that can reproduce faster due to its intrinsic biological characteristics or due to the presence of favourable conditions (abundance of nutrients or temperature, for example). In fact, as the bottom row of Fig. 2 123 92 Page 26 of 40 Journal of Scientific Computing (2022) 92 :92 Ta bl e 1 E xa m pl e 1: R ef er en ce pa ra m et er s of th e m od el to be us ed in Se ct .6 .1 Pa ra m et er U ni ts D es cr ip tio n R ef er en ce s Pa ra m et er U ni ts D es cr ip tio n R ef er en ce s E = 60 kg / cm s2 Y ou ng m od ul us E st im at ed λ lp = 1. 5 1/ d c Ph ag oc yt os is ra te E st im at ed ν = 0. 35 − Po is so n co ef fic ie nt E st im at ed λ̄ lp = λ lp / φ 0 1/ d c R el at iv e ph ag oc yt os is ra te λ s kg / cm s2 Fi rs tL am é pa ra m et er λ pl = 7. 1 1/ d c L eu ko cy te s m ig ra tio n ra te E st im at ed μ s kg / cm s2 Sh ea r m od ul e λ̄ pl = λ pl / φ 0 1/ d c R el at iv e le uk oc yt es m ig ra tio n ρ f = 1 kg / cm 3 Fl ui d ph as e de ns ity E st im at ed π i = 10 m m H g In te rs tit ia lo nc ot ic pr es su re [6 1] ρ s = 2e −3 kg / cm 3 So lid ph as e de ns ity [1 1] π c = 20 m m H g C ap ill ar y on co tic pr es su re [6 1] α = 0. 25 − B io tm od ul us E st im at ed σ 0 = 0. 91 − O sm ot ic re fle ct io n co ef fic ie nt [6 1] φ 0 = 0. 2 − In iti al flu id ph as e [1 3] L bp = 1e 4 1/ c Pa th og en in flu en ce on pe rm ea bi lit y E st im at ed D p = 1e −3 cm 2 / d Pa th og en di ff us io n E st im at ed P c = 20 m m H g C ap ill ar y pr es su re [6 1] D̄ p = D p / φ 0 cm 2 / d R el at iv e pa th og en di ff us io n L p0 = 3. 6e −8 cm / s m m H g H yd ra ul ic pe rm ea bi lit y [6 1] D l = 5e −2 cm 2 / d L eu ko cy te di ff us io n E st im at ed � 0 = 6. 82 e− 5 1/ s N or m al ly m ph at ic flo w E st im at ed D̄ l = D l/ φ 0 cm 2 / d R el at iv e le uk oc yt e di ff us io n k m = 6. 5 m m H g H al f lif e of ly m ph at ic flo w E st im at ed X = 1e −2 cm 2 / d c C he m ot ax is E st im at ed n H il l = 1 − H ill co ef fic ie nt E st im at ed X̄ = X / φ 0 cm 2 / d c R el at iv e ch em ot ax is V m a x = 20 0 − M ax ly m ph at ic flo w E st im at ed γ p = 1. 2e −1 1/ d Pa th og en re pr od uc tio n ra te E st im at ed K = 2. 5e −7 cm 2 / s m m H g Pe rm ea bi lit y [6 1] γ̄ p = γ p / φ 0 1/ d R el at iv e pa th og en re pr od uc tio n ra te (S / V ) = 17 4 1/ cm V es se la re a pe r vo lu m e un it [6 1] 123 Journal of Scientific Computing (2022) 92 :92 Page 27 of 40 92 Fig. 2 Example 1: The interaction of immune response and mechanical parts of the model is evaluated using Young’s modulus (E , first line), phagocytosis rate (λlp , second line) and pathogen reproduction rate (γP , last line) at the time instant in which the peak pathogen concentration is reached. For each of these three parameters, three scenarios are evaluated and plotted as distinct curves in each graphic: E = 60 (scenario 1, solid line), E = 120 (scenario 2, dashed line) and E = 30 (scenario 3, dotted line); λlp = 1.5 (scenario 1, solid line), λlp = 3.0 (scenario 2, solid line) and λlp = 0.75 (scenario 3, solid line); γP = 0.12 (scenario 1, solid line), γP = 0.24 (scenario 2, dashed line), γP = 0.06 (scenario 3, dotted line). The graphics in the left and in the middle columns are composed of two ordinate y−axes. In the first column, the left y−axis represents the concentration of pathogens and the right y−axis the concentration of leukocytes. The second column presents the pressure field and the fraction of fluid phase, respectively, on the left and right y−axes. The third column represents the deformation field on a single y−axis. All figures represent, on the x−axis, the tissue size in centimetres shows, the peak of infection was higher and faster for larger γP values, as expected. How- ever, counter-intuitively, high values of γP are associated with small deformations and small regions of oedema. As a result of more pathogens, the capillary permeability increases to allow more leukocytes to enter into the tissue (Fig. 2, bottom left). The rapid increase in the concentration of leukocytes contained the infection quickly, preventing pathogens from remaining in the tissue long enough to spread. Since the time pathogens stay in the tissue is short, as well as the region in which they spread, the deformation suffered by the tissue also decreases (Fig. 2, bottom right). In the third scenario, with the small value of γP = 0.06, the peak of infection only occurs 24 days after the invasion begins, which leads to a weak inflammatory response. Thus, the concentration of leukocytes is reduced. The low leukocyte concentration results in a long time to remove pathogens that can spread over a larger region. The weak immune response results in more deformation when compared to the cases with faster pathogen dynamics, which highlights the couplings and strong nonlinearities in the model. 123 92 Page 28 of 40 Journal of Scientific Computing (2022) 92 :92 Fig. 3 Example 2: Compression and drainage of a porous structure. Approximate displacement magnitude, Lagrangian porosity, and fluid pressure, on the deformed configuration and computed at t = 0.5 with ν = 0.33 (top) and ν = 0.495 (bottom) For this test we have also used the FEniCS environment along with the direct solver MUMPS. 6.2 Example 2: Compression and Drainage of a Poroelastic Sample With the aim of illustrating the coupling between poroelastic deformation and fluid flow using the formulation (2.2a)–(2.2c) and in the absence of the immune system chemotaxis, we simulate the incremental compression of a box Ω = (0, 1)2 m2. The bottom edge of the boundary constitutes Γ (where according to (2.6a), the deformation is zero and we impose zero fluid flux), whereas Σ is conformed by the top edge (where a sinusoidal-in-time dis- tributed traction tΣ = 2000 sin(π t)n is applied and a constant fluid pressure pin = 0.2MPa is considered) as well as the vertical walls (where we prescribe zero traction and zero fluid pressure). We employ a non-homogeneous initial Lagrangian porosity that is not symmetric with respect to the centre of mass of the domain, given by φ0(x) = 3 5 + 1 10 sin(x1x2), the isotropic power-law porosity-dependent permeability, a compressible neo–Hookean strain energy density yielding the nominal effective stress Seff = μs(Ft − F−1) + λs lnJF−1, and the parameters E = 104 kg/ms2, ν ∈ [0.2, 0.499999], μs = E 2(1 + ν) , λs = Eν (1 + ν)(1 − 2ν) , b = 0, � = 0, ρs = 2 · 10−3 kg/m3, ρ f = 10−3 kg/m3, α = 0.25, μ f = 10−3 m2/s, κ0 = 10−5 m2. 123 Journal of Scientific Computing (2022) 92 :92 Page 29 of 40 92 Fig. 4 Example 2: Compression benchmark. Approximate displacement components, Lagrangian porosity, and fluid pressure, on the deformed configuration and computed at t = 34 (top) and t = 41 (middle), and transients of field variables at points A and B (bottom) The scheme characterised by the finite element spaces in (5.1) is used on a crossed-shaped uniform mesh with 10,000 cells, and the coupled system is simulated until tfinal = 1s using a constant time step Δt = 0.1s. The system is considered initially at rest (u(0) = 0, p(0) = 0, φ f (0) = φ0). This test also serves to assess the robustness of the method. We run a set of simulations varying the Poisson ratio from 0.2 to 0.499999. Regardless of the parameter regime, the numerical results do not show spurious oscillations of pressure, nor unexpectedly small displacements or nonphysical distortions or checker-board patterns in the porosity field. In Fig. 3 we plot deformed geometries and field variables showing the progressive compression of the porous block over three time instants, and using ν = 0.33 and ν = 0.495. As expected, the fluid pore pressure does not have a uniform distribution, and it moves toward Γ on the bottom edge. In addition, we perform a benchmark test where again a time-harmonic loading defined by tΣ = 1 5 (λs + 2μs) sin( 2π15 t)n is applied, now only on part of the top edge of the box Ω = (0, 8) × (0, 5)m2 (on the segment 0 ≤ x1 ≤ 1m, x2 = 5m). The setup of this validation example proceeds similarly as in [49, Sect. 4.2] (see also [84, Sect. 4.2] and [51, 65]). The boundary conditions in these references permit drainage (fixing a zero fluid pressure) on the whole top lid, setting zero fluid flux elsewhere on the boundary, the bottom edge is clamped (imposing zero displacement) while the vertical walls are on rollers (only the horizontal displacement is set to zero), and the remainder of the top edge is traction-free. The parameter values that change with respect to the previous case are E = 3 · 104 N/m2, ν = 0.2, α = 1, κ0 = 10−8 m2, φ0 = 0.33, tfinal = 100 s, Δt = 1 s, and we now use the normalised and isotropic Kozeny–Carman porosity-dependent perme- ability. For this test the mesh is unstructured and graded near the top-left corner. Figure 4 shows snapshots of the approximated field variables at times t = 4, 12, 16s and the bot- 123 92 Page 30 of 40 Journal of Scientific Computing (2022) 92 :92 Fig. 5 Example 3: Error decay with respect to the number of degrees of freedom, computed between approxi- mate solutions computed with a first-order method and manufactured solutions (6.1), at the final adimensional time tfinal = 0.03 tom panels also portray transients of the vertical displacement and porosity recorded at two locations (node A ≈ (0.5,4.5), directly below the centre of the footing strip and node B ≈ (8.0,4.5) approximately at the same height but on the right wall). These plots show the same qualitative behaviour encountered in, e.g., [51]. 6.3 Example 3: Errors with Respect to Manufactured Solutions Since there is no analytical solution currently available for the coupled system (2.2a)–(2.4b), the accuracy of the finite element discretisation is investigated using the following closed- form manufactured solutions cp = t[0.3 exp(x1) + 0.1 cos(πx1) cos(πx2)], cl = t[0.3 exp(x1) + 0.1 sin(πx1) sin(πx2)], u = 0.25t ⎛ ⎝ sin(πx1) cos(πx2) + x21 λs − cos(πx1) sin(πx2) + x22 λs ⎞ ⎠ , p = t sin(πx1x2) cos(πx1x2), φ0 = 0.6 + 0.1 sin(x1x2), (6.1) and φ f = J − 1 + φ0, defined on the unit square domain Ω = (0, 1)2. Together with the strong form of the governing equations, these exact solutions are used to obtain the load and source right-hand side terms. For simplicity we impose Dirichlet boundary conditions for the concentrations cp, cl whereasmixed boundary conditions are considered for the displacement and fluid pressure. On the left side of the boundary Γ , we prescribe the exact displacement from (6.1) and an exact traction Pn on the remainder of the boundary, Σ . The fluid pressure is constrained to match the exact one on Σ and an exact fluid pressure flux is imposed on Γ . 123 Journal of Scientific Computing (2022) 92 :92 Page 31 of 40 92 We use the same neo–Hookean material law as in Example 1, and employ the isotropic Kozeny–Carman permeability κKC iso . The following parameter values are selected λs = 36.4, μs = 22.1, ρs = ρ f = μ f = κ0 = L p0 = Lbp = Lbr = vmax = γp = γl = cmax l = 1, α = 0.5, Dp = 0.9I, Dl = 0.8I, σ0 = λlp = λpl = μl = l0 = p0 = πi = km = πc = χ = pc = S/V = 1, n = 2, they are all regarded adimensional and do not have physical relevance in this case, as we will be simply testing the convergence of the finite element solutions. We generate successively refined simplicial grids and use a sufficiently small (non dimen- sional) time step Δt = 0.01 and simulate a sufficiently short time horizon tfinal = 3Δt , to guarantee that the error produced by the time discretisation does not dominate. Errors between the approximate and exact solutions are plotted against the number of degrees of freedom in Fig. 5. This error history plot confirms the optimal convergence of the finite element scheme (in this case, first order) given by the Remark 4.1, for all variables in their respective norms, where a slightly better rate, of about 1.3, is seen for the porosity. In addition, the average iteration count for the Newton–Raphson method (for each refinement level and for all time steps) is five. 6.4 Example 4: Localisation of Oedema Regions In the next test we present and discuss the coupling dynamics of the model (2.2)–(2.4), considering a two-dimensional domain Ω = (0, 4) × (0, 4) cm2, and a time step size of Δt = 0.1 days. The parameters described in the Table 1 were used. Boundary conditions were applied as follows: u = 0 on the left edge, and p = 0 for the right, top and bottom edges of the domain, and the following initial conditionswere used: cp,0 = 0.001in the region(x1− 2)2 + (x2 −2)2 ≤ 0.03 and cp,0 = 0.0, otherwise. Also, we consider cl,0 = 0.003, φ0 = 0.2 and p0 = 0 in the entire domain. Figure 6 presents the results of a localised oedema formation. Each row presents snapshots of the variables cp , cl , p, u andφ f , respectively, at three selected time instants. The dynamics of the immune response can be observed in Fig. 6 in terms of cp and cl . After the pathogens appear in the tissue, they start to grow (Fig. 6 at time t = 8 days) and reach their peak concentration (Fig. 6 at time t = 13 days). However, in response to these events, leukocytes migrate to the infected site and their concentration also grows (Fig. 6, cl column). As a consequence of the presence of leukocytes, pathogen concentration starts to decrease, as Fig. 6 shows for cl at time t = 18 days. The presence of pathogens triggers a mechanical response due to a sequence of couplings in the model. The endothelium increases its permeability to allow leukocytes to leave the bloodstream and enter the tissue. This increase in the endothelium permeability in turn increases interstitial fluid and pressure (Fig. 6 for φ f and p at time t = 8 days). When the concentration of pathogens reaches its peak (Fig. 6 for cp at time t = 13 days), it is possible to observe that the amount of fluid phase φ f in the tissue also reaches its maximum value. An increase of liquid in the region leads to an increase in both pressure and displacement fields (Fig. 6 for u and p at time t = 13 days). Finally, when pathogens cp are almost eliminated by the leukocytes at time t = 18 days, the dynamics of pressure, fluid phase fraction and displacements, tend to return to their initial values. 123 92 Page 32 of 40 Journal of Scientific Computing (2022) 92 :92 Fig. 6 Example 4. Behaviour of pathogen (first column) and leukocyte (second column) concentrations, pressure field (third column), fraction of fluid phase (fourth column) and displacement field (fifth column) at time t = 8 (top row), t = 13 (middle row) and t = 18 (bottom row), after the solution of the model (2.2a)–(2.2c) and (2.4a)–(2.4b), considering the parameters given by Table 1 and initial conditions given by cp,0 = 0.001, cl,0 = 0.003, φ0 = 0.2, u0 = 0, and p0 = 0 6.5 Example 5: Immune SystemDynamics and Poro-hyperelasticity in a Left Ventricle Next we conduct a series of tests using a ventricular geometry segmented from patient- specific images [78], and where synthetic muscle fibre and collagen sheet directions are constructed using a Laplace–Dirichlet rule-based method [66]. On the basal surface we prescribe zero normal displacement and zero fluid pressure flux, the epicardium is considered stress-free and with a prescribed fluid pressure, and on the endocardium we impose a time- dependent traction tΣ = m0 sin2(π t)n with m0 = 0.1N/cm2, and an endocardial fluid pressure pendo = p0 sin2(π t) having the same time period. No-flux conditions are used for the concentrations on the whole boundary. For this example we use the Holzapfel–Ogden constitutive strain energy stated in (2.1), and employ the anisotropic Kozeny–Carman permeability κKC introduced in (2.3); while the remaining model parameters used for the 3D ventricular test assume the values λs = 27.293 kPa, μs = 3.103 kPa, Dp = 3 × 10−3/φ0 cm 2/h I πi = 10mmHg, n = 5, α = 0.5, Dl = 5 × 10−2/φ0 cm 2/h I vmax = 20, σ0 = 0.91, φ0 = 0.2, p0 = 10.9mmHg, km = 6.5mmHg, πc = pc = 20mmHg, S/V = 174, χ = 10−4 cm3/(h.107cell), γp = 0.13/φ0 cm 3/(h.107cell), λpl = 7.1/φ0, λlp = 1.8/φ0, l0 = 0 cm/s, Lbp = 5000, 123 Journal of Scientific Computing (2022) 92 :92 Page 33 of 40 92 Fig. 7 Example 5. Assigned fibre distribution streamlines, displacementmagnitude, fluid pressure distribution, and porosity on a patient-specific left ventricular geometry at t = 14 (top), and evolution of pathogens and leukocytes concentration (middle and bottom rows, respectively, seen from a slightly different angle) L p0 = 3.6 × 10−8, μ f = 10−3 cm2/s, κ0 = 2.5 × 10−7 cm2, ρs = 2 × 10−3 Kg/cm3, ρ f = 10−3 Kg/cm3, a = 0.496N/cm2, b = 0.041, a f = 0.193N/cm2, b f = 0.176, as = 0.123N/cm2, bs = 0.209, a f s = 0.162N/cm2, b f s = 0.166. An initial concentration of pathogens is considered on a transmural strip, and we simulate the dynamical behaviour of the coupled system over several minutes. The evolution of the pathogens and leukocytes distribution on the deformed ventricular geometry is shown in Fig. 7, where we also show fibre distribution and mechanical fields at half-time. One can observe that the presence of pathogens induces the local entry of leukocytes. Fluid also enters the interstitial space of the tissue, increasing pressure. Although leukocytes can remove part of the pathogens, some of them diffuse through almost the entire domain. As this pathogens’ wave sweeps a significant part of the left ventricle, more parts of the heart are impacted by the inflammatory response and, as a consequence, diffuse oedema is formed. Numerical simulations show that the propagation of the front of pathogens likely depends on the diffusion 123 92 Page 34 of 40 Journal of Scientific Computing (2022) 92 :92 Table 2 Example 5. Iteration counts, total and (average), for fGMRES with the proposed preconditioner, and for either continuous and piecewise linear approximations of displacement, or piecewise linear displacements with bubble enrichment DoFs 1 cpu 2 cpu 4 cpu 8 cpu 16 cpu 32 cpu (a) Displacement approximated with P1 elements 28924 48 (6.86) 88 (12.57) 81 (11.57) 87 (12.43) 80 (11.43) 90 (12.86) 155351 74 (6.73) 139 (12.64) 122 (12.20) 103 (12.88) 89(12.71) 93 (13.29) (b) Displacement approximated with P b 1 elements 72772 74 (6.73) 139 (12.64) 122 (12.20) 103 (12.88) 89 (12.71) 93 (13.29) 438158 50 (6.25) 88 (11.00) 83 (11.86) 85 (10.62) 77 (11.00) 80 (11.43) Note that an increase in the total iterations, for constant average iterations, occurs due to a different number of Newton iterations and replication of pathogens and the wave tail likely depends on the diffusion and efficiency of the leukocytes [52]. The numerical results indicate that, although fluid phase and pressure increase due to the presence of the pathogen, they result in small changes in displacement. In fact, new diagnostic tools have been studied to increase the non-invasive diagnostic accuracy of myocarditis and other myocardial pathologies using for this purpose the detection of increased water contents at some regions instead of an increase in tissue displacement. See, for example, [47, 72, 73]. To conclude this section, for the ventricular geometry we proceed to briefly test the pre- conditioner described in Sect. 5.2. To measure its performance, we look at the total and average number of Krylov iterations in two different meshes corresponding to the left ven- tricle geometry and with respect to P1 and P b 1 (bubble enriched element) conforming finite elements for the displacement to study the impact of the inf-sup stability in the performance. The results are shown in Table 2, and as we have observed that numerically the most difficult time instant is the first one, we report the performance only on it. We highlight that the number of average Krylov iterations is roughly constant with an increasing number of cores. This can be expected due to the use of inexact solvers in the sub-blocks of the Schur complement preconditioner, which yield an adequate approximation of the exact Schur preconditioner, which has at most 3 distinct eigenvalues. It is particularly interesting to observe that this performance is obtained by using only the action of an AMG preconditioner for the chemotaxis block, meaning that (i) the coupling between the chemo- taxis and poroelastic models is not very strong in passive cardiac simulations and that (ii) such block is not necessarily computationally challenging. 7 Final Remarks We have studied a general model capturing the phenomenological features of the interaction between chemotaxis of the immune system in saturated poroelastic media admitting large deformations. The essential properties of the new model include a formulation in terms of displacement, fluid pressure, and Lagrangian porosity, coupled with concentrations for leukocytes and pathogens. In particular, the poro-hyperelastic compartment of the model can be identified, after linearisation and adequately choosing the initial guess, with the three-field formulation for Biot’s poroelasticity from [60, 77].We have proposed a finite elementmethod together with splitting strategies, and have studied their properties in terms of accuracy and 123 Journal of Scientific Computing (2022) 92 :92 Page 35 of 40 92 iteration count. The realisation of the coupling is general enough to accommodate different types of model specifications, including diverse hyperelastic solid laws and showing the interplay of solid deformation, effective stress, and pore fluid pressure build-up. The model also includes species transport through the total velocity, and the system is solved by means of a monolithic finite element method in saddle-point form. Some extensions in the development of this work include the study of long-term behaviour of the oedema formation, as well as performing a thorough exploration of different effects from introducing lymphatic sinks for cp and cl . We also plan to incorporate active stress due to calcium release (in turn generated by the immune reaction-diffusion) and assess quantitative differences in the classical diffusion case. In a more applicative context, a more comprehensive sensitivity analysis is still required to determine the most relevant coupling mechanisms. Also, the validation and verification of the model against patient-specific data are still to be conducted. We have expanded the application discussed in Sect. 6.5 to address further modelling and numerically oriented investigations in the recent work [52]. Regarding numerical schemes, a possible next step is to use mixed formulations for the immune system equations (following, e.g., [43]) to produce mass conservative discretisations. Further investigation is necessary, for instance, regarding the specific form of the anisotropic porosity as well as in designing new coupling mechanisms that will contribute to a better understanding of the formation and termination of myocarditis and myocardial oedema. Funding Open Access funding enabled and organized by CAUL and its Member Institutions. This work has been supported by Universidade Federal de Juiz de Fora (UFJF) through the scholarship “Coordenação de Aperfeiçoamento de Pessoal de Nível Superior” (CAPES) -Brazil-Finance Code 001; by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq)-Brazil Grant numbers 423278/2021-5, 308745/2021-3 and 310722/2021-7; by Fundação de Amparo à Pesquisa do Estado de Minas Gerais (FAPEMIG)-Brazil CEX APQ 02830/17, TEC APQ 03213/17, and TEC APQ 01340/18; by the Monash Mathematics Research Fund S05802-3951284; by the Australian Research Council through the Discovery Project Grant DP220103160; and by the Ministry of Science and Higher Education of the Russian Federation within the framework of state support for the creation and development ofWorld-Class ResearchCentres “Digital biodesign and personalised healthcare” No. 075-15-2020-926. Data Availibility Part of the datasets and finite element implementations generated during and/or anal- ysed during the current study are available from the GitHub repository https://github.com/ruizbaier/ PoroelasticModelForAcuteMyocarditis. Declarations Conflict of interest The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. 123 https://github.com/ruizbaier/PoroelasticModelForAcuteMyocarditis https://github.com/ruizbaier/PoroelasticModelForAcuteMyocarditis http://creativecommons.org/licenses/by/4.0/ 92 Page 36 of 40 Journal of Scientific Computing (2022) 92 :92 References 1. Alnæs, M.S., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., Richardson, C., Ring, J., Rognes, M.E., Wells, G.N.: The FEniCS project version 1.5. Arch. Numer. Softw. 3, 9–23 (2015) 2. Alves, J.R., de Queiroz, R.A.B., Bär, M., dos Santos, R.W.: Simulation of the perfusion of contrast agent used in cardiac magnetic resonance: a step toward non-invasive cardiac perfusion quantification. Front. Physiol. 10, 177 (2019) 3. Alves, J.R., de Queiroz, R.A.B., Weber dos Santos, R.: Simulation of cardiac perfusion by contrast in the myocardium using a formulation of flow in porous media. J. Comput. Appl. Math. 295, 13–24 (2016) 4. Amestoy, P.R., Duff, I.S., L’Excellent, J.-Y., Koster, J.: MUMPS: a general purpose distributed memory sparse solver, In: International Workshop on Applied Parallel Computing, Springer, pp. 121–130, (2000) 5. Anaya, V., Bendahmane, M., Mora, D., Ruiz-Baier, R.: On a vorticity-based formulation for reaction- diffusion-brinkman systems. Netw. Heterog. Media 13, 69–94 (2018) 6. Arnold, D.N., Brezzi, F., Fortin, M.: A stable finite element for the Stokes equations. Calcolo 21, 337–344 (1984) 7. Ateshian, G.A., Weiss, J.A.: Anisotropic hydraulic permeability under finite deformation. J. Biomech. Eng. 132, 111004(7) (2010) 8. Auricchio, F., Beirão da Veiga, L., Lovadina, C., Reali, A.: A stability study of somemixed finite elements for large deformation elasticity problems. Comput. Methods Appl. Mech. Eng. 194, 1075–1092 (2005) 9. Bærland, T., Lee, J.J., Mardal, K.-A., Winther, R.: Weakly imposed symmetry and robust preconditioners for Biot’s consolidation model. Comput. Methods Appl. Math. 17, 377–396 (2017) 10. Barnafi, N., De Oliveira Vilaca, L.M., Milinkovitch, M.C., Ruiz-Baier, R.: Coupling chemotaxis and poromechanics for the modelling of feather primordia patterning, In preparation, pp. 1–28, (2021) 11. Barnafi, N., Di Gregorio, S., Dede’, L., Zunino, P., Vergara, C., Quarteroni, A.M.: A multiscale porome- chanics model integrating myocardial perfusion and systemic circulation, MOX Reports (2021) 12. Barnafi, N., Zunino, P., Dedè, L., Quarteroni, A.: Mathematical analysis and numerical approximation of a general linearized poro-hyperelastic model. Comput. Math. Appl. 91, 202–228 (2021) 13. Basser, P.J.: Interstitial pressure, volume, and flow during infusion into brain tissue. Microvasc. Res. 44, 143–165 (1992) 14. Berger, L., Bordas, R., Burrowes, K., Grau, V., Tavener, S., Kay, D.: A poroelastic model coupled to a fluid network with applications in lung modelling. Int. J. Numer. Methods Biomed. Eng. 32, e02731 (2016) 15. Berger, L., Bordas, R., Kay, D., Tavener, S.: A stabilized finite element method for finite-strain three-field poroelasticity. Comput. Mech. 60, 51–68 (2017) 16. Borregales,M., Radu, F.A., Kumar, K., Nordbotten, J.M.: Robust iterative schemes for non-linear porome- chanics. Comput. Geosci. 22, 1021–1038 (2018) 17. Both, J., Borregales, M., Nordbotten, J., Kumar, K., Radu, F.: Robust fixed stress splitting for Biot’s equations in heterogeneous media. Appl. Math. Lett. 68, 101–108 (2017) 18. Both, J., Kumar, K., Nordbotten, J., Radu, F.: The gradient flow structures of thermo-poro-visco-elastic processes in porous media, arXiv e-prints (2019) 19. Both, J.W., Barnafi, N.A., Radu, F.A., Zunino, P., Quarteroni, A.: Iterative splitting schemes for a soft material poromechanics model. Comput. Methods Appl. Mech. Eng. 388, e114183 (2022) 20. Brenner, S.C., Scott, L.R.: The mathematical theory of finite element methods. Texts in Applied Mathe- matics, Springer, New York (2002) 21. Burger, R.L., Belitz, K.: Measurement of anisotropic hydraulic conductivity in unconsolidated sands: a case study from a shoreface deposit, Oyster, Virginia. Water Resources Res. 33, 1515–1522 (1997) 22. Burtschell, B., Chapelle, D.,Moireau, P.: Effective and energy-preserving time discretization for a general nonlinear poromechanical formulation. Comput. Struct. 182, 313–324 (2017) 23. Chamberland, E., Fortin, A., Fortin,M.: Comparison of the performance of some finite element discretiza- tions for large deformation elasticity problems. Comput. Struct. 88, 664–673 (2010) 24. Chapelle, D., Gerbeau, J.-F., Sainte-Marie, J., Vignon-Clementel, I.E.: A poroelastic model valid in large strains with applications to perfusion in cardiac modeling. Comput. Mech. 46, 91–101 (2010) 25. Chapelle, D., Moireau, P.: General coupling of porous flows and hyperelastic formulations - from thermo- dynamics principles to energy balance and compatible time schemes. Eur. J. Mech., B/Fluids 46, 82–96 (2014) 26. Cherubini, C., Filippi, S., Gizzi, A., Ruiz-Baier, R.: A note on stress-driven anisotropic diffusion and its role in active deformable media. J. Theor. Biol. 430, 221–228 (2017) 27. Choo, J.: Large deformation poromechanics with local mass conservation: an enriched Galerkin finite element framework. Int. J. Numer. Methods Eng. 116, 66–90 (2019) 123 Journal of Scientific Computing (2022) 92 :92 Page 37 of 40 92 28. Cioncolini, A., Boffi, D.: The MINI mixed finite element for the Stokes problem: An experimental investigation. Comput. Math. Appl. 77, 2432–2446 (2019) 29. Colli Franzone, P., Pavarino, L.F., Scacchi, S.:Mathematical Cardiac Electrophysiology, vol. 13. Springer, Berlin (2014) 30. Cookson,A., Lee, J.,Michler, C., Chabiniok, R.,Hyde, E.,Nordsletten,D., Sinclair,M., Siebes,M., Smith, N.: A novel porous mechanical framework for modelling the interaction between coronary perfusion and myocardial mechanics. J. Biomech. 45, 850–855 (2012) 31. Coussy, O.: Poromechanics. Wiley, New York (2004) 32. De Oliveira Vilaca, L.M., Gómez-Vargas, B., Kumar, S., Ruiz-Baier, R., Verma, N.: Stability analysis for a new model of multi-species convection-diffusion-reaction in poroelastic tissue. Appl. Math. Model. 84, 425–446 (2020) 33. Deparis, S., Forti, D., Grandperrin, G., Quarteroni, A.: FaCSI: A block parallel preconditioner for fluid- structure interaction in hemodynamics. J. Comput. Phys. 327, 700–718 (2016) 34. Ehret, A.E., Bircher, K., Stracuzzi, A., Marina, V., Zündel, M., Mazza, E.: Inverse poroelasticity as a fundamental mechanism in biomechanics and mechanobiology. Nature Comm. 10, e1002 (2017) 35. Elman, H., Howle, V.E., Shadid, J., Shuttleworth, R., Tuminaro, R.: A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier–Stokes equations. J. Comput. Phys. 227, 1790–1808 (2008) 36. Farrell, P.E., Gatica, L.F., Lamichhane, B.P., Oyarzúa, R., Ruiz-Baier, R.: Mixed Kirchhoff stress - dis- placement - pressure formulations for incompressible hyperelasticity. Comput. Methods Appl. Mech. Eng. 374, e113562 (2021) 37. Fede