GO_3D_OBS-The Nankai Trough-inspired benchmark geomodel for seismic imaging methods assessment and next generation 3D surveys design (version 1.0)

Detailed reconstruction of deep crustal targets by seismic methods remains a long-standing challenge. One key to address this challenge is the joint development of new seismic acquisition systems and leading-edge processing techniques. In marine environments, controlled-source seismic surveys at regional scale are typically carried out with sparse arrays of ocean bottom seismometers (OBSs), which provide incomplete and down-sampled subsurface illumination. To assess and minimize the acquisition footprint in high-resolution imaging process such as full waveform inversion, realistic crustal-scale benchmark 5 models are clearly required.The deficiency of such models prompts us to build one and release it freely to the geophysical community. Here we introduce GO_3D_OBS a 3D high-resolution geomodel representing a subduction zone, inspired by the geology of the Nankai Trough. The 175 km × 100 km × 30 km model integrates complex geological structures with a visco-elastic isotropic parametrization. It is defined in form of a uniform Cartesian grid containing 33.6e degrees of freedom for a grid interval of 25 m. The size of the model raises significant high-performance computing challenges to tackle large-scale 10 forward propagation simulations and related inverse problems. We describe the workflow designed to implement all the model ingredients including 2D structural segments, their projection into the third dimension, stochastic components and physical parametrisation. Various wavefield simulations we present clearly reflect in the seismograms the structural complexity of the model and the footprint of different physical approximations. This benchmark model shall help to optimize the design of next generation 3D academic surveys in particular but not only long-offset OBS experiments to mitigate the acquisition footprint 15 during high-resolution imaging of the deep crust.


Introduction
To make a step change in our understanding of the geodynamical processes that shape the Earth's crust, we need to improve our ability to build 3D high-resolution multi-parameter models of deep geological targets at the regional scale. To meet this In the following sections we present step by step the workflow that was designed to create the seismologically representative 160 model. The whole procedure was implemented in MATLAB environment.

Geological features
The overall geological setup of our model is mainly (but not only) inspired by the features which were interpreted in the Nankai Trough area. However, these structures can be also found in different margins around the world combined in various configurations. Therefore, our model is not intended to replicate a particular subduction zone and its related geology for 165 geodynamic studies of the targeted region. On the contrary, it was design to have broad features one may encounter in these tectonic environments. Such quantitative subsurface description should challenge performances of high-resolution crustal-scale imaging from the complex waveforms generated in this model.
The skeleton of the initial 2D model is presented in Figure 2(a). It is designed as a cross-section aligned with the direction of subduction and consists of 46 geological units. These units are defined by means of the interfaces creating independent 170 polygons and obtained in the following three steps. First, we pick the spatial position of the nodes in the cross-section (see dots in Figure 2(a)). Second, we interpolate piecewise cubic Hermite polynomials through a subset of nodes to build a boundary (or interface) between the geological units. Such a boundaries represent either lithostratigraphic discontinuities or tectonic features like faults. They can connect and/or intersect each other to create closed regions in the (x,z) plane representing separate geological units. We refer these closed regions to as polygons. With such an approach we can modify the shape of existing or 175 add new features to the initial structure.
To create an overall skeleton of the model, we start from the interfaces that shape the mantle and the crust. In Figure 2(a), the top of the oceanic/continental crust and mantle are marked by solid red and green lines, respectively. The geometry of the subducting slab is inferred from previous seismic imaging results in the Tokai area. The oceanic crust exhibits variability of thickness and local bending consistently with what is proposed by geological studies in the same region (Kodaira et al., 2003;180 Le Pichon et al., 1996;Lallemand et al., 1992;Górszczyk et al., 2019) and addressed as an effect of subducting volcanic ridges and sea mounts. On the right flank of the model (starting at 140 km of model distance) we introduce a large volcanic ridge indicated by the uplift of the bathymetry and the significant thickening of the oceanic crust. Similar ridge (Zenisu Ridge) really exists in the eastern part of the Nankai Trough. It currently approaches the subduction zone and simultaneously undergoes the thrusting process (Lallemant et al., 1989). 185 In the next step we cut the subducting oceanic crust and mantle by series of interfaces (blue dashed lines in Figure 2(a)) to subdivide them into several polygons. These interfaces allow us to project each of the defined block independently in the third dimension (see next section). They are also designed to implement tectonic discontinuities in the oceanic crust and upper mantle. Such discontinuities acting as faults, thrusts or creating horst-graben structures are well documented in subduction zones (Tsuji et al., 2013;Ranero et al., 2005;Azuma et al., 2017;Boston et al., 2014;Vannucchi et al., 2012;von Huene et al., which can be typically found around strike slip faults zones (Huang and Liu, 2017;Ben-Zion and Sammis, 2003;Tsuji et al., 2014). Each of those boundaries can be also further used to implement fluid paths or damage zones.
On top of the oceanic crust we place two layers of subducting sediments (grey dashed and solid lines in Figure 2(a)). They extend starting from the subducting channel up to the right edge of the model. In general the lower layer is thicker and overlay 195 directly the top of the oceanic crust covering its thrusts and faults. The upper layer is thinner, whereas its upper limit delineates the decollement in our model. To provide more complexity, we introduce a small duplexing structure within those two layers (between 65 km and 85 km in Figure 2(a)) similarly to the decollement model presented by Kameda et al. (2017), Hashimoto and Kimura (1999) and Collot et al. (2011).
Just between the edge of the continental crust and the subducting slab, we add a stack of relatively thin sheets corresponding to 200 progressively underplated material (dashed black lines in Figure 2(a)). Typically these kind of structures occur when sediments (which move down on top of the subducting oceanic crust) create a duplex system and are further added to the accretionary wedge (Angiboust et al., 2014(Angiboust et al., , 2016Menant et al., 2018;Agard et al., 2018;Ducea and Chapman, 2018). Including these relatively fine scale structures in the deep part of the model shall significantly affect wavefield propagation.
The accretionary wedge consists of two parts. The land-ward part is formed by large deformed stacked thrusts sheets (solid 205 blue lines in Figure 2(a)) corresponding to the old accretionary prism which extends sea-ward into large out-of-sequence thrusts acting as a backstop now. Similar inner-wedge structures can be found in various subduction zones (Dessa et al., 2004;Raimbourg et al., 2014;Shiraishi et al., 2019;Collot et al., 2008;Cawood et al., 2009;Contreras-Reyes et al., 2010), as well as modelled with analogue sand-box experiments (Gutscher et al., 1998(Gutscher et al., , 1996. The outer-wedge (red shaded polygon in Figure   2(a)) contains the sedimentary prism adapted from Kington and Tobin (2011) and Kington (2012) (Figure 2(b)). It is built out 210 of different types of sediments which underwent deformation and created a sequence of thrusts. The complex-shaped, smallscale, steeply-dipping structures shall impose a challenge for seismic imaging. In our implementation the prism from Figure   2(b) is managed as a single block associated with the respective red shaded polygon in Figure 2(a). Depending on the changes of the shape of this polygon during projection step, we re-sample the prism from Figure 2(b) such that it conforms to the shape of the polygon.

215
Finally, the frontal prism is fed by the layers of incoming sediments (black lines in Figure 2(a)) -relatively thick in the area of trench and thin over the volcanic ridge. These sedimentary layers are associated to those interpreted in Figure 2(b). Similarly, on the land-ward part of the model we implement a sedimentary cover on the backstop including the fore-arc basin over the old part of the prism (between 40 km an 70 km). The whole model is covered with a water layer with strongly variable bathymetry.

220
The next step in our scheme is the geologically guided projection of the initial skeleton (Figure 2(a)) into the third dimension.
We design the set of projection functions which translate the initial inline skeleton in the crossline direction, such that the resulting structure follows concepts related to geology variations along strike direction. For each new inline section, the (x(y 0 ),z(y 0 )) coordinates of the initial node (dot in Figure 2(a)) are firstly projected in the crossline direction with an arbitrarily-chosen deterministic functions of the crossline distance y. In this way for a given inline i we obtain a new node with the coordinates 225 7 https://doi.org/10.5194/gmd-2020-240 Preprint. Discussion started: 6 October 2020 c Author(s) 2020. CC BY 4.0 License. (x(y i ),z(y i )). The general formula of the projection functions can be expressed as a parametric equation of the following form: It is immediate to see that x(y) and z(y) are quadric and linear functions respectively, which control curvature and dipping of a projected structure with the crossline distance y.

230
Once the nodes defining a given interface are projected in the crossline direction, we interpolate the piecewise cubic polynomials between the new nodes in the inline direction. In this way we create a new interface which is slightly shifted with respect to that of the previous inline section. The functions used to project each of the node from the sub-set of nodes defining a given interface are designed such that they are similar but not exactly the same. This means that the interface is not only shifted (as if exactly the same projection function was used for all the nodes defining the interface) but its geometry can be also modified. In It is worth mentioning that some nodes can belong to several interfaces -for example, when they are at the junction of two interfaces. In such a case, these nodes are projected only once with the function assigned to the first interface in order to 245 guarantee a conform space warp of the polygons. Therefore, the projection function applied to the second interface must be partially determined by the fact that one of the nodes (which belongs also to the first interface) was already projected. In such a way, the structural skeleton we create can be viewed as a system of interfaces which are implicitly coupled via the dependencies between the position of the nodes which define them. Therefore, adding new features to the model or modifying existing ones requires careful design of the corresponding projection functions which will honour the previous dependencies between 250 interfaces.
We assumed that the model will be 100-km width, therefore we generate 4001 2D inline sections spaced 25 m apart. To highlight the lateral variability of the structure generated during projection, we show in a perspective view the skeleton of the inline sections extracted every 20 km (Figure 3(b)). First, we design the curved shape of the subduction front. With increasing crossline distance, the oceanic crust and mantle shrinks-back. As a result, the continental part of the mantle and crust elongates 255 seaward. Simultaneously, we compress and extend the respective polygons making their relative size variable. One can track this variability by following the geometry of the interfaces cutting through the oceanic mantle and crust along successive inline sections (blue dashed lines in Figure 3(b)). The absolute depth of the oceanic crust is not only changing along the subduction trend but also gently increasing with crossline distance. Simultaneously, there is an up-growing of the underplating sediments below the edge of the continental crust (black dashed lines in Figure 3(b)) -consistent with the thickening of the subducting 260 sediments as the crossline distance increases. The position of the backstop is also changing from one inline to another. This partially results from the curved shape of the subduction front and was inspired by the 3D interpretation of Tsuji et al. (2017).
At the same time, the length of the inner wedge (solid blue lines) increases. This is coupled with the differences in scale and deformation of the thrusts sheets which are building this unit, as well as with the shape of the forearc basin covering it. Further seaward, the outer wedge (red shaded polygon) becomes shorter and thicker as the crossline distance is increasing. The right 265 flank of the model is marked by the uplift of the bathymetry related to the presence of the incoming volcanic ridge which becomes smaller as the crossline distance increases.
In Figure 3(c-f), the crossline sections extracted at 40 km, 70 km, 100 km, 130 km (vertical red dashed lines in Figure 2(a)) highlight the structural heterogeneity in the crossline direction. This heterogeneity is defined at different scales. Firstly, due to the overall curved shape of the subducting front. Secondly, due to the different projection functions causing warping of the 270 polygons. The structural complexity will be further increased through parametrization and application of stochastic components.

Index and gradient matrix
Once the inline section has been projected in the crossline direction, the newly generated structural skeletons are further 275 mapped into uniform 2D Cartesian grids of dimension 1201 × 7001 (grid size 25 m × 25 m). During this step, we assign an arbitrary index to the grid points representing different geological units. Namely, the grid points inside each of the polygons are now described by a unique integer value. By re-sampling in the horizontal and vertical directions the matrix discretizing the accretionary block shown in Figure 2(b), we reshape the geometry of the outer wedge such that it matches that of the corresponding polygon (red shaded area in Figure 2(a)) in each projected inline section. The deformed sediment layers inside 280 the outer wedge have indexes consistent with those of the sediments layers incoming into the trench since they represent the same geological formations. Filling all the polygons with the corresponding indexes leads to the matrix plotted in Figure 4(a).
We pick random colours to distinguish between the neighbouring units. The matrix contains 46 geological blocks defined by the unique integer index. Such indexation allows us to easily refer to a particular structural unit of the model during the further steps -for example to modify the parametrization.

285
On top of the index matrix, we implement another matrix, referred to as gradient matrix, which allows to introduce spatial variations of the physical parameters in each large-scale unit and around main fault planes (Figure 4(b)). Without this second matrix, the physical parameters in each geological unit would have constant value (related to the integer from matrix of indexes). The spatial variation of the parameters within the same unit can be related to increasing depth in the mantle, layering of the crust, low-velocity zones in the subducting sediments, compaction in the prism or damage zones around the faults etc.

295
In practice, the index matrix and the gradient matrix are combined to implement any physical parameter in each geological formation. To do so, we extract a given geological unit (based on the index matrix) and assign to it a constant parameter value p α . We also extract the corresponding part of the gradient matrix and denote it by G n (normalized gradient matrix). The gradient function associated with G n can be manipulated (scaled, clipped, biased etc.) leading to the matrix G which provides a desired spatial variation of a parameter in a geological unit. The overall formula to compute the physical parameters inside a 300 given model unit is as follow: where p(x, z) is the parameter entry at the position (x, z), while p β scales the entry g(x, z) of G with appropriate physical units.
To better illustrate how the index and gradient matrix can be used to produce a physical model, we consider the following 305 example. Let us imagine that we want to assign V p values inside the part of the oceanic crust delineated by the thick red line in Figure 4(a-b). The integer index assigned to this unit is 11 (Figure 4(a)). We set p α to 4800 m/s and p β to 2900 m/s. Taking into account that the values of G n for this unit span from 0.0 to 1.0 (from the top to the bottom of the crust) the final V p values will be increasing linearly from 4800 m/s to 7700 m/s. One may also desire to implement a gradient discontinuity in the oceanic crust corresponding to the oceanic layer 2-layer 3 boundary. Setting the thickness ratio to 0.3 and 0.7 in the layer 2 and layer 310 3, respectively, would split G n in unit 11 into two sub-matrices given by: As a result we obtain two gradient matrices G 1 and G 2 with normalised values between 0.0 and 1.0. We can now set the p α to 4800 m/s and p β to 1800 m/s for G 1 . As a result we obtain rapid increase of V p between 4800 m/s to 6600 m/s inside layer 2.
For G 2 we use p α equal to 6800 m/s and p β equal to 1000 m/s leading to V p values inside layer 3 ranging between 6800 m/s 315 and 7800 m/s. The resulting boundary between layer 2 and layer 3 is marked with dashed line in Figure 4(c).
We also use the gradient matrix to implement variations of the parameters within the subducting sediments, damage zones and subducting channel which result from a small-scale geological processes (fluid migrations, lithological variations, tectonic compaction, mass wasting etc., (Agudelo et al., 2009;Park et al., 2010). Moreover, we implement lateral parameter discontinuities at the interfaces acting as faults. For example, the stairs in the Moho resulting from the interfaces cutting through the 320 oceanic mantle and crust imply that the gradients on both sides of those interfaces are slightly different. In consequence we create a lateral jump in the parameter values on both side of the fault. This jump is in practice more pronounced near to the top of the mantle (∼5 km to ∼7 km below Moho) and is vanishing with increasing depth to completely disappear at the bottom of the model.
Finally, the gradient matrix describes also two additional types of perturbations. The first type is implemented around the re-325 gions of thickened oceanic crust mimicking subducting volcanic rides (Park et al., 2004;Kodaira, 2000). They are marked by the gradient and velocity variations in Figure 4(b-c) (95 km -105 km and 145 km -155 km). The second type of perturbations is intended to introduce more heterogeneity around the fault planes such as damage zones or fluid paths. Parameter variations within such zones affect the kinematic and dynamic characteristic of the propagating wavefield. They also have potential to generate distinct arrivals -so-called trapped waves (Li and Malin, 2008;Ben-Zion and Sammis, 2003). We implement such 330 anomalies around the faults in the oceanic crust and upper mantle, the major faults at the edge of continental crust, as well as within the large thrusts sheets between inner and outer accretionary wedge. Few small perturbations are also locally added between layers of inner wedge to increase the complexity of this unit.

Physical parameters 335
From the index and gradient matrices described in the previous section, we can now implement the seismic properties in the structural units. In practice, this requires some constraints on the parameters, which might come from field experiments and laboratory measurements. Regional seismic studies often lack resolution to directly map the reconstructed properties into the fine-scale structures of our model. Moreover, they often ignore second-order parameters such as s-wave velocity (V s ), density (ρ), attenuation (Q p and Q p ) and anisotropy to focus on the p-wave velocity (V p ) reconstruction. Therefore, we first develop 340 a V p model from recent FWI case studies in the Nankai trough. Then, we use this V p model as a reference to build the other parameters -namely V s , ρ, Q p , Q s -from empirical relations.
Reconstruction of V p from wide-angle seismic data has long historical records (Christensen and Mooney, 1995;Mooney et al., 1998). Traveltime tomography has produced a vast catalogue of smooth V p models from different regions around the world.
This information (even though values can significantly vary depending on the studied area) is useful to constrain the velocity 345 trend in the main large-scale geological units (such as mantle, continental and oceanic crust). Moreover, recent OBS FWI case studies partially fill this resolution gap (Kamei et al., 2012;Górszczyk et al., 2017) and allow us to refine V p in the short-scale units of our model (Figure 4(c)). We perform this refining of the V p velocities by trial-and-error, until the OBS gathers simulated under acoustic approximation in the V p model exhibit in overall a similar anatomy to the OBS gathers collected in the eastern-Nankai trough (Górszczyk et al., 2017).

350
From the V p model, we build V s and ρ using the empirical polynomial relations of Brocher (2005). These relations were inferred from compiled data on both laboratory measurements, logs, vertical seismic profiling and tomographic studies. They are relevant for V p ranging from 1500 m/s to 8500 m/s, and hence represent the average behaviour of crustal rocks over the depth range covered by our model. The resulting V s and ρ models are presented in Figure 5(a-b). Shear-wave velocities in the shallow sediments are as small as ∼ 530 m/s. Although even lower velocities can be found in the first few meters below seabed, 355 these values already impose significant challenges for wavefield modelling. In the oceanic crust, V s increases from 3200 m/s to 4300 m/s, while it starts at 4600 m/s on top of the upper mantle and gently increases with depth. Density of the oceanic and continental crust varies from 2600 kg/m 3 to 3200 kg/m 3 and from 2500 kg/m 3 to 2900 kg/m 3 respectively, while in the upper mantle ρ starts around 3200 kg/m 3 . Figure 5c shows V p /V s ratio, which varies from 1.6 to 3.0. Initially, the V s model obtained from empirical relations of Brocher (2005) was not producing heterogeneous V p /V s ratio in the subduction channel -as it can 360 occur when fluid overpressure, fluid diffusion or hydro-geologically isolated zones generate variation of the seismic properties in the subduction channel (Kodaira et al., 2002;Collot et al., 2008;Ribodetti et al., 2011). Therefore, we additionally rescale the velocities in this unit using the index and gradient matrices.
We also implement attenuation effects in our model through the Q p and Q s parameters. The choice of consistent approach to constrain plausible Q p and Q s models is quite challenging. This is firstly because the various types of attenuation (intrinsic 365 and scattering-related) can be controlled by countless factors -rock type and mineralogy, porosity, fluid content, saturation etc. Secondly, values of Q factor can differ significantly between studies depending on the used methodology, scale of the attenuation or frequency content of seismic data. Thirdly, accurate reconstruction of Q models from field experiments is not so well documented because the footprint of attenuation in the wavefield remains small in particular in the deep crust. Therefore, to build the Q p and Q s models, we combine few empirical relations between velocity (V p , V s , V p /V s ) and attenuation (Q p , Q s ) 370 which we find consistent for our synthetic model. To build the Q s model we use the following power law Q s =0.0053V 1.25 s (Wiggens et al., 1978). This empirical relation is in good agreement with the Q s estimation of (Olsen et al., 2003) performed in the Los Angeles Basin. The estimated Q s values range from ∼25 in the shallowest sediments to ∼220 in the deep mantle (see Figure 5(e)). To avoid constant ratio between Q p and Q s (Q p =1.5Q s in Olsen et al. (2003)), we generated Q p model ( Figure 5(d)) via another empirical relation Q p =36.8(V p )-49.6 -derived from well-log data by Zhang and Stewart (2008). This 375 formula, even though derived from velocities up to ∼3000 m/s using sonic-log frequency band, produces reasonable Q p variations between ∼35 and ∼260 in our model, and leads to Q p /Q s ratio ranging from ∼1.1 to ∼1.5 ( Figure 5 Additionally, the attenuation of high-frequency wavefield component through the scattering effects is introduced in our model by the means of small scale stochastic perturbations as described in the following section.

Stochastic components 385
In the final step of our model-building procedure, we consider various types of stochastic components, which are intended to introduce perturbations of random character into the model. We design two types of such a random components: (i) small scale parameter perturbations, and (ii) spatial warping of the final 3D cube.
In Figure 4(c) and Figure 5, the parameters are varying smoothly within the units according to the applied gradient functions.
Indeed, such smooth variations are rather idealised vision of the true earth. Therefore, to make it more realistic, we create an another matrix describing small scale perturbations (Figure 6(a)) which shall have a second order impact on the wavefield propagation. Using this matrix, we can scale the models after parametrization and obtain the updated model presented in Figure   6(b) (compare with the smooth velocity variations of the V p model shown in Figure 4(c)). The matrix in Figure 6(a) is obtained 395 by stacking 3D disk-shaped structural elements (SE) (Figure 6(c-f)). The values in the SE (Figure 6(g)) vary from 0 at the edge to 1 at the centre and are further translated into parameter perturbations. We can control the position and size of these SE, as well as their correlation lengths (z,x,y,) -depending on a desired characteristic of the final perturbations.
In Figure 6(c-f), we show four different 3D stacks (10 km × 10 km × 10 km sub-volumes are displayed), corresponding to four different scales of the SE with correlation lengths ratio equal to 0.25 × 1 × 1. The red/blue colours correspond to the 400 positive/negative amplitudes, which are randomly assigned to each of the SE during stacking. These random variations of positive/negative amplitudes guarantee that the final distribution of the stochastic perturbations for each geological unit in the matrix shown in Figure 6(a) has a zero mean. The position of the SE within the stack is controlled by the Cartesian coordinates.
The distances between the centres of neighbouring SE equal (z/2+r z ,x/2+r x ,y/2+r y ), where r z , r x , r y are the spatial shifts randomly picked from the intervals (-z/4,z/4), (-x/4,x/4), (-y/4,y/4). In other words, the neighbouring SE strongly overlap 405 with each other -as can be seen in Figure 6(c).
Looking at the stochastic matrix in Figure 6(a), one can observe that the scale, shape and direction of the perturbations vary from one geological unit to another. For example, perturbations in the prism are finer than those in the crust or mantle. Also their shape in the oceanic crust and prism are more anisotropic than the shape of those implemented in the mantle. In practice, this is implemented by summation of 3D stacks containing different SE. For example, we create stochastic perturbations in the 410 sedimentary layer and the outer prism with only the SE shown in Figure 6(e-f); Perturbations in the inner wedge and underplated units additionally incorporate SE shown in Figure 6(d). Finally, the oceanic crust contains SE at all scales ( Figure 6(c-f)).
This also applies to the mantle and continental crust, however for these units we use SE with different correlation lengths ratio equal to 0.5 × 1 × 1.
In Figure 6(a), we also show that the stochastic perturbations in some of the units are oriented according to the dipping or 415 bending of the structure. Since we know the polynomials defining each of the interfaces in the model (and therefore their position), we shift vertically the columns of the matrix containing stochastic perturbations of a given geological unit such that they follow the shape of the structure. For example, the perturbations which fill the oceanic crust follow the smooth trend generated from the polynomial approximation of the interfaces creating the top of the oceanic crust. Note, that because we stack truly 3D SE within the 3D cube of the same size as the 3D model, our stochastic perturbations follow continuously the geological 420 structures not only along the inline direction but also along the crossline direction.
To analyse what the designed stochastic perturbations mean in terms of the spatial resolution of the model, we compute the 2D wavenumber spectrum (displayed in logarithmic amplitude scale in Figure 6(h-i)) of the V p model with and without stochastic 13 https://doi.org/10.5194/gmd-2020-240 Preprint. Discussion started: 6 October 2020 c Author(s) 2020. CC BY 4.0 License.
components. As shown in Figure 6(h), the spectral amplitudes of the model without stochastic perturbations are focused around the low wavenumber part of the spectrum, which represent structures the size of which is of the order of ∼500 m and larger.

425
In contrast, the high-wavenumber amplitudes are clearly magnified in the spectrum of the model with stochastic components (Figure 6(i)). The overall distribution of the energy added to the background medium (with respect to the spectrum shown in Figure 6(h)) is close to normal, which highlights the randomness of the perturbations. The designed approach led to perturbations of smaller size being in fact of the order of the grid size (25 m). On the other hand, one can also observe in Figure   6(i) a few dipping bands of increased amplitudes. They correspond to the anisotropic shape of the SE following the geological 430 structures -and therefore indicate that our stochastic perturbations are partially predefined and not purely random.

Structure warping
Finally, we apply a second type of stochastic components to the model after projection, parametrization and application of small scale stochastic perturbations. The aim of this step is to perform a small-scale point-wise warping of the input 3D grid 435 using a distribution of random shifts. By small-scale we mean that the structural changes introduced by warping are in general much smaller than those incorporated during the projection step. Their impact on the shape of the wavefield will therefore be weaker than the one induced by the geological structure itself, however more pronounced than that resulting from the small scale perturbations described in previous section. The idea was inspired by Dynamic Image Warping (DIW) (Hale, 2013). DIW is a technique which allows for the estimation of local shifts between two images (2D matrices), assuming that one of the image 440 is a warped version of the other image. This problem is cast as an inverse problem where the unknowns are the shifts that will allow for matching between the two images.
In our case, we use DIW in a forward sense. This means we build a matrix of smoothly varying random spatial shifts of a desired scale and magnitude. Then, we warp each of the inline section from the 3D model grid using those shifts. Figure 7(a-b) shows two 300 km × 300 km matrices that represent distribution of shifts for two different spatial scales. While the scale of 445 the shifts in Figure 7(a) is local (starting from ∼10 km), the spatial extension of the shifts in Figure 7(b) is regional. The blue and red colours indicate the negative and positive shifts, respectively -scaled between ±2km. We take the average of the two distributions to obtain the final distribution shown in Figure 7(c). Our procedure of warping is implemented as follow.
First, for each of the 4001 inline sections, we extract from the distribution shown in Figure 7(c) a 2D sub-matrix with the size of the inline section. for negative (blue colour) and positive (red colour) shifts, respectively. For simplicity, we consider only vertical shifts in this example. However, we also implement warping in the horizontal direction, i.e. m w (x, z)=m i (x + s, z), in the same way as we do in the vertical direction.
Once an inline section has been warped, we extract a sub-matrix of shifts for the next inline section from the distribution of shifts shown in Figure 7(c). We perform this extraction at a slightly modified position with respect to the previous one, typically 460 one grid point further (i.e. 25 m). In this way we maintain the continuity of the warped structures along the crossline direction.
Although this inline by inline procedure is based on the 2D warping, in practice it introduces fully 3D structural changes in the model, at the same time remaining more practical for quality control of the results and computational efficiency. Three depth slices extracted at 5 km, 7.5 km and 10 km depth from the 3D V p models before and after warping highlight how the warping has modified the shape of the structures (Figure 7(e-g)). Importantly, the continuity of the structure is preserved 465 during warping not only along the inline direction but also along the crossline direction. It is also worth remembering that, in this example, the distribution of shifts is random. However, nothing prevent to use a deterministic distribution, for example to produce uplift or bending of certain parts of the model.

Wavefield modelling
In this section, we perform full wavefield modelling in a target of the GO_3D_OBS model to assess the footprint of the 470 structural complexity and the heterogeneity of the physical parameters on the anatomy of these wavefields. Figure  using the TOYXDAC visco-acoustic finite-difference (2 nd order in time and 4 th order in space) and SEM46 visco-elastic spectral-element FWI codes developed within the framework of the SEISCOPE consortium (https://seiscope2.osug.fr).

Acoustic versus visco-acoustic 2D wavefields
The first modelling test demonstrates the influence of the viscous effects on 2D acoustic P-wave modelling. We compare 2D wavefields that are computed without (Q p =10000) and with attenuation. The intrinsic attenuation mechanism we use is based 485 on the generalized Maxwell body including 3 standard linear solid attenuation mechanisms (Yang et al., 2016). For the sake of computational costs, all the models are re-sampled in a uniform Cartesian grid with a 50 m grid interval. Figure 9(a-b) shows the OBS gathers generated with constant and variable Q p models. True-amplitude seismograms are plotted in both panels with the same amplitude scale. Although the two gathers look almost identical, the difference between them -shown in Figure 9(c) -clearly shows the variations of the signal amplitudes. Changes are mostly visible in the first arrivals up to 15 km offset (corresponding to the waves travelling in the shallow parts of the model), as well as in the energetic reflection from the top of the subducting oceanic crust. This reflection is affected by the high attenuation layer within the subduction channel.
Importantly, the variable Q p model introduced not only differences in terms of wavefield amplitude, but also notable phase shifts due to dispersion effects. The inset in Figure 9(c) shows the zoom on interleaved traces extracted within the dashed 495 rectangles in Figure 9(a-b). Arrivals in the blue-shaded traces (data modelled with Q p =10000) are clearly shifted in time with respect to those simulated with variable Q p model. Therefore, incorporating attenuation model during seismic imaging can not only improve the amplitude reconstruction, but it can also have a second order impact on the kinematic correctness of the final image.

2D versus 3D visco-acoustic wavefields 500
We assess now the footprint of 3D effects by comparing the synthetic seismograms computed by 2D and 3D wave modelling.
The data resulting from both simulations are shown in Figure 10 To better understand how the wavefield propagates within the 3D target during the simulation, we extract snapshots every 2.5 s at 10 km depth (see black dotted line in Figure 8 for the depth-slice location). For comparison, we also perform 3D modelling in the 2.5D version of the inline section number 2601. Depth-slices extracted from the V p 3D and 2.5D models are presented in Figure 11(a). The corresponding snapshots are displayed in Figure 11(b-l). The lower-half of each panel (the 515 blue-shaded phases) presents the respective snapshot from the wavefield modelling performed within 2.5D model. It is clear that with increasing time and offset-distance the wavefield simulated in the 3D target becomes more and more complex -as compared to the 2.5D counterpart. In particular, while the blue-shaded wavefield is symmetric with respect to the shooting profile, the wavefield computed in the 3D target shows evidence of significant off-plane propagation. This is indicated both by the slant appearance of the first arrival wave (Figure 11(d-f)), as well as the shape and the spatial location of the reflections.

520
In Figure 11(m) we show analogous comparison as in Figure 10(c) between 3D and 2.5D seismograms. Although the 2.5D data fit better than pure 2D data to their 3D counterparts (probably due to the consistent point-source implementation for 2.5D and 3D modelling) the phase shifts between far-offset traces can still be observed. There is also visible mismatch in therms of amplitude and phase between the traces from the inset containing complex reflections package.
Note, that although our 3D target contains part of the curved subduction front, the acquisition profile is still mainly aligned 525 with the dip direction of the structures within the wedge. Nevertheless, the off-plane wavefield propagation is yet remarkable.
One can imagine that this effect can be further magnified when the geological setting contains more heterogeneities along the crossline direction.

Acoustic versus elastic 3D wavefields
The final modelling example which we present here shows the difference between 3D spectral-element acoustic (in the whole  Figure 12(a-b) shows the two OBS gathers inferred from acoustic and elastic wavefield modelling.
One can observe that the package of wide-angle reflections between 20 and 60 km in Figure 12(a) contains significantly more energetic arrivals as compared to the same arrivals in Figure 12(b). Some of them are difficult to track (or simply not present) in the gather generated using elastic wavefield propagation. This loss of energy can be explained by the P to S conversions.
In contrast, we see in the elastic seismograms some weak arrivals of low apparent velocity between 10 and 20 km, indicating 540 some elastic effects (Figure 12(b)).
The analogous comparison as in Figure 11(c) is presented in Figure 12(c). One can observe that the waveforms from the acoustics modelling (blue-shaded phases) do not create continuous arrivals with their elastic counterpart delineating differences between two types of data.
It is also worth to mention here, that the data simulated with the finite-difference scheme exhibit some evidence of so called 545 "stair-case effect" (see diffraction-like patterns at the short-offset data in Figure 10(a-b)). These artefacts result from the bathymetry projection on the Cartesian grid (50 m grid step). In contrast, the spectral-element mesh can handle accurately the sharp seabed interface which is reflected by lack of similar artefacts in the gathers from Figure 12.

550
Our previous studies oriented on crustal-scale seismic imaging from OBS data were located in the Tokai area of the Nankai Trough (Górszczyk et al., 2017(Górszczyk et al., , 2019. This gave us a true-earth reference to guide the building of a first version of our model.
However, our primary aim was to build a structural model which involves as many as possible geological features observed in subduction zones rather than following rigorously the Nankai Trough geological setting. On this basis, the GO_3D_OBS shall be seen as a generic crustal-scale seismic imaging benchmark model rather than a seismological reference to the Tokai segment. We hope that, based on some feedback from the community, we will be able to incorporate more geological features or modify existing ones.

Acquisition design
In the academic community of marine geosciences, sparse OBS acquisitions remain the primary tool for offshore lithospheric imaging. The acquisition-related issues focus mainly on the shooting strategy and the optimal OBS spacing (Brenders and 560 Pratt, 2007a). This is the consequence of the fact that the number of available instruments is limited ( ∼ 100 maximum). The spacing between OBS must not be chosen at the expense of the maximum offset sampled by the acquisition as this maximum offset should be enough to record diving waves that reach the deepest structures (namely the upper mantle). On the other hand, the sparse acquisition design shall minimise model distortions resulting from an incomplete structure illumination associated with limited coverage of the surface acquisition. These distortions need to be assessed since they can generate significant bias 565 during geological interpretation.
To investigate the footprint of various acquisition setting on the imaged target, realistic synthetic models are necessary for the survey design. Processing of different datasets generated in such models would allow for optimization of the acquisition setup to find the best compromise between deployment effort, image resolution and target sampling. This was our main motivation beyond the design of the GO_3D_OBS model. One of the ongoing projects now focuses on the investigations of the off-570 plane wavefield propagation during 2D crustal-scale FWI. From the modelling tests shown in this paper (Figure 10), we can conclude that these effects may have a detrimental impact on the imaging results. This shall stimulate the extension of routinely performed 2D OBS surveys toward the optimally designed future 3D deployments. We can evoke here the 2001 Seize France Japan (SFJ) 2D OBS experiment (Operto et al., 2006) which was conducted in the Tokai area of the Nankai Trough using 100 OBS deployed with a dense 1 km spacing. Today, 3D sparse OBS deployments are conducted with a similar number of 575 stations with an aim to apply FWI. Therefore which of those two different acquisition setups are closer to the optimal setting?
Was the 2D SFJ profile oversampled or maybe the today's 3D experiments are undersampled? We believe that for a new OBS experiments we shall look for justification of a proposed acquisition geometry rather than use an ad hoc configuration.

Usefulness for tomography
Although the processing method we pair the most with our model is FWI, there is no limits in terms of application of other 580 tomographic techniques, for example ray-theory based modelling techniques such as ray-tracing or eikonal solvers. While FAT is typically producing rather smooth velocity models, the other variants of tomography have potential to build models with improved resolution. Evaluation of those tomographic methods against GO_3D_OBS model can lead to their further development. Traveltime tomography techniques are in general computationally much less intensive than waveform inversion methods. It is therefore worth to push their efficiency to the resolution-limits and validate them with complex benchmarks.

585
This kind of development of tomographic methods indirectly contributes also to the FWI popularisation since the tomographic models are usually used as initial models for FWI.

Uncertainty estimation
Together with development of different imaging techniques, methods for quantitative estimation of their uncertainty and resolution limits shall be proposed. Indeed, this subject is not trivial when facing real-data processing and therefore the issue is 590 often skipped by authors. The obvious problem is that the true geological structure is never exactly known, and therefore the model of this structure -derived under given approximation -has no direct reference point to compare. On the other hand, the resolution and accuracy of the seismic imaging methods are always limited. Therefore, it shall bring us to reflection about the potential pitfalls related to over-interpretations during real-data case studies. Evaluating the uncertainty of the inversion methods with realistic synthetic tests can bring us closer to the estimation of the error during real-data processing. Not only 595 because we exactly know the underlying structure, but also because of the possibility to investigate the influence of different approximations, experiment geometries or noise.

Further development
At the current stage of development, we believe that GO_3D_OBS is sufficiently realistic to explore the validity of various tomographic and inversion methods on synthetic datasets with different acquisition designs. Future branches of development 600 of this benchmark can cover different aspects including both structural components and parametrisation. For instance, one of the interesting topics is extension toward anisotropy. Building such an anisotropy model consistent with geological studies of deep crust could be challenging due to general difficulty of precise anisotropy estimation. Nevertheless, it would significantly increase authenticity of the corresponding seismic data. On the other hand, generating accompanying dataset for the current models can further broaden its impact. Possible dataset could include sparse 3D and dense 2D OBS deployments, as well as 605 corresponding streamer data. Such dataset could be directly use for a given type of processing. In particular, generating the seismic data without and with realistic type and level of noise would help to further evaluate the impact of noise on the accuracy of the seismic imaging methods.