Uy Cong

Nguyen

Adnan Ibrahimbegovic
email: adnan.ibrahimbegovic@utc.fr

Visco-plasticity stress-based solid dynamics formulation and time-stepping algorithms for stiff case

Keywords: Helligner-Reissner formulation, visco-plasticity, dynamic response, stiff solids, energy conserving/decaying schemes

In this paper, the visco-plastic behavior of solids is formulated based on the Hellinger-Reissner variational principle. The latter provides an optimal basis to obtain an increased accuracy of computed stress by using independent space-discretizations of both displacement and stress fields. The formulation is further adapted to solid dynamics and in particular to so-called stiff case that is characterized with a large difference in stiffness and/or evolution time between different deformation modes, which presents a great challenge when integrating the dynamic response. For this class of problem, we here propose the algorithmic modifications leading to energy conserving or decaying schemes that can either preserve the total energy or dissipate the energy of higher modes, and thus control the overall stability of numerical computations with large steps over long period of time.

Several numerical examples are presented to illustrate a very satisfying performance of the proposed approach and algorithm.

Introduction

This work seeks to contribute on formulating the time integration scheme for nonlinear dynamics of stiff systems undergoing visco-plasticity phenomena. The proposed solution concerns two research fields in computational mechanics that have been receiving a great attention over the years. The classical algorithm for visco-plasticity can be found in literature, see [Zienkiewicz andCormeau, 1974, Hughes andTaylor, 1978]. In these studies, the partial differential governing equations of visco-plasticity are discretized and solved numerically via the traditional finite element method. Only displacement is considered as independent field in the discrete approximation of the weak form of governing equations, which is then examined under static or dynamic loading conditions by using operator split method for plasticity [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF] that always lead to the first-order time integration schemes [START_REF] Chorin | Product formulas and numerical algorithms[END_REF]. There are several other researches relevant to visco-plasticity problem, especially in metal forming and cutting [START_REF] Pipard | Elastovisco-plastic modeling of mild steels for sheet forming applications over a large range of strain rates[END_REF], Mahnken et al., 2013].

Here, the conventional visco-plasticity problem is formulated by using the Helligner-Reissner approach and the finite element method. This approach enlarges the solution space including displacement and stress fields, which increases the stress accuracy, see [START_REF] Auricchio | The dimensional reduction modelling approach for 3d beams: Differential equations and finite-element solutions based on hellinger-reissner principle[END_REF], Auricchio et al., 2015]. Less computational cost can be obtained by using this method since it can yield an accurate solution even in a coarse mesh and especially in incompressible condition, see [START_REF] Viebahn | A simple and efficient hellinger-reissner type mixed finite element for nearly incompressible elasticity[END_REF]. This advantage can be observed in the numerical examples with both displacement and stress fields considered as independent fields. This visco-plasticity problem is further generalized to dynamics, in which the inertia force plays role in the governing equations of a system.

In order to deal with dynamics applications of interest for this work, proper time-stepping scheme should be proposed to compute structural responses over long period of time. Here, the stability of the time integration scheme must be maintained over a large number of time steps. Several time-stepping schemes for dynamic problems are successfully developed, which are proved to remain stable for long-term computation, see [START_REF] Bathe | On a composite implicit time integration procedure for nonlinear dynamics[END_REF], Hilber et al., 1977, Erlicher et al., 2002]. In the visco-plasticity problem, the dissipative property of energy does not guarantee the numerical stability of time-stepping schemes, which is very important quality for numerical computation. For that reason, the energy conserving scheme is introduced so that the total energy of system, consisting of the kinetic and potential energies, is always preserved in the elastic regime unless the visco-plastic dissipation starts to evolve due to visco-plastic strain in a system when the stresses satisfy the yield criteria.

The proposed scheme can support the computation for a long period of time for several particular physic problems, providing the estimate of plastic evolution in long-time. The underlying idea of this scheme is the combination of mid-point time integration scheme and algorithmic equations for the computation of displacement and stress within an interval of time step. In the view of computational accuracy, the mid-point time integration scheme is well-known for the higher accuracy compared to other schemes such as backward or forward Euler method [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF]. For higher accuracy of stress with coarse finite element mesh, we also introduce a modified energy conserving scheme that can lead to the energy decaying scheme, since it is desirable to dissipate the contribution of higher frequencies and corresponding modes to the computed nonlinear response of a stiff system. We follow the approach [Ibrahimbegovic andMamouri, 2002, Armero, 2006], the variation of velocities and stress computation within a time interval is modified so that it reduces the participation of high-frequency modes in controllable manner via dissipation coefficients α and β.

The paper outline is as follows. The theoretical formulation of visco-plasticity problem using Hellinger-Reissner functional is presented in the Section 2. It is noted that this study is restricted to structural behavior undergoing only small strain. Section 3 is devoted to numerical implementation of visco-plasticity problem including the discretization of the physical domain and time-dependent variables. In addition, the energy conserving and decaying schemes are formulated to control the overall energy of system and dissipate higher frequency modes. In Section 4, several selective simulations of the visco-plasticity problem in a thick-walled pipe under dynamic loading scenarios are simulated. In Section 5, our finding is finally summarized by the conclusions and discussions.

Theoretical formulation of visco-plasticity problem

In this section, the Hellinger-Reissner variational principle is employed for visco-plasticity at small strain. Let Ω be a bounded domain in R 3 with piece-wise smooth boundary ∂Ω, satisfying

∂Ω = ∂Ω 1 ∪ ∂Ω 2 and ∂Ω 1 ∩ ∂Ω 2 = ∅ (2.1)
where ∂Ω 1 denotes the displacement boundary and ∂Ω 2 the stress boundary. The linear total strain (u) tensor in small strain theory is defined by the symmetric gradient of displacement displacement u field as in Eq. (2.2).

(u) := ∇ s u = 1 2 ∇u + (∇u) T (2.2)
The symbol ∇u denotes the gradient of a displacement u vector, e.g. ∇u = [u i,j ]. The above total strain tensor is assumed as a combination of elastic strain e tensor and visco-plastic strain vp tensor, hence we have the relation e = -vp . The problem of visco-plasticity at small strains is characterized by the strain energy functional Π int (u) as follows:

Π int (u) = Ω ψ e ( -vp ) dV (2.3)
where the elastic free energy density function ψ e ( -vp ) is a function of elastic strain tensor:

ψ e ( -vp ) = 1 2 ( -vp ) : C : ( -vp ) (2.4)
with C as the fourth-order elasticity tensor. The symbol ':' stands for the double contraction operation, e.g. A : B = A ij B ij (two second order tensor) and A : B = A ijkl B kl (fourth and second order tensors). From such elastic free energy function, the relationship between stress σ tensor and elastic strain e tensor can be derived as

σ := ∂ψ e (•) ∂ = C : ( -vp )
In visco-plasticity problem, the yield function φ(σ) governs the evolution of visco-plastic strain vp along with time. The negative value of the yield function indicates the elastic regime without the development of visco-plastic strain ( ˙ vp = 0). Any non-negative value of the yield function implies the visco-plastic regime, with the corresponding change of visco-plastic strain ( ˙ vp = 0) according to:

φ(σ) := sgn F |F | m ; F = 3J 2 /σ 2 y -1; φ(σ) ≤ 0 ⇒ ˙ vp = 0 φ(σ) > 0 ⇒ ˙ vp = 0 (2.5)
where σ y is the yield stress and J 2 is the stress invariant. The dot in ˙ vp denotes the derivative of visco-plastic strain with respect to time, or the rate of visco-plastic strain. In this study, we choose the parameter value m = 1 in Eq. (2.5); we thus recover the classical viscoplasticity model (e.g. [START_REF] Ibrahimbegovic | Classical plasticity and viscoplasticity models reformulated: theoretical basis and numerical implementation[END_REF]), where the corresponding yield function accepts all admissible values in R as follows:

φ ≡ F = 3J 2 /σ 2 y -1; J 2 = 1 2 dev(σ) : dev(σ) (2.6)
where the operator dev(σ) stands for the deviatoric part of stress tensor σ. The external load potential Π ext (u) and kinetic energy K( u) are written as follows:

Π ext (u) = Ω f • udV + ∂Ω 2 t • udA ; K( u) = 1 2 Ω ρ u • udV (2.7)
with the body forces f in Ω and traction vector t on ∂Ω 2 , ρ as density and the overhead dot in u denotes the derivative with respect to time. The dot between two given vectors denotes the dot product, e.g. a • b = a i b j .

The stress tensor σ is considered as an additional independent variable by appealing to the Legendre transformation for the free energy function in Eq. (2.8). The complementary stored energy χ(σ) is a function of stress field, in which C -1 is the inverse of the fourth-order elasticity tensor.

ψ e ( -vp ) = σ : ( -vp ) -χ(σ); χ(σ) = 1 2 σ : C -1 : σ (2.8)
From Eq. (2.8) and Eq. ( 2.3), the Hellinger-Reissner functional E HR (σ, u) is now written as functional of the stress σ and displacement u fields.

E HR (σ, u) = Ω [σ : ( -vp ) -χ(σ)] dV -Π ext (u) + K( u) (2.9)
The least action principle can be used to obtain corresponding variational equations of the functional E HR (σ, u) above, which can be written as equivalent d'Alembert equations (e.g. [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF]) in terms of the virtual displacement δu, the virtual stress δσ as follows:

G u (u, σ; δu) := δ u E HR = Ω δ : σdV + Ω δu • ρüdV - Ω δu • f dV - ∂Ω 2 δu • tdA = 0 G σ (u, σ; δσ) := δ σ E HR = Ω (δσ : ( -vp ) -δσ : ∂ σ χ(σ)) dV = 0
(2.10) where the virtual strain δ tensor is defined as ∇ s δu in the same manner of Eq. (2.2) and the symbol ∂ σ χ(σ) denotes the derivative of complementary stored energy with respect to stress field. The corresponding Euler-Lagrange equations can be obtained directly from the variations in Eq. (2.10).

G u (u, σ; δu) = - Ω δu • (div σ + f -ρü)dV = 0 ⇒ div σ + f = ρü G σ (u, σ; δσ) = Ω δσ : ( -vp -∂ σ χ(σ)) dV = 0 ⇒ -vp = ∂ σ χ(σ) (2.11)
We can appeal to the second law of thermodynamics by means of the Clausius-Duhem inequality to guarantee the thermodynamical admissibility of the constitutive relations

D vp = σ : ˙ -ψe = σ : ˙ vp ≥ 0 (2.12)
where D vp denotes the visco-plastic dissipation. In principle, all the stress values are admissible, but those outside the elastic domain are penalized by an additional term P (•) directly proportional to the penalty factor 1/η (with η as the viscosity coefficient depending on material property). The modified visco-plastic dissipation D vp η is now written as follows

-D vp η [σ] = -σ : ˙ vp + 1 η P (φ(σ)) (2.13)
The simplest choice of the penalty term is a quadratic functional P (•) placing higher penalty on the stress states which are further outside the elastic domain.

P (φ(σ)) =    1 2 φ(σ) 2 ; φ(σ) ≥ 0 0; φ(σ) < 0 ⇒ d dφ P (φ(σ)) =< φ(σ) > (2.14)
where < • > is the Macauley parenthesis, which is defined as:

< φ >:= (φ(σ) + |φ(σ)|)/2 = φ(σ); φ(σ) ≥ 0 0; φ(σ) < 0 (2.15)
The evolution of visco-plasticity can be derived by means of the principle of maximum viscoplastic dissipation. From all possible states σ * satisfying the yield criterion φ (σ * ) ≥ 0, the one σ that maximizes visco-plastic dissipation D vp η is the solution. This problem can also be written as constrained minimization problem.

D vp η (σ) = max φ(σ * )≥0 D vp η (σ * ) ⇔ -D vp η (σ) = min φ(σ * )≥0 -D vp η (σ * ) (2.16)
The Kuhn-Tucker optimality conditions for the minimization problem defined in the above equation lead to the evolution equation of visco-plastic strain.

0 = ∂D vp η (σ) ∂σ := ∂[-σ : ˙ vp + 1 η P (φ(σ))] ∂σ = -˙ vp + 1 η < φ > ∂φ ∂σ (2.17)
Hence, the rate of visco-plastic strain is written as

˙ vp = 1 η < φ > ∂φ ∂σ =: β(σ); ∂φ ∂σ = 3 σ 2 y dev(σ) (2.18)
3 Numerical implementation of visco-plasticity problem

In order to provide the solution for the system of differential equations in Eqs. (2.11) and (2.18) above, the physical domain is discretized via the finite element approximation, while the nonlinear dynamic response is obtained by means of a proposed time-stepping scheme. The main computational task of such a procedure is to calculate the values of the nodal variables and interpolation parameters over a single time interval, which is discussed further in this section.

Finite element approach

The stress and displacement fields are interpolated independently. At time t n+1/2 , the displacement field u n+1/2 and acceleration field ün+1/2 are interpolated from nodal displacement d n+1/2 and nodal acceleration a n+1/2 vectors via shape function N matrix, respectively.

u n+1/2 = N d n+1/2 ; ün+1/2 = N a n+1/2
The corresponding strain approximation can be derived from displacement field as follows n+1/2 = Bd n+1/2 ; B = ∇ s N Meanwhile, the stress field in each element is interpolated directly using a complete linear polynomial S matrix, which is similar to the stress interpolation [Pian andSumihara, 1984, Wilson andIbrahimbegovic, 1990] in terms of stress interpolation parameters τ n+1/2 matrix.

σ n+1/2 = Sτ n+1/2
For convenience, the index n+1/2 for the obvious terms is only shown where it is needed. From Eq. (2.10), the discretized weak forms for a typical element Ω e are written in the matrix forms as follows.

G e u = δd T Ω e B T Sτ dV + Ω e ρN T NadV - Ω e N T f dV - ∂Ω e 2 N T tdA = 0 G e σ = δτ T Ω e S T (Bd -vp -∂ σ χ(σ)) dV = 0 (3.1)
The (3.

3)

The linearization of Eq. (3.1b) in discretized weak form leads to

Lin[G e σ ] = G e σ (τ, u) + ∆G e σ (∆τ , ∆u) (3.4)
In general, the set of governing equations can be computed with any selected value of α ∈ [0, 1].

However, here we choose the second-order mid-point scheme, by selecting α = 1/2.

σ n+α = (1 -α)σ n + ασ n+1 ; α = 1/2 ˙ vp n+α = β(σ n+α ) ⇒ vp n+1/2 = vp n + ∆t 2 β(σ n+1/2 ) ∂ σ χ (σ) = C -1 : σ (3.5)
From the above auxiliary results, the variation ∆G e σ (•) are derived as follows

∆G e σ (∆τ , ∆u) = δτ T Ω e S T BdV ∆d n+1/2 - Ω e S T D vp n+1/2 SdV ∆τ n+1/2 D vp n+1/2 := ∆t 2 β (σ n+1/2 ) + C -1 β (σ n+1/2 ) = ∂β(σ n+1/2 )/∂σ n+1/2 (3.6)
Similarly we also notice that δτ T is arbitrary virtual stress state so the local governing equation inside a discrete element Ω e is derived as follows

h e,(i+1) n+1/2 (τ , u) = h e,(i) n+1/2 (τ , u) + ∆h e,(i) n+1/2 (∆τ , ∆u) = 0 h e,(i) n+1/2 (•) := Ω e S T BdV d (i) n+1/2 - Ω e S T vp,(i) n+1/2 + ∂ σ χ(σ (i) n+1/2 ) dV ∆h e,(i) n+1/2 (•) := F e ∆d e,(i) n+1/2 -H e,(i) n+1/2 ∆τ e,(i) n+1/2 H e n+1/2 := Ω e S T D vp n+1/2 SdV (3.7)
At this point, one can establish the full system of governing equations, in which both nodal displacement and stress are unknown variables.

A nel e=1 F e,T A nel e=1 ηM e -H e n+1/2 F e ∆τ e,(i) n+1/2 ∆d e,(i) n+1/2 = - A nel e=1 r e,(i) u,n+1/2 (•) h e,(i) n+1/2 (•) (3.8)
In the way of programming this problem in an open development platform (such as FEAP-Finite Element Analysis Program) [Taylor, 2014], it is necessary to condense the previous system so that only nodal displacement is the unknown field. The stress field would be computed correspondingly from the displacement fields. From Eq. (3.8b), the stress state in each element can be rewritten as

∆τ e n+1/2 = H e,-1 n+1/2 h e,(i) n+1/2 (•) + H e,-1 n+1/2 F e ∆d e n+1/2
(3.9)

Substituting the above equation into Eq. (3.3a) yields the condensed system of equations, which depends only on the incremental displacements at time step t n+1/2 . Assembling each element contributions over entire discrete domain yields the governing equation of system in matrix form.

A nel e=1 F e,T n+1/2 H e,-1 n+1/2 F e + ηM e ∆d e,(i) n+1/2 = -A nel e=1 r e,(i) u,n+1/2 + F e,T H e,-1 n+1/2 h e,(i) n+1/2 ⇔ K ef f,(i) n+1/2 ∆d (i) n+1/2 = -R (i) n+1/2 (3.11)
In each iteration, the step of displacement is computed via Eq. (3.11). Subsequently, the step of stress is computed at the level of each element from equation Eq. (3.9). The stress and displacement fields are correspondingly updated after each computational iteration of Newton-Raphson method.

τ (i+1) n+1/2 = τ (i) n+1/2 + ∆τ (i) n+1/2 d (i+1) n+1/2 = d (i) n+1/2 + ∆d (i) n+1/2
(3.12)

Energy conserving scheme

An energy conserving (EC) scheme is developed to guarantee the stability of the computation over long period for visco-plasticity problem. To achieve the second-order accurate O(∆t 3 ) solution [START_REF] Artioli | Second-order accurate integration algorithms for von-mises plasticity with a nonlinear kinematic hardening mechanism[END_REF], the mid-point time integration scheme is employed. For enforcing the second-order accuracy in visco-plasticity, we must simultaneously integrate visco-plastic deformation with both displacement and stress, as opposed to the operator split algorithm ( [START_REF] Ibrahimbegovic | Classical plasticity and viscoplasticity models reformulated: theoretical basis and numerical implementation[END_REF]], [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF]) that separates the integration of these two and thus reduces the time-stepping scheme to first order (e.g. [START_REF] Chorin | Product formulas and numerical algorithms[END_REF]). Accordingly, the nodal velocity v n+1/2 and acceleration a n+1/2 at the time step t n+1/2 can be formulated in the form of nodal displacement and velocity.

v n+1/2 = (d n+1 -d n )/∆t a n+1/2 = (v n+1 -v n )/∆t (3.13)
The increment of displacement field within an interval of time step is computed by using the mid-point approximation.

d n+1 -d n = 1 2 ∆t(v n+1 + v n ) (3.14)
This equation can be formed in an alternative expression showing the simple update of displacement vector.

d n+1 = d n + u; u := 1 2 ∆t(v n+1 + v n ) (3.15)
The stress state at t n+1/2 is selected in the following algorithmic form. This proposal, which is left open at the moment, will be explained shortly.

τ alg n+1/2 := 1 2 (τ n+1 + τ n ) (3.16)
By choosing a test displacement vector δd = d n+1 -d n , a work done both the internal f int and external f ext can be elaborated from the weak form Eq. (3.3b) of balance equation.

(d n+1 -d n ) T f int n+1/2 = (d n+1 -d n ) T f ext n+1/2
The internal work includes two obvious components including kinetic and potential energy.

(d n+1 -d n ) T f int n+1/2 = (d n+1 -d n ) T Ω B T Sτ n+1/2 dV + (d n+1 -d n ) T Ω ρN T Na n+1/2 dV
The work done by inertia force can be simplified thanks to results obtained in Eq. (3.13b) and (3.14).

(d n+1 -d n ) T Ma n+1/2 = 1 2 (v n+1 + v n ) T M (v n+1 -v n ) = 1 2 v T n+1 Mv n+1 - 1 2 v T n Mv n = K n+1 -K n (3.17)
where K n+1 is the kinetic energy of system at time t n+1 . It is noted that the mass matrix M is symmetric, so this identity holds v T n+1 Mv n -v T n Mv n+1 = 0. In the followings, the subscript n + 1/2 for matrix F and H is dropped for shorter expression. The work done by potential force can be also expressed as a combination of 2 groups.

(d n+1 -d n ) T F T τ n+1/2 = 1 2 (d n+1 -d n ) T F T (τ n+1 + τ n ) = 1 2 d T n+1 F T τ n+1 - 1 2 d T n F T τ n + 1 2 d T n+1 F T τ n - 1 2 d T n F T τ n+1 L (3.18)
Since the scalar L has a mix variables from both time step t n and t n+1 , further simplification of L should be elaborated. Via exploiting the equation Eq. (3.7b) with a test stress function δτ T = τ n+1 -τ n , the following identity is derived based on the fact that the matrix H is always symmetric, even when the system goes into visco-plastic regime.

(τ n+1 -τ n ) T Fd n+1/2 -Hτ n+1/2 =0 ⇒ 1 2 d T n+1 F T τ n+1 - 1 2 τ T n+1 Hτ n+1 - 1 2 d T n F T τ n + 1 2 τ T n Hτ n = 1 2 d T n+1 F T τ n - 1 2 d T n F T τ n+1 L (3.19)
The matrix H can be decomposed in two matrices, the first one H vp n+1/2 is a function of viscoplastic strain rate ˙ vp n+1/2 and the other H e is a function of material elastic modulus. It is clear that the matrix H vp n+1/2 vanishes under the elastic regime, where the visco-plastic strain does not evolve ( ˙ vp n+1/2 = 0). With the last identity, Eq. (3.18) can be rewritten.

(d n+1 -d n ) T F T τ n+1/2 = 1 2 d T n+1 F T τ n+1 - 1 2 d T n F T τ n + 1 2 d T n+1 F T τ n+1 - 1 2 τ n+1 Hτ n+1 - 1 2 d T n F T τ n + 1 2 τ n Hτ n = d T n+1 F T τ n+1 - 1 2 τ T n+1 H e τ n+1 P n+1 -(d T n F T τ n - 1 2 τ T n H e τ n ) Pn + ( 1 2 τ T n H vp n+1/2 τ n - 1 2 τ T n+1 H vp n+1/2 τ n+1 D vp n+1/2 ) = P n+1 -P n + D vp n+1/2 (3.20)
where P n+1 is the potential energy of system at time t n+1 and D vp n+1/2 is the visco-plastic dissipation within an interval of time step. From the above derivation, one can express the balance of the internal energy and external work in following equation.

(K n+1 + P n+1 ) E n+1 -(K n + P n ) En +D vp n+1/2 = (d n+1 -d n ) T f ext n+1/2 ∆W n+1/2 E n+1 -E n + D vp n+1/2 =∆W n+1/2 (3.21)
The symbols E n and E n+1 denote the total energy of system at time step t n and t n+1 , respectively. We can conclude that the proposed algorithmic constitutive equation in Eq. (3.16) would conserve the total energy E n+1 = E n +∆W n+1/2 for any bounded external loading as long as the system undergoing elastic regime (D vp n+1/2 = 0). This property is viewed as "unconditionally stable" [START_REF] Zienkiewicz | The finite element method: its basis and fundamentals[END_REF]. Otherwise, once the system undergoes the visco-plastic regime, the total energy can not be conserved anymore due to the existence of visco-plastic dissipation (D vp n+1/2 = 0). However, this non-negative visco-plastic dissipation is always bounded via the controllable penalty factor as shown in Section 2.

Energy decaying scheme

An energy decaying (ED) scheme is an alternative approach to control potential instability problems of time-stepping algorithms, with an additional advantage of a providing superior accuracy in stress computations by eliminating the high-frequency noise especially for a set of stiff equations (e.g. see [START_REF] Ibrahimbegovic | Energy conserving/decaying implicit time-stepping scheme for nonlinear dynamics of three-dimensional beams undergoing finite rotations[END_REF]). The latter implies a large difference between max/min eigenvalues in tangent stiffness, which can also occur in the present case when one part of the domain has turned visco-plastic and the rest still remains elastic. In this case, we will enforce so-called a modified constitutive algorithm with the ability of decaying contribution of higher frequencies over each interval of time step. Such algorithmic updates for of displacement and stress at time step t n+1/2 are represented by modifying the conservative terms, as defined in Eqs. (3.14) and (3.16), with the corresponding dissipative terms:

d n+1 -d n := ∆t 2 (v n+1 + v n ) + α∆t(v n+1 -v n ) τ alg n+1/2 := 1 2 (τ n+1 + τ n ) + β(τ n+1 -τ n ) (3.22)
where α and β are the dissipation coefficients, which are typically chosen as α, β ∈ [0, 0.1]. The corresponding dissipative correction terms α∆t(v n+1 -v n ) and β(τ n+1 -τ n ) remain typically small, or practically vanish, for low frequency modes with a small difference of successive values within a typical time step; the latter establishes the consistency of the present modification in the spirit of the EC scheme, for the case where high frequency modes are not triggered. Meanwhile, the correction terms are less likely to vanish in high frequency modes, where the corresponding successive values within a time step can be quite significant. Therefore, the energy decaying scheme is only introducing a slight perturbation, which is not likely to fully reduce the scheme performance to the first-order accuracy. With such a choice of algorithmic update of displacement and stress updates, the velocity and stress at time step t n+1 take a modified form:

v n+1 = 1 (1/2 + α)∆t (d n+1 -d n ) - 1/2 -α 1/2 + α v n τ n+1 = 1 1/2 + β τ alg n+1/2 - 1/2 -β 1/2 + β τ n (3.23)
One can apply the same implementation as presented in Section 3.2. The work done by inertia force takes a new form.

(d n+1 -d n ) T Ma n+1/2 = 1 2 v T n+1 Mv n+1 - 1 2 v T n Mv n + α(v n+1 -v n ) T M (v n+1 -v n ) D K,n+1/2 = K n+1 -K n + D K,n+1/2 (3.24)
where D K,n+1/2 is the numerical dissipation of kinetic energy within an interval of time step.

In the followings, the subscript n + 1/2 for matrix F and H is dropped for shorter expression.

The work done by potential force takes a combination form of 3 groups.

(d n+1 -d n ) T F T τ n+1/2 = 1 2 (d n+1 -d n ) T F T (τ n+1 + τ n ) = 1 2 d T n+1 F T τ n+1 - 1 2 d T n F T τ n + 1 2 d T n+1 F T τ n - 1 2 d T n F T τ n+1 L + β(τ n+1 -τ n ) T H(τ n+1 -τ n ) D P,n+1/2 (3.25)
where D P,n+1/2 is the numerical dissipation of potential energy within an interval of time step. The similar identity to the one in Eq. (3.19) can be written as

(τ n+1 -τ n ) T Fd n+1/2 -Hτ n+1/2 =0 ⇒ 1 2 d T n+1 F T τ n+1 - 1 2 τ T n+1 Hτ n+1 - 1 2 d T n F T τ n + 1 2 τ T n Hτ n + β(τ n+1 -τ n ) T H(τ n+1 -τ n ) = 1 2 d T n+1 F T τ n - 1 2 d T n F T τ n+1 L (3.26)
With this identity on hand, the work done by potential force takes a new form.

(d n+1 -d n ) T F T τ n+1/2 = d T n+1 F T τ n+1 - 1 2 τ T n+1 H e τ n+1 P n+1 -(d T n F T τ n - 1 2 τ T n H e τ n ) Pn + ( 1 2 τ n H vp n+1/2 τ n - 1 2 τ n+1 H vp n+1/2 τ n+1 ) D vp n+1/2 + β(τ n+1 -τ n ) T H(τ n+1 -τ n ) D P,n+1/2 = P n+1 -P n + D vp n+1/2 + D P,n+1/2
(3.27) From the above derivation, the balance of the internal energy and external work can be formed in following equation.

(K n+1 + P n+1 ) E n+1 -(K n + P n ) En +D vp n+1/2 + D K + D P D KP,n+1/2 = (d n+1 -d n ) T f ext n+1/2 ∆W n+1/2 E n+1 -E n + D vp n+1/2 + D KP,n+1/2 =∆W n+1/2 (3.28)
where D KP,n+1/2 is the numerical dissipation of total energy due to the dissipative constitutive algorithm in Eq. 3.23. It is obvious that switching between energy conserving/decaying schemes can be executed flexibly without changing the main code of user-defined element since the update of displacement and velocity updates can be programmed in a separate subroutine for time-stepping integration algorithm. Additionally, the amount of numerical dissipation is fully controllable via an appropriate selection of dissipation coefficients α and β.

Representative numerical examples

In order to illustrate the favorable properties of the newly proposed EC and ED schemes for visco-plastic problem, several numerical simulations are presented and compared against more conventional visco-plasticity approach. All the computations are programmed as a userdefined element in computer code FEAP v8.4, see [Taylor, 2014]. The output data is then post-processed via Matlab scripts for plotting. The computations are performed for both quasi-static and dynamic pressure loads applied to thick-walled pipe by using the proposed finite element formulation based on Hellinger-Reissner variational pinciple, hybrid-stress interpolations and time-stepping schemes developed for viscoplasticity in the previous section. The computed results are compared against either the exact (analytic) solution in the elastic regime or the available results in visco-plasticity regime [Lubliner, 2008], as well as against the previous study results [START_REF] Hughes | Unconditionally stable algorithms for quasi-static elasto/visco-plastic finite element analysis[END_REF]. Additionally, the performance of the proposed element is compared against the quadrilateral rectangle element embedded B-bar method (Q4/P1).

Quasi-static loadings

The elasticity and visco-plasticity problem of a pipe under a quasi-static internal pressure is considered. In quasi-static analysis, the inertial term, proportional to mass matrix, is neglected in Eq. (3.8). The pipe is assumed very long compared to the outer radius, so that it is subjected to plane strain constraint with zero deformation along the pipe axis. Thanks to the axisymmetry of this 3D problem, the domain can be discretized by finite elements using the shape functions varying only in r direction. The inner and outer radius of pipe are chosen as r 1 = 1 and r 2 = 2. The properties of pipe material are listed as: elastic modulus E = 30.10 6 , Poisson's ratio ν = 0.3, yield stress σ y = 30.10 3 and 1/η = 10 -8 . Regarding the boundary conditions, on the inner end we apply pressure, while the outer end of the pipe is free of stress. The domain is discretized into 8 such elements, each with 3 nodes and 2 Gauss points as in Fig. 5a in r direction, which provide the second order displacement accuracy. This interpolation be handled with isoparametric elements using quadratic shape functions, see [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF]. Thus, the displacement u matrix and strain matrix are functions of nodal displacement d matrix, see [START_REF] Zienkiewicz | The finite element method: its basis and fundamentals[END_REF].

u = 3 i=1 N i u i → u = Nd; r θ z = ∂u ∂r u r 0 → = Bd (4.1)
The elasticity matrix in Eq. (3.5c) is rewritten for axisymmetric and plane strain conditions.

C = E (1 + v)(1 -2v)   1 -v v v v 1 -v v v v 1 -v   (4.2)
The interpolation of stress field is selected based on the equilibrium equations. The latter, expressed in polar coordinates with zero body forces and plane strain constraint, reduces to:

∂σ rr ∂r + 1 r (σ rr -σ θθ ) = 0; ∂σ θθ ∂θ = 0; ∂σ zz ∂z = 0 (4.3)
Since the quadratic shape functions for displacement field in this element result with linear strain field, we choose similarly the corresponding stress approximations. In particular, the stress σ rr should be a linear function to avoid the trivial case that it is a constant and equal to σ θθ . If the stress σ zz and σ θθ are constants in order to satisfy the Eqs. (4.3b,c), they can not satisfy the condition for plane strain (σ rr + σ θθ )ν = σ zz . As the result, all of three stress components are selected as linear functions. When the element gets smaller, Eq. ( 4.3) would be satisfied in a "weak" sense. The chosen approximation for stress field can be interpolated from 6 stress parameters, which become part of unknowns in this problem:

  σ rr σ θθ σ zz   =   1 r 0 0 0 0 0 0 1 r 0 0 0 0 0 0 1 r           τ 1 τ 2 τ 3 τ 4 τ 5 τ 6         ⇒ σ = Sτ (4.4)
All the ingredients of chosen discrete approximations are readily replaced into the theoretical formulation presented in Section 3. The analytical solution of stress field for a pipe under internal pressure p, see [Lubliner, 2008], is given as follows. As shown in Fig. 9a, the pipe under an elastic regime is analyzed by quasi-static analysis, in which the internal pressure p 0 load varies from 0 to 5.10 3 within 100 time steps. A coarse mesh with 8 by 8 elements is required by using enhanced quadrilateral element (Q4/P1) and the result is extracted at Gauss points laying on radial line AB, as shown in Fig. 4a. Meanwhile, a mesh with only 8 elements is adequate by using proposed element, as shown in Fig. 5a. The solution from the proposed element matches exactly the analytic solution while the stress computed by Q4/P1 element show a staggered trend, as shown in Figs. 2b-d. The total time of CPU (by a standard portable computer) to solve this problem is 0.38s by proposed element and 0.86s by Q4/P1 element. It is clear that the proposed element has advanced computational cost with respect to mesh generation and also computational time. To alleviate the high incompressibility constraint and enhance the accuracy of two-dimensional elements in this type of problem, one can also employ other methods such as enhancing displacement gradient by F-bar method (see [START_REF] Neto | Design of simple low order finite elements for large strain analysis of nearly incompressible solids[END_REF], Neto et al., 2005]) and/or generate a finer mesh. As shown in Fig. 3, the pipe under a visco-plastic regime is analyzed by quasi-static analysis, in which the internal pressure p 0 load stay constant at the value of 21.10 3 . After 15 seconds, the visco-plastic zone develops to the middle of the pipe. This current result is simulated with time step ∆t = 10 -4 . The zone with F + 1 > 1 equivalent to φ > 0 , which has non-zero visco-plastic deformation, corresponds to the region with r 1 ≤ r ≤ r p = 1.4r 1 , as shown in Fig. 3a. Meanwhile the elastic zone is located at r p ≤ r ≤ r 2 . The current solution matches with previous numerical study by Hughes and Taylor, see [Hughes and Taylor, 1978], which is in agreement with the exact solution given by Prager and Hodge with the ratio r p /r 1 = 1.4, see page 271 of [Lubliner, 2008]. Compared to elastic analysis, the similar trend of radial stress can be observed, as shown in Fig. 3b. Meanwhile, the circumferential and longitudinal stresses show clearly two different zones due to the existence of visco-plasticity, as shown in Fig. 3c,d. The stress field is also compared to that from Q4/P1 element under perfectly plastic regime. The location r p of separation line is the same as that from the proposed element, as shown in Fig. 4b. It is also noticed that the solution from Q4/P1 element show a staggered trend, as shown in Fig. 3b-d. For the proposed element, the visco-plastic problem is solved within 3 iterations under quadratic convergence rate (energy residual up-to 10 -15 ). The proposed element shows an advanced smooth and accurate stress field compared to Q4/P1 element in the pipe.

σ r = -p (r 2 /r 1 ) 2 -1 r 2 2 r 2 -1 ; σ θ = p (r 2 /r 1 ) 2 -1 r 2 2 r 2 + 1 ; σ z = 2νp (r 2 /r 1 ) 2 -1 (4.5) (

Dynamic loadings

The visco-plastic problem inside selected thick-walled pipes is simulated under dynamic loads. There are two models with different radius for the simulation. The first one has r 1 = 1 and r 2 = 2 referred as "soft" model, while the other has r 1 = 5 and r 2 = 10 referred as "stiff" model. The difference between two concerns the corresponding ratio between the smallest and largest eigenvalues of the tangent stiffness. The same number of elements is used to perform the computations for both models. The results presented subsequently for cases 1-3 belong to the "soft" model and for cases 4-6 belong to the "stiff" model. The chosen time step in all simulations is ∆t = 10 -5 , which is able to capture high frequency response.

Regarding the energy conserving (EC) scheme, only one triangle pulse with a peak of pressure p 0 is introduced to both models as shown in Fig. 5a, which subsequently leads to the free vibration of system. Several cases with different pulse loading duration are examined, as listed in Table 1. Regarding the energy decaying (ED) scheme, an additional pulse of short duration is later introduced into the system. This pulse is applied on end-nodes of each element with opposite direction, as shown in Fig. 5b, in order to activate high frequency modes of the system. The results are extracted at nodes 1, 9 and 17 as shown in Fig. 5a and 1,2,3 21e2 8e-4 30e-4 31e-4 0.01,0.05,0.1 2 4,5,6 21e2 8e-4 30e-4 31e-4 0.01,0.05,0.1 Visco-plastic regime (ED scheme) 1 1,2,3 42e3 8e-4 30e-4 31e-4 0.01,0.05,0.1 2 4,5,6 42e4 8e-4 30e-4 31e-4 0.01,0.05,0.1 Table 2: Load groups for ED scheme The results in Fig. 6 show that the total internal energy of system is conserved exactly under free vibration. The case 1-3 is from soft model and case 4-5 is from stiff model. After releasing the external load, the total energy of system, including kinetic and potential energies, remains constant and coincides with the external energy introduced to system by the pulse. For example, the total energy of system is around 0.09 in case 1 and 0.02 in case 4, as shown in Figs. 6a,b. Longer duration of loading leads to higher total energy gained by the system. It is noted that these cases are fully under elastic regime without any evolution of visco-plasticity strain. The dynamic response of system including displacement, velocity and stress are recorded at node number 1, 9 and 17, which are, respectively, at the inner radius, in the middle and at the outer radius of the mesh, as shown in Fig. 5a. As shown in Fig. 7, the amplitude of displacement at node 1 is higher than those at other nodes since the pulse is applied at this node. It is noted that the displacement field in soft model is dominated by one frequency, whereas this is no longer the case in stiff model. Fig. 8: Velocity under elastic regime (EC scheme)

A similar observation is made regarding the velocity and stress fields, as shown in Figs. 8 and9, where higher values of absolute displacement leads to higher stress. During the elastic regime, the norm of deviatoric part of largest stress || dev(σ)|| always remains lower than the critical value 2 3 σ y ≈ 24.5 × 10 3 (see Eq. (2.6)), which is equivalent to φ < 0. The same numerical examples are further examined under visco-plastic regime. Here, higher values of pressure pulse have to be introduced in order to activate visco-plastic regime. In particular, the stiff model requires a higher pulse amplitude compared to the soft model. During visco-plastic regime, the total energy of system is not conserved due to the evolution of viscoplastic dissipation, as shown in Fig. 10. However, since the visco-plastic dissipation remains rather small, the total energy of system is approximately comparable to the total external work after releasing the pulse. For example, the total external work is around 38 in Case 1 and 800 in Case 4, as shown in Figs. 10a,b. The visco-plastic dissipation increases overall; however it stays constant sometimes when there is no excursion in visco-plastic regime. This trend can be explained clearly by corresponding variations of the stress field, as in Fig. 11.

The stresses are extracted from only first Gauss point of elements 1, 4 and 8 (numbered starting from the inner radius), as shown in Fig. 5b. During the period that every element has the value || dev(σ)|| < 2 3 σ y (or equivalent to φ < 0), the whole system is under elastic regime and there is no evolution of visco-plastic strain, which keeps the corresponding viscoplastic dissipation constant. By contrast, the visco-plastic dissipation evolves when the yielding criteria is satisfied, precisely || dev(σ)|| > 2 3 σ y . This energy dissipation depends on type of material in which the viscosity coefficient η is taken into account. In Fig. 12, we show the total energy of system in the elastic regime, computed by ED scheme. The only dissipation sources are introduced by energy-decaying scheme, where higher frequency dissipation is controlled by dissipation factors α, β that regulate potential and kinetic dissipation energies. In both soft and stiff models, there are 3 cases with different dissipation factors (α, β) ∈ {(0.01, 0.01), (0.05, 0.05), (0.1, 0.1)}. It is obvious that the dissipative energy increases proportionally with respect to the chosen values of dissipation factors (α, β) as shown in Figs. 12a,b. In case 3 of soft model, the total energy decreases slightly due to the dissipative energy over time as shown in 12c. The same pattern can be better recognized after the injection of high frequency modes contribution produced by the second pulse. This high frequency pulse does not dissipate neither kinetic nor potential energy of the soft model as much as it dissipates in case 6 of the stiff model. The noise generated by the second pulse is dissipated out quickly for both kinetic and potential energies, as shown in Fig. 12d. When dissipation factors (α, β) approach zeros, Figs. 12a,b will agree with Figs. 6e,f in the period before introducing the second pulse at ∆t 2 = 30.10 -3 .

In Fig. 13, the stresses are taken again in elements 1, 4 and 8 at the Gauss points as in previous simulation. As shown in Figs. 13a,c,e, the second pulse in soft model does not bring much noise to the stress, except for the an obvious effect that it increases the magnitude of stresses. However, the perturbation of high frequency modes can be viewed from cases 4-6 in the stiff model. As shown in Fig. 13f, high frequency stresses are quickly filtered out of the system with α = β = 0.1, which is faster than the other cases with smaller values of α and β, as shown in Fig. 13b,d As shown in Fig. 14, the total energy of system under visco-plastic regime are computed by ED scheme. There are in total 3 sources of dissipation, namely visco-plastic dissipation, numerical dissipation of energy-decaying scheme for kinetic and for potential energy. The direct correlation between the higher values of α and β versus resulting numerical dissipation is obviously captured in Figs. 14a and14b. The visco-plastic dissipation increases in time whenever the system satisfies the yield criteria, as shown in Figs. 14c and14d. It is recognized that numerical dissipation is much higher than for the visco-plasticity. From the results in Fig. 12a,b and Fig. 14a,b, the trend and relative dissipation of system due to ED scheme is comparable between elastic and visco-plastic regimes.

The stresses are taken again in elements 1, 4 and 8 at the first Gauss points, as presented in Fig. 15. Similar to the previous simulation, the second pulse bring only little noise to the stress in soft model, as shown in Figs. 15a,c,e. The same observation on the energy-decaying on high frequency mode's contribution no longer applies for the stiff model, as shown in Figs. 15b,d,f The stiff model is discretized via 16 elements as in Figs. 16a,c and via 32 elements as in Figs. 16b,d, respectively. The energy-decaying scheme results obtained with dissipation factors (α, β) ∈ {(0.01, 0.01), (0.1, 0.1)} are then evaluated with the two above discrete models, where only elastic regime is considered. In Figs. 16a,c, the higher values of α and β are chosen, the faster higher frequency modes are filtered out for the stiff system with coarse mesh. In contrast, the lower values of α and β are chosen, the slower higher frequency modes are filtered out from the system with fine mesh, as in Figs. 16b,d. Moreover, any finer mesh can better resolve high frequencies, and thus produces smaller numerical dissipation and less significant on amplitude decay due to higher-frequency modes. From the above observations, the most proper choice for the numerical dissipation factors α and β depends both on the type of system (soft vs. stiff), and on the element size in the mesh. This best choice would be left to the users depending upon the kind of results that are the most acceptable for them. Practically, it is recommended that the higher values of dissipation factors (α, β) should be chosen for a system with coarser mesh than those for a system with finer mesh.

The visco-plastic problem itself dissipates energy through the visco-plastic strain. The energy decaying scheme can be employed if the user expects that the external forces can lead to a significant high frequency modes, in which the contribution from these modes should be dissipated especially on a coarse mesh where they cannot be properly resolved. Meanwhile, the energy conserving scheme can be selected within any case, regardless of elastic or visco-plastic regimes. It is recorded that the dynamic problem coupled with visco-plastic regime is solved within 3 iterations under quadratic convergence rate (energy residual up-to 10 -15 ). The solution method for visco-plasticity problem proposed herein validates the choice of the variational principle and the hybrid-stress finite element approximation to deliver a superior accuracy of stress field. Namely, by considering both displacement and stress fields as independent variables does increase the problem size with respect to conventional displacement-based finite element approach, but it also provides the possibility for using coarse mesh with high stress accuracy and a very efficient solution employing consistent tangent elastoviscoplastic modulus.

The main thrust in this work is further directed to dynamics problem by developing the energy conserving (EC) and the energy decaying (ED) schemes. Both of these schemes are capable of controlling stability of numerical computations. The proposed algorithm shows very satisfying results for a number of numerical simulations, where the computation's stability is enforced by preserving the total energy in elastic regime under EC scheme, or by eventual numerical dissipation that is proved capable of decaying of high frequencies modes which cannot be resolved properly on a coarse finite element mesh under ED scheme. The dissipation parameters α and β in ED scheme can be chosen with respect to mesh size. Practically, the higher values of α and β yield faster dissipation of inappropriate high-frequency modes for a system with coarse mesh. By contrast, these parameters can be selected relatively lower for a system with fine mesh which has fewer unresolved high frequencies and thus automatically reduce the numerical dissipation.

Thanks to the fast convergence and stability of the time-stepping scheme, the proposed formulation can have a potential for other non-linear dynamic problems. In particular, the hybrid-stress interpolation can provide not only an enhanced accuracy of stress field, but also the desirable smoothness of stress approximation, which is of special interest for capturing the wave propagation phenomena.

nFig. 1 :

 1 Fig. 1: The domain Ω and its boundary ∂Ω with indicated boundary conditions

Fig. 2 :

 2 Fig. 2: Computed stress under elastic regime at p 0 = 5.10 3

Fig. 3 :

 3 Fig. 3: Computed stress under visco-plastic regime by constant p 0 = 21.10 3

  1

Fig. 6 :

 6 Fig. 6: Energy conserving under elastic regime (EC scheme)

Fig. 7 :

 7 Fig. 7: Displacement under elastic regime (EC scheme)

Fig. 9 :

 9 Fig. 9: Stress under elastic regime (EC scheme)

Fig. 11 :

 11 Fig. 11: Stress under visco-plastic regime (EC scheme)

Fig. 13 :

 13 Fig. 13: Stress under elastic regime (ED scheme)

Fig. 15 :

 15 Fig. 15: Stress under visco-plastic regime (ED scheme)

Fig. 16 :

 16 Fig. 16: Stress under elastic regime (stiff model-ED scheme)

  at elements 1, 4 and 8 as shown in Fig. 5b. Several cases with different values for dissipation factors (α, β) are shown in Table2.

	N1	...	N9	...	N17
							E1	...	E4	...	E8
	(a) One pulse: EC scheme		(b) Two pulses: ED scheme
				Fig. 5: Mesh and pressure pulses
		Simulation	Model Case	p 0	∆t 1
	Elastic regime	1		1,2,3 21e2 4e-4,6e-4,8e-4
	(EC scheme)		2		4,5,6 21e2 4e-4,6e-4,8e-4
	Visco-plastic regime	1		1,2,3 42e3 4e-4,6e-4,8e-4
	(EC scheme)		2		4,5,6 42e4 4e-4,6e-4,8e-4
			Table 1: Load groups for EC scheme
	Simulation		Model Case	p 0	∆t 1 ∆t 2	∆t 3	α = β
	Elastic regime
	(ED scheme)

  .

		5000						400
				|| || 1.1	|| || 4.1	|| || 8.1				|| || 1.1	|| || 4.1	|| || 8.1
		4000						300
	||dev( )||	2000 3000					||dev( )||	200
		1000						100
		0	0	2	4	6	8	0	0	2	4	6	8
					Time (s)		10 -3				Time (s)		10 -3
				(a) Case 1					(b) Case 4
		5000						400
				|| || 1.1	|| || 4.1	|| || 8.1				|| || 1.1	|| || 4.1	|| || 8.1
		4000						300
	||dev( )||	2000 3000					||dev( )||	200
		1000						100
		0						0
			0	2	4	6	8		0	2	4	6	8
					Time (s)		10 -3				Time (s)		10 -3
				(c) Case 2					(d) Case 5
		5000						400
				|| || 1.1	|| || 4.1	|| || 8.1				|| || 1.1	|| || 4.1	|| || 8.1
		4000						300
	||dev( )||	2000 3000					||dev( )||	200
		1000						100
		0						0
			0	2	4	6	8		0	2	4	6	8
					Time (s)		-3				Time (s)		10 -3
				(e) Case 3					(f) Case 6

  .

		10	10 4					8	10 4
				|| || 1.1	|| || 4.1	|| || 8.1				|| || 1.1	|| || 4.1	|| || 8.1
		8						6
	||dev( )||	4 6					||dev( )||	4
		2						2
		0	0	2	4	6	8	0	0	2	4	6	8
					Time (s)		10 -3				Time (s)		10 -3
				(a) Case 1					(b) Case 4
		10	10 4					8	4
		8		|| || 1.1	|| || 4.1	|| || 8.1		6		|| || 1.1	|| || 4.1	|| || 8.1
	||dev( )||	4 6					||dev( )||	4
		2						2
		0	0	2	4	6	8	0	0	2	4	6	8
					Time (s)		10 -3				Time (s)		10 -3
				(c) Case 2					(d) Case 5
		10	10 4					8	10 4
				|| || 1.1	|| || 4.1	|| || 8.1				|| || 1.1	|| || 4.1	|| || 8.1
		8						6
	||dev( )||	4 6					||dev( )||	4
		2						2
		0						0
			0	2	4	6	8		0	2	4	6	8
					Time (s)		10 -3				Time (s)		-3
				(e) Case 3					(f) Case 6

Acknowlegdement

This work is financially supported under the research project SELT-TUM by Agence Nationale de la Recherche (ANR). Moreover, Prof. Adnan Ibrahimbegovic is also supported by Institut Universitaire de France (IUF). These sources of funding are gratefully acknowledged.


Radhouan El Zrelli
email: radhouan.elzrelli@gmail.com

Lotfi Rabaoui

Rubén H Roa-Ureta

Nicola Gallai

Sylvie Castet

Michel Grégoire

Nejla Bejaoui

Pierre Courjault-Radé

Economic impact of human-induced shrinkage of Posidonia oceanica meadows on coastal fisheries in the Gabes Gulf (Tunisia, Southern Mediterranean Sea)

Keywords: Gabes Gulf, Posidonia oceanica, Ecosystem services, Coastal fisheries, Economic valuation, Phosphogypsum

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Over 40% of the oceans are impacted by pollution [START_REF] Halpern | Spatial and temporal changes in cumulative human impacts on the world's ocean[END_REF]. In the last fifty years, the number of marine dead zones has doubled with each passing decade [START_REF] Diaz | Spreading dead zones and consequences for marine ecosystems[END_REF][START_REF] Limburg | Ocean Deoxygenation: A Primer[END_REF]. Marine ecosystems contaminated by pollutants such as toxic runoff, dredged materials, oil spills, trash and debris, industrial wastes and sewage sludge, are persistently disrupted [START_REF] Rabaoui | Anthropogenicinduced acceleration of elemental burial rates in blue carbon repositories of the Arabian Gulf[END_REF]. The Mediterranean basin is among the most impacted regions in the world together with the North Sea, South and East China Seas, east coast of North America, Red Sea and Arabian Gulf [START_REF] Halpern | A global map of human impact on marine ecosystems[END_REF][START_REF] Halpern | Spatial and temporal changes in cumulative human impacts on the world's ocean[END_REF][START_REF] Halpern | Recent pace of change in human impact on the world's ocean[END_REF]. Consequently, the socio-economic status of some coastal fisheries has become threatened by habitat degradation [START_REF] Pedro | Shifting baselines in a Mediterranean small-scale fishery[END_REF]. Such an environmental and socioeconomic status regarding coastal fisheries is the rule for all countries bordering the Mediterranean basin, where more than 250 million people live on the coasts (UNEP, 2012) along with more than 230 million tourists per year [START_REF] Pons | Culture of mass tourism -Doing the Mediterranean in the age of banal mobilities[END_REF]. The United Nations Environment Program (UNEP) has estimated that 650 million tons of sewages, 129 thousand tons of mineral oil, 60 thousand tons of mercury, 36 thousand tons of phosphates and 380 tons of lead are dumped into the Mediterranean Sea each year (UNEP, 2012).

The good health status of essential fish habitats like seagrass ecosystems is an essential factor for the maintenance of Mediterranean coastal fishing economic activity because of its ecosystemic functions for benthic and demersal commercial species [START_REF] Unsworth | Seagrass meadows support global fisheries production[END_REF].

Seagrass meadows are among the most valuable coastal ecosystems worldwide in terms of goods and services they provide [START_REF] Costanza | The value of the world's ecosystem services and natural capital[END_REF][START_REF] Hemminga | Seagrass ecology[END_REF].

According to [START_REF] Waycott | Accelerating loss of seagrasses across the globe threatens coastal ecosystems[END_REF], almost one third of the areal extent of seagrass has disappeared globally since the end of the XIX th century. In the Mediterranean basin, Posidonia oceanica is the most important endemic seagrass species and it forms meadows or beds covering an area between 21 and 55 thousand km 2 [START_REF] Pasqualini | Mapping of Posidonia oceanica using aerial photograph and side-scan sonar: Application off the Island of Corsica (France)[END_REF]. Among the ecosystemic services provided by P. oceanica is the fishing income [START_REF] Costanza | The value of the world's ecosystem services and natural capital[END_REF][START_REF] Campagne | The seagrass Posidonia oceanica: Ecosystem services identification and economic evaluation of goods and benefits[END_REF]Vasallo et al., 2015;[START_REF] Boudouresque | The necromass of the Posidonia oceanica seagrass meadow: fate, role, ecosystem services and vulnerability[END_REF]. Direct monetary assessment of coastal fisheries through transformation of benthic and demersal species volumes into financial value testify of the concrete human benefits derived from ecological functions supplied by the seagrass meadows. Consequently, they can be considered as an ecosystem service sensu [START_REF] Liquete | Integrated valuation of a nature-based solution for water pollution control: Highlighting hidden benefits[END_REF] and as "supporting services" sensu MEA (2005). A first global valuation of ecosystemic services of all types of seagrass at the worldwide scale has been estimated at about 19,000 US$ ha -1 y -1 [START_REF] Costanza | The value of the world's ecosystem services and natural capital[END_REF]. After this founding work, economic valuations of these systems have been sparse (Dewsburry et al., 2016). This is all the truer in the case of the Mediterranean Sea, where few papers have addressed the value of ecosystem services provided by P. oceanica seagrass (Vasallo et al. 2013;[START_REF] Campagne | The seagrass Posidonia oceanica: Ecosystem services identification and economic evaluation of goods and benefits[END_REF].

The seagrass meadows of the Gulf of Gabes (SE Tunisia; Fig. 1) used to be considered as the most extensive beds in the Mediterranean Sea (De Gaillande, 1970a). But since the beginning of XX th century, they declined drastically because of multiple anthropogenic pressures (e.g. [START_REF] Brahim | Variability in the structure of epiphytic assemblages of Posidonia oceanica in relation to human interferences in the Gulf of Gabes, Tunisia[END_REF][START_REF] El Zrelli | Biomonitoring of coastal pollution in the Gulf of Gabes (SE, Tunisia): use of Posidonia oceanica seagrass as a bioindicator and its mat as an archive of coastal metallic contamination[END_REF]. According to Zaouali et al. (1993), the decline of P. oceanica meadows from the Gulf of Gabes has occurred in three main phases: the first phase started with the development of the destructive dredge (or gangava), used to collect sponges, which had induced a first local degradation of the meadows leading to the initiation of a slow silting process in the central area of the Gulf [START_REF] Poizat | Hydrodynamisme et sédimentation dans le golfe de Gabès (Tunisie)[END_REF]. The second phase happened between the mid and late sixties when a major shrinkage of Posidonia beds took place following the development of bottom trawling techniques for shrimping. This led to limit the Posidonia-covered areas to coastal strips only and to a depth rarely exceeding 20-25 m, and consequently has triggered the silting process from the centre of the Gulf to the coast [START_REF] Poizat | Hydrodynamisme et sédimentation dans le golfe de Gabès (Tunisie)[END_REF][START_REF] Ktari-Chakroun | Les fonds chalutables de la région sud-est de la Tunisie (golfe de Gabès)[END_REF]. As for the third phase, it began in the mid-seventies when the remaining coastal strip of P. oceanica beds gradually disappeared around the first unit of Gabes phosphate fertilizer plant, put into operation in 1972 [START_REF] Darmoul | Effets des rejets industriels de la région de Gabès sur le milieu récepteur[END_REF][START_REF] Darmoul | Pollution dans le Golfe de Gabès (Tunisie) : bilan de six années de surveillance[END_REF]. By the second half of the seventies, the large quantities of dumped phosphogypsum wastes (≈14 400 T. d -1 of dry PG;

El Zrelli et al., 2018a) in the surrounding marine environment were responsible for degradation of Posidonia beds in front of the industrial complex due to the increase of turbidity coupled with the radioactive and chemical contamination of seawater [START_REF] Darmoul | Effets des rejets industriels de la région de Gabès sur le milieu récepteur[END_REF][START_REF] Darmoul | Pollution dans le Golfe de Gabès (Tunisie) : bilan de six années de surveillance[END_REF]El Zrelli et al., 2018b;El Zrelli et al., 2019a). Several studies carried out across the central area of Gabes have confirmed a continuous degradation of the coastal belt of P. oceanica triggered by phosphogypsum releases [START_REF] Pergent | L'environnement marin côtier en Tunisie[END_REF]Zaouali, 1993;[START_REF] El Zrelli | Biomonitoring of coastal pollution in the Gulf of Gabes (SE, Tunisia): use of Posidonia oceanica seagrass as a bioindicator and its mat as an archive of coastal metallic contamination[END_REF]. This continuous anthropogenic-induced regression of P. oceanica meadows has led to significant changes in the biodiversity of the impacted coastal area over the last few decades, causing thus significant changes in the catch composition and incomes of the coastal fisheries sector. In this context, several studies reported that the fisheries resources of the Gulf of Gabes are highly exploited [START_REF] Hattab | An ecosystem model of an exploited southern Mediterranean shelf region (Gulf of Gabes, Tunisia) and a comparison with other Mediterranean ecosystem model properties[END_REF][START_REF] Halouani | A spatio-temporal ecosystem model to simulate fishing management plans: A case of study in the Gulf of Gabes (Tunisia)[END_REF][START_REF] Halouani | An end-to-end model to evaluate the sensitivity of ecosystem indicators to track fishing impacts[END_REF], thus threatening the socioeconomic situation of Gabes coastal fisheries (FAO ArtFiMed, 2011a). In 1993, ~39 thousand people were directly employed in the coastal fisheries segment and only ~20 thousand in 2009 (FAO ArtFiMed, 2011a). Production varied from ~36 thousand tons in [START_REF] Sakamoto | Akaike information criterion statistics[END_REF]-1990[START_REF] Pergent | L'environnement marin côtier en Tunisie[END_REF] to less than 14 thousand tons in 2009 (FAO ArtFiMed, 2011a). No studies have been conducted so far on the changes on coastal fisheries neither on the valuation of economic changes related to the losses in P. oceanica seagrass habitats. This work was conducted with the aims i) to describe the decadal changes in the catch and composition of coastal fisheries in relation to the decline of P. oceanica beds and increase of phosphogypsum wastes discharges, and ii) to assess the losses of economic value proceedings from goods and benefits provided by coastal fisheries in connection with the decline of P. oceanica meadows in the Gulf of Gabes.

Materials and Methods

Study area

The Gulf of Gabes extends between Ras Kaboudia in the north and the Libyan/Tunisian boundary in the south (Fig. 1a). It is the widest area of continental shelf in the Mediterranean Sea together with the North-Central Adriatic shelf. It covers about 36 thousand km 2 [START_REF] Hattab | An ecosystem model of an exploited southern Mediterranean shelf region (Gulf of Gabes, Tunisia) and a comparison with other Mediterranean ecosystem model properties[END_REF]. The study area concerns the central part of Gulf of Gabes which extends between the north-western coast of Djerba Island and El Akarit covering more than 100 km of coastline (Fig. 1). This area hosts three main landing sites: Ghannouch, Gabes and Zaratt (Fig. 1a). Among these, Gabes landing site is the most important since its commercial landings represent more than 80% of the total commercial landings from the 3 sites. They are subdivided into "small pelagics" and "coastal fisheries". Since 1984, the small pelagics group represents an average of ≈ 80% of the commercial landings of Gabes governorate.

Fisheries Data collection

Annual volumes of landed catches of coastal and small pelagic fisheries were collected, over 36 years extending from 1980 to 2016. The 1980-1993 data were collected from the General Fisheries Commission in Tunisia (Commission Générale des Pêches de Tunisie; CGP, 1996) ; while the 1994-2015 data were provided by the Tunisian Directorate-General for Fisheries and Aquaculture (Direction Générale de la Pêche et de l'Aquaculture: DGPA; http://www.onagri.nat.tn). Data collected were checked from various reports and since the catch volumes of coastal fisheries were not solid enough for the two years 1984 and 1986, these two records were excluded from the economic analyses (Appendix I).

Unfortunately, no official records of coastal and small pelagic fisheries catches were available before the eighties. According to FAO ArtFiMed (2011b), the coastal benthic/demersal landed volumes, before the eighties, were low and mostly composed of cephalopods (e.g. the common octopus, Octopus vulgaris) sold mainly in local markets (M'Rabet et al., 2011). At that time, there were no estimations for volumes sold in local fish markets or indirectly by fishermen.

Economic data collection

Official economic data were obtained from the CGP, between 1990 and 2014. The starting year of 1990 was determined by (i) uncertain reliability of data before this date, (ii) its coincidence with the beginning of the second phase of landed volumes evolution (Appendix I), and (iii) with the end of the quick degradation phase of Posidonia oceanica meadows. The market prices of fish sold in Gabes market were in Tunisian Dinars (TND). These prices were updated (http://mon-convertisseur.fr) and converted to Euros (www.xe.com/fr/currencyconverter/).

Prediction of the Area of Posidonia meadows over the period 1970 to 2014

The quantification of marine area covered by Posidonia meadows since the 60s has been established for four years using cartographic data of previous works (Appendix II).

Posidonia maps have been digitized and their areas has been estimated using ArcMap software (version 10.2.2, Fig. 1b-e). Although this is a small sample, it is informative because it covers a long period of time and it has a clear trend of linear shrinking (Fig. 2a). Taking advantage of this informative sample, we fitted a linear regression model of Posidonia meadows area versus the year to predict the extension of these meadows in the years without measurements. These predictions and their prediction standard errors were then used for determining the effect of the Posidonia meadows shrinkage on fisheries output.

Relationship between P. oceanica meadow extension and discharges of phosphogypsum

We built a database with the annual meadow extension prediction and discharge of phosphogypsum in megatons into the Gulf of Gabes spanning the period of 1972 to 2016.

Assuming that discharges of phosphogypsum have an additive effect, we examined the relation between cumulative annual discharges and the extension of Posidonia meadows with generalized linear models (glm). Taking into account that Posidonia meadow extension is a right-skewed positive random variable, we fit the glm assuming a gamma distribution with the log link function. The Gamma distribution model is appropriate for positive, continuous and right-skewed data with a variance that is nearly constant on the log scale.

Since meadow extension was not directly measured data but rather a prediction from the linear model described in §2.4, it was necessary to account for the prediction error when fitting the generalized linear model. To that effect we fitted the glm to each one of 1000 realizations of the predicted meadow extension, with each realization coming from a truncated normal distribution defined by a mean (μ) equal to the predicted extension, a standard deviation (σ) equal to the predicted standard error of prediction, and bounds of [μ-0.25σ, μ+0.25σ]. This marginalization of the distribution of the meadow extension data accounted for some, though not all, the uncertainty from replacing actual meadow extension data with the predictions from the linear regression between year and the four available direct measurements of meadow extension.

Economic analysis

Our goal is to assess the economic impact on commercial coastal fisheries in the central part of Gulf of Gabes over the last twenty-five years in relation to the gradual shrinkage of P. oceanica seagrass meadows. We value the provisioning services by P.

oceanica meadows sold on official Tunisian market using the market price-based approach (Silvis and Van der Heid, 2013), where the main output of the service offered by P. oceanica is the production of benthic and demersal commercial fishery stocks. This method consisted in valuation of the part of benthic/demersal species value dependent on Posidonia. We assume that the production in tons of benthic/demersal species at a time t (Qt) is as follows:

t Q = a t p t = b t z t (1)
where p is the area of P. oceanica per hectare (ha) at time t, z the set of ecological factors that contribute to the production of benthic/demersal species per ha, a and b are the corresponding yield coefficients. The total value of Posidonia corresponds to Qt multiplied by the market price of benthic/demersal at time t (mt), =

From the equation 2, we can determine the hectare value of Posidonia at time t:

= (3)
We also measured the economic losses due to the degradation of Posidonia (ELDP) meadows in the Gulf of Gabes. For that, we compared the maximum potential quantity of Posidonia, P, with the quantity pt, which is the quantity of Posidonia at time t. Then the loss value (lt) can be measured by the following formula:

= ( -) (4)

Fishing effort modelling (2000-2016)

Fishing effort data, including the number of boats and number of fishermen, were collected on yearly basis from the annual reports of fisheries statistics published by DGPA.

To relate fisheries output to P. oceanica meadows extension, we built a database consisting of annual fisheries output in tons of landings, fishing effort in number of boats operating per year, the prediction of meadows extension and the prediction standard error. The fishery effort and output were split between two fleets operating in the Gulf of Gabes, the coastal fleet and the small pelagic fish fleet. This database covered the years between 2000 and 2016. The limiting factor was the unavailability of fishing effort of both fleet components before 2000.

Generally, the raw fishing yield data show differences between the coastal and the small pelagic fish fleets (Fig. 2b). The former shows a fairly constant annual yield around 2-3 thousand tons and a rise in effort over the last few years, while the latter shows a decreasing yield and constant effort (Fig. 2c).

To test the effect of meadows extension, we fitted a glm with the fisheries output as response variable and the fishing effort and the meadows extension as potential predictor variables, each fleet component separately. The fishing output data was assumed to have a gamma probability distribution with a log link function. As in the previous analyses, we accounted for the uncertainty in the prediction of meadows cover by marginalizing over 1000 realizations of the annual series of meadow shrinkage with the same parameters used to generate the data for the modelling of meadows extension versus cumulative phosphogypsum discharge.

To test the effect of the meadow extension on the fishing output, we fitted a glm where both the fishing effort and the meadow extension were included as predictors versus a model where the fishing effort was the single predictor. We then calculated the Akaike Information Criterion (AIC) difference between the model with both predictors versus the model with fishing effort as single predictor. Under the hypothesis that the extension of meadows did not contribute to the fishing output in addition to the effect of fishing effort, the distribution of the AIC difference across the 1000 realizations would be symmetrical and centered on 2 [START_REF] Sakamoto | Akaike information criterion statistics[END_REF]. Conversely, under the hypothesis of the effect of the meadow extension on fishing output the distribution of the AIC difference would be skewed to the left, and the stronger the effect the stronger the skewness.

Results

3.1. Historical decline of P. oceanica P. oceanica meadows have shrunk drastically along the last 5 decades. In spite of having only four data points the decline was well parameterized by the linear regression prediction, with rather narrow confidence bands (Fig. 2a). The reasons for this success with just four data points were that these points are well spaced in time and they are falling neatly on the linear trend (Fig. 2a). At the estimated rate of shrinkage, the meadows will disappear any year now, although the precision of this prediction allows for residual patches well into the future.

Economic valuation

Table 1 summarizes the main results of the economic valuation of Posidonia in the Gulf of Gabes in 1990 and 2014. First, we observed that the value has increased by 263% in the period, from €2 million to more than €7.5 million. This increase is due to the high increase in market price of benthic species combined with the rise of benthic species production. In the same periods, the yield of fish has varied from 0.024 to 0.17 per hectare of Posidonia. Thus, the value of Posidonia per hectare was €48 in 1990 and more than €500 in 2014. Between 1990 and 2014, the loss of cover of Posidonia meadows in the Gulf of Gabes has resulted in a loss of value from €4 million initially (1990) to almost €60 million currently (2014).

The phosphogypsum-Posidonia-fisheries connection

The hypothesis of the effect of cumulative phosphogypsum discharge on predicted Posidonia cover is overwhelmingly more supported than the alternative of no effect (Fig. 3a). This is noted in much lower AIC values for the hypothesis that included the cumulative phosphogypsum discharge effect versus the null model that did not include that effect (Fig. 3a). In addition, the largest p-value of the null hypothesis of no cumulative phosphogypsum discharge effect across 1000 realization of the predicted Posidonia cover was very close to zero. Examination of the coefficients showed that every new megaton of phosphogypsum discharged into the Gulf of Gabes reduces Posidonia meadow cover by 1 hectare on average (Fig. 3b). Further, in the period from 2000 to 2016, the hypothesis of the effect of Posidonia meadows shrinkage is not supported for the coastal fishery but it is strongly supported for the small pelagics fishery (Fig. 3c &d). In the latter, the loss of 1 hectare of Posidonia meadows is related on average the loss of 1 tonne of fisheries yield.

Discussion

Since the early 70s, the central part of the Gulf of Gabes has been severely impacted by high anthropogenic pressures of industrial origin [START_REF] Mosbahi | Response of benthic macrofauna to multiple anthropogenic pressures in the shallow coastal zone south of Sfax (Tunisia, central Mediterranean Sea)[END_REF][START_REF] El Kateb | Foraminiferal-based biotic indices to assess the ecological quality status of the Gulf of Gabes (Tunisia): Present limitations and future perspectives[END_REF][START_REF] Fourati | Sources and spatial distribution of dissolved aliphatic and polycyclic aromatic hydrocarbons in surface coastal waters of the Gulf of Gabès (Tunisia, Southern Mediterranean Sea)[END_REF], represented mainly by the huge untreated phosphogypsum quantities released directly in the coastal environment since 1972 [START_REF] El Zrelli | Heavy metal contamination and ecological risk assessment in the surface sediments of the coastal area surrounding the industrial complex of Gabes city Gulf of Gabes, SE Tunisia[END_REF], 2018a). We have shown in this work that over the same period seagrass cover loss occurred linearly and connected to the cumulative effect of phosphogypsum discharges. The most likely mechanism of decline was a widespread process of siltation of the sea bottom. Further, Cyanophyceae proliferation increased markedly because of nutrients enrichment discharged by the phosphate processing industries directly into sea water (Zaouali, 1993). The benthic ecosystem became dominated by invertebrates of muddy substrata enriched with organic matter (Zaouali, 1993).

This transformation of the benthic/demersal environment can be considered as a shift in the biocenosis of the central Gulf of Gabes which occurred between [START_REF] Berman | Optical transparency, chlorophyll and primary productivity in the Eastern Mediterranean near the Israeli coast[END_REF][START_REF] Estrada | Life and productivity of the open sea[END_REF]and 1990(Zaouali, 1993;;[START_REF] Guillaumont | Pollution impact study in Gabes Gulf (Tunisia) using remote sensing data[END_REF].

The results showed that Posidonia meadows have an important contribution to the local economy with a per hectare value of €511 and a total loss of value of ~ €60 million per year (Table 1). As previously mentioned, very few papers addressed this value in the scientific literature especially for the endemic Mediterranean Posidonia. For example, Campagne et al.

(2015) valued the Posidonia seagrass at €513 per hectare and a total of €55 million. But these values have been evaluated for several services only one of which was the contribution to the coastal fisheries. In addition, it is valid only in the socioeconomic French context [START_REF] Campagne | The seagrass Posidonia oceanica: Ecosystem services identification and economic evaluation of goods and benefits[END_REF]. Nevertheless, an attempt to determine the economic value of Cymodocea nodosa seagrass meadows for local fisheries at the Island of Gran Canaria using standard market prices has been carried out giving a monetary value at €866 per ha [START_REF] Tuya | Economic assessment of ecosystem services: Monetary value of seagrass meadows for coastal fisheries[END_REF] which is comparable to our valuation.

The loss estimate is an indicator that demonstrates the potential of economic gains in developing the Posidonia meadows. Indeed, the degradation of Posidonia meadows seems constant at least since 1990. If we assume that this degradation is linear and so is the increase in loss, a cumulated loss from 1990 to 2014 would give a total loss of more than €750 million.

Such a monetary loss should be considered by stakeholders that depend essentially on the fishing activity.

Comparing the economic losses due to the degradation of P. oceanica from the Gulf of Gabes to the added values of Gabes Industrial Complex (represented mainly by the GCT (Groupe Chimique Tunisien) factories), during the 1990-2014 period, we found that Posidonia meadows of Gabes Gulf are economically more important (Table 2). In fact, the economic losses related to the shrinkage of P. oceanica represented, in 2014, 39% and 64% of the added value of all Tunisian GCT factories and those of Gabes only, respectively. The cumulative economic losses, over the 1990-2014 period (15 years), reached ~ €780 million, so that the quarter of added values of Tunisian GCT factories and ~ 43% of the added values of Gabes GCT factories (Table 2). According to the study conducted by the European Community in 2017 on the economic losses related to the industrial pollution in Gabes (EC, 2017), the estimated annual cost of environmental degradation in Gabes Gulf was found to reach a maximum of €34 million (≈76 million Tunisian Dinar), which represents nearly 27%

of the added values of Tunisian GCT factories. However, it is worth noting that this latter economic valuation did not take into consideration the impact of phosphogypsum pollution on tourism and fisheries sector in the Gulf of Gabes in isolation, but also the impact of air pollution on other agricultural sectors, tourism, and human health. The fisheries sector per se represents ~47 % (~ €16 million) of this latter estimated amount (EC, 2017). Since the economic loss estimated during this study was found to be close to €59 million in 2014 (Table 2), it appeared that the EC ( 2017) study largely underestimated the economic losses in Gabes area (€16 million versus €59 million), although it took into account a wider range of economic sectors.

Ecological aspects of the Posidonia meadows shrinkage concern the life cycle of benthic/demersal coastal species, the role of Posidonia as primary producer and the morphological particularities of the central Gulf of Gabes. The strictly linear relationship applied between landing volumes of coastal fisheries and Posidonia meadows area is the most critical approximation used to carry out the evaluation. An estimation of the degree of residence of benthic commercial species has been proposed to improve the economic evaluation accuracy of Posidonia ecosystem services [START_REF] Jackson | Use of a seagrass residency index to apportion commercial fishery landing values and recreation fisheries expenditure to seagrass habitat service[END_REF]. According to the available literature, the residence degree can be evaluated for the principal commercial landed species in the Gabes governorate (FAO, 2016). More than ~ 90% of the benthic and demersal commercial species landed in the Gabes governorate basically depend on the presence of Posidonia meadows (FAO ArtFiMed, 2011b;FAO, 2016). Consequently, the linear approximation retained to qualify the link between landed volumes and the seagrass beds is not far off with regards of the economic evaluation.

Phytoplankton constitutes the first step of the marine food chain sustaining all the aquatic life (e.g. [START_REF] Field | Primary production of the biosphere: integrating terrestrial and oceanic components[END_REF] and especially oceanic life (Duarte and Cebrian, 1996). It feeds on mineral elements, mostly nitrogen and phosphorus. They are used to produce organic material via the photosynthesis process. The organic matter is reusable by zooplankton and then consumed by a wide variety of filter-feeding organisms. Nevertheless, Mediterranean waters are meso-oligotrophic to oligotrophic i.e. more or less depleted of mineral nutrients [START_REF] Estrada | Life and productivity of the open sea[END_REF][START_REF] Herut | Response of East Mediterranean surface water to Saharan dust: On-board microcosm experiment and field observations[END_REF]. In this situation of overall nutrient deficiency, the eastern Mediterranean basin is the most oligotrophic region [START_REF] Berman | Optical transparency, chlorophyll and primary productivity in the Eastern Mediterranean near the Israeli coast[END_REF] and the Gulf of Gabes is known as the most oligotrophic region of the eastern basin [START_REF] Drira | Dynamics of dinoflagellates and environmental factors during the summer in the Gulf of Gabès (Tunisia, Eastern Mediterranean Sea)[END_REF] because of weakness of mineral stock inherited from arid hinterlands and limited intrusions of Atlantic waters into the Gulf [START_REF] Bel Hassen | Summer phytoplankton pigments and community composition related to water mass properties in the Gulf of Gabes[END_REF]. Therefore, coastal waters of the Gulf offer very low levels of nutrients and little to sustain organisms which should result in low levels of biological productivity. Nevertheless, the Gulf of Gabes has the peculiarity of being highly productive. This contrast is a matter of debate. Some authors emphasise the role played by certain types of phytoplankton and zooplankton to explain the productivity of the Gabes Gulf [START_REF] Hamdi | The characterisation and summer distribution of ultraphytoplankton in the Gulf of Gabès (Eastern Mediterranean Sea, Tunisia) by using flow cytometry[END_REF][START_REF] Drira | Spatial and temporal variations of microphytoplankton composition related to hydrographic conditions in the Gulf of Gabès[END_REF]Béjaoui et al., 2020).

An alternative view could be formulated in terms of the food web system. The combination of the well-established role of primary producer of Posidonia coupled with the great extent of the continental shelf covered by these seagrass meadows may be a pertinent approach to the paradoxical situation -ultra-oligotrophy vs. high productivity -of the Gulf of Gabes. The hypothesis herein proposed is to view the shallow-water (<50m) of the central part of Gabes Gulf food web as benthic-dominated rather than plankton-dominated.

Firstly, Posidonia seagrass are nursery habitats because the leafy underwater canopy they create provides shelter for small invertebrates, small fish and juveniles of larger fish species [START_REF] Jehly | Phytoplankton composition off the coast of Point Revellata, Corsica, France[END_REF]. Other invertebrates grow nestled between the blades or in the sediments. The accumulation of smaller organisms amongst and on the Posidonia blades, as well as the Posidonia itself, attracts bigger animals. Thus, Posidonia can be home to many types of fish and invertebrates.

Secondly, Posidonia meadows contribute to global nitrogen cycling [START_REF] Costanza | The value of the world's ecosystem services and natural capital[END_REF][START_REF] Duarte | Seagrass community metabolism: Assessing the carbon sink capacity of seagrass meadows[END_REF][START_REF] Unsworth | Recognising the necessity for Indo-Pacific seagrass conservation[END_REF]. The role of Posidonia is essential because it retains mineral nutrients by reducing the water flow under the canopy of leaves, thus enabling nutrients to be trapped [START_REF] Jehly | Phytoplankton composition off the coast of Point Revellata, Corsica, France[END_REF]. Consequently, the Posidonia biocenosis favoured expansion of phytoplankton populations as a permanent source of nutrients like nitrogen, phosphorous, and organic carbon used by them [START_REF] Jehly | Phytoplankton composition off the coast of Point Revellata, Corsica, France[END_REF].

Nevertheless, the primary productivity of Posidonia is not the major direct source of carbon and nitrogen for higher trophic levels of the food web and most Posidonia meadow invertebrate herbivores utilize epiphytic species [START_REF] Unsworth | Recognising the necessity for Indo-Pacific seagrass conservation[END_REF]. Posidonia is an adequate support for epiphytic species [START_REF] Boudouresque | The outstanding traits of the functioning of the Posidonia oceanica seagrass ecosystem[END_REF][START_REF] Duarte | Seagrass ecosystems: their global status and prospects[END_REF]) and many species of algae and microalgae, bacteria and invertebrates that grow as "epiphytes" directly on living Posidonia leaves. Small invertebrates, such as crustaceans and snails, feed on epiphytes. They are in turn consumed by larger crustaceans and fish and are important links in the coastal food web [START_REF] Duarte | Seagrass ecosystems: their global status and prospects[END_REF][START_REF] Vizzini | Analysis of the trophic role of Mediterranean seagrasses in marine coastal ecosystems: A review[END_REF][START_REF] Vizzini | The role and contribution of the seagrass Posidonia oceanica (L.) Delile organic matter for secondary consumers as revealed by carbon and nitrogen stable isotope analysis[END_REF]. Further, some epiphytic bacteria can extract nitrogen from the environment and make it available to larger animals [START_REF] Jehly | Phytoplankton composition off the coast of Point Revellata, Corsica, France[END_REF]. Epiphytes are commonly associated to Posidonia leaves in the Gulf of Gabes [START_REF] Brahim | Variability in the structure of epiphytic assemblages of Posidonia oceanica in relation to human interferences in the Gulf of Gabes, Tunisia[END_REF].

Lastly, dead Posidonia leaves also play an important role in coastal ecosystems. When the leaves die, they decay on the sediment or are washed onto the beach, supporting a diverse community of decomposers that prosper on rotting material.

Taking into account the situation of small pelagic fisheries, it is worth noting that the efficiency of the purse seine technique coupled with the benthic to pelagic biocenosis shift hypothesis has strongly favoured the development of small pelagic fishing activity.

Nevertheless, its decline at the end of the eighties (Fig. 2b &c) most probably corresponded to the effects of pollution load of the marine environments which reached a larger area in the central part of Gabes Gulf [START_REF] Rabaoui | On the relationship between the diversity and structure of benthic macroinvertebrate communities and sediment enrichment with heavy metals in Gabes gulf Tunisia[END_REF][START_REF] El Zrelli | Biomonitoring of coastal pollution in the Gulf of Gabes (SE, Tunisia): use of Posidonia oceanica seagrass as a bioindicator and its mat as an archive of coastal metallic contamination[END_REF], 2019b & c). The high proportion of toxic phytoplankton indicates that quality of marine waters is still poor [START_REF] Feki | What are the potential drivers of blooms of the toxic dinoflagellate Karenia selliformis? A 10-year study in the Gulf of Gabes, Tunisia, southwestern Mediterranean Sea[END_REF][START_REF] Abdennadher | Factors determining the dynamics of toxic blooms of Alexandrium minutum during a 10-year study along the shallow southwestern Mediterranean coasts[END_REF]. Nowadays, part of small pelagic catch by coastal fisheries techniques (as beach seine for example) has increased, especially the sardinella volume (Fig. 3a &b). This could be the consequence of the high nutrients amounts in coastal sea water. It is probably the reason why the coastal fisheries volumes are irregular from one year to the next and can be, from times to times, relatively high (Fig. 2b &c).

A first step towards a better understanding of the vulnerability of the Gabes Gulf could consist in introducing an economic measure of this vulnerability in terms of producer surplus.

The producer surplus loss measures the difference between the amounts a producer of a good receives, the market price, and the minimum amount the producer is willing to accept for the good. This minimum amount is called the marginal cost. The producer surplus is represented in Fig. 4, where we assume that the demand for fish of the Gulf of Gabes is perfectly elastic.

The degradation of the Posidonia seagrass increases the marginal cost for the fishermen. This is translated into a move of the supply curve of fish from S0 to S1 or S2 (Fig. 4). The producer surplus loss is estimated through the shaded area in Fig. 4: the more the marginal costs increase, the more the producer surplus losses increase.

Our indicator value and loss can be integrated in this figure. The value is represented by the rectangle pOQ1S1, when considering a loss of Posidonia from Q0 to Q1, and the rectangle pOQ2S2, when considering a loss of Posidonia from Q0 to Q2. The loss is represented by the rectangle S1Q1Q0S0, when considering a loss of Posidonia from Q0 to Q1, and the rectangle S2Q2Q0S0, when considering a loss of Posidonia from Q0 to Q2. The producer surplus can be considered as much accurate to value the Posidonia seagrass when the loss is important.

Indeed, if we consider to change from S0 to S2 and use our value indicator, we can observe graphically in Fig. 4 that the value will be very low compared to the real producer surplus loss. However, the producer surplus method is not easy to apply in a large territory because of the need of data on marginal cost or production from fishermen which is difficult to obtain.

Furthermore, we focus our analysis on the marketed direct use value of the seagrass.

Thus, we studied the supply side and the impact of Posidonia seagrass on the fishermen benefits. However, the economic value of Posidonia seagrass is not limited to this as seen in [START_REF] Campagne | The seagrass Posidonia oceanica: Ecosystem services identification and economic evaluation of goods and benefits[END_REF]. Indeed, it provides benefits through the regulation and maintenance services, water purification value, water oxygenation and cultural services. In order to complete our economic valuation of Posidonia seagrass in the Gulf of Gabes and give more accurate information to decision makers, a future line of research would be to realize these valuations, including interviewing various stakeholders, other than fishermen, living around the Gulf of Gabes and valuate the benefits of the cultural service. factories at the national (Tunisia) and regional (Gabes) scales, during the 1990-2014 period. (1) estimates made based on the economic valuation carried out during the present study; (2) 1965(De Gaillande, 1970b), c. 1981[START_REF] Darmoul | Pollution dans le Golfe de Gabès (Tunisie) : bilan de six années de surveillance[END_REF]), d. 1990(Zaouali, 1993), e. 2014[START_REF] El Zrelli | Biomonitoring of coastal pollution in the Gulf of Gabes (SE, Tunisia): use of Posidonia oceanica seagrass as a bioindicator and its mat as an archive of coastal metallic contamination[END_REF]).

Conclusion

A

  gradual intensification of sponge, then shrimp fishing activities during the first half of the 20 th century has triggered the process of Posidonia meadows area reduction in the centralGulf of Gabes. Subsequently, between 1972 and[START_REF] Estrada | Life and productivity of the open sea[END_REF], successive commissioning of phosphate fertilizer plant units in Gabes has induced a high environmental stress via phosphogypsum releases into coastal marine environments. Furthermore, the gradual disappearance of the seagrass beds caused the gradual silting up of the central part of the Gulf which progressed from centre to coast. This silting process contributed to the waters transparency decrease and thus to gradual disappearance of the seagrass beds. From an economic point of view, the loss of economic value provided by coastal fisheries income correlated to the decline of the Posidonia meadows has been valued at about €60 million in 2014 and a cumulated loss of more than €750 million from 1990 to 2014.

Figures CaptionsFig. 1 .

 1 Figures CaptionsFig. 1. Location of the Gabes Gulf (a) and the historical decline of P. oceanica (right panel: b.

Fig. 2 .

 2 Fig. 2. Historical trends of observed (asterisks) and predicted Posidonia cover in the Gulf of Gabes, from a linear model of annual decay plus two standard errors (a), and historical variations of "coastal" (b) and "small pelagic" (c) landing volumes and fishing efforts as well as the Posidonia area, during 1970-2015.

Fig. 3 .

 3 Fig. 3. Hypotheses testing using the Akaike Information Criterion (AIC) for the effect of fishing effort and predicted Posidonia cover on the (a & b) fishing yield by two fleets in the Gulf of Gabes and (c & d) coastal and pelagic fishing fleets.

Fig. 4 .

 4 Fig. 4. Estimate of the producer surplus loss after the Posidonia oceanica seagrass service degradation.

Table 1 .

 1 Economic value of Posidonia oceanica meadows.

	Year	1990	2014
	Market price (€)	1 965	2 960
	Quantity (tons)	1 062	2 560
	Posidonia oceanica meadows (ha) and {in percent of total land area} 43 563 {33.51} 14 826 {11.41}
	Total value of Posidonia oceanica (€)	2 086 640	7 576 473
	Yield of ton of benthic species per hectare of Posidonia oceanica	0.024	0.173
	Per hectare value of Posidonia (€)	48	511
	Economic losses due to the degradation of Posidonia oceanica (€)	4 140 278	58 854 676

Table 2 .

 2 Economic losses linked to the Posidonia oceanica degradation versus the added values of the GCT (Groupe Chimique Tunisien)

  data collected from the European Community study conducted in 2017(EC, 2017).

	Years / Period		1990	2014	1990-2014
	Economic losses due to the Posidonia oceanica degradation (ELDPO, M€) (1)	4	59	780
	Added values of all GCT factories in Tunisia (AVTF) (2)	Added Values (M€) Ratio ELDPO/AVTF (%)	-24 117	152 39	3316 24
	Added values of Gabes GCT factories (AVGF) (2)	Added Values (M€) Ratio ELDPO/AVGF (%)	-12 132	92 64	1824 43

Acknowledgments

The authors would like to thank the editor and the reviewer who contributed to improve the quality of the manuscript through their constructive comments and suggestions. We are also grateful to all those who helped in data collection and analysis. This study is dedicated to the memory of Habib Chairet (09.11.1960Chairet (09.11. /29.10.2018) ) for his environmental protection efforts in the coastal oasis of Chatt Sidi Abd Essalam.