New continuous strain‐based description of concrete's damage‐permeability coupling

This contribution aims at improving the finite element modelling of concrete's damage‐permeability coupling using continuous and elastic‐based damage approaches. With this aim in view, a review of existing local and nonlocal damage formulations is performed showing how they differ in terms of damage estimation, crack propagation, and spatial discretization restrictions (when the energy‐based Hillerborg criterion is applied). Another critical review of damage‐based permeability laws is achieved to pinpoint their limitations. Indeed, the use of a bounded and strictly increasing damage variable does not allow (1) an accurate description of permeability in the presence of macrocracks that are physically unbounded and (2) a proper description of the permeability decrease as those macrocracks are closed. Accordingly, a new strain‐based permeability law is suggested. Its analytical validation is based on the BIPEDE tensile test for one loading‐unloading cycle. This new formulation intends to enhance the predictiveness of damage‐permeability coupling models leading to a better assessment of large structure's durability and serviceability.


Summary
This contribution aims at improving the finite element modelling of concrete's damage-permeability coupling using continuous and elastic-based damage approaches. With this aim in view, a review of existing local and nonlocal damage formulations is performed showing how they differ in terms of damage estimation, crack propagation, and spatial discretization restrictions (when the energy-based Hillerborg criterion is applied). Another critical review of damage-based permeability laws is achieved to pinpoint their limitations.
Indeed, the use of a bounded and strictly increasing damage variable does not allow (1) an accurate description of permeability in the presence of macrocracks that are physically unbounded and (2) a proper description of the permeability decrease as those macrocracks are closed. Accordingly, a new strain-based permeability law is suggested. Its analytical validation is based on the BIPEDE tensile test for one loading-unloading cycle. This new formulation intends to enhance the predictiveness of damage-permeability coupling models leading to a better assessment of large structure's durability and serviceability.

| INTRODUCTION
The tightness of strategic and large structures such as dams and nuclear containment buildings is due to the low permeability of concrete and, mostly, to the continuous control and monitoring of their cracking state throughout their operational lifespan. With that regard, the ability to predict the evolution of damage and cracking states in time, in addition to the resultant permeability, are key issues for the durability and serviceability assessment process.
As it has been demonstrated by numerous experimental works, the permeability of concrete is strongly related to its damage state: 1. Darcy's mode: When damage is low and diffuse, the increase of concrete's permeability is due to the deployment of microcracks within its volume. This leads to a generalized increase of the porous network's connectivity. In the absence of strain localization, concrete is seen as a damage-dependent porous media that justifies the hypothesis of a transfer mode governed by Darcy's law. [1][2][3][4] 1. Fracture-based crack models: These models are based on elastic (or plastic) fracture mechanics to explicitly describe discontinuities within the concrete volume. Amongst this class, one can enumerate cohesive crack models, 12 lattice and network models, [13][14][15] and LSM-based* inclusive models using either nodal enrichment XFEM † or elemental enrichment EFEM ‡ . 16 For such approaches, both, the location and opening of the cracks are known and quantifiable, which facilitate the use of Poiseuille's law. However, for the precracked state, a continuous analysis of the domain is required, which raises the question of the transition from continuous (noncracked) to discrete (cracked) modelling. In Grassl and Bolander 15 and Larsson et al, 17 it is suggested to use a continuous strategy until crack initiation before pursuing calculations in the cracked area using a discrete approach. The threshold to switch control from continuous to cohesive modelling can be set arbitrarily based on strain, stress, or energetic quantities 15,18,19 or defined according to an error analysis between strong discontinuity (SD) gradient-based approach (as a valid approximation of cohesive models) and weak discontinuity (WD) integral-based approach (continuous). 20 It is worth mentioning, however, that the SD approach is a particular case of the WD methods using a sharper and singular weighting function (Green's function). 21 This shows that the WD method and the continuous damage approach can be used to describe equivalently the discontinuous nature of cracks. 2. Continuous damage-based models: These models are based on damage mechanics, which regards cracking of concrete as a continuous loss of rigidity. 22 In the particular case of strain-based damage criterion, the damage variable (usually denoted d) evolves from 0 to 1 and is computed based on the gap between a certain strain damage threshold and the developed equivalent strain. The way this so-called equivalent strain is computed distinguishes local models from gradient-enhanced ones and integral nonlocal-based ones. 20,21 For the first class (local models), the equivalent strain is computed for each integration point and is mesh sensitive if not regularized. For the second type of models (implicit gradient and nonlocal approaches), mesh sensitivity is alleviated as an internal length is introduced to account for a constant and material-dependent microstructural scale (usually related to the aggregates' size 23 ) over which energy dissipation occurs. Nonlocal strategies are generally more time consuming compared with local ones and require more refined mesh to preserve their nonlocal nature. This might lead to considerable computational memory and time cost when dealing with large concrete structural volumes. One can also notice the constant bridging with fracture mechanics even though the continuum media are used: Local damage models are usually regularized using the Hillerborg energy-based approach, 24 and the equivalent crack openings 20,25 and paths 26 are postprocessed using adapted algorithms to overcome the incompatibility between the crack's discrete nature and the used continuous approach. On that issue, the dissipated energy in a continuum domain is supposed to fully contribute to crack propagation and is, consequently, converted into an equivalent crack length. 27 This hypothesis remains strong since, experimentally, only a part of the dissipated fracture energy contributes to crack propagation ahead of the crack tip; this amount being intrinsic to the fracture process zone (FPZ) development is estimated at 40% of the total fracture energy. [28][29][30] The remaining 60% is associated with frictional energy dissipation. Eventually, despite particular shortcomings related to spurious boundary effects for nonlocal models and mesh dependency for local ones, the ability of continuous elasticity-based damage models to describe initiation and propagation of mode I cracks is not altered when proper regularization and identification techniques are used. Moreover, continuous finite element (FE) models remain, within the currently available computational capacities, the most efficient support for large concrete structures' (~10 m) modelling.
In this contribution, only the permeability associated with mode I cracking is of interest. Accordingly, mode II and mode III cracking modes should be addressed differently as the physics of their associated permeability differ. To serve our purpose, continuous damage-based models are retained and discussed according to their validity domains, their evaluation of damage around the FPZ, their evaluation of the crack opening in the localization point, and their effect on the resulting permeability. So this paper is twofold: 1. The first part of this paper recalls the mathematical framework of the used local, SD, and WD nonlocal modelling strategies (using a strain-based damage criterion) prior to permeability estimation. It also demonstrates the limits of existing damage-based permeability laws that lack a state variable and mathematical formwork able to describe accurately the permeability's evolution from a low damage state to a macrocracked one. 2. The second part introduces a new strain-based permeability law, addressing those drawbacks. It is more physically representative of transfer properties during the transition phase from lowly to highly damaged state. It also describes better the experimental results from the BIPEDE tensile test (one loading-unloading cycle) compared with existing damage-based permeability laws.

| DAMAGE-BASED DESCRIPTION OF DAMAGE-PERMEABILITY COUPLING
In the framework of a staggered approach, the mechanical calculation (ie, damage and crack opening estimation) represents the first step of the damage-permeability coupling. In this section, the damage formulation principles are recalled along with the postprocessing steps required to define an equivalent crack opening in the continuous domain.
In terms of the dissipated fracture energy, energetic criteria are also addressed to ensure mesh independency of local models and softening laws shape independency of both local and nonlocal damage formulations.

| Isotropic damage modelling
In the case of isotropic damage description, the general behaviour law relating the stress σ to the strain ε variable writes where E is the nondamaged stiffness and d is the isotropic damage variable increasing from 0 to 1 depending on the following loading function: where ε eq is an equivalent strain (see, for instance, the μ-Mazars damage criterion presented in Mazars et al 31 ) that can only increase with respect to the second law of thermodynamics starting from an initial value of ε d0 , ε d0 being the strain threshold starting from which the loss of stiffness is obtained. || f (t)|| t = max ( f (τ ≤ t)) is a norm that ensures the monotonous increase over time of a time-dependent function f .
). D c is, accordingly and by definition, representative of the area over which energy dissipation occurs for a given crack propagation. And ϕ defines the spatial distribution of such dissipation around the point of interest (ie, the point of localization- Figure 2).
where ε loc eq estimation involves the first and second strain tensor's invariants. 31 Depending on the used weighting function ϕ (Table 1), different models are obtained. By setting ϕ equal to the Dirac function δ, local damage models are obtained, the Green's function G leads to the second-order gradient-enhanced models, which are descriptive of an SD gradient-based approach, 20,21 and the classical Gaussian function ρ leads to the WD formulation 20,21 and stress-based (SWD) integral models. 33 It is also possible to derive a stress-based secondorder gradient-enhanced models (SSD) when the Green's function is used instead if the Gaussian one. In these 2 last cases, the effective internal length l ′ c depends on, both, the internal length l c and the stress distribution around the point of interest. It is worth underlining that all damage models, herein, depend on a strain-based damage criterion (ε loc eq in Equation 6); the "stress-based" terms refer rather to the energetic regularization (preventing mesh dependency-ε eq in Equation 6). Compared with SD and WD formulations, stress-based approaches SSD and SWD are more representative of the stress concentration ahead of the crack tip and unloading behind it. Stress-based formulations lead to a more physical regularization close to stress-free edges (geometrical boundaries or fully developed cracks). Also, they can Dufour et al 20 Peerlings et al 21 Dufour et al 20 Peerlings et al 21 φ and θ are angles in the principal stress basis (σ I , σ II , σ III ) α = 2 Giry et al 33 Stress-based second-order gradient enhanced (SSD) l c σ f t φ and θ are angles in the principal basis (σ I , σ II , σ III ) -prevent damage attraction by boundaries due to the truncation of the interaction domain of integration points without requiring a preconditioning of the denominator in Equation 6. Consequently, as the stress state is reduced due to damage, the effective internal length tends to zero and the weighting function tends to become singular and equal to a Dirac impulsion for the considered position. In other words, and for a cracked zone, the evolution of SWD and SSD damage models is equivalent to the one using local models. Generally speaking, these weighting functions (Table 1) define the shape of the FPZ over an influence domain sized Influence domain equivalent to a segment sized D c : Influence domain equivalent to a circle with a radius D c : Influence domain equivalent to a sphere with a radius D c : Naturally, this influence domain should be the same as the designated FPZ size (D c = l FPZ ) over which the shape is defined. This means that, for an objective description of energy dissipation along the FPZ (Figure 2A), the influence domain D c = l FPZ should be independent from the used weighting function to prevent a biased computation of the damage variable d and the equivalent strain ε eq x ð Þ ( Figure 2B). However, the internal length l c should be a seen as a fitting parameter and needs to be defined by solving Equation 7. l c can be considered as a subjective internal length of the nonlocal models whereas D c = l FPZ is a physical parameter that is objective regarding the nonlocal formulations. For the weighting functions in Table 1, the objective formulations in Table 2 are obtained.
In the case of local formulations and in the absence of an internal length, the influence domain D c is brought to the unit of the spatial discretization h (also called the mesh characteristic length). This requires an additional regularization step to prevent mesh dependency. Yet it remains problematic in terms of the maximal cracks' number estimation, which is more related to the FPZ size than the spatial discretization h.

| Crack opening estimation
As damage is described according to a continuous approach, a postprocessing step is required to compute the equivalent crack opening in the considered continuous domain. Herein, the 1D strategy defined in Mattalah et al 25 and accounting for elastic unloading is retained (Equation 8). The equivalent cracking strain ε ck then writes ; 2D ffiffiffi π 6 3 r ; 3D φ and θ are angles in the principal stress basis (σ I , σ II , σ III ) SSD where ε I is the elastic strain in the first principal direction (supposed to be the mode I cracking direction) and <x> + = max (x; 0). And for all weighting functions in Table 2 (in addition to local models), the equivalent crack opening w ck is computed according to Equation 9. 20 Equation 9 remains also true for local modelling strategies considering an influence domain equal to the spatial discretization unit D c = h. One can notice that, for nonlocal formulations, by using Equation 9 and Table 2, such definition does not depend on the used weighting function nor the subjective internal length l c (D c = l FPZ = c ste ). Also, it is worth underlining that the crack opening estimation is based on the integrated strain over the influence domain D c . This leads to the implicit hypothesis of one crack per influence domain. Physically, such hypothesis is only convenient if the influence domain is related to the size of the FPZ. Consequently, for nonlocal models, it suffices to choose D c = l FPZ to be physically representative. However, for local models, the influence domain is equal to the mesh characteristic size h (user-defined), which can be bigger or smaller than l FPZ . From a conceptual point of view and following a 1D tensile load, this can be problematic in the presence of nonnull residual stresses leading to multicracking: For h ≪ l FPZ , the number of cracks per l FPZ might be overestimated as more than one FE within the FPZ length can be damaged. On the contrary, for h ≫ l FPZ , and following an inverse reasoning, the number of cracks can be underestimated. Eventually, this means that local models can only be predictive regarding the number of cracks if the mesh discretization is equal to the FPZ size. Otherwise, an additional hypothesis has to be considered requiring a prior knowledge of the cracks' number n fiss per FE-in which case each crack opening writes w ck ¼ ε ck n fiss l FPZ -.

| Hillerborg energy-based criterion
Another issue related to the definition of the crack opening is its dependence on the used softening law (and damage evolution law) with respect to the considered fracture energy. In the case of nonlocal formulations, the weighting functions insure a mesh independency of the dissipated energy. However, for local models, mesh independency is not guaranteed by default and requires the application of regularization techniques such as the Hillerborg energy-based criterion. 24 It allows the implementation of a material descriptive behaviour law that energy is constant (σ, ε) regardless of the used spatial discretization (case of local approaches) and regardless of the postpeak softening laws (case of local and nonlocal models). Practically, the criterion (Equation 10) leads to the identification of the postpeak parameters for a given influence domain D c and a given fracture energy G F . In the case of local models, the Hillerborg criterion has a true regularization meaning since D c = h is mesh dependent. Whereas for nonlocal formulations, Equation 10 becomes mesh independent as D c = l FPZ regardless of the used spatial discretization and should be looked at as an identification formula of the postpeak parameters. This writes for all damage models in Table 2 and softening laws (Equations 3-5): For non-local models: where G F , f t , and E are known parameters identified at the specimen scale. One can notice that the residual stress when damage is equal to 1 needs to be null for the dissipated energy to be finite and defined. The application of Equation 10 leads to different validity domains of the postpeak softening laws depending on the damage formulation (local or nonlocal). For the sake of conciseness, the theoretical developments and their analyses are given in Appendix A. Eventually, one can refer to Tables 3 and 4 (summarizing the validity domain of the softening laws described above [ ) to select the convenient softening law and damage modelling strategy in Table 2. For other softening laws, the same analysis should be performed. An insight on the effects of regularization and damage formulation on the local and global mechanical responses is highlighted in Appendix B based on an FE analysis of a 1D beam tensile test. It shows how the energy-based regularization restrains spatial discretization and improves damage modelling objectivity. However, one should keep in mind that energy-based regularization does not ensure local objectivity in terms of damage evaluation within the FPZ (better objectiveness is achieved for the crack opening evaluation at the point of localization) nor the global response objectivity in the presence of structural effects such as snap backs. This issue also affects the permeability evaluation if based on the damage variable as it will be detailed in the next section.

| Limitations of damage-based descriptions of permeability
As mentioned in the introductory part, 2 transfer modes have been extensively studied: a first one through a low and diffuse damaged medium and a second through a localized crack.
Requires numerical solving Requires numerical solving  3,4 given that the performed analysis covered the whole specimen instead of the central part where damage increases. In that regard, results from Picandet et al 5 are more reliable since they cover and characterize only the damaged zone. Assuming the equivalence between compressive and tensile damage at low damage values, tests in previous studies [3][4][5] led to the definition of a damage-based exponential permeability law for 0 ≤ d ≤ 0.15. As d increases, the exponential law shows a considerable overestimation of the apparent permeability. Consequently, the Taylor expansion of the exponential law at d = 0 is preferred (up to the third order) herein as suggested in Choinska et al. 9 Another issue raised by experimental results in Choinska et al 3,4 is the existence of irreversible k i D and reversible k r D parts of damage-induced permeability at low levels (Equation 11). Such irreversible tendencies have also been observed in Bary 35 Choinska et al 3,4 ; based on a global analysis of damage in a compressive test (loaded specimen)   and Gérard 36 through a BIPEDE tensile test (concrete disc subjected to uniaxial tension- Figure 8) where advanced tensile loading has been applied until the development of macrocracks. This justifies the higher values of damage (up to 80%) that are considered in the models suggested in Gawin et al 37 and Dal Pont et al 38 based on the BIPEDE test. 36 Such models are debatable since at this high damage level, strain localization leads to the development of macrocracks that are governed rather by the Poiseuille's mode than Darcy's. As for the observed residual permeability and in the presence of macrocracks, it can be due to a given residual crack opening or, more likely, to a residual diffuse damage in the cracked domain. Indeed, for a complete unloading, the macrocrack can be considered fully closed whereas the network of microcracks around remains active. Such configuration favours the use of Darcy's law rather than Poiseuille's, which is valid for crack openings that are at least higher than 20 μm. 6 Accordingly, the irreversibility of permeability is herein associated with concrete damage and not with a residual macrocrack opening. Eventually, for whatever damage-permeability coupling model, one can notice that the variation of their parameter has more effect on the computed permeability compared with the effect of the damage softening law variation (Figure 3).
For the purpose of further developments, the reference permeability law retained herein is the following: with values issued from experimental results in Picandet et al 5 (α = 11.3, β = 1.6) and based on experimental observations in Choinska et al 3,4 It is important to underline that the previous permeability law has been derived in a local damage modelling sense despite the fact that only the mean response (mean damage of the influence domain) is accessible experimentally. Consequently, and from a numerical point of view, the computed permeability over the influence domain (Equation 12) is strongly dependent on the damage formulation (LOC/SW/SD/SSW/SSD in Figure B4). This has a nonnegligible effect on k D estimation ( Figure 4): For nonlocal formulations and in the case of an unloaded concrete at d = 15%, k D is 10 times lower than the reference value (LOC descriptive of the FPZ mean response). Gaps are 15 times higher for d = 100%. This simply means that the fitting coefficients in Table 5 for Equation 11 are only valid for local models and need to be redefined in the case of nonlocal modelling of damage. Another alternative, more physically representative of the FPZ's mean response, is to associate Darcy's transfer mode with a mean quantity that is less damage dependent (at least from a local point of view), for instance, the mean damage value that would require an additional Another drawback of Equation 11, also related to its damage dependence, is the absence of a path between loaded specimen and unloaded one that hinders the description of permeability evolution during a loading-unloading cycle.

| Permeability through a macrocracked medium
Poiseuille's law is usually used to describe permeability through a macrocrack 5,6,8,11,39 (Equation 13). It consists of writing the apparent permeability k P (P in k P stands for Poiseuille) as a function of the squared crack opening value § (which is unbounded in opposition with the damage variable): where ξ ≤ 1 is a factor to account for the roughness and shape effects of the through crack domain. 5,6 When ξ = 1, the domain is supposed to be perfectly smooth. The expression suggested in Rastiello et al 6 based on real-time water permeability evolution through a localized crack during a splitting test (Equation 14) seems more appropriate compared to others where ξ is considered constant (between 0.01 and 0.1 in Picandet et al 5 ). Indeed, as the crack opening becomes larger, the roughness effects are supposed to be relatively less and ξ should tend towards 1 ( Figure 5). It is worth mentioning that experimental results cover only crack openings up to 160 μm whereas the condition of ξ ≤ 1 is verified until w ck = 268 μm. So an extrapolation is performed for crack opening values between 160 and 268 μm.
Numerically, a macrocracked FE domain sized D c contains a macrocrack of an opening w ck and also uncracked domain around sized D c − w ck . The permeability k F (F stands for fracture in k F ) when considering the cracked domain only writes 40   (Table 6). A comparison between the different matching laws is depicted in Figure 6A-C: 1. In literature, 6,9,10,42 one can find a linear matching between the 2 modes where the increase of permeability, when damage exceeds 15%, is only due to the increase of the crack opening value. However, for values of damage between 15% and 90%, the effective crack opening remains quasinull, which leads to a constant permeability within that domain. 2. Other linear models suggest a damage-weighted addition of the 2 modes. 3,9 This has the advantage of not having to limit the damage value to 15% since at high damage values, only the Poiseuille's mode is activated. Still, this solution seems to underestimate permeability at low damage values when the crack's roughness effect is considered (timorous lines in Figure 6A-C). On the other hand, when roughness is not considered (bold lines in Figure 6A-C), the damage-pondered linear law overestimates Darcy's permeability at high damage values leading to a nonphysical decrease of permeability as strain localizes and Poiseuille's mode is activated.

Lastly, the damage-pondered logarithmic law 3,41 (Equation 16) offers a smoother transition from Darcy's to
Poiseuille's mode. However, when the crack's roughness effects are accounted for, it shows a nonphysical decrease of permeability as damage is initiated and also underestimates the resulting permeability for 0 ≤ d ≤ 0.15. In the

Rastiello et al 6 and Asali et al 41
Ezzedine El Dandachy 42 case of the BIL softening law, the decrease also occurs at d~80% as Poiseuille's mode has more weight even if its contribution remains less than Darcy's.
Herein, the log-type matching is retained for further developments (Equation 16). Two main issues need to be addressed: The first one is the weighting of Darcy's and Poiseuille's modes so as to prevent permeability decrease with damage. And the second is the damage dependence of the matching law preventing an accurate description of permeability as concrete is unloaded. Indeed, for a fully damaged element (d = 1), the only active mode is Poiseuille's. During an unloading cycle, the damage variable due to tensile loads remains equal to one (one is reminded that the damage variable can only increase) whereas the crack opening and Poiseuille's flow decrease. Such configuration leads to an unrealistic residual permeability that is either overestimated (k = k D (d = 1)) in the case of linear matching 3 or null in the case of damage-pondered laws (k = k F (d = 1) = 0 in the absence of any residual crack opening). 42

| NEW STRAIN-BASED DESCRIPTION OF DAMAGE-PERMEABILITY COUPLING
The previous section has presented different damage-based matching laws and the related drawback that can be identified. Several enhancements can be suggested aiming at a more physical description of permeability evolution from a low diffuse damage to a highly localized damage state. Particularly, the foreseen law should • limit damage model dependence on the permeability law (especially for Darcy's transfer mode), • ensure a monotonously increasing permeability with damage, • distinguish reversible and irreversible parts of permeability, and • account for loading-unloading cycles.

| Matching law
1. Strain-based permeability for Darcy's transfer mode: Let us first express the permeability as the sum of a reversible and an irreversible part in line with experimental results in literature. 3,4,35,36 The reversible permeability due to damage k r D would be associated with the transfer properties of the reversible closure and opening of microcracks or the elastic behaviour of the damaged porous network subjected to loading. The irreversible part k i D would be descriptive of the fluid's flow through the residual and irreversible openings of the mechanically microcracked domain. In opposition to the reversible part, the irreversible one can only increase in time depending on the loading level. Its maximum value is obtained when the macrocracked state is reached. At this stage, it is supposed that, in the vicinity of the crack lips, concrete is unloaded and the maximum damage level has been achieved. In other words, the microcracks network around the macrocrack is within its residual mechanical state. From that point on, any evolution of the transfer properties is directly associated with the variation of the macrocrack opening. As a consequence, and within the suggested framework, the residual opening of macrocracks is supposed null; irreversibility is rather associated with the microcracked network. In a hydraulic sense, it is possible to define an equivalent residual macrocrack opening leading to the same permeability evaluation. However, from a physical point of view, a residual opening of the macrocrack would be more representative of the crack lips' mismatch during the closure or the surface alteration (during a compressive cycle) and should be addressed separately. In Equation 11, the total permeability k D and irreversible one k i D have already been defined. So what lacks is the definition of the reversible part k r D ¼ k D − k i D , a loading path between the total and irreversible permeability values and a state variable to prevent damage dependence of k D . This last issue is particularly of interest to limit the dependence of Darcy's transfer mode modelling on the damage formulation (local or nonlocal- Figure B4).
with 0 ≤ d r ≤ 1 a "reversible" and strain-based variable to prevent damage dependence of Darcy's law and d lim a fitting parameter depending on the measured residual permeability. γ is a loading path parameter varying between 0 and 1 chosen herein arbitrarily as 0 ≤ γ ¼ ε eq ε eq ≤ 1: when ε eq ¼ ε eq , γ = 1, and d lim ¼ jd r j t -during the unloading phase: Accordingly, the reversible part of Darcy's contribution writes 2. Strain-based permeability for Poiseuille's transfer mode: Poiseuille's mode in Equation 15 is already expressed in terms of the crack opening rather than damage. Such formulation does not require any additional modification. The challenge regarding Equation 15 remains related to the experimental quantification of the crack's roughness and tortuosity effects on the measured Poiseuille's flow. 3. Log-type matching: As demonstrated in Section 2.2.3, the use of damage as a weighting function does not allow a proper description of permeability at low damage values, as it does not ensure an increasing permeability with damage and, being monotonously increasing, does not allow an accurate description of the permeability of an unloaded concrete. Accordingly, the foreseen weighting function needs to be reversible and only active within the intermediate zone 0.15 ≤ d < 0.99. Such pondering can be achieved using the previously defined "reversible" damage variable (using ε eq instead of ε eq in Equations 3-5), and its associated weighting function has to be decreasing only when 0.15 ≤ d. This leads to a new strain-based (and damage independent) permeability law (Equation 20): with δ a fitting parameter defining the rate of pondering between the 2 modes. In the absence of further experimental results, δ should be adjusted to get an increasing permeability in the loading direction. Finally, and as shown in Figure 7, the new matching law (Equation 20) is more accurate and physically acceptable for both the loading and unloading phases compared with the original logarithmic one ( Figure 6): The permeability transition from Darcy's to Poiseuille's mode is achieved successfully with respect to the monotonously increasing permeability with damage-the permeability transition from the loaded to the residual (unloaded) configuration is properly defined so as to insure a nonnull residual permeability.
This new formulation remains valid whether unilateral effects and residual strains are present in the damage model or not. 31 In a local sense (for the FPZ characterization), one should keep in mind that the fitting of Equation 20 depends on the used damage softening law. Particularly, the condition of increasing permeability leads to a smaller value of δ (δ = 0.5) in the case of BIL law compared with LIN/EXP softening laws (δ = 0.8). This is due to the slow evolution of damage after initiation when using the BIL law to meet with energetic criteria. However, still in a local sense, its dependence on the damage formulation remains limited given the strain-dependent nature of Equation 20. Indeed, from a theoretical point of view, the estimation of crack opening values is sensibly the same using local or nonlocal models when the energetic criteria are respected (see Appendix B). The effect of damage formulation is expected to be more pronounced in the presence of structural effects (snap backs). In that case, different fittings would be required depending on the used damage formulation.

| Analytical validation based on the BIPEDE tensile test
Experimental results covering the evolution of concrete's permeability under loading-unloading cycles for low and high damage values are not numerous in the literature. In this contribution, the BIPEDE tensile tests developed in Gérard 36 is considered. It consists of applying uniaxial tensile loads to a disc glued to steel plates ( Figure 8). Tension is applied to the plates and is transmitted to the disc until cracks develop.
The disc is 40 mm thick and has a diameter of 110 mm. The plates are holed so as to allow fluid (water) transfer through the concrete disc (hole's diameter is 55 mm). Test results are provided in terms of water flow through the concrete specimen q (m/s). 36 The specimen's equivalent permeability can be estimated using Equation 21 supposing that the stationary state is constantly achieved.
with ΔP ¼ 0:25 MPa; e ¼ 40 mm; μ ¼ μ water ¼ 10 − 3 Pa:s; and S ¼ π The purpose of this paragraph is not to show the validity of a given mechanical model but to demonstrate that the suggested matching law in Equation 20 can be fitted to realistic experimental results. It is worth mentioning that, experimentally, one or 2 cracks were obtained (not forcibly centred). Herein, one crack is supposed to develop and local strains at the point of localization are directly used to define-according to a given softening law-the resultant damage and crack opening values (explicit analytical scheme). Also, in the framework of this contribution, only homogeneous FE modelling strategies of concrete are mentioned. However, nothing prevents the use of the suggested strain-based matching law (Equation 20) in the case of mesoscopic modelling strategies where concrete is considered as a biphasic material (aggregates and cement paste). 43 For the numerical application, the considered concrete properties are the following: f t = 4.5 MPa, G f = 80 N/m, and E = 30 GPa. The influence domain D c is estimated at 2 d max = 6 cm. In Figure 9, the adjustment of Equation 20 to experimental results is performed using the least squared method where only roughness effects (α w ) and irreversibility parameter d lim are adjusted. In terms of crack opening evolution, analytical developments lead to the same results for LOC/ SD/WD/SSD/SWD approaches when the energetic Hillerborg criterion is properly applied (Appendix B). By comparing the obtained strain-based curve (Equation 20) with damage-based one (Equation 16) in Figure 9, one can notice the accuracy of the newly suggested formula to describe permeability evolution during both the loading and unloading phases using LIN/BIL and EXP softening laws. For this test, the roughness effect on the computed Poiseuille's permeability is less pronounced compared with the one in Rastiello et al. 6 For that reason, the decrease of permeability with damage using Equation 16 is limited at low damage values. As shown in Figure 6, this drawback of damage-based approaches is more pronounced as the roughness and tortuosity effects increase. As for the unloading phase, it is clear that Equation 20 leads to a more physically admissible results as Equation 16 provides a null residual permeability. The strong scattering of experimental data as concrete is fully unloaded does not allow a clear understanding of whether the loading and unloading permeability paths are superimposed or not. With that regard, an in-depth experimental investigation of permeability evolution of an unloaded concrete is still lacking. For this contribution, the path is considered reversible and no hysteresis is included. This can be improved, if clearly demonstrated by experimental results, by modifying the shape of Equation 19 for the unloading phase accordingly. Nevertheless, our main interest is served since the suggested strain-based permeability law in Equation 20 allows appropriate fitting of realistic experimental results regardless of the used damage modelling strategy (in the absence of structural effects).

| CONCLUSIONS
In this contribution, a new strain-based permeability log-type matching law is suggested instead of damage-based ones. Being exclusively strain dependent, the new law does not depend on unilateral effects or damage formulation (as far as the local behaviour of the FPZ is concerned). As demonstrated through the BIPEDE tensile test, it allows better description of permeability decrease during the unloading phase and ensures a monotonously increasing permeability as the loading advances. As convincing experimental works of permeability evolution during the unloading phase are lacking, no hysteresis is considered herein for the reversible permeability part. Our future work shall investigate concrete hydric behaviour more thoroughly for more than one loading-unloading cycle to better characterize the permeability during unloading phases.
• Linear softening law: In Equation 3 (LIN), the unknowns are the damage threshold ε d0 and the failure strain ε u ≥ ε d0 ( Figure 1A). After simple calculations (Equation 10), one can derive : the computed behaviour is perfectly brittle: ε u = ε d0 and the fracture energy per an element sized D c is equal to the elastic one.
• Bi-linear softening law: In Equation 4 (BIL), the unknowns are the damage threshold ε d0 , the failure strains ε ′ u and ε u , and the n parameter ( Figure 1B). The number of parameters exceeds the number of available equations. For that reason, additional hypotheses need to considered 29,30 : The first slope is defined using the intrinsic fracture energy that represents 40% of the total fracture energy. Eventually, the second slope is deduced to verify the global Hillerborg criterion in Equation 10.
Compared to the linear softening law, values of D c are more restrained and cannot exceed 4 5 G F E f 2 t given the use of the intrinsic energy ( 2 5 G F ) rather than the total one G F .
• Exponential softening law: In Equation 5 (EXP), the unknowns are the damage threshold ε d0 and the shape parameter B t ( Figure 1C). Calculations are more complex than the 2 previous ones but eventually lead to the following system of equations: To prevent the damage threshold f t E from being inferior to f t E , another regularization strategy is suggested in Mazars et al. 31 It consists of considering the linearized postpeak behaviour law (EXPL-tangent line at ε ¼ instead of the exponential one ( Figure A1).

This simplifies calculations (Equation A6
) since by definition, the yield strength becomes equal to the ultimate tensile strength of concrete. However, it does not reduce the number of unknowns. Indeed, even if the damage threshold is known, the amount of the fracture energy needed for regularization remains to be identified. This amount, denoted αG F , should allow the the area under the behaviour law to approaches the best G F /D c .
where G EXPL F stands for the dissipated energy using the linearized law regularization.
A simple functional analysis of the condition G EXPL leads to the conclusion that G EXPL F →G F only when α → + ∞. Hence, the computed fracture energy will always be underestimated compared to the desired one leading to the computation of a more brittle behaviour. Eventually, the choice of α = 1 is the only one acceptable with The regularization of the exponential law is controversial for small mesh sizes (D c ≤ 2 5 G F E f  solution to be acceptable, the damage threshold should remain strictly positive, hence the existence of a minimal value D c, min . In opposition to the LIN and BIL softening laws, the computed behaviour for EXP law shows a hardening phase for small D c values: To dissipate the same fracture energy and not overestimate the tensile strength, the damage threshold needs to be reduced, which induces a hardening behaviour until the tensile strength is reached ( Figure B2A,B). The main drawback, though, is the underestimation of the damage threshold (at least from a numerical point of view) Table A1) and overestimation of the failure strain at low damage values. This result is controversial, though, from a physical point of view since the damage strain threshold is physically constant and cannot be equal to zero nor depend on the mesh size (for local models). This restrains even more the use of the EXP law when both energetic and damage strain threshold are considered.

APPENDIX B EFFECT OF LOC/SD/WD/SSD/SWD FORMULATIONS ON THE GLOBAL AND LOCAL BEHAVIOUR
To illustrate the utility of the regularization steps presented above, a 1D beam under pure tensile loads is considered. It has a length of L = 1 m, is fixed on one side, and undergoes a tensile load on the other ( Figure B1). To prevent a diffuse damage state and be able to characterize strain and damage distribution around a localized crack, a centred default is considered. The selected weak point has a tensile strength that is reduced by 1% compared to the rest of the beam. The following mechanical properties are considered: f t = 4.5 MPa, G f = 80 N/m, and E = 30 GPa. The FPZ width is supposed to be D c = l FPZ ≈ 3 * d max = 6 cm (d max = 2 cm).
As for the spatial discretization, values of h are issued from Table 4 defining the validity domain of each postpeak law and local/nonlocal damage modelling: The identified parameters for LIN, BIL, and EXP laws are summarized in Table A1.
The implemented behaviour laws are depicted in Figure B2. One can notice, in the case of local models, how the computed brittleness increases with the size of the considered FE. Figure B2A,B shows the hardening behaviour of the EXP softening law for small mesh sizes. In the case of stress-based nonlocal models, one can notice how the area under the (σ, ε) curve is higher than the one of strain-based nonlocal models given the decrease of their internal length with stress concentration at the crack tip (Table A1 and Figure B2D,F). Nevertheless, the dissipated energy around the FPZ remains the same.

B.1 | Local response
By local response, it is meant the damage distribution and energy dissipation over the influence domain D c (FPZ), and crack opening evolution at the weak point (WP) in Figure B1. Their estimation can be performed using explicit and analytical developments. In Figure B3, the maximum damage value at the weak point d WP is depicted in the case of local and nonlocal formulations, whereas in Figure B4, damage distribution over the influence domain d FPZ is shown for nonlocal configurations (the damage variable is homogeneous over the mesh size h in the case of local models). In Figure B5, the crack strain evolutions at the weak point are provided. At this position, one can see the dependence of the damage and crack opening evolutions on both the used softening law ( Figure 1) and damage modelling strategy (local Figure B5A-E and nonlocal Figure B5D,F/ Figure B4). In terms of the maximum damage value (at the weakest point- Figure B3), more gaps between local and nonlocal approaches are observed when the mesh size h recedes from the FPZ size l FPZ . This result underlines again the previous remark mentioned in Section 2.1.3 for crack opening estimation that local damage models are only comparable with nonlocal ones when the considered spatial discretization is equal to the FPZ size (according to the hypothesis of one crack per FPZ width). For D c = h = l FPZ = 6 cm, absolute differences between the damage evolutions at the weakest point ( Figures B3D,F and B4) due to the use of different softening laws and different weighting functions remain less than 0.2. As for the equivalent crack opening values, absolute differences due to the shape of the softening laws are less than 3 μm ( Figure B5D,F). Within these margins, all damage modelling strategies and softening laws can be considered equivalent in terms of stress-strain laws ( Figure B4A) and energy dissipation over the FPZ ( Figure B6). But even for the same dissipated fracture energy G FPZ F , it is worth mentioning that nonlocal modelling of damage around the FPZ is more realistic (particularly stress-based approach), as it decreases away from the localization point ( Figure B4A), compared with local models where damage is constant over the influence domain.

B.2 | Global response
As one can see in Figure B2C-F, the LIN and EXP behaviour laws are about the same within their respective validity domains in Table 4. One should expect similar response at the structural scale as well. As for the BIL law, more gaps are expected given the different energy dissipation paths. Herein, the study is limited to the EXP law and LOC, WD, and SWD damage formulations (similar tendencies are expected for LOC, SSD, and SSD). The obtained structural responses are shown in Figure B7A-C. They are expressed in terms of the normalized global stress (force over the beam's section) vs the displacement over the beam's length. Figure B7D-F shows the postprocessed crack opening as a function of the imposed displacement. The numerical problem is solved using Cast3m FE code 45 in which an arc-length FIGURE B1 1D beam under tensile test with a predefined position of default for damage localization technique is implemented to trace the equilibrium path when the global response exhibits eventual snap backs. From Figure B7A,D, one can observe that the stress-based nonlocal formulations SWD is more sensitive to the spatial discretization compared with strain-based ones WD; the more refined the mesh is, the more brittle the global structural behaviour gets. This result is expected given that the performed regularization (Table 3) is performed based on a continuous solving where the influence domain tends to zero as damage localizes. In FE analysis such configuration is never reached since, by definition, the influence domain cannot be less than the mesh characteristic length h. Consequently, the computed brittleness is underestimated as the mesh size gets bigger. In other words, the use of refined mesh in the case of stress-based models is recommended in line with the regularization in Table 3 to obtain realistic responses. Generally and for the same mesh size h, the SWD result remains bounded by LOC and WD responses ( Figure B7B,E) as previously underlined. As for LOC models ( Figure B7C,F), the hardening behaviour of the EXP law obtained for h< 2 5 G F E f 2 t induces a mesh-dependent global response (such results are not obtained for LIN and BIL laws, which remain brittle to quasibrittle). However, for h> 2 5 G F E f 2 t and ε d0 = f t /E, the Hillerborg regularization ensures the objectiveness of the global response.

B.3 | Discussion
As it has been explicitly shown in this section, the use of continuous damage-based models is able to describe concrete cracking whether local or nonlocal approaches are selected. However, the obtained results can largely differ in the absence of additional energetic and usage restraints. With that regard, one can define 3 criteria according to which a given modelling strategy can be selected or discarded: • Local response objectivity: It refers to the strain, stress, damage, and crack opening estimation within the considered domain of influence. At this scale, and in terms of local quantities' (particularly strain and damage) distribution away from the localization point, nonlocal models are more physically representative compared with local ones. However, the stress-based damage models are the only ones able to reproduce the decrease of the influence domain Table 2) due to the stress unloading in the vicinity of crack lips. 33 Such differences are expected to FIGURE B2 Behaviour law using LIN/BIL/EXP softening laws with respect to the Hillerborg energy criterion. For local models: ; and E ð Þ h ¼ ; and E ð Þ h ¼ affect the modelling of permeability as well especially if it is damage based. When the mean mechanical behaviour over the influence domain is considered (represented by the averaged crack opening or mean damage and strains), and regardless of the used postpeak law, local models show approximately the same response as nonlocal ones when the mesh size is equal to the FPZ size. In other words, local and nonlocal models are locally comparable when they both (1) have the same influence domains and (2) verify the same energetic criterion. This is true also in terms of the cracks' number as only one crack can develop within a given influence domain; physically, the influence domain of a crack has the same size as the FPZ over which dissipation occurs. • Semiglobal response objectivity: By semiglobal is meant the comparison of local and global quantities as usual experiments provide a partial description of concrete behaviour at the local scale (FPZ). For instance, in the case of splitting or flexion tests, the maximal crack opening (local quantity) is monitored along with the applied force (global . For nonlocal models with l FPZ = 6 cm : (D) SD/SW and (F) SSD/SSW quantity) whereas the stress state around the FPZ is (experimentally) uncharacterized. Since structural effects are considered, the objectiveness of the semiglobal response is conditioned by the global one and vice versa. Obviously, such aim can only be achieved for the mean mechanical response when local and nonlocal models have the same domain of influence (local objectiveness). • Global response objectivity: In the case of local models, the objectiveness of the global (or structural) response is ensured when the Hillerborg energy-based and damage threshold criteria are applied. Within the validity domain of each postpeak law, the global response scarcely varies with the mesh size. For nonlocal models, the stress-based ones seem particularly mesh dependent as their influence domain tends to localize within a defined mesh size even if, locally at the FPZ level, the dissipated energy is constant. Their mesh dependency seems to decrease as the spatial discretization gets smaller, which can become heavily time consuming for large structures. Generally speaking, at the structural scale, the local, nonlocal, and stress-based nonlocal models can hardly provide the same response even though the energy criterion is applied. Nevertheless, when local objectiveness is ensured at least for the mean behaviour of the FPZ, LOC and WD/SD remain distant from the SWD/SSD configurations that show naturally more brittleness at the structural scale.