A unified modelling and simulation for coupled anomalous transport in porous media and its finite element implementation

This paper presents an unified mathematical and computational framework for mechanics-coupled “anomalous” transport phenomena in porous media. The anomalous diffusion is mainly due to variable fluid flow rates caused by spatially and temporally varying permeability. This type of behaviour is described by a fractional pore pressure diffusion equation. The diffusion transient phenomena is significantly affected by the order of the fractional operators. In order to solve 3D consolidation problems of large scale structures, the fractional pore pressure diffusion equation is implemented in a finite element framework adopting the discretised formulation of fractional derivatives given by Grunwald–Letnikov (GL). Here the fractional pore pressure diffusion equation is implemented in the commercial software Abaqus through an open-source UMATHT subroutine. The similarity between pore pressure, heat and hydrogen transport is also discussed in order to show that it is possible to use the coupled thermal-stress analysis to solve fractional consolidation problems.


Introduction
This paper introduces a unified formalism for coupled anomalous diffusion processes in porous media as well as a simple open-source implementation in FE commercial code. It was noticed that in soft tissue hydrated with water such as articular cartilage [1] and meniscus [2][3][4][5][6] the interstitial water pressurizes when the tissue is loaded. The fluid pressurization has been hypothesized to be a major factor in the load-support mechanism and in the response to friction of these types of tissues [1]. When the physiological load is applied, this load is shared between the solid skeleton and the interstitial fluid. Theoretical studies have demonstrated that the interstitial fluid (i.e. pore pressure) may sustain up B O. Barrera obarrera@brookes.ac.uk 1 to 90% or more of the total applied load particularly in the initial phases of loading [1,7,8].
Experimental studies have revealed that the time-varying response of pore pressure diffusion is anomalous and it is mainly due to variable fluid flow rates which are, in turn, caused by spatially and temporally varying permeability [2]. This is mainly due to the fact that the permeability, hence the rate of fluid flow, often depends on the porosity [9] which varies throughout the tissue [2]. Variations in permeability occur, for example, when the fluid flow impacts the geometry or the micro-structural features such as the configuration of the pores. Experiments on water flow in building materials, sand and zeolite among other minerals, highlighted that the permeability changes during the water flow process as a result of the microstructural rearrangement of grains/pores. Iaffaldano at al. [10] hinted that during compaction of sand, permeability might decrease due to the fact that the fluid carries solid particles which then close some of the pores. Essentially the configuration of the medium, in particular the ratio between closed/open pores, changes during the process. Fluid might be trapped in the medium leading to a slower fluid flow rate. On the contrary, if during the water diffusion process some of the pores open creating conductive microchannels, permeability might increase. Therefore, water can be transported for large distances in a reduced time determining a faster diffusion process. These types of phenomena are well captured by a form of Darcy's law in which fractional operators are involved in the rate of fluid flow [11,12].
It is then fundamental to be able to understand and model the time-varying interstitial fluid pressurization, the order of the fractional operators and its role in load bearing of soft hydrated tissues. It is shown that a common mathematical formulation is able to describe both stress-driven anomalous diffusion, biomechanical processes in soft tissues and a number of multiphysics problems. Coupled models of transient flux (temperature, mass, pore pressure) with deformable solids rely on appropriately considering the interaction among separate physical fields. The physics refers to common types of physical processes, e.g., heat transfer (thermo-), pore water or pore pressure movement (hydro/poro-), concentration field (concentro or diffuso/convecto/advecto-), stress and strain (mechano-), dynamics (dyno-), chemical reactions (chemoor chemico-), electrostatics (electro-), and magnetostatics (magneto-). The mathematics of these phenomena is similar. The differences lay mainly on the coupling terms, and how they are modelled, as they are dictated by the physics of the problem. It is clear that in the case of thermoelasticity and poroelasticity the coupling stress-diffusion in the transport equation is confined to a source/sink term. Different is the case of hydrogen diffusion-mechanics in which the coupling is both in the flux (as there is a term which is related to the gradient of hydrostatic stress) and in a sink term (related to the role of plastic strain).
Here the focus of the paper is to present an unified mathematical and computational framework for mechanicscoupled "anomalous" transport phenomena in porous media and provide the link how the procedure can be applied in a number of multiphysics problems. The anomalous pore pressure diffusion which contain a fractional Darcy's law is derived and implemented in the finite element framework using the discretised formulation given by Grunwald-Letnikov (GL). The coupled poro-mechanics problem is then, for the first time, implemented in Abaqus, through a UMATHT routine, in order to adopt the coupled temperaturedisplacement procedure available in Abaqus to effectively solve diffusion-mechanics problems. The same approach has been used when coupling hydrogen diffusion mechanics problems to highlight some material deterioration phenomena. For example the coupled elasto-plastic response with models of diffusional hydrogen transport provided new insights into a number of hydrogen embrittlement mechanisms [13][14][15].
It is also important to know when a given coupled problem can be solved in an uncoupled manner. The strength of coupling is discussed in terms of deriving dimensionless parameters for thermoelasticity/poroelasticity in the same fashion as [16]. Dimensionless parameters are also here derived for the first time for hydrogen diffusion-elastoplasticity. Summarizing the main novelties of the paper are: • Numerical implementation of the anomalous pore pressure diffusion equation in a finite element framework. • Details of the Abaqus implementation in UMATHT subroutine and comparison with analytical solutions of a fractional boundary value problems. • Derivation of dimensioneless parameters to evaluate the strength of coupling for diffusion mechanics problems The paper is structured as follows. First the coupled transportmechanics problem is porous solid is introduced within the context of poroelasticity. The anomalous pore pressure diffusion equation, needed to solve the poroelastic set of equations, is subsequently derived. The strength of coupling in diffusion-mechanics problem is then discussed in Sect. 3. Section 4 focuses on the numerical implementation of the pore pressure diffusion equation, the similarities with heat and hydrogen transport are clearly identified. Numerical results including uncoupled and fully coupled poromechanics such as consolidation problems are presented in Sect. 5.

Fully coupled transport and mechanical behaviour in deformable porous solids
In the following we recap the main equations needed to solve coupled pore pressure diffusion with mechanics in the context of the theory of poroelasticity. For the sake of simplicity we limit the discussion to linear elastic solid. In appendix ("Appendix A") parallelisms with thermo-elasticity and hydrogen diffusion mechanics are reported. The classical linear model of transient flow and deformation of a homogeneous fully saturated elastic porous medium depends on an appropriate coupling of the fluid pressure and solid stress [17,18]. A change in applied stress produces a change in fluid pressure or fluid mass and a change in fluid pressure or fluid mass is responsible for a change in the volume of the porous material. The coupling term affects only the hydrostatic part of the stress tensor. Considering that stresses are positive when they are tensile and pressure is positive when it is compressive (i.e. tr(σ ) = −3 p, where p is the pore pressure), the stress tensor is written as follows: where λ = K − 2 3G is the Lamé constant, G, K are the shear and bulk modulus, respectively, p is the pore pressure and α is the Biot coefficient. It is useful to express the linear poroelastic constitutive equations in terms of strain: With ν being the Poisson's ratio. It is also important for the following to derive an expression for the volumetric strain or dilatation d by taking the trace of both sides of Eq. 2: Equation 3 can also be written as follows: d being the volumetric/dilatation strain and σ H = σ ii /3 the hydrostatic component of stress.
Equations 1, A.1 and B.1 in "Appendices A and B" respectively have the same structure. These equations govern coupled stress-diffusion phenomena. The diffusion equations of the pore pressure p, temperature θ and hydrogen concentration c L are derived from the same set of equations (i.e. conservation of pore fluid, energy, mass along with the constitutive equations relating fluxes with gradient of p, θ, c L (i.e. Darcy's, Fourier's or modified Fick's laws).

Anomalous pore pressure diffusion equation
The anomalous pore pressure is controlled by a governing diffusion equation. This equation is found by considering conservation of mass for the pore fluid ζ , along with a constitutive equation that relates the fluid flux j p to the pore pressure gradient through non integer order operators. Mass conservation states that the rate of change of the total mass in a volume is equal to the flux j p through the surface ∂ which, in absence of sources or sinks, reads as follows: where n is the outward normal to ∂ . In the following the linear kinematic/small deformation theory is taken into account. Nonlinear kinematic terms in small deformation regime is often associated with partially saturated conditions (Skempton coefficient B<1), which the theory presented in the paper is not suited for. Considering that the theory presented here is limited to the linear kinematic case and applying the divergence theorem on the surface integral, the following equation is derived: For Eq. 5 to be true for any region , the integrand must vanish, i.e. ∂ζ ∂t + ∇ · j p = 0 .
Equation 7 represents conservation of mass for the pore fluid.
The relationship between the fluid-flux vector j p and the pore pressure p is given by a modified version of Darcy's law involving fractional operators, which in the isotropic case and without inertia's term takes the form [11,12]: where k is the permeability and μ the dynamic viscosity. In the following k μ will be indicated as λ β . D β 0 indicates the Caputo's fractional time derivative [19]. The Caputo's fractional derivative of order β of the function ∇ p is defined as: The expression is valid for n − 1 < β < n, and is the Euler's Gamma function. One of the advantage of using the definition of Caputo fractional derivative is that initial conditions in terms of integer order derivatives are obtained. This is important for applications in real physical problems. If α → n the Caputo's derivative is reduced to classical integer order derivatives; moreover the Caputo fractional derivative of a constant is zero as integer order derivative. Given that in the cases considered in this paper the lower bond of integration a = 0, Caputo's fractional derivative of order β a D β t is indicated as D β 0 . The fractional derivative method offers the possibility to model with reduced number of parameters all of the anomalous diffusion behaviours by changing the order of the derivative. The drawback is that it is difficult to link the order of the derivatives with microstructural features. In order to understand if the choice of a macroscopic fractional model is sensible,recently approaches based on Bayesian model selection have been developed [20][21][22] and the error in the parameters estimation quantified [20][21][22][23][24].
It is generally assumed that the variation of the pore fluid content ζ depends linearly on the hydrostatic component of stress σ H (shear stresses have no influence) and on the pore pressure. Hence ζ can be expressed as follows [16]: where B is the Skempton coefficient which essentially gives indication of the increase in fluid pressure of an elastic isotropic porous material subjected to undrained loading. B is related to the drain/undrained bulk moduli: K and K u and the Biot coefficient α as follows: The value of B is always between 0 and 1. When B is 1, the hydrostatic stress or the dilatation strain is completely transferred into changing pore pressure. When B equals to 0 indicates no change in pore pressure after applying the volumetric stress or strain field. Laboratory studies indicate the value of B depends upon the fluid saturated pore volume of the sample [25]. A value of B < 1 indicates a partially saturated media. The theory here considers saturated solid for which B is close to 1. The variation of fluid content with time takes the following form: Substituting Eqs. 12 and 8 into Eq. 7 and we obtain the following fractional partial differential equation: In here the following rule of fractional derivatives is implied [19]: It is important to note that, in case of β = 0, λ β = k μ with dimension of [L] 4 [F] [T ] . In the case of β = 0, λ β has dimension of [L] 4 [F][T ] 1−β . The pore pressure diffusion equation can also be written in terms of strains. In order to do so we start by combining Eqs. 10, 11 and 4 in order to have an expression relating the variation of fluid content depending on dilatation and pore pressure: Differentiating Eq. 15 with respect to time we obtain: Substituting Eqs. 16 and 8 into Eq. 7 we obtain the following fractional partial differential equation in terms of strain:

The strength of coupling in diffusion mechanics problems
When dealing with coupled models of transient flux (temperature, mass, pore pressure) with deformable solids it is important to highlight when these problems can be treated in an uncoupled manner. In other words, it is important to quantify when is appropriate to compute stress/strain fields ignoring the transient flux and viceversa. For this purpose in the following dimensionless parameters are identified for the case of poroleasticity, thermoelasticity and hydrogen diffusion-elasto plasticity.

Poroelastic coupling
Equations 1 and 2 show that the influence of the pore pressure on the volumetric part of the stress or strain tensor is governed by the Biot coefficient α. If α = 0 pore pressure would have no influence on the stress/strain tensor. The Biot coefficient is the fluid volume change induced by bulk volume changes in the drained condition. It is shown that, for a large class of porous media, the value of the Biot coefficient α is significantly greater than zero therefore stresses/strain cannot be computed without computing the pore pressure [16]. Equations 13 and 17 show that the pore pressure is governed by a diffusion-type equation containing an additional coupling term that serves as a source-sink: with r p being a source term. From Eq. 17: (19) or from Eq. 13 These source terms in Eqs. 19 and 20 are the coupling terms relating the pore pressure increment to the increment of the volumetric part of either the stress or strain tensor through the Skempton coefficient B. This effect can be viewed as a transient Skempton-type effect, in the sense that either the change in time of hydrostatic stress or of dilatation strain (i.e. changing in volume) gives rise to a local increase in the pore pressure. Equation 4 relates the dilation strain with hydrostatic stress and pore pressure: Now let us consider two cases: (a) Constant hydrostatic stress: ∂σ H ∂t = 0. This reduces Eq. 21 to: In the case of constant hydrostatic state of stress, a positive change in pore pressure contributes to lowering the dilation strain rate. Substituting this into Eq. 19 we obtain that the magnitude of source terms becomes: MagnitudeSource : This reduces Eq. 21 to: Equation 24 shows that a positive increment of pore pressure causes an increase of the rate of hydrostatic stress. Substituting this into Eq. 20 we obtain that the magnitude of source terms becomes: From Eqs. 23 and 25 the strength of coupling is given by Bα.
If Bα << 1 then it is possible to ignore the source terms in the pore pressure diffusion equation written in term of stress Eq. 13 or strain Eq. 17.
It is important to note that substituting Eq. 25 into Eq. 13 we obtain Therefore if α B << 1 then is it possible to ignore mechanical effects. The pore pressure field can then be calculated by solving the standard uncoupled diffusion equation:

Thermoelastic coupling
The temperature diffusion is governed by a diffusion-type equation given below which contains an additional coupling term that serves as a sink (more details in "Appendix C", i.e. Eqs. C.5 and C.6) : with r θ being a sink term and It is clear that the influence of the temperature on volumetric part of the stress or strain tensor is governed by the thermal expansion coefficient γ as stated by Eqs. A.1) and A.2. If γ = 0 temperature would have no influence on the stress/strain tensor. The coefficient of thermal expansion γ describes how the size of an object changes with a change in temperature. Specifically, it measures the fractional change in size per degree change in temperature at a constant pressure. It is shown that for a wide range of material γ is significantly greater than zero, therefore stresses/strain cannot be computed without computing the temperature [16].
It is possible to relate the dilation strain with hydrostatic stress and temperature (see Eq. A.4 in "Appendix A"): As done in the poroelastic coupling section we can consider two cases: (a) Constant hydrostatic stress: ∂σ H ∂t = 0. This would mean from Eq. 31 that: Equation 32 indicates that a positive increment in temperature generates a positive increment in dilation strain, therefore a change in volume of the solid, as it is expected. Substituting this into Eq. 29 we obtain that the magnitude of source terms becomes: This would mean from Eq. 31 that: Equation 34 shows that, when the volumetric strain is constant, an increase in temperature is compensated by a decrease in hydrostatic stress. Substituting this into Eq. 30 we obtain that the magnitude of source terms becomes: From Eqs. 33 and 35 the strength of coupling is given by ρC v << 1 then it is possible to ignore the sink terms in the temperature diffusion equation.
It is important to note that substituting Eq. 33 into Eq. C.5 we obtain: Therefore, if 9K β 2 θ ρC v << 1 then the temperature field can be calculated by solving the uncoupled diffusion equation for temperature:

Hydrogen diffusion-mechanics coupling
Hydrogen atoms in metals reside either at normal interstitial sites (NILs) or at trapping sites such as dislocations, grain boundaries, carbide/matrix interfaces, microvoids and other defects. We denote the hydrogen concentration in the lattice as c L (number of H atoms per unit volume) and use c X to denote the concentration associated with the hydrogen in the traps. The total concentration of hydrogen is given by: c T = c L + c X . The hydrogen diffusion equation reads as follows ("Appendix D", Eq. D.8): with r c being a source term, c L hydrogen concentration at the lattice sites, j c = − D L c L RT ∇μ hydrogen flux. μ = μ 0 + RT ln c L −σ H V H is the chemical potential. The source terms is expressed as below: where N X is the density of traps which increases with increasing plastic strain. Hydrogen concentration influences both (a) the volumetric part of the stress/strain tensor and (b) plastic deformation.

Analogy between coupled transport equations and numerical implementation
The form of the transport equations in Eqs. 18, 28 and 38 is similar and they can be summarized here below in a generic form: where X is the degree of freedom which is θ for the the heat diffusion, c L for hydrogen transport and p for pore pressure diffusion equations. The similarities of these equations are summarized in Table 1.
The transport equation in Eq. 41 in the frame of finite element method (FEM) is solved step by step with numerical integration following an implicit scheme (Newton-Raphson algorithm). Considering that: Heat flux: J = j`= −K t ∇ · θ Hydrogen flux: Hydrogen source: Applying the divergence theorem and writing the weak form Eq. 43 becomes: δ X is an arbitrary variational field. Considering that: and q = −J · n (46) Equation 44 becomes: Substituting Eq. 42 into Eq. 47 the following relation is derived: Now we indicate with U : : S i state variables.
and source/sink

Numerical Implementation of the fractional pore pressure diffusion equation in the UMATHT subroutine
The pore pressure diffusion equation (Eq. 18) can be solved by substituting in Eq. 48 to obtain: A list of the terms to be coded in UMATHT routine are given below: U p : U p : where: • Source: where: • Flux: The flux term in Eq. 8 contains a fractional derivative. The fractional derivatives must be discretized. To this purpose the Grunwald-Letnikov (GL) fractional derivative is used [19]: For sufficiently small t the GL fractional derivative coalesces with the Caputo's fractional derivative.
• Jacobian: Usually FE codes that use an implicit Newton-Raphson integration scheme allow the time increment to be determined automatically to optimize the run time. The GL formula for evaluating of the fractional derivatives has been derived assuming a constant increment (i.e. the time) and, to the best of our knowledge, a corresponding formulation for a variable increment is not available in the literature; furthermore, the automatic time increment requires the definition of a tolerance criterion, that is difficult to define without knowledge of the elastic and inelastic parts of the strain. For these two reasons we have currently limited ourselves to using this model with a fixed time increment. In order to evaluate the GL derivative the history of strain at each Gauss Point must be stored leading possibly to a considerable amount of memory when analysing large FE models. A number of strategies to overcome this problem have been explored and are discussed in [26].

Uncoupled problem and comparison with analytical solutions
The numerical solution of a simple 1D pore pressure diffusion problem can compare with the analytical solution. A bar made of porous material of length l = 50 mm subjected to an initial pore pressure p = 0 Pa throughout the bar. At time t = 0, the pressure at one end is increased to p = 100 Pa. The governing fractional pore pressure diffusion equation Eq. 13 (σ H = 0) and the boundary conditions are written below: The analytical solution was also found (see "Appendix E") and it reads as follow: and Note that in Eq. 65 in the case of β = 0 the Mittag Leffler function E 1−β,1 become E 1,1 = exp, the exponential function. In this case the solution in Eq. 65 is identical to the classical Fick's diffusion solution [13]:  This phenomenon can be appreciated by considering the characteristic time for diffusion i.e. time necessary for the pore pressure to diffuse along the bar. This is calculated from the third term in Eq. 65 being equal to zero: ]c n sin nπ x l = 0. This equality is reached when the characteristic timet is equal tō t = ( l 2 λ ) 1 1−β . The steady state solution reads: p(x, t) = p 0 + x l ( p l − p 0 ). Figure 1 shows the evolution of the characteristic timē t = ( l 2 λ ) 1 1−β with β. In the case of β = 0 the characteristic time of the standard diffusion case is recovered, i.e. t = ( l 2 λ ). Increasing β (with 0 < β < 1) the characteristic time decreases which means that the transient diffusion phenomenon is faster, therefore is it important to consider the time step to adopt when solving the mathematical problem. In order to capture the transient phenomenon when increasing the value of β, the time step will need to be reduced. Figure 2 shows the evolution of the pore pressure along the bar for β = 0, 0.1, 0.3, 0.5. The Abaqus simulations were conducted using quadratic continuum coupled temperaturedisplacements hybrid element CPE8HT. It can be seen that the full transient solutions for the case of β = 0, 0.1, 0.3, 0.5 at different times are in perfect agreement. The input files and the UMATHT subroutine are attached as electronic supplement material.

Fractional consolidation: comparison with analytical solutions
We solve here a fully coupled poromechanics problem. Let us consider a confined consolidation experiment in which we apply a compressive stress on a cylinder made of porous material saturated with water as pictured in Fig. 3. The cylinder is insulated except at the base through which water can flow out. During the deformation process the applied load is sustained by both the pore pressure and the solid skeleton. As for the Terzaghi's solution, the applied load is firstly sustained by the pore pressure. As fluid flows out of the sample the pore pressure decreases and the solid skeleton starts deforming. We are interested here to model the transient pore pressure diffusion phenomenon and to understand the effect of the order of the fractional derivative β. The test is modelled through the 1D uniaxial strain poroelastic problem. Considering that there is variation of field quantities only in z-direction, the equilibrium equation is given by ∂σ zz ∂z = 0, recalling Eq. 1 adapted for the 1D case in which the only non zero component of strain is zz , we obtain: Equation 68 is then coupled with the pore pressure diffusion equation in Eq. 17 adapted for the 1D uniaxial strain version below: Combining Eqs. 68 and 69 we obtain the following pore pressure diffusion equation: . The boundary conditions are those shown in Fig. 3: A constant compressive stress in the z-direction is applied to the cylinder at z = 0: where −P A is the applied compressive stress. The initial pore pressure is derived for undrained conditions i.e.: The analytical solution in terms of pore pressure reads as follows: where: where E 1−β,1 is the Mittag-Leffler function. In the case of β = 0 the solution identical to the classical Terzaghi's solution in which E 1−β,1 = exp. As done in the previous section, it is interesting to observe that the characteristic time for diffusion for this problem is equal tō t = ( 4h 2 λ ) 1 1−β . This is calculated as the time needed for the transient diffusion process to finish i.e. the third term in Eq. 76 It is possible to understand the role of the order of the fractional derivative β in the transient solution by observing Fig. 4 in which the variation of the characteristic timet with β is shown. It can be noted that, as in the pure diffusion case, the transient phenomenon is faster when the order of the fractional derivative β increases, i.e. the characteristic time decreases. However, compared to the pure diffusion case in Fig. 1, the rate of the decrease in the duration of the transient phenomenon with increasing β is slower in the consolidation problem. Figure 5 shows the evolution of the pore pressure along the cylinder for β = 0, 0.1, 0.5. The Abaqus simulations were conducted using quadratic continuum coupled temperature-displacements hybrid element CPE8HT. The material parameters used in the simulations are summarized in Table 2. In this case we keep the time step and the total time constant (dt = 4.8 × 10 −6 , total time t = 4.8 × 10 −4 ) for the three simulations with the value of β = 0, 0.1, 0.5. Figure 5a-c show that the rate of pore pressure diffusion increases with increasing order of fractional derivative β. This also implies that the rate at which the fluid flows out of the sample during the consolidation problem is higher when the value of β increases.

Conclusions
This paper presents a unified formalism for coupled anomalous diffusion processes in porous media as well as a simple open-source implementation within FE software. Here the anomalous pore pressure diffusion, which contains a Darcy's law involving non-integer order operators, is derived and implemented in the finite element framework using the discretised formulation given by Grunwald-Letnikov (GL). The coupled poro-mechanics problem is then, for the first time, implemented in a UMATHT routine in the commercial software ABAQUS. The similarity of the formalism between coupled temperature-displacement and pore pressure-displacement procedures is used. The unified framework presented here is adaptable for other multiphysics problems ruled by the same set of PDEs such as thermoelasticity and hydrogen diffusion mechanics. The differences among multiphysics problems lay on how the coupling terms are modelled as these are associated with the specific physics of the phenomena. Here we have summarized the fundamental equations for coupling, implementing and solving diffusion equations involved in fractional poroelasticity/thermoelasticty and hydrogen diffusion-mechanics highlighting how the coupling between stress-strain and transport differs in these three cases and the strength of the coupling terms. We have derived dimensionless parameters that allow to understand when is appropriate to simplify and uncouple the problem, i.e. solving the diffusion equation without consider the stress state and viceversa. Numerical simulations of uncoupled and coupled poromechanics problems are successfully compared with analytical solutions of fractional boundary value problems whose details are also given in the appendix.
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://creativecomm ons.org/licenses/by/4.0/.

Appendix A: Thermoelasticity
The classical linear model of heat transfer and deformation of a homogeneous elastic material depends on an appropriate coupling of the temperature and solid stress. A change in applied stress produces a change in temperature field and change in temperature is responsible for deforming the solid. The coupling term affects only the hydrostatic part of the stress tensor. We indicate with θ = [K ] the temperature rise from T 0 to T , θ = T 0 − T .
is the coefficient of (linear) thermal expansion. It is also useful to express the linear thermoelastic constitutive equations in terms of strain: It is also important for the following to derive an expression for the volumetric strain or dilatation d by taking the trace of both side of Eq. A.2: .3 can also be written as follows: d being the volumetric/dialation strain.

Appendix B: Hydrogen diffusion-elasticplastic model
It is quite well known that hydrogen being the smallest element diuses more easily through metals than any other element. Even at low concentration hydrogen often leads to the embrittlement of metals for reasons that are not yet well understood but are certainly related to the speed at which it can diuse to highly stressed regions. Hydrogen atoms move through the metals by normal interstitial site (NILS) diusion or dislocation transport. A build up of hydrogen is responsible for volumetric change of the material (swelling) and a gradient of hydrostatic stress drives hydrogen diffusion. Moreover, the plastic properties of metals (i.e. yield strength) are affected by the presence of hydrogen. Here we consider the hydrogen induced material softening which has been linked to the experimental observations of the effect of hydrogen in enhancing dislocation mobility. We indicate with c L = [mol]/[m 3 ] the hydrogen concentration in the normal interstitial site (NILS) and assume that this is responsible for the swelling of the material. The hydrogen diffusion-elasticity equation can be expressed as: It is also important for the following to derive an expression for the volumetric strain or dilatation d by taking the trace of both side of Eq. A.2: .3 can also be written as follows: d being the volumetric/dilation strain.
• Hydrogen diffusion-plasticity. Following [27] we consider an isotropic Von Mises plasticity model in which the flow stress is a function of the hydrogen content. In particular we model the effect of hydrogen-induced material softening, this is to be viewed as an attempt to describe the experimental observations of the effect of hydrogen on dislocation mobility. The yield condition can be expressed as: σ e being the Von Mises stress and σ y = f ( p , c L ) the yield stress which is function of hydrogen concentration as: where σ H 0 is the initial yield strength in the presence of hydrogen that decreases with increasing hydrogen con- is a monotonically decreasing function of hydrogen concentration at NILS.
shown that the following partial differential equation holds [16]: where c v [J /K gK ] is the specific heat at constant strain. Equation C.5 can be written in terms of dilation strain [16]:

Appendix D: Hydrogen transport equation
Hydrogen atoms move through a metal by normal interstitial site (NILS) diffusion or dislocation transport. Hydrogen atoms reside either at NILS or at trapping sites such as dislocations, grain boundaries, carbide/matrix interfaces, microvoids and other defects. The vast majority of sites are the NILS and the minor fraction of sites are the traps. We denote the hydrogen concentration in the lattice as c L (number of H atoms per unit volume) and use c X to denote the concentration associated with the hydrogen in the traps. The total concentration of hydrogen is given by: c T = c L + c X . Mass conservation states that the rate of change of the total hydrogen concentration in a volume is equal to the flux j c (mol m −2 s −1 ) through the surface ∂ : where n is the outward normal to ∂ . The driving force for diffusion is the chemical potential gradient, the flux j c can be expressed as follows: where D L is the diffusion coefficient for hydrogen. μ is the chemical potential defined as follows: In Eq. (D.5) there are two unknowns: c L and c X . However following Oriani's theory the H concentration in the lattice is in equilibrium with the concentration in the traps. This means that once we know the concentration in the lattice we can calculate the concentration in the traps. Also the concentration in the traps c X is a function of the density of trapping sites N X which increases with ε p : Substituting (D.6) in (D.5) we obtain: where D e f f is the effective diffusivity. Equation D.7 can then be arranged as follows:

Appendix E: Analytical solution of the fractional diffusion equation
The boundary value problem Eq. 64 discussed in Sec. 5.1 with the governing fractional pore pressure diffusion Eq. 13 in the case of σ H = 0 together with Dirichlet boundary conditions is reported below: p(0, t) = p 0 ; p(l, t) = p l (E.3) withλ = K B α λ β Substitute the variable p = v + P E , the solution of Eq. E.3 is of the type: t); (E.5) Which can be summarized in two sets of problems: and ∇ 2 P E = 0 ( E . 1 0 ) The solution of Eq. E.12 is of the type: Regarding Eq. E.9, let's substitute v = D −β 0 u in Eq. E.9: Hence we can write that: We can now solve Eq. E.15 with separation of variables: u = u(x, t) = (t) (x). Equation E.15 then becomes: We need then to find a solution for the following equations: Let's first analyzeλ( ∂ 2 ∂ x 2 ) = −λ 2 which can be written as: From which we derive that: The solution to BVP in Eq. E.3 considering Eq. E.6, i.e. p = v + P E is as follows: v = p 0 + p 0 − p l l x + ∞ n=1 E 1−β,1 − n 2 π 2 l 2λ t 1−β c n sin nπ l x (E.31)