H S Aghamiry
email: h.aghamiry@ut.ac.ir

A Gholami
email: agholami@ut.ac.ir

S Operto
email: operto@geoazur.unice.fr

ADMM-based multi-parameter wavefield reconstruction inversion in VTI acoustic media with TV regularization

Keywords: Waveform inversion, Inverse theory, Seismic anisotropy, Controlled source seismology

Full waveform inversion (FWI) is a nonlinear waveform matching procedure, which suffers from cycle skipping when the initial model is not kinematically-accurate enough. To mitigate cycle skipping, wavefield reconstruction inversion (WRI) extends the inversion search space by computing wavefields with a relaxation of the wave equation in order to fit the data from the first iteration. Then, the subsurface parameters are updated by minimizing the source residuals the relaxation generated. Capitalizing on the wave-equation bilinearity, performing wavefield reconstruction and parameter estimation in alternating mode decomposes WRI into two linear subproblems, which can solved efficiently with the alternating-direction method of multiplier (ADMM), leading to the so-called iteratively refined wavefield reconstruction inversion (IR-WRI). Moreover, ADMM provides a suitable framework to implement bound constraints and different types of regularizations and their mixture in IR-WRI. Here, IR-WRI is extended to multiparameter reconstruction for VTI acoustic media. To achieve this goal, we first propose different forms of bilinear VTI acoustic wave equation. We develop more specifically IR-WRI for the one that relies on a parametrisation involving vertical wavespeed and Thomsen's parameters δ and . With a toy numerical example, we first show that the radiation patterns of the virtual sources generate similar wavenumber filtering and parameter cross-talks in classical FWI and IR-WRI. Bound constraints and TV regularization in IR-WRI fully remove these undesired effects for an idealized piecewise constant target. We show with a more realistic long-offset case study representative of the North Sea that anisotropic IR-WRI successfully reconstruct the vertical wavespeed starting from a laterally homogeneous model and update the long-wavelengths of the starting model, while a smooth δ model is used as a passive background model. VTI acoustic IR-WRI can be alternatively performed with subsurface parametrisations involving stiffness or compliance coefficients or normal moveout velocities and η parameter (or horizontal velocity).

INTRODUCTION

Full waveform inversion (FWI) is a waveform matching procedure which provides subsurface model with a wavelength resolution. However, it suffers from cycle skipping when the initial model is not accurate enough according to the lowest available frequency. To mitigate cycle skipping, the search space of frequency-domain FWI can be extended by wavefield reconstruction inversion (WRI) [START_REF] Van Leeuwen | Mitigating local minima in fullwaveform inversion by expanding the search space[END_REF], 2016). In WRI, the search space is extended with a penalty method to relax the wave-equation constraint at the benefit of the observation equation (i.e., the data fit) during wavefield reconstruction. Then, the subsurface parameters are estimated from the reconstructed wavefields by minimizing the source residuals the relaxation generated. If these two subproblems (wavefield reconstruction and parameter estimation) are solved in alternating mode [START_REF] Van Leeuwen | Mitigating local minima in fullwaveform inversion by expanding the search space[END_REF] rather than by variable projection [START_REF] Van Leeuwen | A penalty method for PDEconstrained optimization in inverse problems[END_REF], WRI can be recast as a sequence of two linear subproblems capitalizing on the bilinearity of the scalar Helmholtz equation with respect to the wavefield and the squared slownesses (the Helmholtz equation is linear with respect to the wavefield for a given squared slowness model and is linear with respect to the squared slownesses for a given wavefield). Aghamiry et al. (2019c) improved WRI by replacing the penalty method with the augmented Lagrangian method [START_REF] Nocedal | Numerical Optimization[END_REF]) and solve it using the alternating-direction method of multipliers (ADMM) [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction of multipliers[END_REF], leading to iterativelyrefined WRI (IR-WRI). Although ADMM was originally developed for separable convex problems, the bilinearity of the wave equation makes IR-WRI biconvex, which allows for the use of ADMM as is (Boyd et al. 2010, Section 3.1.1). Moreover, a scaled form of ADMM draws clear connection between WRI and IR-WRI in the sense that it shows that IR-WRI reduces to a penalty method in which the right-hand sides (RHSs) in the quadratic objective functions associated with the observation equation and the wave equation are iteratively updated with the running sum of the data and source residuals in iteration (namely, the scaled Lagrange multipliers). This RHS updating, the lacking feature in WRI, efficiently refines the solution of the two linear subproblems when a fixed penalty parameter is used (Aghamiry et al. 2018b(Aghamiry et al. , 2019c)). Later, Aghamiry et al. (2018cAghamiry et al. ( , 2019b) ) interfaced bound constraints and edge-preserving regularizations with ADMM to manage large-contrast media. Also, to preserve the smooth components of the subsurface when edge preserving regularizations are used, Aghamiry et al. (2018aAghamiry et al. ( , 2019a) ) combine blocky and smooth promoting regularization in the framework of IR-WRI by using infimal convolution of Tikhonov and Total Variation (TV) regularization functions.

IR-WRI was mainly assessed for wavespeed estimation from the scalar Helmholtz equation. The objective of this study is to develop and assess the extension of IR-WRI to multi-parameter reconstruction in VTI acoustic media. To achieve this goal, we first need to review different formulations of the VTI acoustic wave equation and different subsurface parametrisations for which bilinearity of the wave equation is fulfilled, in order to keep the parameter estimation subproblem linear. Since the wavefield reconstruction requires to solve a large and sparse system of linear equations, second-order or fourth-order wave equation for pressure will be favored at the expense of first-order velocity-stress formulations to mitigate the number of unknowns during wavefield reconstruction. However, we stress that different forms of the wave equation can be used to perform wavefield reconstruction and parameter estimation, provided they provide consistent solutions (Gholami et al. 2013b, Their Appendix A and B). Also, we will favor subsurface parametrisation involving the vertical wavespeed v0 and the Thomsen's parameter and δ at the expense of that involving stiffness coefficients according to the tradeoff analysis of Gholami et al. (2013b, Their Appendix A and B).

When bilinearity of the wave equation is fulfilled, IR-WRI can be extended to multi-parameter estimation following the procedure promoted by Aghamiry et al. (2019b), where TV regularization and bound constraints are efficiently implemented in the parameter-estimation subproblem using variable splitting schemes [START_REF] Glowinski | Splitting Methods in Communication[END_REF]. The splitting procedure allows us to break down the non-differentiable TV regularization problem into two easy-to-solve subproblems: a least-squares quadratic subproblem and a proximity subproblem [START_REF] Goldstein | The split Bregman method for L1regularized problems[END_REF].

In this study, we perform a first assessment of multi-parameter IR-WRI for v0 and using a toy inclusion example and the more realistic synthetic North Sea case study tackled by Gholami et al. (2013b) and Gholami et al. (2013a). With the toy example, we first show that the radiation patterns and the parameter cross-talks have the same effects as in classical FWI when TV regularization is not applied. Then, we show how TV regularization fully removes bandpass filtering and cross-talk effects generated by these radiation patterns for this idealized piecewise constant target. With the North Sea example, we first show the resilience of IR-WRI against cycle skipping when using a crude initial v0 model. The reconstruction of the v0 model is accurate except in the deep smooth part of the subsurface which suffers from a deficit of wide-angle illumination, while the reconstruction of is more challenging and requires additional regularization to keep its update smooth and close to the starting model. However, we manage to update significantly the long-to-intermediate wavelengths of by IR-WRI, unlike Gholami et al. (2013a). Also, comparison between mono-parameter IR-WRI for v0 and multi-parameter IR-WRI for v0 and allows us to gain qualitative insights on the sensitivity of the IR-WRI to in terms of subsurface model quality and data fit.

This paper is organized as follow. We first discuss the bilinearity of the acoustic VTI wave equation as well as its implication on the gradient and the Gauss-Newton Hessian of the parameterestimation subproblem. From the selected bilinear formulation of the wave equation, we develop bound-constrained TV-regularized IR-WRI for VTI acoustic media parametrized by 1/v 2 0 , 1 + 2 and √ 1 + 2δ. Third, we assess IR-WRI for VTI acoustic media against the inclusion and North Sea case studies. Finally, we discuss the perspectives of this work.

THEORY

In this section, we first show that the VTI acoustic wave equations is bilinear with respect to the wavefield and model parameters. Then, we rely on this bilinearity to formulate multi-parameter IR-WRI in VTI acoustic media with bounding constraints and TV regularization.

Bilinearity of the wave equation: preliminaries

If we write the wave equation in a generic matrix form as

A(m)u = s, (2.1)
where A ∈ C (n×nc)×(n×nc) is the wave-equation operator, u ∈ C (n×nc)×1 is the wavefield vector, s ∈ C (n×nc)×1 is the source vector, m ∈ R (n×nm)×1 is the subsurface parameter vector, n is the number of degrees of freedom in the spatial computational mesh, nc is the number of wavefield components, nm is the number of parameter classes, then the wave equation is bilinear if there exists a linear operator

L(u) ∈ C (n×nm)×(n×nm) such that L(u)m = y(u), (2.2)
where y(u) ∈ C (n×nm)×1 . Bilinearity is verified when the lefthand side of the wave equation can be decomposed as:

A(m)u = B M(m) Cu + Du, (2.3)
where M(m) is a block matrix, whose blocks of dimension n × n are either 0 or of the form diag(mi), and matrices B, C and D don't depend on m. The operator diag(•) denotes a diagonal matrix of coefficients • and mi is the subsurface parameter vector of class i. By noting that diag(x)z = diag(z)x, the block diagonal structure of M allows one to rewrite the term B M(m) Cu as B C (u)m, where C (u) is a block matrix, whose blocks of dimension n × n are either 0 or diagonal with coefficients depending on u. It follows from equation 2.3 and the above permutation between u and m that the wave equation can be re-written as

A(m)u -s = B C (u)m + Du -s = L(u)m -y(u). (2.4)
Moreover, in the framework of multi-parameter analysis, it is worth noting that (2.5) where the left-hand side is the so-called virtual source associated with m k [START_REF] Pratt | Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[END_REF]) and e k ∈ C (n×nm)×1 denotes a column vector whose kth component is one while all the others are zeros. Accordingly, the normal operator L T L, i.e., the Gauss-Newton Hessian of the parameter estimation subproblem in IR-WRI, is formed by the zero-lag correlation of the virtual sources and, hence is extremely sparse. We also point that L and its associated normal operator L T L are block diagonals if B is diagonal. This means that, if the model parameters are first sorted according to their position in the mesh (fast index) and second according to the parameter class they belong to (slow index), then the diagonal coefficients of the off-diagonal blocks describe the inter-parameter coupling. In the following section, we show the bilinearity of the acoustic VTI wave equation based upon the above matrix manipulations.

∂A(m) ∂m k u = L(u)e k ,

Bilinearity of the acoustic VTI wave equation

First-order velocity-stress wave equation

We first consider the frequency-domain first-order velocity-stress wave equation in 2D VTI acoustic media [START_REF] Duveneck | Acoustic VTI wave equations and their application for anisotropic reverse-time migration[END_REF][START_REF] Duveneck | Stable P-wave modeling for reversetime migration in tilted TI media[END_REF][START_REF] Operto | Computationally-efficient three-dimensional visco-acoustic finite-difference frequency-domain seismic modeling in vertical transversely isotropic media with sparse direct solver[END_REF])

-îωv x,l = diag(b)∇xu x,l , -îωv z,l = diag(b)∇zu z,l , -îωu x,l = diag(c11)∇xv x,l + diag(c13)∇zv z,l -îωs l , -îωu z,l = diag(c13)∇xv x,l + diag(c33)∇zv z,l -îωs l , (2.6)
where î = √ -1, ω is the angular frequency, v x,l ∈ C n×1 and v z,l ∈ C n×1 are the horizontal and vertical particle velocity wavefields, s l ∈ C n×1 denote the pressure sources, and u x,l ∈ C n×1 and u z,l ∈ C n×1 are the so-called horizontal and vertical pressure wavefields [START_REF] Plessix | A parametrization study for surface seismic full waveform inversion in an acoustic vertical transversely isotropic medium[END_REF]. The subscript l ∈ {1, 2, ..., ns} is the source index, where ns denotes the number of sources. The subsurface properties are parametrized by the buoyancy b ∈ R n×1 (inverse of density) and the stiffness coefficients cij ∈ R n×1 . Operators ∇x and ∇z are finite difference approximation of first order derivative operators with absorbing perfectly matched layer (PML) coefficients [START_REF] Bérenger | A perfectly matched layer for absorption of electromagnetic waves[END_REF].

Gathering equation 2.6 for all sources leads to the following matrix equation:

M b 0 0 Mc 0 ∇ ∇ 0 + îωI V U = îω 0 S ,
(2.7) where I is the identity matrix,

∇ = ∇x 0 0 ∇z , U = U h Uv , V = Vx Vz , S = Sx Sz , (2.8) with Vx = vx,1 vx,2 ... vx,n s ∈ C n×ns ,
and analogously for Vz, U h , Uv, and

Sx = Sz = s1 s2 ... sn s ∈ C n×ns .
Furthermore, When U and V are known, this system can be also recast as a new linear system in which the unknowns are the model parameters, thus highlighting the bilinearity of the wave equation. For the lth source, the new equations become

M b = diag(b) 0 0 diag(b) , Mc = diag(c11) diag(c13) diag(c13) diag(c33) .
    L11 0 0 0 L21 0 0 0 0 L32 L33 0 0 0 L43 L44         b c11 c13 c33     = îω     0 -v x,l 0 -v z,l s l -u x,l s l -u z,l     , (2.9)
where

         L11 = diag(∇xu x,l ), L21 = diag(∇zu z,l ), L32 = L43 = diag(∇xv x,l ), L33 = L44 = diag(∇zv z,l ).
Equations 2.7 and 2.9 are equivalent forms of the original equation 2.6. The former expresses the discretized wavefields as the unknowns of a linear system, whose coefficients depend on the known subsurface parameters, while the latter expresses the model parameters as the unknowns of an another linear system, whose coefficients depend on the known wavefields. In the framework of WRI, this bilinearity allows one to recast the waveform inversion problem as two linear subproblems for wavefield reconstruction and parameter estimation, which can be solved efficiently in alternating mode with ADMM (Aghamiry et al. 2019c,b). In the next section, we show the bilinearity of the second-order frequency-domain wave equation, which may be more convenient to solve with linear algebra methods than the first-order counterpart, since it involves fewer unknowns for a computational domain of given size.

Second-order wave equation

Following a parsimonious approach (e.g., [START_REF] Operto | 3D finite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver: A feasibility study[END_REF]), we eliminate v x,l and v z,l from equation 2.6 to derive a system of two second-order partial differential equations as

-ω 2 u x,l = diag(c11)∇xdiag(b)∇xu x,l + diag(c13)∇zdiag(b)∇zu z,l -ω 2 s l , -ω 2 u z,l = diag(c13)∇xdiag(b)∇xu x,l + diag(c33)∇zdiag(b)∇zu z,l -ω 2 s l .
(2.10) Equation 2.10 defines a tri-linear equation with respect to buoyancy, stiffness parameters and pressure wavefields. A tri-linear function is a function of three variables which is linear in one variable when the other two variables are fixed. In this study, we will assume that density is constant and equal to 1 to focus on the estimation of the anisotropic parameters. If heterogeneous density needs to be considered, the second-order wave equation can be recast as a bilinear system if the first-order wave equation is parametrized with compliance coefficients instead of stiffness coefficients (see Appendix A and [START_REF] Vigh | Elastic full-waveform inversion application using multicomponent measurements of seismic data collection[END_REF][START_REF] Yang | A review on the systematic formulation of 3D multiparameter full waveform inversion in viscoelastic medium[END_REF]. Alternatively, the parameter estimation can be performed with the bilinear firstorder wave equation, equation 2.9, while wavefield reconstruction is performed with the second-order wave equation, equation 2.10, or the fourth-order wave equation reviewed in Appendix C for sake of computational efficiency. We continue by assuming that the density is constant and equal to 1 and parametrize the VTI equation in terms of vertical wavespeed v0 and Thomsen's parameters and δ [START_REF] Thomsen | Weak elastic anisotropy[END_REF]). Accordingly, we rewrite equation 2.10 as (2.11) where ∇xx = ∇x∇x and ∇zz = ∇z∇z. We write this 2n × 2n linear system in a more compact form as

-ω 2 diag 1/v 2 0 u x,l = diag(1 + 2 )∇xxu x,l + diag( √ 1 + 2δ)∇zzu z,l -diag(1/v 2 0 )s l , -ω 2 diag 1/v 2 0 u z,l = diag( √ 1 + 2δ)∇xxu x,l + ∇zzu z,l -diag 1/v 2 0 s l ,
A(m)U = S, (2.12)
where S T = ω 2 S T x S T z , U is defined as in equation 2.8, and the matrix A is given by

ω 2 diag(mv 0 ) + diag(m )∇xx diag(m δ )∇zz diag(m δ )∇xx ω 2 diag(mv 0 ) + ∇zz ,
where the model parameters are

m =   mv 0 m m δ   =   1/v 2 0 1 + 2 √ 1 + 2δ   .
Equation 2.12 is linear in wavefields U when the model parameters m are known. When U is known, this system can be recast as a linear system in which the unknowns are the model parameters.

        L1 . . . L l . . . Ln s         m =         y1 . . . y l . . . yn s         , (2.13)
where L l is given by

ω 2 diag(u x,l ) diag(∇xxu x,l ) diag(∇zzu z,l ) ω 2 diag(u z,l ) 0 diag(∇xxu x,l )
and

y l = ω 2 s l ω 2 s l -∇zzu z,l .
Note that each block of L l is diagonal. In the next section, we develop multi-parameter acoustic VTI IR-WRI with bound constraints and TV regularization. We give the most general formulation in which all the tree parameter classes mv 0 , m , m δ are optimization parameters (updated by the inversion). However, one may process some of them as passive parameters or update the parameter classes in sequence rather than jointly. In this case, the linear system associated with the parameter estimation subproblem, equation 2.13, will change. Table 1 presents this system for the different possible configurations.

ADMM-based acoustic VTI wavefield reconstruction inversion

We consider the following bound-constrained TV-regularized multivariate optimization problem associated with the wave equation described by equation 2.12: (2.14) where ∂x and ∂z are, respectively, first-order finite-difference operators in the horizontal and vertical directions with appropriate boundary conditions, C = {x ∈ R 3n×1 | m l ≤ x ≤ mu} is the set of all feasible models bounded by a priori lower bound m l and upper bound mu, D = d1 d2 ... dn s ∈ C nr ×ns with d l denoting the recorded data for the lth source, each including nr samples (the number of receivers), P ∈ R nr ×2n is a linear observation operator which samples the wavefields at the receiver positions. Here, we assume that the sampling operator is identical across all sources (stationary-recording acquisitions). However, one may used a specific operator for each source.

min U,m∈C |∂xm| 2 + |∂zm| 2 , subject to PU = D, A(m)U = S,
We solve this constrained optimization problem with ADMM [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction of multipliers[END_REF]Aghamiry et al. 2019b), an augmented Lagrangian method with operator splitting, leading to the following saddle point problem

min U,m∈C max D, S |∂xm| 2 + |∂zm| 2 + D, PU -D + λ0 PU -D 2 F + S, A(m)U -S + λ1 A(m)U -S 2 F , (2.15)
where • 2 F denotes the Frobenius norm of •, λ0, λ1 > 0 are penalty parameters and D ∈ C nr ×ns and S ∈ C n×ns are the Lagrange multipliers (dual variables).

The last two lines of Equation 2.15 shows that the augmented Lagrangian function combines a Lagrangian function (left terms) and a penalty function (right terms). Also, scaling the Lagrange multipliers by the penalty parameters, D = -D/λ0 and S = -S/λ1, allows us to recast the augmented Lagrangian function in a more convenient form (Boyd et al. 2010, Section 3

.1.1) min U,m∈C max D, S |∂xm| 2 + |∂zm| 2 + λ0 PU -D -D 2 F -D 2 F + λ1 A(m)U -S -S 2 F -S 2 F , (2.16)
where the scaled dual variables have been injected in the penalty functions.

In the WRI framework, the augmented Lagrangian method provides an efficient and easy-to-tune optimization scheme that extends the parameter search space by introducing a significant relaxation of the wave equation at the benefit of the observation equation during the early iterations, while satisfying the two equations at the convergence point. We solve the saddle point problem, equation 2.16, with the method of multiplier, in which the primal variables, U and m, and the dual variables, D and S, are updated in alternating mode. The dual problem is iteratively solved with basic gradient ascent steps. Accordingly, we immediately deduce from equation 2.16 that the scaled dual variables D and S are formed by the running sum of the constraint violations (the data and source residuals) in iterations. They update the RHSs (the data and the sources) of the quadratic penalty functions in equation 2.16 to refine the primal variables U and m at a given iteration from the residual source and data errors (this RHS updating is a well known procedure to iteratively refine solutions of ill-posed linear inverse problems). The Table 1. The linear systems corresponding to the update of mv 0 , m and m δ when they are active or passive during the inversion. In the first three columns, "a." denotes an active parameter and "p." denotes a passive parameter. mv 0 m m δ The corresponding linear system for updating the model a. p. p.

ω 2 diag(u x,l ) ω 2 diag(u z,l ) L l ∈C 2n×n mv 0 m∈R n×1 = ω 2 s l -diag(m )∇ xxux,l -diag(m δ )∇ zz u z,l ω 2 s l -∇ zz u z,l -diag(m δ )∇ xxux,l y l ∈C 2n×1 p. a. p. diag(∇ xxux,l) L l ∈C n×n m m∈R n×1 = ω 2 s l -ω 2 diag(mv 0 )u x,l -diag(m δ )∇ zz u z,l y l ∈C n×1 p. p. a. diag(∇ zz u z,l ) diag(∇ xxux,l) L l ∈C 2n×n m δ m∈R n×1 = ω 2 s l -ω 2 diag(mv 0 )u x,l -diag(m )∇ xxux,l ω 2 s l -∇ zz u z,l -ω 2 diag(mv 0 )u z,l y l ∈C 2n×1 p. a. a. diag(∇ xxux,l) diag(∇ zz u z,l ) 0 diag(∇ xxux,l) L l ∈C 2n×2n m m δ m∈R 2n×1 = ω 2 s l -ω 2 diag(mv 0 )u x,l ω 2 s l -∇ zz u z,l -ω 2 diag(mv 0 )u z,l y l ∈C 2n×1 a. a. p. ω 2 diag(u x,l ) diag(∇ xxux,l) ω 2 diag(u z,l ) 0 L l ∈C 2n×2n mv 0 m m∈R 2n×1 = ω 2 s l -diag(m δ )∇ zz u z,l ω 2 s l -∇ zz u z,l -diag(m δ )∇ xxux,l y l ∈C 2n×1 a. p. a. ω 2 diag(u x,l ) diag(∇ zz u z,l ) ω 2 diag(u z,l ) diag(∇ xxux,l) L l ∈C 2n×2n mv 0 m δ m∈R 2n×1 = ω 2 s l -diag(m )∇ xxux,l ω 2 s l -∇ zz u z,l y l ∈C 2n×1
bi-variate primal problem, equation 2.14, is biconvex due to the bilinearity of the wave equation highlighted in the previous section. Therefore, it can be broken down into two linear subproblems for U and m, which can be solved efficiently in alternating mode with ADMM after noting that the TV regularizer is convex (Aghamiry et al. 2019b).

As pointed out by Aghamiry et al. (2019c), a key advantage of augmented Lagrangian methods compared to penalty methods is that fixed penalty parameters λ0 and λ1 can be used in iterations, because the Lagrange multipliers progressively correct for the constraint violations generated by the penalty terms through the above mentioned RHS updating.

Starting from an initial model m and zero-valued dual variables, the kth ADMM iteration embeds the following steps (see Appendix B for the complete development):

2.5.1

Step 1: The primal wavefield reconstructions.

Build regularized wavefields by solving the following multi-RHS system of linear equations with direct or iterative methods suitable for sparse matrices:

λ0P T P + λ1A(m k ) T A(m k ) U = λ0P T [D + Dk ] + λ1A(m k ) T [S + Sk ] .
(2.17) By choosing a small value of λ0/λ1, the reconstructed wavefields closely fit the observations during the early iterations, while only weakly satisfying the wave equation. Problem 2.17 can be also interpreted as an extrapolation problem to reconstruct U, when the observation equation (i.e. PU = D) is augmented with the wave equation.

To mitigate the computational burden of the wavefield reconstruction, we solve equation 2.17 with a fourth-order wave equation operator following the parsimonious approach of [START_REF] Operto | Computationally-efficient three-dimensional visco-acoustic finite-difference frequency-domain seismic modeling in vertical transversely isotropic media with sparse direct solver[END_REF], while the subsequent model estimation subproblem relies on the bilinear wave equation provided in equation 2.11. The elimination procedure allowing to transform the system of two secondorder wave equations for ux and uz, equation 2.11, into a fourthorder wave equation for ux coupled with the closed-form expression of uz as a function of ux is reviewed in Appendix C.

2.5.2

Step 2: The primal model estimation.

The reconstructed wavefields, equation 2.17, are injected in the linearized equation 2.13 to update the subsurface parameters by solving the following system of linear equations:

λ1 ns l=1 L T l L l + ∇ T Γ∇ + ZI m = λ1 ns l=1 L T l (y k l + sk l ) + ∇ T Γ(p k + pk ) + Z(q k + qk ), (2.18)
where

∇ = ∂x ∂z . (2.19)
In equation 2.18, p and q are auxiliary primal variables, which have been introduced to solve the bound-constrained TV-regularized parameter estimation subproblem with the split Bregman method (Appendix B). The vectors p and q are the corresponding dual variables. These auxiliary primal and dual variables are initialized to 0 during the first ADMM iteration.

The operator Z is a diagonal weighting matrix defined as

Z =   ζv 0 I 0 0 0 ζ I 0 0 0 ζ δ I   ∈ R 3n×3n + , (2.20)
where ζv 0 , ζ , ζ δ > 0 control the relative weights assigned to the bound constraints applied on the three parameter classes. Note that the bound constraints introduce also a damping (DMP) or zeroorder Tikhonov regularization in the Hessian of equation 2.18.

In the same way, Γ is a diagonal matrix defined as

Γ = Γ11 0 0 Γ22 , (2.21)
where

Γ11 = Γ22 =   γv 0 I 0 0 0 γ I 0 0 0 γ δ I   ∈ R 3n×3n + ,
and γv 0 , γ , γ δ > 0 control the soft thresholding that is performed by the TV regularizer, equation 2.23. We remind that augmented Lagrangian methods seek to strictly satisfy the constraints at the convergence point only, not at each iteration. Therefore, the relative values of these penalty parameters have a significant impact upon the path followed by the inversion to converge toward this convergence point.

2.5.3

Step 3: The TV primal update.

Update the TV primal variable p via a TV proximity operator. Set

z ← ∇m k+1 -pk = zx zz , then p k+1 ← prox Γ -1 (z) = ξ • zx ξ • zz , (2.22)
where

ξ = max(1 - 1 Γ √ z 2 x + z 2 z , 0) (2.23)
and x • y denotes the Hadamard (component-wise) product of x and y. Also, the power of 2 indicates the Hadamard product of x with itself, i.e. x 2 = x • x.

2.5.4

Step 4: The bounding constraint primal update.

Update the primal variable q via a projection operator, which has the following component-wise form

q k+1 ← proj C (m k+1 -qk ), (2.24)
where the projection operator is given by proj C (x) = min(max(x, m l ), mu).

2.5.5

Step 5: Dual updates.

Update the scaled dual variables with gradient ascent steps

         Sk+1 ← Sk + S -A(m k+1 )U k+1 , Dk+1 ← Dk + D -PU k+1 , pk+1 ← pk + p k+1 -∇m k+1 , qk+1 ← qk + q k+1 -m k+1 ,
(2.25)

2.5.6
Step 6: Check the stopping condition.

Exit if the preset stopping conditions are satisfied else go to step 1. We will describe the stopping criteria of iterations in the following "Numerical example" section for each numerical example

Hyperparameter tuning

We tune the different penalty parameters by extending the procedure reviewed by Aghamiry et al. (2019b,a) to multiparameter reconstruction.

We start from the last subproblem of the splitting procedure and set the penalty parameters contained in Γ. These hyperparameters control the soft thresholding performed by the TV regularization, equation 2.22. We set γi = 2% max z 2 ix + z 2 iz , where the subscript i ∈ {v0, , δ} denotes the parameter class (mv 0 , m , m δ ). This tuning can be refined by using a different thresholding percentage for each parameter class adaptively during iterations or according to prior knowledge of the geological structure, coming from well logs for example. Also, we use the same weight for the damping regularization associated with the bound constraints and the TV regularization: ζi = γi. Then, we select λ1 as a percentage of the mean absolute value of the diagonal coefficients of ns l=1 L T l L l during the parameter estimation subproblem, equation 2.18. This percentage is set according to the weight that we want to assign to the TV regularization and the bound constraints relative to the wave equation constraint during the parameter estimation. Parameter λ1 may be increased during iterations to reduce the weight of TV regularization and bound constraints near the convergence point. We found this adaptation useful when we start from very crude initial models. Finally, we set λ0 such that λ = λ1/λ0 is a small fraction of the highest eigenvalue ξ of the normal operator A(m) -T P T PA(m) -1 during the wavefield reconstruction subproblem, equation 2.17, according to the criterion proposed by van [START_REF] Van Leeuwen | A penalty method for PDEconstrained optimization in inverse problems[END_REF]. In all the numerical tests, we use λ = 1e-2ξ and λ = 1e-0ξ for noiseless and noisy data, respectively. This tuning of λ is indeed important because it controls the extension of the search space. A too high value of λ reduces the weight of PU -D 2 2 during the wavefield reconstruction and makes IR-WRI behave like a reduced approach. Conversely, using a small value for λ fosters data fitting and expends the search space accordingly. However, a too small value can lead to a prohibitively high number of iterations of the augmented Lagrangian method before the wave equation constraint is fulfilled with sufficient accuracy. Moreover, when data are contaminated by noise, a too small value for λ will make the wavefield reconstruction over-fit the data and drive the algorithm to be a poor minimizer. We always use λ as a fixed percentage of ξ in iterations.

NUMERICAL EXAMPLES

Inclusion test

We first assess multi-parameter IR-WRI with a simple inclusion example for noiseless data. The experimental setup in terms of model, acquisition geometry and frequency selection is identical to that used by Gholami et al. (2013b). The vertical velocity v0 and the Thomsen's parameters δ and are 3 km/s, 0.05 and 0.05, respectively, in the homogeneous background model and 3.3 km/s, 0.1 and 0.2, respectively, in the inclusion of radius 100 m. Nine frequencies between 4.8 Hz and 19.5 Hz are processed simultaneously during IR-WRI and a maximum of 25 iterations is used as stopping criterion for iterations. An ideal fixed-spread acquisition is used, where 64 sources and 320 receivers surround the inclusion, hence providing a complete angular illumination of the target.

Although we use (1/v 2 0 , √ 1 + 2δ, 1 + 2 ) as optimization parameters during our inversions, we show the reconstructed models under the form of v0, δ and for comparison with the results of Gholami et al. (2013b). Note that the radiation patterns of the (1/v 2 0 , √ 1 + 2δ, 1 + 2 ) parametrisation are scaled versions of those of the (v0, δ, ) parametrisation: namely, they exhibit the same amplitude variation with scattering angle, with however different amplitudes. This means that, for an equivalent regularization and parameter scaling, we expect similar resolution and trade-off effects with these two parametrisations. Note that v0 is provided in km/s in our inversion such that the order of magnitude of 1/v 2 0 is of the order of δ and .

We start with bound-constrained IR-WRI with damping (DMP) regularization only (γi = 0, i ∈ {v0, , δ} in equation 2.18) and perform three independent mono-parameter reconstructions for mv 0 , m δ and m , respectively. For each mono-parameter inversion, the true model associated with the optimization parameter contains the inclusion, while the true models associated with the two passive parameters are homogeneous (Fig. 1). For all three inversions, the starting models are the true homogeneous background models. In other words, the data residuals contain only the footprint of the mono-parameter inclusion to be reconstructed. This test is used to assess the intrinsic resolution of IR-WRI for each parameter reconstruction, independently from the cross-talk issue (Gholami et al. 2013b). It is reminded that this intrinsic resolution is controlled by the frequency bandwidth, the angular illumination provided by the acquisition geometry and the radiation pattern of the optimization parameter in the chosen subsurface parametrisation. As this test fits the linear regime of classical FWI, we obtain results (Fig. 1) very similar to those obtained by Gholami et al. (2013b, Their Fig. 9), where the shape of the reconstructed inclusions is controlled by the radiation pattern of the associated parameter. The reader is referred to Gholami et al. (2013b) for a detailed analysis of these radiation patterns. Then, we perform the joint reconstruction of mv 0 , m δ and m , when the true model contains an inclusion for each parameter class (Fig. 2). With our parameter scaling and subsurface parametrisation, v0, and δ are reconstructed with well balanced amplitudes compared to the results of Gholami et al. (2013b, Their Fig. 10), where v0 has a dominant imprint in the inversion. Comparing the models reconstructed by the mono-parameter and multi-parameter inversions highlights however the wavenumber leakage generated by parameter cross-talks (Figs. 1 and2). Then, we complement DMP regularization with TV regularization during bound-constrained IR-WRI for the above monoparameter and multi-parameter experiments (Figs. 3 and4). The results show how TV regularization contributes to remove the wavenumber filtering performed by radiation patterns (compare for example Figs. 1 and3) as well as the cross-talk artifacts during the multi-parameter inversion (compare Figs. 2 and4). Although this toy example has been designed with a piecewise constant model which is well suited for TV regularization, yet it highlights the potential role of TV regularization to mitigate the ill-posedness of FWI resulting from incomplete wavenumber coverage and parameter trade-off.

Synthetic North Sea case study

Experimental setup

We consider now a more realistic 16 km × 5.2 km shallow-water model representative of the North Sea [START_REF] Munns | The Valhall field: a geological overview[END_REF]. The reader is also referred to Gholami et al. (2013a) for an application of acoustic VTI FWI on this model. The true model and the initial (starting) models for v0, δ and are shown in Fig. 5. The subsurface model is formed by soft sediments in the upper part, a pile of low-velocity gas layers above a chalk reservoir, the top of which is indicated by a sharp positive velocity contrast at around 2.5 km depth, and a flat reflector at 5 km depth (Fig. 5a). The initial v0 model is laterally homogeneous with velocity linearly increasing with depth between 1.5 to 3.2 km/s (Fig. 5b), while the δ and initial models are Gaussian filtered version of the true models. Note that our initial/background v0, δ and models are cruder than those used by Gholami et al. (2013a). The fixed-spread surface acquisition consists of 80 (reciprocal) explosive sources spaced 200 m apart at 75 m depth on the sea bottom and 320 (reciprocal) hydrophone receivers spaced 50 m apart at 25 m depth. Accordingly, the pressure wavefield is considered for the acoustic VTI inversion.

A free-surface boundary condition is used on top of the grid and the source signature is a Ricker wavelet with a 10 Hz dominant frequency. We compute the recorded data in the true models (Figs. 5a,c,e) with the forward modelling engine described in Appendix C. During IR-WRI, we use the same forward engine to compute the modelled data according to an inverse crime procedure. Commonshot gathers computed in the true model and in the initial model for a shot located at 14 km are shown in Fig. 6. The seismograms computed in the true model are dominated by the direct wave, the diving waves from the sedimentary overburden, complex packages of pre-and post-critical reflections from the gas layers, the top of the reservoir and the deep reflector. The refracted wave from the deep reflector is recorded at a secondary arrival between -14km and -10km offset in Fig. 6a. Also, energetic reverberating P-wave reflections are generated by the wave guide formed by the shallow water layer and the weathering layer. They take the form of leaking modes with phase velocities higher than the water wave speed [START_REF] Operto | On the role of density and attenuation in 3D multi-parameter visco-acoustic VTI frequency-domain FWI: an OBC case study from the North Sea[END_REF]. The seismograms computed in the starting model mainly show the direct wave and the diving waves, these latter being highly cycle skipped relative to those computed in the true model. We perform both mono-parameter IR-WRI for mv 0 and multi- The background models used for the passive anisotropic parameters (either and δ or δ alone) are the smooth models (Fig. 5d,f), which means that the IR-WRI results will be impacted upon by the smoothness of the passive parameters.

We perform IR-WRI with small batches of two frequencies with one frequency overlap between two consecutive batches, moving from the low frequencies to the higher ones according to a classical frequency continuation strategy. The starting and final frequencies are 3 Hz and 15 Hz and the sampling interval in one batch is 0.5 Hz. The stopping criterion for iterations and for each batch is given by kmax = 15 or

A(m k+1 )U k+1 -S F ≤ ε b & PU k+1 -D F ≤ ε d , (3.1)
where kmax denotes the maximum iteration count. For noiseless data, ε b and ε d are set to 1e-3 and 1e-5, respectively. For noisy data, they are set to 1e-3 and the noise level of the batch, respectively.

We perform three paths through the frequency batches to improve the IR-WRI results, using the final model of one path as the initial model of the next one (these cycles can be viewed as outer iterations of IR-WRI). The starting and finishing frequencies of the paths are [3,6], [4, 8.5], [6,15] Hz respectively, where the first element of each pair shows the starting frequency and the second one is the finishing frequency.

Convexity and sensitivity analysis

Before discussing the IR-WRI results, we illustrate how WRI extends the search space of FWI for the North Sea case study. For this purpose, we compare the shape of the FWI misfit function with that of the parameter-estimation WRI subproblem for the 3 Hz frequency and for a series of v0 and models that are generated according to

v0(α) = v0 true + |α|[v0 init -v0 true ],
(3.2a)

(β) = true + |β|[ init -true], (3.2b)
where v0 true and v0 init denote the true and the initial v0 models, respectively (Fig. 5a,b) and -1 ≤ α ≤ 1. Similarly, true and init are the true and the initial models, respectively (Fig. 5e,f) and -1 ≤ β ≤ 1. Finally, we use the true δ model (Fig. 5c) to generate the recorded data and the smoothed version (Fig. 5d) as a passive parameter to evaluate the misfit function. The misfit functions of the classical reduced-approach FWI as well as that of WRI are shown in Fig. 7. The WRI misfit function is convex, while that of FWI exhibits spurious local minima along both the α and β dimensions. Also, the sensitivity of the misfit function to v0 is much higher than that of for the considered range of models as already pointed out by Gholami et al. (2013b), Gholami et al. (2013a) and Cheng et al. (2016, Their Fig. 2). The weaker sensitivity of the misfit function to makes the joint update of v0 and challenging.

Mono-parameter IR-WRI

We start with mono-parameter bound-constrained IR-WRI when v0 is the optimization parameter and δ and are passive parameters. We first consider noiseless data. The final v0 models inferred from bound-constrained IR-WRI with DMP regularization only and DMP+TV regularizations are shown in Fig. 8. Also, direct comparisons between the logs extracted from the true model, the initial model, and the IR-WRI models at x = 3.5 km, x = 8.0 km and x = 12.5 km are shown in Fig. 9. Although the crude initial v0 model and the smooth δ and passive models, the shallow sedimentary part and the gas layers are fairly well reconstructed with the two regularization settings. The main differences are shown at the reservoir level and below. Without TV regularization, the reconstruction at the reservoir level is quite noisy and the inversion fails to reconstruct the smoothly-decreasing velocity below the reservoir due to the lack of diving wave illumination at these depths. This in turn prevents the focusing of the deep reflector at 5 km depth by migration of the associated short-spread reflections. When TV regularization is used, IR-WRI provides a more accurate and cleaner image of the reservoir and better reconstructs the sharp contrast on top of it. It also reconstructs the deep reflector at the correct depth in the central part of the model, while the TV regularization has replaced the smoothly-decreasing velocities below the reservoir by a piecewise constant layer with a mean velocity (Fig. 9).

To generate more realistic test, we add random noise with Gaussian distribution to the data with SNR=10 db and re-run the mono-parameter bound-constrained IR-WRI with DMP and TV regularizations (Fig. 10). The direct comparison between the true model, the initial model and the IR-WRI models of Fig. 10 at distances x = 3.5 km, x = 8.0 km and x = 12.5 km are shown in Fig. 11. The noise degrades the reconstruction of the gas layers both in terms of velocity amplitudes and positioning in depth when only DMP regularization is used. Also, the sharp reflector on top of the reservoir is now unfocused and mis-positioned in depth accordingly (Fig. 11). The TV regularization significantly reduces these amplitude and mis-positioning errors (Figs. 10b and11) and hence produces v0 models which are much more consistent with those obtained for noiseless data (Figs. 8b and9).

To assess how the differences between the velocity models shown in Fig. 8 impact waveform match, we compute time-domain seismograms in these models as well as the differences with those computed in the true model (Figs. 12 and13). The time-domain seismograms and the residuals shown in Fig. 12 give an overall vi- sion of the achieved data fit, while the direct comparison between the recorded and modelled seismograms shown in Fig. 13 allow for a more detailed assessment of the waveform match for a specific arrival.

A first conclusion is that, for all of the models shown in Fig. 8, the main arrivals, namely those which have a leading role in the reconstruction of the subsurface model (diving waves, pre-and postcritical reflections), are not cycle skipped relative to those computed in the true model. We note however more significant residuals for the reverberating guided waves when TV regularization is used. These mismatches were generated by small wavespeed errors generated in the shallow part of the model by the TV regularization. These artifacts can be probably corrected by deactivating or decreasing the weight of the TV regularization locally. For noiseless data, the direct comparison between the seismograms computed in the true model and in the reconstructed ones show how the DMP+TV regularization improves the waveform match both at pre-and post-critical incidences relative to the DMP regularization alone (Fig. 13a,b). For noisy data, the data fit is slightly degraded by noise when DMP regularization is used, while the TV regularization produces a data fit which is more consistent with that obtained with noiseless data (Compare Figs. 13a and13c for DMP regularization, and Figs. 13b and13d for DMP+TV regularization). It is striking to see the strong impact of noise on the quality of the v0 reconstruction when DMP regularization is used (compare Figs. 8a and 10a, and Figs. 9 and 11, green curves), compared to its more moderate impact on the data fit (Compare Figs. 13a and13c). This highlights the ill-posedness of the FWI, which is nicely mitigated by the prior injected by the TV regularization as illustrated by the consistency of the v0 models inferred from noiseless and noisy data (Compare Figs. 13b and13b).

Multi-parameter IR-WRI

We continue with multi-parameter bound-constrained IR-WRI when v0 and are optimization parameters and δ is a passive Figure 11. North Sea case study with noisy data (SNR=10 db). Direct comparison along the logs at x = 3.5 (left), x = 8.0 (center) and x = 12.5 km (right) between the true velocity model (black), the initial model (dashed line) and the bound constrained IR-WRI estimated models with DMP, Fig. 8a, (green) and TV regularization, Fig. 8b,(red). optimization parameters allows us to better conciliate the fit of the pre-and post-critical reflections (Compare Figs. 13a and 19a for DMP regularization,and Figs. 13b and 19b for DMP+TV regularization). With noiseless data, DMP regularization is enough to achieve a high-fidelity data fit, which looks better than that obtained with DMP+TV regularization (Compare Figs. 19a and 19b).

We repeat now the joint inversion when data are contaminated by Gaussian random noise with a SNR=10 db. The final v0 and models inferred from bound-constrained IR-WRI with DMP and DMP+TV regularization are shown in Fig. 16. Direct comparisons along logs extracted from the true models, the initial models, and the models inferred from the bound constrained multi-parameter IR-WRI with DMP and DMP+TV regularization at x = 3.5 km, x = 8.0 km and x = 12.5 km are shown in Fig. 17. For the v0 reconstruction, a trend similar to that shown for the mono-parameter inversion is shown, with a more significant impact of the noise on the velocity model reconstructed with DMP regularization compared to the one reconstructed with DMP+TV regularization (Compare Figs. 14 and16). As expected, the impact of noise is more significant on the second-order model, even when TV regularization is used, in the sense that the estimated perturbations of the initial model have smaller amplitudes compared to the noiseless case (Fig. 17b).

The time-domain seismograms are shown in Figs. 18 and19c,d. The data fit obtained with DMP regularization has been significantly degraded compared to that obtained with noiseless data (compare Figs. 19a and 19c), while the data fit obtained with DMP+TV regularization is more consistent with noiseless and noisy data (compare Figs. 19b and 19d). This is consistent with the previous conclusions drawn from the mono-parameter inversion.

DISCUSSION

We have extended the ADMM-based wavefield reconstruction inversion method (IR-WRI), originally developed for monoparameter wavespeed reconstruction (Aghamiry et al. 2019c,b), to multi-parameter inversion in VTI acoustic media. We have first discussed which formulations of the VTI acoustic wave equation are bilinear in wavefield and subsurface parameters. First-order velocity stress form is often more convenient than the second-order counterpart to fulfill bilinearity, in particular if density (or buoy- ancy) is an optimization parameter. However, it may be not the most convenient one for frequency-domain wavefield reconstruction as the size of the linear system to be solved scales to the number of wavefield components. To bypass this issue, wavefield reconstruction and parameter estimation can be performed with different wave equations, provided they give consistent solutions (Gholami et al. (2013a) and this study).

Bilinearity allows us to recast the parameter estimation subproblem as a linear subproblem and hence the waveform inversion as a biconvex problem. Although ADMM has been originally developed to solve distributed convex problems, it can be used as is to solve biconvex problem [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction of multipliers[END_REF]. From the mathematical viewpoint, biconvex problems should have superior convergence properties compared to fully nonconvex problem (this discussion is out of the scope of this study but we refer the interested reader to [START_REF] Benning | Preconditioned admm with nonlinear operator constraint[END_REF]. Alternatively, ADMMlike optimization can be used to perform IR-WRI heuristically without bilinear wave equation, hence keeping the parameter estimation subproblem, equation B.4a, nonlinear. Accordingly, equation B.4a would be solved with a Newton algorithm rather than with a Gauss-Newton one. This nonlinear updating of the parameters may however require several inner Newton iterations per IR-WRI cycle, while Aghamiry et al. (2019c) showed that one inner Gauss-Newton iteration without any line search was providing the most efficient convergence of IR-WRI when bilinearity is fulfilled.

Indeed, the bilinearity specification limits the choice of subsurface parametrisation for parameter estimation. In the general case of triclinic elastodynamic equations, a subsurface parametri- sation involving buoyancy and stiffness or compliance coefficients will be the most natural ones, as they correspond to the coefficients of the equation of motion and the Hooke's law. In the particular case of the VTI acoustic wave equation, we have developed a bilinear wave equation whose coefficients depend on the vertical wavespeed v0 and the Thomsen's parameters δ and . Although v0 and are coupled at wide scattering angles, the (v0, , δ) parametrisation was promoted by Gholami et al. (2013b) and Gholami et al. (2013a) because the dominant parameter v0 has a radiation pattern which doesn't depend on the scattering angle, and hence can be reconstructed with a high resolution from wide-azimuth long-offset data. The counterpart is that updating the long wavelengths of the secondary parameter is challenging and requires so far a crude initial guess of its long wavelengths (Fig. 5e,f), which can be used as prior to regularize the update. Comparing the results of IR-WRI when is used as a passive parameter and as an optimization parameter shows that the sensitivity of the inversion to remains small provided that a reasonable guess of its long wavelengths are provided in the starting model (for the models, compare Figs. [8][9][10] for the data fit, compare Fig. 13 and Fig. 19). This limited sensitivity of FWI to the short-to-intermediate wavelengths of in the (v0, , δ) prompted for example Debens et al. (2015) to estimate a crude background model with a coarse parametrisation by global optimization.

Among the alternative parametrisations proposed for VTI acoustic FWI, [START_REF] Plessix | A parametrization study for surface seismic full waveform inversion in an acoustic vertical transversely isotropic medium[END_REF] proposed the (vn, v h , δ) or the (vn, η, δ) parametrisations for long-offset acquisition, while [START_REF] Alkhalifah | A recipe for practical full-waveform inversion in anisotropic media: An analytical parameter resolution study[END_REF] promoted the (v h , η, ) parametrisation, where vn = v0 √ 1 + 2δ is the so-called NMO velocity, v h = v0 √ 1 + 2 is the horizontal velocity and η = ( -δ)/(1+2δ) represents the anellipticity of the anisotropy. For the parametrisation promoted by [START_REF] Plessix | A parametrization study for surface seismic full waveform inversion in an acoustic vertical transversely isotropic medium[END_REF], the VTI equation developed by [START_REF] Zhou | An anisotropic acoustic wave equation for VTI media[END_REF] given by to be smooth, the above equation can be approximated as

ω 2 diag 1 v 2 n uq + 2diag(η)∇xx(up + uq) = sq, ω 2 diag 1 v 2 n up + ∇xx(up + uq) + diag 1 √ 1 + 2δ ∇zzdiag 1 √ 1 + 2δ up = sp, ( 4
ω 2 diag 1 v 2 n uq + 2diag(η)∇xx(up + uq) = sq, ω 2 diag 1 v 2 n up + ∇xx(up + uq) + diag 1 1 + 2δ ∇zzup = sp (4.2)
which is bilinear in wavefields and parameters (1/v 2 n , η, 1/(1 + 2δ). This implies that δ can be involved as an optimization parameter if necessary. Note that, if this smoothness approximation is used to update the parameters during the primal problem, the modelled data and the source residuals can be solved with the exact equation to update the dual variables with more accuracy. For the (v h , η, ) parametrisations [START_REF] Alkhalifah | A recipe for practical full-waveform inversion in anisotropic media: An analytical parameter resolution study[END_REF], according to Zhou H. ( 2006) the VTI equations with smooth δ can be written as

ω 2 diag 1 v 2 h uq + ∇xx(up + uq) -diag 1 1 + 2η ∇xx(up + uq) = sq, ω 2 diag 1 v 2 h up + diag 1 1 + 2η ∇xx(up + uq) + diag 1 1 + 2 ∇zzup = sp, (4.3)
which is bilinear in wavefields and parameters (1/v 2 h , 1/(1 + 2η), 1/(1 + 2 )).

CONCLUSION

We have shown that ADMM-based IR-WRI can be extended to multi-parameter reconstruction for VTI acoustic media. The gradient of the misfit function for the parameter estimation subproblem involves the so-called virtual sources, which carry out the effect of the parameter-dependent radiation patterns. This suggests that, although IR-WRI extends the search-space of FWI to mitigate cycle skipping, it is impacted by ill-posedness associated with parameter cross-talk and incomplete angular illumination as classical FWI. We have verified this statement with a toy numerical example for which we have reproduced the same pathologies in terms of resolution and parameter cross-talks as those produced by classical FWI during a former study. We have illustrated how equipping IR-WRI with bound constraints and TV regularization fully remove the ill-posedness effects for this idealized numerical example. We have provided some guidelines to design bilinear wave equation of different order for different subsurface parametrisations. Although bilinearity puts some limitations on the choice of the subsurface parametrisation, it recasts the parameter estimation subproblem as a quadratic optimization problem, which can be solved efficiently with Gauss-Newton algorithm. Application on the long-offset synthetic case study representative of the North Sea has shown how IR-WRI can be started from a crude laterally-homogeneous vertical velocity model without impacting the inversion with cycle skipping, when a smooth δ parameter is used as a passive background model. However, a smooth initial model, albeit quite crude, is necessary to guarantee the convergence of the method to a good solution, either when is used as a passive or as an optimization parameter. For this case study where the low velocity gas layers and the smooth medium below the reservoir suffer from a deficit illumination of diving waves, the TV regularization plays a key role to mitigate the ill-posedness. parameters, and hence the bilinearity of the wave equation. For the lth source, the new equation becomes

    L11 0 0 0 L21 0 0 0 0 L32 L33 0 0 0 L43 L44         b ∫11 ∫13 ∫33     =     0 -îωv x,l 0 -îωv z,l îωs l -∇xv x,l îωs l -∇zv z,l     , (A.3) where          L11 = diag(∇xu x,l ), L21 = diag(∇zu z,l ), L32 = L43 = îωdiag(u x,l ), L33 = L44 = îωdiag(u z,l ).
To develop the second-order wave equation, we eliminate v x,l and v z,l from equation A.1. We obtain the following equation

1 ω 2 ∇xdiag(b)∇xu x,l + diag(∫11)u x,l + diag(∫13)u z,l = s l , 1 ω 2 ∇zdiag(b)∇zu z,l + diag(∫13)u x,l + diag(∫33)u z,l = s l , (A.4)
which is bilinear with respect to buoyancy, compliance parameters and pressure wavefields. With known buoyancy and compliance parameters, we get the following 2n × 2n linear system to estimate wavefields for all sources A(m)U = S, (A.5)

where A is given by

1 ω 2 ∇xdiag(b)∇x + diag(∫11) diag(∫13) diag(∫13) 1 ω 2 ∇zdiag(b)∇z + diag(∫33) . and m =     b ∫11 ∫13 ∫33     ,
and S is defined in equation 2.8. When U is known, this system can also be recast as a new linear system in which the unknowns are the model parameters as

        L1 . . . L l . . . Ln s         m =         y1 . . . y l . . . yn s         , (A.6)
where L l is given by where the diagonal weighting matrix Ξ is defined as

Ξ =   λ1I 0 0 0 Γ 0 0 0 Z   ∈ R [3+ns]n×[3+ns]n .
Equation B.6 is now quadratic and admits a closed form solution as given in equation 2.18 (step 2 of the algorithm).

The only remaining tasks consist in determining the auxiliary primal variables (p, q), equations B.4b,c, and the auxiliary dual variables (p, q). They are initialized to 0 and are updated as follows: The primal variables p is updated through a TV proximity operator, which admits a closed form solution given by equation 2.22 (see Combettes & Pesquet 2011) (step 3 of the algorithm). The primal variable q is updated by a projection operator, which also admits a closed form solution given by equation 2.24 (step 4 of the algorithm). Finally, the duals are updated according to gradient ascend steps (step 5 of the algorithm) pk+1 = pk + p k+1 -∇m k+1 qk+1 = qk + q k+1 -m k+1 .

(B.7)

APPENDIX C: USING FOURTH-ORDER EQUATION FOR WAVEFIELD RECONSTRUCTION (STEP 1 OF THE ALGORITHM)

The wavefield reconstruction subproblem, equation 2.17, can be written as the following over-determined system where P is the sampling operator of each component of the wavefield at receiver positions and P = [ 1 2 P 1 2 P] (the coefficient 1 2 results because the (isotropic) pressure wavefield u recorded in the water is given by u = 1 2 (ux + uz)). We can eliminate Uz from the equation C.1 to develop a fourth-order partial-differential equation for Ux and then update Uz from its explicit expression as a function of Ux without any computational burden. By multiplying the second row of equation C.1 by diag(m δ ) and taking the difference with the first row, we find ). The former can be accurately discretized with the 9-point finite-difference method of [START_REF] Chen | An optimal 9-point finite difference scheme for the Helmholtz equation with PML[END_REF], while the anelliptic term can be discretized with a basic second-order accurate 5-point stencil without generating significant inaccuracies in the modelling [START_REF] Operto | Computationally-efficient three-dimensional visco-acoustic finite-difference frequency-domain seismic modeling in vertical transversely isotropic media with sparse direct solver[END_REF]).

  I, C = ∇ and D = îωI. Equation 2.7 is linear in U and V when the model parameters embedded in M b and Mc are known.

Figure 1 .

 1 Figure 1. Inclusion test: Mono-parameter IR-WRI results. For v 0 (first column), δ (second column), and (third column). The initial models are the true homogeneous background models. The vertical and horizontal profiles in the true model (blue) and estimated model (red) are extracted across the center of the inclusions.

Figure 2 .

 2 Figure 2. Inclusion test: Multiparameter IR-WRI results for joint update of v (first column), δ (second column), and (third column). The initial models are the homogeneous background models.

Figure 3 .Figure 4 .

 34 Figure 3. Inclusion test: Same as Fig. 1, but with TV regularization.

Figure 5 .Figure 6 .

 56 Figure 5. North Sea case study. (a) True v 0 model. (b) Initial v 0 model. (c) True δ model. (d) Smoothed δ model, which is used as passive background model during inversion tests. (e) True model. (f) Smoothed model, which is used as a passive background model during the mono-parameter inversion and as an initial model during the joint reconstruction of v 0 and . The vertical dashed lines in (b) and (f) indicate the location of vertical logs.

Figure 7 .

 7 Figure 7. North Sea case study. 3-Hz misfit function for the first iteration when v 0 and are active parameters and δ is kept fixed to the smooth background model (Fig. 5f). (a) Classical reduced approach, (b) WRI. The variables α and β parametrize the v 0 and models, respectively, for which the misfit function is computed (see text for details).

Figure 8 .

 8 Figure 8. North Sea case study with noiseless data. Mono-parameter IR-WRI. Vertical wavespeed v 0 models inferred from IR-WRI with (a) DMP regularization, (b) TV regularization.

Figure 9 .

 9 Figure 9. North Sea case study with noiseless data. Direct comparison along the logs at x = 3.5 (left), x = 8.0 (center) and x = 12.5 km (right) between the true velocity model (black), the initial model (dashed line) and the bound constrained IR-WRI models shown in Fig. 8. The logs of the IR-WRI models obtained with DMP regularization and DMP+TV regularization are the green lines and the red lines, respectively.

Figure 10 .

 10 Figure 10. North Sea case study with noisy data (SNR=10 db). Monoparameter bound constrained IR-WRI with (a) DMP (b) TV regularization.

Figure 12 .Figure 13 .

 1213 Figure 12. North Sea case study. Time-domain seismograms computed in the subsurface models obtained by (v 0 ) mono-parameter IR-WRI. (a-b) Noiseless data: (a) DMP regularization, (b) DMP+TV regularization. (cd) Same as (a-b) for noisy data. (e-h) Residuals between the seismograms computed in the true model (Fig. 6a) and those shown in (a-d). The seismograms are plotted with a reduction velocity of 2.5 km/s.

Figure 14 .

 14 Figure 14. North Sea case study with noiseless data. Multi-parameter IR-WRI results. (a-b) v 0 models obtained with (a) DMP (a) and (b) DMP+TV regularizations. (c-d) Same as (a-b) for .

Figure 15 .

 15 Figure 15. North Sea case study with noiseless data. Multi-parameter bound-constrained IR-WRI results. Direct comparison along logs at x = 3.5 (left), x = 8.0 (center) and x = 12.5 km (right) between the true model (black), the initial model (dashed line) and the models inferred from IR-WRI with DMP (green) and DMP+TV regularization (red) for (a) v 0 and (b) (Fig. 14).

Figure 16 .

 16 Figure 16. North Sea case study with noisy data (SNR=10 db). Multiparameter bound-constrained IR-WRI results. (a-b) v 0 models obtained with (a) DMP (a) and (b) DMP+TV regularizations. (c-d) Same as (a-b) for .

  .1) is bilinear in wavefields and parameters (1/v 2 n , η), where up = √ 1 + 2δuz and uq = ux -√ 1 + 2δuz. Note that if δ is assumed

Figure 17 .

 17 Figure 17. North Sea case study with noisy data. Multi-parameter boundconstrained IR-WRI results. Direct comparison along logs at x = 3.5 (left), x = 8.0 (center) and x = 12.5 km (right) between the true model (black), the initial model (dashed line) and the models inferred from IR-WRI with DMP (green) and DMP+TV regularization (red) for (a) v 0 and (b) (Fig. 16).

Figure 18 .

 18 Figure 18. North Sea case study. Same as Fig. 12 for (v 0 , ) multi-parameter IR-WRI.

Figure 19 .

 19 Figure19. North Sea case study. Same as Fig.13for (v 0 , ) multi-parameter IR-WRI.

1 ω 2 F+

 2 ∇xdiag(∇xu x,l ) diag(u x,l ) diag(u z,l ) 01 ω 2 ∇zdiag(∇zu z,l ) 0 diag(u x,l ) diag(u z,l )and y l = s l s l .APPENDIX B: SOLVING THE OPTIMIZATION PROBLEM, EQUATION 2.16, WITH ADMMStarting from an initial model m 0 , and setting the dual variables D and S equal to zero, the kth ADMM iteration for solving equation 2.16 reads as (seeAghamiry et al. 2019c,b,d, for more details)U k+1 ← arg min U C m k , Dk , Sk (U), (B.1a) m k+1 ← arg min m∈C C U k+1 , Sk (m), (B.1b) Sk+1 ← Sk + S -A(m k+1 )U k+1 , (B.1c) Dk+1 ← Dk + D -PU k+1 , (B.1d)whereC m k , Dk , Sk (U) = λ0 PU -D -Dk 2 λ1 A(m k )U -S -Sk 2 F (B.2)andC U k+1 , Sk (m) = |∂xm| 2 + |∂zm| 2 + λ1 A(m)U k+1 -S -Sk 2 F . (B.3)The regularized wavefields U are the minimizers of the quadratic cost function C m k , Dk , Sk (U), equation B.2, where m k , Dk and Sk are kept fixed. Zeroing the derivative of C m k , Dk , Sk (U) gives the wavefields as the solution of a linear system of equations defined by equation 2.17 (step 1 of the algorithm). The regularized wavefields are then introduced as passive quantities in the cost function C U k+1 , Sk (m), equation B.3, which is minimized to estimate m over the desired set C. We solve this minimization subproblem with the splitting techniques. Accordingly, we introduce the auxiliary primal variables p and q to decouple the 1 and 2 terms and split the parameter estimation subproblem into three sub-steps for m, p and q[START_REF] Goldstein | The split Bregman method for L1regularized problems[END_REF] Aghamiry et al. 2019b):m k+1 ← arg min m λ1 A(m)U k+1 -S -Sk 2 F , (B.4a) + ∇m -p kpk 2 Γ + mq kqk 2 Z , p k+1 ← arg min p |px| 2 + |pz| 2 (B.4b) + ∇m k+1 -ppk 2 Γ , q k+1 ← arg min q∈C m k+1 -qqk 2 Z , T •x, ZandΓ are diagonal matrices defined in equation 2.20 and 2.21, respectively.From the linearized equation 2.13, the subproblem for m, equation B.4a, can be written as m k+1 ← arg min m

  (mv 0 ) + diag(m )∇xx diag(m δ )∇zz diag(m δ )∇xx ω 2 diag(mv 0 ) + ∇zz

  (mv 0 ) + diag(m )∇xx diag(m δ )∇zz ω 2 diag(mv 0 ) + diag(m -m 2 δ )∇xx -ω 2 diag(m δ • mv 0 ) δ = m δ • m δ .The second equation of C.2 provides us the closed-form expression of Uz as a function of UxUz = AzUx + Bz, -m 2 δ )∇ xx+ω 2 diag(mv 0 ) ω 2 diag(mv 0 •m δ ) , Bz = diag(m δ )Sz -Sx ω 2 diag(mv 0 •m δ ) .(C.4) Injecting equation C.4 into the first equation of C.2 leads to the over-determined system satisfied by Ux ω 2 diag(mv 0 ) + diag(m )∇xx + diag(m δ )∇zzAz of solving the (nr + 2n) × 2n linear system C.1 to update U, Ux is updated by solving the (nr + n) × n linear system C.5 and then Uz is updated using C.4 without significant computational overhead. The wave equation operator in equation C.4 has been broken down into an elliptic wave equation operator and anelliptic correction term (the term of Az related to diag(m -m 2 δ ), equation C.4

ACKNOWLEDGMENTS

We would like to thank the Editor Huajian Yao, Xin Cheng and Daniel Köhn for their comments. This study was partially funded by the SEISCOPE consortium (http://seiscope2.osug.fr), sponsored by AKERBP, CGG, CHEVRON, EQUINOR, EXXON-MOBIL, JGI, PETROBRAS, SCHLUMBERGER, SHELL, SINOPEC and TOTAL. This study was granted access to the HPC resources of SIGAMM infrastructure (http://crimson.oca.eu), hosted by Observatoire de la Côte d'Azur and which is supported by the Provence-Alpes Côte d'Azur region, and the HPC resources of CINES/IDRIS/TGCC under the allocation 0596 made by GENCI.

We consider the frequency-domain first-order velocity-stress equation in 2D VTI acoustic media with compliance notation as

where ∫ij ∈ R n×1 are the compliance coefficients, and the other notations are defined after equation 2.6. Gathering equation A.1 for all sources results in the following matrix equation:

where

and the other notations are defined in equation 2.7. Equation A.2 is linear in U and V when the model parameters embedded in B and S are known. When U and V are known, this system can be recast as a new linear system in which the unknowns are the model


E Ender

M A Subirana

A Raab

E M Krupp

D Schaumlöffel

J Feldmann

Why is NanoSIMS elemental imaging of arsenic in seaweed 3 (Laminaria digitata) important for understanding of arsenic 4 biochemistry in addition of speciation information? 5 6

Brown seaweed such as Laminaria digitata is known to accumulate arsenic to more than 100 16 mg/kg. How the algae can tolerate such high level of arsenic has traditionally been studied by 17 arsenic speciation analysis using HPLC-ICPMS, but the knowledge of its molecular forms 18

has not yet given any answers. Here we demonstrate for the first time that the combination of 19 speciation analysis with high resolution imaging by NanoSIMS and TEM identifies not only 20 the molecular structures of arsenic but also the location of arsenic in cells and cell 21 substructures in a brown seaweed species.

22

The majority of 117 mg/kg arsenic in L. digitata fronds was in the form of inorganic arsenic 23 (53%) and arsenosugars (32%) and only 1.5% of total arsenic as arsenolipids (mainly as 24 AsHC and AsPL). The lateral resolution of 300 nm and the concentration of arsenic was high 25 enough for the localization of arsenic in the cells of the seaweed using NanoSIMS. The 26 majority of arsenic was found in the cell walls and cell membrane, while the inside of the cell 27 was almost arsenic free, which is not expected if the majority of arsenic species are 28 hydrophilic. The NanoSIMS images questions the integrity of the arsenic species during 29 extraction for the speciation analysis and that inorganic arsenic is unlikely to occur freely in The use of analytical atomic spectrometry and in particular ICPMS has been instrumental to 37 identify the molecular forms of arsenic in biological samples over nearly three decades. The 38 complex arsenic biochemistry of more than 100 naturally occurring organoarsenicals in 39 addition to arsenite and arsenate made it necessary to use various HPLC methods coupled to 40 ICPMS and ESI-MS detection in off-line or on-line mode to identify the full spectrum of 41 arsenic species in biological samples even without standards. i However, neither the 42 biosynthetic pathways nor the classification whether these metabolic transformations are 43 detoxification pathways for inorganic arsenic or the formation of beneficial arsenic 44 containing biomolecules have been elucidated; although credible suggestions have recently 45 been made about the arsenosugars. ii What is missing is the localisation of the arsenic inside 46 the cells in order to answer the above mentioned questions about the biochemistry of arsenic.

48

It is well known that seaweeds bioaccumulate large concentrations of arsenic from seawater 49 with bioconcentration factor between 1,000 and 100,000. iii Seaweeds and in particular brown 50 macroalgae can reach concentrations up to 200 mg As/kg. iv Most of the arsenic in seaweeds 51 is present as organoarsenicals in the form of arsenosugars (As-sugars), dimethylarsinic acid 52 (DMA) and also in the lipid soluble forms of arsenic containing hydrocarbons (AsHC), fatty 53 acids (AsFA) and phospholipids (AsPL) v,vi . Only Hijiki spp vii and Laminaria digitata viii,ix 54 seem to contain large amounts of inorganic arsenic.

56

Here we would like to illustrate the complementary information of using detailed arsenic 57 speciation analysis and mapping of arsenic at sub-cellular level using the brown seaweeds 58 Laminaria digitata as examples to gain more detailed insight into the bioaccumulative 59 behaviour of arsenic especially when the arsenic species are known. Special attention will be 60 given to the differences in the sample preparation procedures. While often the samples are 61 freeze-dried before speciation analysis by HPLC-ICPMS is performed in order to report on 62 concentration relative to dry mass, is it important to have fresh samples available for 63 transmission electron microscope (TEM) and nanoscale secondary ion mass spectrometry 64 (NanoSIMS) analysis which need to be prepared by cryofixation using high-pressure freezing 65 followed by cryo-substitution and resin embedding. This is necessary in order to conserve the 66 cellular ultrastructure and to fix the elements and not distribute them at a sub-cellular level. diameter were prepared with a biopsy punch and transferred into gold-plated copper specimen 97 carriers filled with hexadecane. High Pressure Freezing was then carried out using a Leica EM 98 ICE (Leica Microsystems, Milton Keynes, UK). For freeze substitution samples were post-99 fixed with 2% osmium tetroxide (OsO 4 ) in acetone and transferred into an automatic freeze 100 substitution system (AFS 2, Leica Mircosystems) following the programme outlined in through a stepwise infiltration of the samples with a graded resin and acetone series (10, 30, 103 50, 70, and 90% resin in acetone). Finally, embedded in 100% resin the samples were 104 polymerized at 60 o C for at least 24 hours. For TEM analysis, 90 nm sections were cut from 105 the resin blocks using a diamond knife (Diatome Ltd, Switzerland) on an ultra-microtome 106 (Leica UC6, Leica Microsystems) and placed onto copper grids (TAAB, UK). For NanoSIMS 107 analysis, 300 nm sections were cut by the same method and placed on silicon wafers (Wafer 108 Solution, Le Bourget du lac, France). The arsenic concentration in the freeze-dried L. digitata was 117 +/-30 mg As/kg (Table 1),

160
which is in the normal range as reported previously. viii The accuracy and precision of the total 161 arsenic was evaluated using a seaweed CRM and the recovery was quantitative. The speciation of the hydrophilic fraction of arsenic gave the expected distribution of arsenic 165 species, with 53% of the total arsenic in the form of inorganic arsenic (62 +/-19 mg As/kg).

166

Arsenosugars in the form of As-sugars (-glycerol, phosphate ester and sulfonate) contributed 167 32% (37 mg As/kg), while DMA was only a minor compound with less than 1 % of total As.

The arsenolipid fraction contributed also only 1.5% of total arsenic of which the majority 169 65% was in the form of arsenosugar containing phospholipids (AsPL) (structures see Figure 170 1). The identification and quantification was done using RP-HPLC-ICPMS/MS coupled 171 simultaneously to ESI-MS. The assignment of the arsenic peaks was achieved by using the 172 retention time as the first criteria, which should coexist with the corresponding protonated 173 mass of the molecular mass (ESI-MS) as can be seen in Figure S1a The NanoSIMS analysis method has been optimized with regard to maximal spatial/lateral 225 resolution for sub-cellular imaging and high sensitivity for detection of arsenic at trace level.

226

For high lateral resolution the Cs + primary ion beam was trimmed using a diaphragm (D1-4)

242

An advantage of the NanoSIMS is that it combines high lateral resolution with high sensitivity.

243

In order to accumulate a sufficient number of counts for the detection of trace amounts of 244 arsenic a compromise had to be found between the dwell time of the primary beam per pixel, In Figure 3ii sulphur was mainly found in the thick cell walls. This intense sulphur signal can 264 be related to the presence of fucoidans, the sulphated polysaccharides typically found in seaweed cell walls of brown algae such as L. digitata xix . Moreover, sulphur was also detected 266 in the chloroplast which can be explained by sulphur-containing proteins involved in 267 photosynthesis, such as ferredoxin and iron-sulphur clusters.

268

The phosphorous image is an accumulation of 5 scans. P was detected in concentrated small 269 spots inside the cell which can be related to in 0.5-1 µm sized vacuoles, which can be seen as 270 electronically dense sub-cellular structures in the TEM image (Figure 2). In green algae such In Figure 3iii arsenic was found in a pattern consistent with the cell wall/intercellular material.

278

The image is an accumulation of 5 scans. The accumulation of arsenosugars in the cell wall is however more difficult to understand.

317

Since no arsenosugars have been found in seawater so far, the accumulation of these arsenic 324 digitata is the AsSugar-SO 3 ). These compounds again would be rather unstable in water and 325 hydrolyse quickly. The fact that no hotspot inside the cell which could be assigned to the Golgi 326 apparatus has been found, indicates that the arsenic transformation inside the cell needs to be 327 fast. Alginate generation is taking place within 20 minutes. This would explain that the cytosol 328 and the sub-cellular structures inside the cell are virtually arsenic free.

30

  the seaweed. Whether inorganic arsenic and the arsenosugars are bound directly to the 31 polymeric carbohydrates alginates or fucoidans in the seaweed is unclear and needs further

  reagents were purchased from Sigma-Aldrich (Saint-Quentin Fallavier, 72 France) unless stated otherwise. All solutions, dilutions and preparations were made with water 73 (18.2 MΩ cm) obtained from a Milli-Q system (Millipore, Bedford, Ma, USA) unless stated 74 otherwise. A 10 µL drop of an arsenic ICP Standard, 1000 µg/mL (PlasmaCAL, Villebon-sur-75 Yvette, France), was dried on a silica wafer and used as standard for mass calibration and mass 76 resolution tuning for As in NanoSIMS. 77 Chemicals used for extraction were of lab-grade quality (Fisher, UK). Chemicals used for total 78 As determination were of trace-element grade (VWR, UK). HPLC grade methanol (Rathburn, 79 UK), ammonium carbonate (Sigma, UK) and formic acid (mass-spec quality, Fluka, UK) were 80 used for the preparation of HPLC eluents. Nitric acid (70 % p.a. Fisher, UK) and hydrogen 81 peroxide (32 % Fisher, UK) were used for extraction and digestion. For quantification of As-82 species sodium dimethylarsinic acid (DMA, Chemservice, USA) was dissolved in deionized 83 water and diluted as necessary. Identification of anionic As-species was aided by species 84 standards of DMA, MA (Chemservice, USA), arsenite and arsenate (BDH, UK). Seaweed 85 CRM (ERM-CD-200) a brown algae Fucus vesiculosus was purchased from JRC-Institute for 86 Reference Materials and Measurements (Geel, Belgium).

  at a beach south of Aberdeen, on the east coast of Scotland 90 (57.139856 N, -2.051430 W) in November 2018. Freshly washed ashore whole thalli of 91 Laminaria digitata were harvested.

  conducted at the Microscopy and Histology Core Facility at the 95 University of Aberdeen. Two replicates of L. digitata were prepared. Algae discs of 4 mm 96

  totals arsenic analysis and speciation analysis 111 Laminaria digitata fronds were cleaned by scraping and washing off epiphytes and sediment, freeze 112 dried and ground with mortar and pestle to a fine powder. Only the young fronts were analysed. 113 114 Extraction of lipophilic arsenic 115 Algae (1 g) were extracted with hexane (2 × 5 mL), followed by CH 2 Cl 2 /MeOH (2:1 v/v, 2 × 116 5 mL), and both extracts were evaporated to dryness. For speciation analysis only the MeOH 117 soluble part of the CH 2 Cl 2 /MeOH extract was used. Petursdottir et al. x was used for the extraction. In brief 0.1 g algae were mixed 121 with 10 mL of extraction solution (1 % (v/v) nitric acid and 2 % (v/v) hydrogen peroxide) and 122 heated in a microwave (Mars5, CEM, UK) to 95ºC for 30 min. The solution was centrifuged 123and the supernatant used for determination of inorganic arsenic.

  determined after standard microwave digestion with HNO 3 /H 2 O 2 using a Mars5 127 system (CEM, UK) in solid L. digitata. An ICPMS/MS (Agilent 8800, UK) was used in mass 128 shift mode to determine As on m/z 91 as AsO + . 129 130 Determination of hydrophilic and lipophilic arsenic species 131 Hydrophilic arsenic species were determined by using anion exchange HPLC-ICPMS/MS. The 132 parameters are given in the supplementary material. The lipophilic arsenic speciation was 133 performed by coupling reverse-phase HPLC simultaneously to ICPMS/MS and ESI-qTOF-MS 134 as described previously xi . The details are given in the supplementary material. 136 TEM analysis 137 Samples were viewed on the transmission electron microscopes JEM 1400 plus (JEOL, UK) at 138 an accelerating voltage of 100 kV using an AMT UltraVUE camera. The TEM instrument was 139 used to determine quality of sample preparation and for high resolution imaging of sub-cellular 140 structures. Copper grids have been used for sample sections. (Cameca, Genneviliers, France) was used for high resolution secondary ion 144 mass spectrometry analysis. A primary Cs + ion source served for the mapping of 145 electronegative elements. Detectors were tuned for carbon ( 12 C 2 -), nitrogen ( 12 C 14 N -), 146 phosphorus ( 31 P -), sulphur ( 32 S -) and arsenic ( 75 As -). As the concentration of As in the sample 147 was relatively low, an As standard was used for mass calibration and mass resolution tuning 148 process. Cs + ions were implanted onto the sample surface in order to increase sensitivity. It 149 should be noted that N cannot be directly measured using NanoSIMS, therefore carbon-150 nitrogen (CN -) cluster ions are measured to detect N in the sample. Images were acquired with 151 a raster size ranging of 35 μm, divided into 512 x 512 pixels with a dwell time of 10 ms per 152 pixel. Individual sample sections were located using the NanoSIMS Charged Coupled Device 153 (CCD) camera. The field of view of the CCD camera is about 500 × 600 μm and the optical 154 resolution around 1 μm, which allows for an overview of the sample and a selection of regions 155 of interest for NanoSIMS analysis.

Figure 1 :

 1 Figure 1: generic structures of major organoarsenic compounds in L. digitata: arsenosugar

218 Figure 2 .

 2182 Figure 2. TEM image obtained from 90 nm section of resin embedded L. digitata cell. Image

227

  with a circled aperture of 150 µm and finally focused with the primary focusing lens (EOP) to 228 its final spot size of about 50 nm at the sample surface. This lateral resolution was confirmed 229 in our previous work by applying the knife-edge method on a standard sample xviii and is in 230 agreement with the instrument specifications. In this work, the actual lateral resolution in the image of the seaweed sample was determined by line scans at a sharp edge in the CN -and S - 232 images (Figures3ib and 4iia) using the 16-84% criterion (knife-edge method). Depending on 233 the location in the image the actual resolution was about 300 to 400 nm. Note, that in biological 234 samples not only the NanoSIMS primary beam size, but also the structure of the sample as well 235 as the quality of sample preparation which influences the element distribution in the sample 236 structure contributes to the measured lateral resolution. Therefore in biological samples usually 237 a lower resolution is measured than in standards for the knife-edge method. The selected pixel 238 numbers (512 x 512) for a 35 µm x 35 µm image results in a theoretical pixel size of 68 nm 239 which is slightly higher than the smallest possible probe size of 50 nm but largely sufficient in 240 view of the measured actual lateral resolution for mapping subcellular structures in seaweed 241 cells.

245

  Figures 3i-iii show high resolution imaging of a L. digitata cell by NanoSIMS. Both, carbon

271

  small vacuoles are known to be acidocalcisomes containing polyphosphate and calcium and 272 are involved in osmoregulation.xx,xxi,xxii Although those acidocalcisomes have never been 273 identified in brown algae, they have been found in bacteria to human cells, which make them 274 likely to occur in L. digitata. Moreover, phosphorus was clearly detected in the phospholipid 275 containing cell membranes adjacent to the cell walls.276 277

  Figure 3iii (b and c) are a colour merge images of the 279 relative distribution of arsenic (white) and nitrogen (red) or phosphorus (blue). It further 280 demonstrates that the arsenic is mainly present in the cell wall with phosphorus outlining the 281 cell membranes and subcellular structures.

(

  speciation in L. digitata, it is not expected that all arsenic would be in the 284 cell membrane, since only a small proportion of arsenic is lipophilic (1.5%). Almost all arsenic 285 is in the form of either inorganic arsenic or in the form of hydrophilic arsenosugars (together 286 approx. 82%). No arsenic was detected inside the cells by NanoSIMS in either the cytosol, the 287 vacuoles or acidocalcisomes (their location was indicated by the phosphorus map). Hence, the 288 localisation of arsenic in L. digitata is significantly different to the accumulation of arsenic in 289 rice root cells. xxiii Inorganic arsenic in rice roots is bound mainly to phytochelatins while DMA 290 seemed unbound. xxiv In contrast brown seaweed with a large proportion of inorganic arsenic 291 Hijiki spp.) showed no significant amounts of arsenic phytochelatin complexes. xxv Hence, the 292 lack of arsenic in the cytosol and in particular in the vacuoles of L. digitata is not surprising.

318

  species need to come from biotransformation reactions either directly within the macroalgae or 319 from bacteria on the surface of the seaweeds. xxx If the arsenosugars are transported through the 320 cell membrane into the cell wall, then certain transporters which excrete those compounds need 321 to exist. Otherwise, the arsenosugars could be linked to the alginates which are produced in the 322 Golgi apparatus and then directly incorporated into the carbohydrate structure. They might bind 323 directly to the carbohydrates via the alcohol groups or the sulfonate (major arsenosugar in L.

Figure 3i . 340 Figure 3ii .

 3i3403ii Figure 3i. High resolution imaging of carbon (via C 2 -) (a) and nitrogen (via CN -) (b) in L.

Table 1 :

 1 and S1b (supplementary 174 material). Further evidence for the identification were the accurate masses and MS/MS 175 pattern as described by previously byRaab et al. xii , where the arsenolipid profile of a brown 176 seaweed Saccharina latissima was studied in detail. Total arsenic and hydrophilic arsenic species quantified; total arsenic (iAs),

	177
	178
	179
	182
			tAs	tAs hydro	iAs	DMA	As-Sugar
				(hydrophilic)			(sum)
			mg As/kg	mg As/kg	mg As/kg	mg As/kg	mg As/kg
		CRM CD-	55 ± 1.0	44 ± 0.67	0.090 ± 0.010	8.6 ± 0.3	32
		200 *
		L. digitata	117 ± 30	100 ± 22	62 ± 19	0.93 ± 0.11	37 ± 3.9
				(85.5%)	(53.0%)	(0.8%)	(31.6%)
	183	*certified reference value for tAs: 55 ± 4.0 mg/kg
	184
	185

180

inorganic arsenic determined as arsenate (iAs), Sum of all arsenosugars (As-Sugar). Fractions 181 are given as mass % of tAs. Errors are given as +/-SD for n=3.

Table 2 :

 2 Total arsenic and arsenolipids quantified; arsenic containing hydrocarbons (AsHC),

	186	arsenosugar containing phospholipids (AsPL).
			tAs lipid	AsHC	AsPL
			(lipophilic)	(sum)	(sum)
		as mg As/kg	1.7 ± 0.22	0.14	1.1
				(8.2%) a	(65%) a
		Mass % of
		total As	1.5%	0.12%	0.94%

Table 3 :

 3 identified As-lipids and their concentration using the ICPMS signal and compound-

	196	independent quantification (n=3).
		detected arsenolipid *	detected accurate mass **	µg As / kg algae dry mass
		(AsFA422)	C 22 H 36 AsO 3 (m/z 423)	4.2 ± 0.12
		(AsHC360)	C 19 H 42 AsO (m/z 361)	71 ± 6.7
		(AsHC374)	C 20 H 44 AsO (m/z 375)	23 ± 5.3
		(AsHC388	C 21 H 46 AsO (m/z 389)
		+ AsHC402)	+ C 22 H 48 AsO (m/z 403)	18 ± 2.4
		(AsPL930)	C 43 H 84 AsO 14 P (m/z 931)	109 ± 14
		(AsPL944)	C 44 H 86 AsO 14 P (m/z 945)	69 ± 11
		(AsPL958)	C 45 H 88 AsO 14 P (m/z 959)	636 ± 93
		(AsPL972)	C 46 H 90 AsO 14 P (m/z 973)	30 ± 4.4
		(AsPL986)	C 47 H 92 AsO 14 P (m/z 987)	97 ± 13
		(AsPL1000)	C 48 H 94 AsO 14 P (m/z 1001)	14 ± 1.7
		(AsPL1014)	C 49 H 96 AsO 14 P (m/z 1015)	119 ± 18
	197

*number in parentheses gives the non-protonated molecular mass 198 **For more information of the As-lipid structures see

Raab et al. 2013 vi

  such as Na + , K + , Ca 2+ , but also Zn 2+ and Cu 2+ . xiv Arsenic however in its 297 inorganic forms is either neutral (As(OH) 3 ) or anionic (HAsO 4 2-) and will not bind to the or fucoidans, which are abundant in the cell walls in brown seaweed. The fucoidans however 300 bind sulphate as sulfuric acid ester. If arsenic binds as arsenate to the carbohydrate in a similar 301 way is unclear. The resulting arsenic acid esters would be very unstable and will quickly 302 hydrolyse and form again inorganic arsenic. xxvi It has been noted that not only AsHC and AsPL Although the binding of arsenite or arsenate to the constituents of the cell wall and membrane 309 is unclear, the accumulation of arsenate into the cell wall is however mechanistically 310 conceivable. The arsenic accumulation into brown seaweed depends on the salinity and hence 311 also on the arsenic concentration in the seawater. xxviii Arsenate might directly accumulate from

	303	can be extracted, some unstable arsenic species have also been extracted with DCM/MeOH. A
	304	significant part of the extract eluted however in the void of a reverse phase HPLC method. This
	305	was likely to be attributed to unstable non-polar arsenolipids as shown before. xxvii Arsenic acid
	306	esters would quickly hydrolyse to inorganic arsenic during the sample preparation step.
	307
	308
	293
	294	It seems that most arsenic is located in the cell walls containing polysaccharides such as
	295	alginate, fucoidan and cellulose. Those alginates are also known to accumulate mono and

296

divalent cations 312 seawater or indirectly when it is taken up through phosphate transporters and then excreted as 313 arsenite. Efflux of arsenite is well known, while the efflux of arsenate has only recently been 314 suggested to occur as well. xxix 315 316

obtained with an Cs + plasma primary ion source for detection of negative secondary ions: 35 x