Hicham Badri
email: hicham.badri@inria.fr

Hussein Yahia
email: hussein.yahia@inria.fr

Khalid Daoudi
email: khalid.daoudi@inria.fr

Fast and Accurate Texture Recognition with Multilayer Convolution and Multifractal Analysis

A fast and accurate texture recognition system is presented. The new approach consists in extracting locally and globally invariant representations. The locally invariant representation is built on a multi-resolution convolutional network with a local pooling operator to improve robustness to local orientation and scale changes. This representation is mapped into a globally invariant descriptor using multifractal analysis. We propose a new multifractal descriptor that captures rich texture information and is mathematically invariant to various complex transformations. In addition, two more techniques are presented to further improve the robustness of our system. The first technique consists in combining the generative PCA classifier with multiclass SVMs. The second technique consists of two simple strategies to boost classification results by synthetically augmenting the training set. Experiments show that the proposed solution outperforms existing methods on three challenging public benchmark datasets, while being computationally efficient.

Introduction

Texture classification is one of the most challenging computer vision and pattern recognition problems. A powerful texture descriptor should be invariant to scale, illumination, occlusions, perspective/affine transformations and even non-rigid surface deformations, while being computationally efficient. Modeling textures via statistics of spatial local textons is probably the most popular approach to build a texture classification system [START_REF] Lazebnik | A sparse texture representation using local affine regions[END_REF][START_REF] Zhang | Local features and kernels for classification of texture and object categories: A comprehensive study[END_REF][START_REF] Varma | A statistical approach to material classification using image patch exemplars[END_REF][START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF][START_REF] Crosier | Texture classification with a dictionary of basic image features[END_REF][START_REF] Liu | Texture classification from random features[END_REF][START_REF] Liu | Sorted random projections for robust texture classification[END_REF]. Based on this Bag-of-Words architecture, these methods try to design a robust local descriptor. Distributions over these textons are then compared using a proper distance and a nearest neighbor or kernel SVMs classifier [START_REF] Scholkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond[END_REF]. Another alternative to regular histograms consists in using multifractal analysis [START_REF] Varma | Locally invariant fractal features for statistical texture classification[END_REF][START_REF] Xu | A projective invariant for textures[END_REF][START_REF] Xu | Combining powerful local and global statistics for texture description[END_REF][START_REF] Xu | A new texture descriptor using multifractal analysis in multi-orientation wavelet pyramid[END_REF][START_REF] Ji | Wavelet domain multifractal analysis for static and dynamic texture classification[END_REF]. The VG-fractal method [START_REF] Varma | Locally invariant fractal features for statistical texture classification[END_REF] statistically represents the textures with the full PDF of the local fractal dimensions or lengths, while the methods in [START_REF] Xu | A projective invariant for textures[END_REF][START_REF] Xu | Combining powerful local and global statistics for texture description[END_REF][START_REF] Xu | A new texture descriptor using multifractal analysis in multi-orientation wavelet pyramid[END_REF][START_REF] Ji | Wavelet domain multifractal analysis for static and dynamic texture classification[END_REF]] make use of the boxcounting method to estimate the multifractal spectrum. Multifractal-based descriptors are theoretically globally invariant to bi-Lipschitz transforms that include perspective transforms and texture deformations. A different approach recently presented in [START_REF] Sifre | Rotation, scaling and deformation invariant scattering for texture discrimination[END_REF] consists in building a powerful local descriptor by cascading wavelet scattering transformations of image patches and using a generative PCA classifier [START_REF] Bruna | Invariant scattering convolution networks[END_REF]. Unfortunately, while these methods achieve high accuracy on some standard benchmark datasets, little attention is given to the computational efficiency, which is crucial in a real-world system.

We present in this paper a new texture classification system which is both accurate and computationally efficient. The motivation behind the proposed work comes from the success of multifractal analysis [START_REF] Xu | A projective invariant for textures[END_REF][START_REF] Varma | Locally invariant fractal features for statistical texture classification[END_REF][START_REF] Xu | Combining powerful local and global statistics for texture description[END_REF][START_REF] Xu | A new texture descriptor using multifractal analysis in multi-orientation wavelet pyramid[END_REF][START_REF] Ji | Wavelet domain multifractal analysis for static and dynamic texture classification[END_REF]. Given an input texture, the image is filtered with a small filter bank for various filter orientations. A pooling operator is then applied to improve robustness to local orientation change. This process is repeated for different resolutions for a richer representation. This first step generates various low-pass and high-pass responses that form a locally invariant representation. The mapping towards the final descriptor is done via multifractal analysis. It is well known that the multifractal spectrum encodes rich texture information. The methods in [START_REF] Xu | A projective invariant for textures[END_REF][START_REF] Xu | Combining powerful local and global statistics for texture description[END_REF][START_REF] Xu | A new texture descriptor using multifractal analysis in multi-orientation wavelet pyramid[END_REF][START_REF] Ji | Wavelet domain multifractal analysis for static and dynamic texture classification[END_REF] use the box-counting method to estimate the multifractal spectrum. However, this method is unstable due the limited resolution of real-world images. We present a new multifractal descriptor that is more stable and improves invariance to bi-Lipschitz transformations. This improvement is validated by extensive experiments on public benchmark datasets. The second part of our work concerns training strategies to improve classification rates. We propose to combine the generative PCA classifier [START_REF] Sifre | Rotation, scaling and deformation invariant scattering for texture discrimination[END_REF][START_REF] Bruna | Invariant scattering convolution networks[END_REF] with kernel SVMs [START_REF] Scholkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond[END_REF] for classification. We also introduce two strategies called "synthetic training" to artificially add more training data based on illumination and scale change. Results outperforming the state-of-the-art are obtained over challenging public datasets, with high computational efficiency.

The paper is organized as follows : section 2 describes the proposed descriptor, section 3 presents the proposed training strategies, section 4 presents classification results conducted on 3 public datasets as well as a comparison with 9 state-of-the-art methods.

Robust Invariant Texture Representation

The main goal of a texture recognition system is to build an invariant representation, a mapping which reduces the large intra-class variability. This is a very challenging problem because the invariance must include various complex transformations such as translation, rotation, occlusion, illumination change, non-rigid deformations, perspective view, among others. As a result, two similar textures with different transformation parameters must have similar descriptors. An example is given in Figure 1. Not only the system should be accurate, but it should be also computationally efficient. Otherwise, its use in a real-world system would be limited due to the long processing time to extract the descriptor. Our goal in this paper is to build both an accurate and fast texture recognition system. Our Matlab non-optimized implementation takes around 0.7 second to extract the descriptor on a medium size image (480 × 640) using a modern laptop. The processing time can be further decreased by reducing the resolution of the image without sacrificing much the accuracy. This is due to the strong robustness of our descriptor to scale changes via accurate multifractal statistics that encode rich multi-scale texture information. We explain in this section how we build the proposed descriptor, the motivation behind the approach and the connection with previous work.

Overview of the Proposed Approach

The proposed descriptor is based on two main steps : 1. Building a locally invariant representation : using multiple high-pass filters, we generate different sparse representations for different filter orientations. A pooling operator is applied on the orientation to increase the local invariance to orientation change. The process is repeated for multiple image resolutions for a richer representation.

2. Building a globally invariant representation : the first step generates various images that encode different texture information. We also include the multi-resolution versions of the input to provide low-pass information. We need a mapping that transforms this set of images into a stable, fixed-size descriptor. We use multifractal analysis to statistically describe each one of these images. We present a new method that extracts rich information directly from local singularity exponents. The local exponents encode rich multi-scale texture information. Their log-normalized distribution represents a stable mapping which is invariant to complex bi-Lipschitz transforms. As a result, the proposed multifractal descriptor is proven mathematically to be robust to strong environmental changes.

Locally Invariant Representation

A locally invariant representation aims at increasing the similarity of local statistics between textures of the same class. To build this representation, we construct a simple convolutional network where the input image is convolved with a filter bank for various orientations, and then pooled to reduce local orientation change. The multilayer extension consists in repeating the same process for various image resolutions on the low-pass output of the previous resolution, which offers a richer representation.

Given an input texture I, the image is first low-pass filtered with a filter ψ l to reduce small image domain perturbations and produce an image J 1,0 . This image is then filtered with multiple zero-mean high-pass filters ψ k,θ , where k denotes the filter number and θ its orientation. High-pass responses encode higher-order statistics that are not present in the low-pass response J 1,0 . A more stable approach consists in applying the modulus on the high-pass responses, which imposes symmetric statistics and improves invariance of the local statistics. Applying multiple filtering with multiple different filters naturally increases the amount of texture information that are going to be extracted further via multifractal analysis. In order to increase the local invariance to orientation, we apply a pooling operator φ θ : R i×j×n → R i×j on the oriented outputs for each filter :

J 1,k = φ θ (|J 1,0 ψ k,θ | , θ = θ 1 , ..., θ n ) , k = 1, ..., K, (1)
where n is the number of orientations and i × j is the size of the low-pass image.

As a result, we obtain 1 low-pass response and K high-pass responses, each image is encoding different statistics. For a richer representation, we repeat the same operation for different resolutions s = 2 0,...,-L , where s = 1 is the finest resolution and s = 2 -L is the coarsest resolution. The image generation process is then generalized as follows :

J s,k =      I ψ l k = 0 , s = 1 ↓ (J 2s,0 ψ l ) k = 0 , s = 1 φ θ (|J s,0 ψ k,θ | , θ = θ 1 , ..., θ n ) k = 1, ..., K, (2)
where ↓ denotes the downsampling operator. We found that calculating statistics on multiple resolutions instead of a single one increases significantly the robustness of the descriptor. This can be expected because two textures may seem "more similar" at a lower resolution. As a result, the intra-class variability decreases as the resolution decreases, but keeping higher resolution images is important to ensure extra-class decorrelation.

Dimensionality Reduction with Pooling

Using multiple filters ψ k,θ increases dramatically the size of the image set. Knowing that each image J s,k will be used to extract statistics using multifractal analysis, this will result in a very large descriptor. One resulting issue is the high dimensionality of the training set. Another one is the processing time as the statistics should be applied on each image. We propose to merge different high-pass responses J s,k together to reduce the number of images. A straightforward approach would be to gather various images {J s,k , k = t, .., u} and then apply a pooling operator φ r that is going to merge each image subset into one single image J s,kt,..,u :

J s,kt,..,u = φ r ( J s,k , k = t, .., u). (3)
As a result, the number of high-pass responses will be decreased ; this leads to a reduced size descriptor. The pooling operator φ r can be either the mean or the min/max functions. We take φ r as a maximum function in this paper. An example is given in Figure 2 for one resolution s = 0 using 6 high-pass filters and one low-pass filter. The number of images is reduced from 7 to 3. For 5 resolutions (s = 2 0,...,-4 ), the total number of images goes from 35 to 15, which is an important reduction.

I J s,k ↓ J0,0 J0,1 J0,2 J0,3 J0,4 J0,5 J0,6 ↓ φr ↓ φr ↓ J0,0 J 0,k 1,...,3 J 0,k 4,...,6
Fig. 2: Image generation example applied on the texture input I for one resolution using 6 high-pass filters. The images J 0,1...6 are a result of the orientation pooling (eq. 2). The 6 images are reduced to 2 images using a pooling operator φ r on similar responses to reduce the dimensionality. The same process is repeated for multiple resolutions.

Globally Invariant Representation

Once the set of low-pass and high-pass images is generated, we need to extract global statistics, a mapping into a fixed-size descriptor, which is globally invariant to the complex physical transformations. We propose to use a new multifractal approach to statistically describe textures suffering from strong environmental changes. To understand the difference between the proposed method and the previous work, we first present the standard fractal and multifractal analysis framework used by the previous methods, we then introduce the proposed approach.

Multifractal Analysis In a nutshell, a fractal object E is self-similar across scales. One characteristic of its irregularity is the so-called box fractal dimension. By measuring a fractal object on multiple scales r, the box fractal dimension is defined as a power-law relashionship between the scale r and the smallest number of sets of length r covering E [START_REF] Falconer | Techniques in Fractal Geometry[END_REF]:

dim(E) = lim r→0 log N (r, E) -log r , (4)
Using squared boxes of size r, this dimension can be estimated numerically, known as the box-counting method. Multifractal analysis is an extension of this important notion. A multifractal object F is composed of many fractal components F 1,...,f . In this case, a single fractal dimension is not sufficient to describe this object. The multifractal spectrum is the collection of all the associated fractal dimensions that describe the multifractal object.

It is easy to show mathematically that the fractal dimension is invariant to bi-Lipschitz transformations [START_REF] Xu | Viewpoint invariant texture description using fractal analysis[END_REF], which includes various transformations such as nonrigid transformations, view-point change, translation, rotation, etc.. As a result, the multifractal spectrum is also invariant to these transformations. This makes the multifractal spectrum an attractive tool to globally describe textures. However, the box-counting method gives a rather crude estimation of the real fractal dimension. The fractal dimension is estimated for each fractal set using a log-log regression. As the resolution r is supposed to be very small (r → 0), using small-sized boxes on a relatively lowresolution image results in a biased estimation due to the relatively low-resolution of real-world images [START_REF] Arneodo | The thermodynamics of fractals revisited with wavelets[END_REF]. It has been used as the core of various recent multifractal texture descriptors [START_REF] Xu | A projective invariant for textures[END_REF][START_REF] Xu | Combining powerful local and global statistics for texture description[END_REF][START_REF] Xu | A new texture descriptor using multifractal analysis in multi-orientation wavelet pyramid[END_REF][START_REF] Ji | Wavelet domain multifractal analysis for static and dynamic texture classification[END_REF] that use the same box-counting method to build the final descriptor. We present a different method to statistically describe textures using multifractal analysis. Contrary to previous methods, we use a new measure which is based on the distribution of local singularity exponents. It can be shown in fact that this measure is related to the true multifractal spectrum, and its precision is proven by the high-accuracy of the proposed descriptor. Moreover, this approach is computationally efficient, which permits to achieve high accuracy at reduced processing time.

Proposed Multifractal Descriptor

The proposed method first estimates the local singularity exponents h(x) on each pixel x, and then applies the empirical histogram followed by log operator to extract the global statistics ϕ h = log(ρ h + ). This operation is performed on all the resulting images of the first step, which results in multiple histograms ϕ hi . The concatenation of all these histograms forms the final descriptor.

Let J be an image, and µ ψ (B(x, r)) = B(x,r) (J ψ r )(y)dy a positive measure, where ψ r is an appropriate wavelet at scale r (Gaussian in our case) and B(x, r) a closed disc of radius r > 0 centered at x. Multifractal analysis states that the wavelet projections scale as power laws in r [START_REF] Turiel | Reconstructing images from their most singular fractal manifold[END_REF][START_REF] Yahia | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. Part I: singularity analysis[END_REF][START_REF] Turiel | The multifractal structure of contrast changes in natural images: From sharp edges to textures[END_REF]. We use a microcanonical evaluation [START_REF] Yahia | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. Part I: singularity analysis[END_REF] which consists in assessing an exponent h(x) for each pixel x :

µ ψ (B(x, r)) ≈ α(x)r h(x) , r → 0. ( 5
)
The validity of equation ( 5) has been tested on a large dataset [START_REF] Turiel | The multifractal structure of contrast changes in natural images: From sharp edges to textures[END_REF], which proves that natural images exhibit a strong multifractal behavior. Introducing the log, the formula is expressed as a linear fit :

log(µ ψ (B(x, r))) ≈ log(α(x)) + h(x) log(r) , r → 0. (6)
Rewriting the equation in the matrix form permits to calculate all the exponents at once by solving the following linear system :

    1 log(r1) . . . . . . 1 log(r l )     A log(α(x1)) • • • log(α(x N )) h(x1) • • • h(x N ) η =     log(µ ψ (B(x1, r1))) . . . log(µ ψ (B(x N , r1))) . . . . . . . . . log(µ ψ (B(x1, r l ))) . . . log(µ ψ (B(x N , r l )))     b , (7)
argmin η ||Aη -b|| 2 2 , h(x i ) = η(2, i), (8)
where N is the number of pixels of the image J, l is the number of scales used in the log-log regression. This matrix formulation is computationally efficient and plays an important role in the speed of the proposed method. Given the local exponents h(x), which is an image of the same size of J that describes the local irregularities at each pixel, we need to extract now a fixed-size measure that globally describes the statistics of h(x). Using the box-counting method, this would require extracting all the fractal fractal sets F h = {x | h(x) ≈ h}, and then calculating the box-counting dimension for each set F h . As discussed before, this approach leads to a crude estimation of the true multifractal spectrum due to the actual low-resolution of real-world images. Moreover, a log-log regression should be performed on each fractal set. Instead, we propose to use the empirical histogram ρ h followed by a log operator :

ϕ h = log(ρ h + ), (9)
where ≥ 1 is set to provide stability. The distribution of the local exponents is an invariant representation which encodes the multi-scale properties of the texture. The log acts as a normalization operator that nearly linearizes histogram scaling and makes the descriptor more robust to small perturbations. This way, we have access to reliable statistics 1 . This log-histogram is calculated on each image generated in the first step, which results in a set of histograms ϕ h 1,...,M , where M is the total number of generated images. The final descriptor ϕ is constructed by concatenating ( ) all the generated histograms :

ϕ = M m ϕ hm ; (10)
A descriptor example is given in Figure 3. This descriptor ϕ is the result of the concatenation of 14 log exponents histograms calculated on the images generated with the first step of the method presented in section 2.2 and further explained in Figure 2. Three images are generated for each scale s ; a low-pass response is presented in red, and two high-pass responses are presented in black and gray in the figure 2 .

Analysis

The basic multifractal framework consists in generating multiple images and then extracting statistics using multifractal analysis. Multifractal descriptors are mathematically invariant to bi-Lipschitz transforms, which even includes non-rigid transformation and view-point change. The proposed method follows the same strategy, but is substantially different from the previous methods. The differences lie in both the image generation step and the statistical description. For instance, the WMFS method [START_REF] Ji | Wavelet domain multifractal analysis for static and dynamic texture classification[END_REF] low first high range second high range

s = 2 0 s = 2 -1 s = 2 -2 s = 2 -3 s = 2 -4 s = 2 0 s = 2 -1 s = 2 -2 s = 2 -3 s = 2 -4 s = 2 0 s = 2 -1 s = 2 -2 s = 2 -3
Fig. 3: A descriptor example using a low-pass response and two high-pass responses for 5 resolutions s = 2 0,...,-4 . The exponents log-histogram is calculated for each response and for multiple image resolutions s.

generates multiple images for multiple orientations, each oriented image is then analyzed using Daubechies discrete wavelet transform as well as using the wavelet leaders [START_REF] Wendt | Wavelet leaders and bootstrap for multifractal analysis of images[END_REF]. The multifractal spectrum (MFS) is then estimated for each image, for a given orientation using the box-counting method. Each MFS is then concatenated for a given orientation and the final descriptor is defined as the mean of all the descriptors over the orientation. Contrary to this method, we use different high-pass filters instead of one single analyzing wavelet, which permits to extract different statistics. Generating multiple descriptors for multiple orientations is computationally expensive. In contrast, we generate only one descriptor. To ensure local robustness to orientation, we apply a pooling operator on the filtered responses. This approach is much more computationally efficient. Finally, the core of our method is the new multifractal descriptor which permits to extract accurate statistics, contrary to the popular box-counting method as explained in the previous section. The proposed method takes about 0.7 second to extract the whole descriptor on an image of size 480 × 640, compared to 37 seconds as reported in the state-of-the-art multifractal method [START_REF] Ji | Wavelet domain multifractal analysis for static and dynamic texture classification[END_REF]. Experiments show that the proposed descriptor permits also to achieve higher accuracy, especially in large-scale situations when the extra-class decorrelation is a challenging issue.

Pre and Post Processing

Pre-processing and post-processing can improve the robustness of a texture recognition system. For instance, the method in [START_REF] Xu | A new texture descriptor using multifractal analysis in multi-orientation wavelet pyramid[END_REF] performs a scale normalization step on each input texture using blob detection. This step first estimates the scale of the texture and then a normalization is applied, which aims at increasing the robustness to scale change. Other texture classification methods such as [START_REF] Varma | Locally invariant fractal features for statistical texture classification[END_REF] use Weber's law normalization to improve robustness to illumination. We do not use any scale normalization step such as [START_REF] Xu | A new texture descriptor using multifractal analysis in multi-orientation wavelet pyramid[END_REF][START_REF] Ji | Wavelet domain multifractal analysis for static and dynamic texture classification[END_REF], we rather use sometimes histogram equalization to improve robustness to illumination change. We also use a post-processing on features vector ϕ using wavelet domain soft-thresholding [?]. This step aims at increasing the intra-class correlation by reducing small histogram perturbations (for more details, please refer to the supplementary material).

Classification and Training Strategies

The second part of our work concerns the training aspect of the texture recognition problem. The globally invariant representation offers a theoretically stable invariant representation via accurate multifractal statistics. However, there are other small transformations and perturbations that may occur in real-world images and this is where a good training strategy will help us to take advantage of the proposed descriptor in practice. We work on two ideas :

1. The choice of the classifier can improve recognition rates : we introduce a simple combination between the Generative PCA classifier [START_REF] Sifre | Rotation, scaling and deformation invariant scattering for texture discrimination[END_REF] and SVMs [START_REF] Scholkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond[END_REF]. 2. The lack of data is an issue, how to get more data? : Given an input training texture image, we synthetically generate more images by changing its illumination and scale. We call this strategy "synthetic training".

Experiments on challenging public benchmark datasets, including a large-scale dataset with 250 classes, validates the robustness of the proposed solution.

Classification

Support Vector Machines SVMs [START_REF] Scholkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond[END_REF] are widely used in texture classification [START_REF] Xu | A projective invariant for textures[END_REF][START_REF] Xu | A new texture descriptor using multifractal analysis in multi-orientation wavelet pyramid[END_REF][START_REF] Xu | A new texture descriptor using multifractal analysis in multi-orientation wavelet pyramid[END_REF][START_REF] Ji | Wavelet domain multifractal analysis for static and dynamic texture classification[END_REF][START_REF] Liu | Sorted random projections for robust texture classification[END_REF][START_REF] Liu | Texture classification from random features[END_REF]. Commonly used kernels are mainly RBF Gaussian kernel, polynomials and χ 2 kernel. Extension to multiclass can be done via strategies such as one-vs-one and one-vs-all. In this paper, we use the one-vs-all strategy with an RBF-kernel. It consists in building a binary classifier for each class as follows : for each class, a positive label is assigned to the corresponding instances and a negative label is affected to all the remaining instances. The winning class c svm can be chosen based on probability estimates [START_REF] Chang | LIBSVM: A library for support vector machines[END_REF] or a simple score maximization :

c svm = argmax 1≤c≤Nc {f svm (x, c)} , f svm (x, c) = Mc i=1 α c i y c i K(x c i , x) + b c , (11)
where α c i are the optimal Lagrange multipliers of the classifier representing the class c, x c i are the support vectors of the class c, y c i are the corresponding ±1 labels, N c is the number of classes and x is the instance to classify.

Generative PCA Classifier The generative PCA (GPCA) classifier is a simple PCAbased classifier recently used in [START_REF] Bruna | Invariant scattering convolution networks[END_REF][START_REF] Sifre | Rotation, scaling and deformation invariant scattering for texture discrimination[END_REF]. Given a test descriptor x, GPCA finds the closest class centroid E({x c }) to x, after ignoring the first D principal variability directions. Let V c be the linear space generated by the D eigenvectors of the covariance matrix of largest eigenvalues, and V ⊥ c its orthogonal complement. The generative PCA classifier uses the projection distance associated to P V ⊥ c :

c pca = argmin 1≤c≤Nc ||P V ⊥ c (x -E({x c })) || 2 . ( 12
)
Classification consists in choosing the class c pca with the minimum projection distance.

GPCA-SVM Classifier

We propose to combine GPCA and SVMs in one single classifier. The idea behind this combination comes from the observation that SVMs and GPCA often fail on different instances. As a result, a well-established combination of these classifiers should theoretically lead to improved performance. We propose a combination based on the distance between the score separation of each classifier output

c f inal = c svm if f svm (x, c svm ) -f svm (x, c pca ) ≥ th svm c pca otherwise, (13)
where th svm is a threshold parameter. The score separation gives an idea of SVMs' accuracy to classify a given instance. Another similar approach would be using probability estimates [START_REF] Chang | LIBSVM: A library for support vector machines[END_REF] instead of the score. If the measure f svm (x, c svm ) -f svm (x, c pca ) is relatively important, this means that SVMs are quite "confident" about the result. Otherwise, the classifier selects the GPCA result. Determining the best threshold th svm for each instance is an open problem. In this paper, we rather fix a threshold value for each experiment. We generally select a small threshold for small training sets and larger thresholds for larger sets. Even if this strategy is not optimal, experiments show that the combination improves the classification rates as expected.

Synthetic Training

One important problem in training is coping with the low amount of examples. We propose a simple strategy to artificially add more data to the training set by changing illumination and scale of each instance of the training set. While this idea seems simple, it can have a dramatic impact on the performance as we will see in the next section.

Multi-Illumination

Training Given an input image I, multi-illumination training consists in generating other images of the same content of I but with different illumination. There are two illumination cases ; the first one consists in uniform changing by image scaling of the form aI, where a is a given scalar. The second case consists in nonuniform changing using histogram matching with a set of histograms. The histograms can come from external images, or even from the training set itself (for example by transforming or combining a set of histograms).

Multi-Scale

Training Given an input image I, multi-scale training consists simply in generating other images of the same size as I by zooming-in and out. In this paper, we use around 4 generated images, 2 by zooming-in and 2 others by zooming-out.

Texture Classification Experiments

We present in this section texture classiffication results conducted on standard public datasets UIUC [24,1], UMD [25] and ALOT [START_REF] Burghouts | Material-specific adaptation of color invariant features[END_REF]27], as well as a comparison with 9 state-of-the-art methods. Implementation details In order to build a fast texture classification system, we use only two high-pass filtering responses, which results in 3 histograms per image resolution 3 . The number of the image scales is fixed to 5. The filter bank consists in high-pass wavelet filters (Daubechies, Symlets and Gabor). A more robust descriptor can be built by increasing the number of filters and orientations. Filtering can be parallelized for faster processing. While augmenting the number of filters slightly improves classification results, the minimalist setup presented above, coupled with the training strategies introduced in this paper, permits to outperform existing techniques while offering in addition computational efficiency.

Evaluation

We evaluate the proposed system and compare it with state-of-the-art methods for 50 random splits between training and testing. The evaluation consists in three steps :

1. log-histogram vs. box-counting : We evaluate the precision of our log-histogram method and compare it with the box-counting method used in previous methods. 2. Learning efficiency : We compare the proposed GPCA-SVM combination with single GPCA and SVM results and see how the proposed synthetic training strategy improves classification rates. 3. We compare our main results with 9 state-of-the-art results. log-histogram vs. box-counting In this experiment, we replace the log-histogram step of our approach with the box-counting method widely used in the previous multifractal methods to see if the proposed log-histogram leads to a more accurate bi-Lipschitz invariance. The results are presented in Figure 5. As can be seen, the log-histogram approach leads to higher performance, especially when more data is available. This clearly shows that indeed, the log-histogram leads to a better bi-Lipschitz invariance, as theoretically discussed before. The log-histogram is a simple operation that permits our system to achieve high computational efficiency. Learning Efficiency In this experiment, we first compare the proposed GPCA-SVM combination with single GPCA and SVM classifiers using the proposed descriptor. Each dataset is presented in the form D x (y) where x is the name of the dataset and y is the training size in number of images. The best results are in bold. As can be seen in Table 1, the GPCA-SVM does indeed improve classification rates. We expect to get even better results with a better strategy to set the threshold parameters th svm as in the proposed experiments, the threshold is fixed for all the instances. Now we compare the results with and without the proposed synthetic training strategy. As can be seen, synthetic training leads to a dramatic improvement. This is a very interesting approach as it increases only the training time. The system can achieve higher recognition accuracy for almost the same computational effiency. For the UMD and ALOT datasets, we use uniform illumination change with the multiplicative parameter a in the range [0.9, 0.95, 1.05, 1.1]. For the UIUC dataset, we use the nonuniform illumination change with two histograms. For the multi-scale training, we use only four generated images (two by zooming-in and two other by zooming-out), which increases the training set 9 times in the UMD and UIUC datasets (no mutli-scale training is used for the ALOT dataset). 4 . As can be seen, the proposed method outperforms the published results on the 3 datasets. Compared to the leading method [START_REF] Sifre | Rotation, scaling and deformation invariant scattering for texture discrimination[END_REF], our system seems to better handle viewpoint change and non-rigid deformations. This is clearly shown in the results on the UIUC dataset that exhibits strong enviromental changes. This result can be expected as the scattering method builds invariants on translation, rotation and scale changes, which does not include viewpoint change and non-rigid deformations. Contrary to this, using accurate multifractal statistics, our solution produces descriptors that are invariant to these complex transformations. The proposed system maintains a high performance on the UMD dataset. It is worth noting that on this dataset, the images are of high resolution (1280 × 960), which gives an advantage over the UIUC dataset. However, we did not use the original resolution, we rather rescale the images to half-size for faster processing. The high accuracy shows that the proposed multifractal method is able to extract robust invariant statistics even on low-resolution images.

On the large-scale dataset ALOT, the proposed method maintains high performance.

Recall that this dataset contains 250 classes with 100 samples each. This is a very challenging dataset that evaluates the extra-class decorrelation of the produced descriptors.

A robust descriptor should increase the intra-class correlation, but should also decrease the extra-class correlation and this has be evaluated on a large-scale data set which contains as many different classes as possible. The results on the ALOT dataset clearly show a significant performance drop of the leading multifractal method WMFS. The proposed solution in fact outperforms the WMFS method even without synthetic train as can be seen in Table 1. This proves that the proposed descriptor is able to extract a robust invariant representation.

DescriptorsFig. 1 :

 1 Fig.1: Intra-class variability demonstration. The three textures 1, 2 and 3 exhibit strong changes in scale and orientation as well as non-rigid deformations. As can be seen, the proposed descriptor is nearly invariant to these transformations (see section 2).

Datasets

  Description The UIUC dataset [24,1] is one of the most challenging texture datasets presented so far. It is composed of 25 classes, each class contains 40 grayscale images of size 480 × 640 with strong scale, rotation and viewpoint changes in uncontrolled illumination environment. Some images exhibit also strong non-rigid deformations. Some samples are presented in Figure4. The UMD dataset [25] is similar to UIUC with higher resolution images (1280 × 960) but exhibits less non-rigid deformations and stronger illumination changes compared to UIUC. To evaluate the proposed method on a large-scale dataset, we choose the ALOT dataset[START_REF] Burghouts | Material-specific adaptation of color invariant features[END_REF]27]. It consists of 250 classes, 100 samples each. We use the same setup as the previous multifractal methods[START_REF] Ji | Wavelet domain multifractal analysis for static and dynamic texture classification[END_REF]: grayscale version with half resolution (768 × 512). The ALOT dataset is very challenging as it reprensents a significantly larger number of classes (250) compared to UIUC and UMD (25) and very strong illumination change (8 levels of illumination). The viewpoint change is however less dramatic compared to UIUC and UMD.

Fig. 4 :

 4 Fig. 4: Texture samples from the UIUC dataset [24,1]. Each row represents images from the same class with strong enviromental changes.

Fig. 5 :

 5 Fig.5: Comparison between the box-counting method and the proposed log-histogram approach for various dataset training sizes(5, 10 and 20). The proposed approach leads to a more accurate descriptor.

Table 1 :

 1 Classification rates comparison using GPCA-SVM and synthetic training. Discussions We compare the proposed method MCMA (Multilayer Convolution -Multifractal Analysis) with 9 state-of-the-art methods for 50 random splits between training and testing, for different training sizes. Results are presented in Table 2. The best results are in bold

	D ALOT (50)

A mathematical relationship between the log exponents histogram and the multifractal spectrum is presented in the supplementary material.

[START_REF] Zhang | Local features and kernels for classification of texture and object categories: A comprehensive study[END_REF] A histogram was discarded for s = 2 -4 in the second high response (in gray) due to the large size of the filter which is larger than the actual size of the input image at resolution s = 2 -4 .

Except for ALOT dataset, we use 3 high-pass responses for a more robust representation.

Detailed results with standard deviation can be found in the supplementary material.

Acknowledgements

Hicham Badri's PhD is funded by an INRIA (Direction of Research) CORDI-S grant. He is making a PhD in co-supervision with INRIA and Mohammed V-Agdal University -LRIT, Associated Unit to CNRST (URAC 29).

Conclusion

This paper presents a fast and accurate texture classification system. The proposed solution builds a locally invariant representation using a multilayer convolution architecture that performs convolutions with a filter bank, applies a pooling operator to increase the local invariance and repeats the process for various image resolutions. The resulting images are mapped into a stable descriptor via multifractal analysis. We present a new multifractal descriptor that extracts rich texture information from the local singularity exponents. The descriptor is mathematically validated to be invariant to bi-Lipschitz transformations, which includes complex environmental changes. The second part of paper tackles the training part of the recognition system. We propose the GPCA-SVM classifier that combines the generative PCA classifier with the popular kernel SVMs to achieve higher accuracy. In addition, a simple and efficient "synthetic training" strategy is proposed that consists in synthetically generating more training data by changing illumination and scale of the training instances. Results outperforming the state-of-the-art are obtained and compared with 9 recent methods on 3 challenging public benchmark datasets, while ensuring high computational efficiency.


Pietro Sternai
email: pietro.sternai@univ-orleans.fr

Laurent Jolivet

Armel Menant

Taras Gerya

Driving the upper plate surface deformation by slab rollback and mantle flow

Keywords: slab tearing, slab rollback, mantle flow, surface deformation

The relative contribution of crustal and mantle processes to surface deformation at convergent plate margins is still controversial. Conflicting models involving either extrusion mechanisms or slab rollback, in particular, were proposed to explain the surface strain and kinematics across the Tethyan convergent domain. Here, we present new high-resolution 3D thermo-mechanical numerical joint models of continental collision, oceanic subduction and slab tearing, which for the first time allow self-consistent reproduction of first-order Tethyan tectonic structures such as back-arc rifting and large-scale strike-slip faults accommodating continental escape. These models suggest that mantle flow due to slab rollback and tearing can modulate the surface strain and kinematics by locally enhancing trench retreat and dragging the upper plate from below. These results highlight the active role of the asthenospheric flow in driving the surface strain, not only by modulating the vertical stresses and producing dynamic topography but also through sub-horizontal motion. We discuss the implications of these findings based on observations across the Aegean-Anatolian and eastern Indian-Eurasian domains, though similar considerations may as well apply to other settings.

Introduction

Plate fragmentation, mountain building, formation of extensional basins, and major strike-slip fault zones characterized the long-term evolution of the Tethyan convergent domains (e.g., [START_REF] Ricou | Tethys Reconstructed-Plates, Continental Fragments and Their Boundaries since 260-Ma from Central-America to South-Eastern Asia[END_REF]. Surface deformation is classically attributed to either collision-related processes such as crustal shortening, extrusion and gravitational spreading (e.g., [START_REF] Tapponnier | Slip-line field theory and large scale continental tectonics[END_REF]; Le [START_REF] Pichon | Kinematics, topography, shortening and extrusion in the India-Eurasia collision[END_REF] or subduction-related mechanisms such as slab-pull, rollback and trench retreat (e.g., [START_REF] Dewey | Lithospheric stress, deformation, and tectonic cycles: the disruption of Pangaea and the closure of Tethys[END_REF][START_REF] Royden | The tectonic expression slab pull at continental convergent boundaries[END_REF]. However, subduction of the Tethyan lithosphere and collision along the African, Arabian, Indian and Eurasian margins often coexisted, interacting with each other to set jointly the surface kinematics and strain. For instance, oceanic subduction below the Hellenic domain coexisted with active convergence across the neighbouring Bitlis-Zagros region since approximately the late Oligocene (e.g. [START_REF] Sengör | The North Anatolian transform fault: its age, offset and tectonic significance[END_REF][START_REF] Jolivet | Mediterranean extension and the Africa-Eurasia collision[END_REF][START_REF] Allen | Late Cenozoic reorganization of the Arabia-Eurasia collision and the comparison of short-term and long-term deformation rates[END_REF]. GPS measurements show a westward motion of Anatolia and a southward propagation of the Hellenic trench with respect to Eurasia [START_REF] Reilinger | Geodetic constraints on the tectonic evolution of the Aegean region and strain accumulation along the Hellenic subduction zone[END_REF], resulting in an overall counter-clockwise rotation (Figure 1a). Most of this rotation is accommodated by the north-Anatolian fault (NAF), a strike-slip fault zone joining the Bitlis-Zagros domain to the east and the Aegean/Hellenic extensional back-arc region to the west (e.g., [START_REF] Mckenzie | Active tectonics of the Mediterranean region[END_REF][START_REF] Sengör | The North Anatolian transform fault: its age, offset and tectonic significance[END_REF][START_REF] Armijo | Westward propagation of the north Anatolian into the northern Aegean: timing and kinematics[END_REF]. Tomographic models (Li et al. 2008) show a low-velocity anomaly below western Turkey, which was interpreted as a major tear in the Hellenic slab [START_REF] Boorder | Late Cenozoic mineralization, orogenic collapse and slab detachment in the European Alpine Belt[END_REF][START_REF] Govers | Lithosphere tearing at STEP faults: response to edges of subduction zones[END_REF] (Figure 1b). Recent studies foster the hypothesis that the mantle flow induced by slab rollback codetermines the surface deformation [START_REF] Faccenna | Slab detachment beneath eastern Anatolia: A possible cause for the formation of the North Anatolian fault[END_REF][START_REF] Jolivet | From mantle to crust: Stretching the Mediterranean[END_REF][START_REF] Pérouse | Bridging onshore and offshore present-day kinematics of central and eastern Mediterranean: Implications for crustal dynamics and mantle flow[END_REF]Faccenna et al., 2013;[START_REF] Jolivet | Aegean Tectonics: strain localization, slab tearing and trench retreat[END_REF]. [START_REF] Zhang | Continuous deformation of the Tibetan Plateau from global positioning system data[END_REF][START_REF] Gan | Present day crustal motion within the Tibetan Plateau inferred from GPS measurements[END_REF][START_REF] Reilinger | Geodetic constraints on the tectonic evolution of the Aegean region and strain accumulation along the Hellenic subduction zone[END_REF]). b,d) Tomographic models (Li et al. 2008) showing the Vp anomalies at the depths indicated. Stars show the locations of slab tears. Red and blues lines in all panels represent the major strike-slip fault zones and thrusts and subduction fronts, respectively.

Another example spanning larger scales and a longer temporal evolution is southeastern Asia in which collision of India with Eurasia and the Andaman-Sumatra subduction contribute jointly to surface deformation across the Himalaya, Indochina and Indonesia since ~45 Ma (e.g., [START_REF] Royden | Surface Deformation and Lower Crustal Flow in Eastern Tibet[END_REF][START_REF] Tapponnier | Oblique stepwise Rise and growth of the Tibet plateau[END_REF][START_REF] Royden | The geological evolution of the Tibetan Plateau[END_REF][START_REF] Replumaz | Amount of Asian lithospheric mantle subducted during the India/Asia collision[END_REF]. Geodetic measurements relative to stable Eurasia show a prominent clockwise rotation around the Eastern Himalayan Syntaxes (EHS) characterized by eastward motion of eastern Tibet and western Sichuan, south-eastward motion in northern Yunan and south to south-eastward motion in southern Yunan [START_REF] Zhang | Continuous deformation of the Tibetan Plateau from global positioning system data[END_REF][START_REF] Gan | Present day crustal motion within the Tibetan Plateau inferred from GPS measurements[END_REF][START_REF] Arrajehi | Geodetic constraints on present day motion of the Arabian Plate: implications for Red Sea and Gulf of Aden rifting[END_REF]) (Figure 1c). Several strike-slip fault zones, which arrange from the convergent domain north of Tibet, such as the Altyn-Tagh Fault, to the extensional domain east and south of it, such as the Red River Fault (RRF), accommodate the majority of this motion (e.g., [START_REF] Molnar | Cenozoic tectonics of Asia: Effects of a continental collision[END_REF][START_REF] Tapponnier | Slip-line field theory and large scale continental tectonics[END_REF]. Several basins, such as the South China Sea, have developed east of Tibet and offshore the Indochina peninsula in the Oligocene and Miocene [START_REF] Tapponnier | Propagating extrusion tectonics in Asia: New insights from simple experiments with plasticine[END_REF][START_REF] Taylor | Origin and history of the South China Sea basin, in The tectonic and geologic evolution of southeast Asian seas and islands, part 2[END_REF][START_REF] Tapponnier | On the mechanics of the collision between India and Asia, in Collision tectonics[END_REF]). The main direction of extension across these regions was N-S and the relationships between the formation of these basins and large-scale shear zones such as the RRF have been widely debated: these topographic lows formed either as a pull-apart basin at the south-easternmost extremity of the RRF (i.e., extrusion: e.g., [START_REF] Tapponnier | Propagating extrusion tectonics in Asia: New insights from simple experiments with plasticine[END_REF][START_REF] Tapponnier | On the mechanics of the collision between India and Asia, in Collision tectonics[END_REF] or as more classical back-arc basins behind the retreating subduction (i.e., slab retreat: e.g., [START_REF] Taylor | Origin and history of the South China Sea basin, in The tectonic and geologic evolution of southeast Asian seas and islands, part 2[END_REF][START_REF] Jolivet | Japan Sea, opening history and mechanism, a synthesis[END_REF][START_REF] Fournier | Backarc extension and collision: an experimental approach to the tectonics of Asia[END_REF]). The deep mantle structures associated with subduction of ancient Tethyan lithosphere vary significantly along the collision boundary (Figure 1d). In particular, while deep and shallow slabs may be still (spatially) connected in the central Himalayas (Li et al. 2008), there is no evidence for such a connection beneath the eastern Himalayas. Tomographic images show a low velocity zone below the EHS at depth larger than ~200km separating the Indian and Andaman-Sumatra slabs (Li et al. 2008). This slab tear might have facilitated southward to westward rollback of the Andaman-Sumatra subduction after the onset of collision between India and Eurasia, and the mantle flow through the slab window and around the EHS might have codetermined the surface strain (Figure 1c,d) [START_REF] Holt | Correlated crust and mantle strain fields in Tibet[END_REF].

Conflicting models involving extrusion tectonics originated by the Arabia/India-Eurasia collision or rollback of the Tethyan torn slabs were proposed to explain the surface strain and kinematics across both the Aegean-Anatolian and eastern Himalayan regions (e.g., [START_REF] Tapponnier | Slip-line field theory and large scale continental tectonics[END_REF][START_REF] Tapponnier | On the mechanics of the collision between India and Asia, in Collision tectonics[END_REF][START_REF] Jolivet | Right-lateral shear along the Northwest Pacific Margin and the India-Eurasia Collision[END_REF][START_REF] Royden | Surface Deformation and Lower Crustal Flow in Eastern Tibet[END_REF][START_REF] Armijo | Westward propagation of the north Anatolian into the northern Aegean: timing and kinematics[END_REF][START_REF] Fournier | Backarc extension and collision: an experimental approach to the tectonics of Asia[END_REF][START_REF] Becker | Mantle conveyor beneath the Tethyan collisional belt[END_REF]. In addition, slab rollback is often proposed as a passive source of space to accommodate extrusion (e.g., [START_REF] Tapponnier | Propagating extrusion tectonics in Asia: New insights from simple experiments with plasticine[END_REF][START_REF] Armijo | Westward propagation of the north Anatolian into the northern Aegean: timing and kinematics[END_REF]). However, the southward migration of the Hellenic trench is sensibly faster than the westward motion of Anatolia [START_REF] Reilinger | Geodetic constraints on the tectonic evolution of the Aegean region and strain accumulation along the Hellenic subduction zone[END_REF]) (Figure 1a), and classical extrusion models fail in predicting widespread extension such as that across the Ordos, Sichuan, Yunan regions in absence of a "free boundary" (e.g., [START_REF] Jolivet | Right-lateral shear along the Northwest Pacific Margin and the India-Eurasia Collision[END_REF][START_REF] Jolivet | Japan Sea, opening history and mechanism, a synthesis[END_REF][START_REF] Royden | Surface Deformation and Lower Crustal Flow in Eastern Tibet[END_REF]; [START_REF] Wang | Present-Day Crustal Deformation in China Constrained by Global Positioning System Measurements[END_REF]. These flaws suggest that slab rollback is in fact an active player and not just a passive source of space. These geological, geophysical and geodetic observations, which provide insights on the long-term and recent geologic history, can constrain self-consistent three-dimensional (3D) numerical thermo-mechanical geodynamic models in order to achieve a first-order quantification of the relative contribution of the deep and shallow dynamics in setting the surface deformation across convergent plate boundaries. Previous numerical studies demonstrated along trench variations and complex mantle and crustal processes and surface deformation in subductioncollision systems (e.g., [START_REF] Van Hunen | Continental collision and slab break-off: A comparison of 3-D numerical models with observations[END_REF][START_REF] Capitanio | Subduction and slab breakoff controls on Asian indentation tectonics and Himalayan western syntaxis formation[END_REF][START_REF] Li | Collision of continental corner from 3-D numerical modeling[END_REF][START_REF] Duretz | Slab detachment in laterally varying subduction zones: 3D numerical modeling[END_REF]. Although some implications for the large-scale tectonics across the Tethyan domain have been suggested [START_REF] Li | Collision of continental corner from 3-D numerical modeling[END_REF][START_REF] Duretz | Slab detachment in laterally varying subduction zones: 3D numerical modeling[END_REF], the relative contributions of crustal and mantle flow processes to surface deformation remained widely unaddressed.

Methods

The purpose of this paper is to assess the relative contributions of the lithospheric and mantle dynamics in setting the surface deformation at the transition between collisional and subduction domains. To this aim, we use 3D Cartesian thermo-mechanical modeling to resolve self-consistently the joint effects of oceanic subduction, slab tearing and continental collision on mantle flow and surface deformation at high resolution and accounting for realistic crustal and mantle rheologies. In addition to the reference model setup and boundary conditions, we outline below some of the principal features of the numerical model, but more detailed information can be found in [START_REF] Gerya | Robust characteristics method for modelling multiphase visco-elasto-plastic thermo-mechanical problems[END_REF], [START_REF] Zhu | 3-D dynamics of hydrous thermal-chemical plumes in oceanic subduction zones[END_REF], [START_REF] Gerya | Introduction to Numerical Geodynamic Modelling[END_REF] and [START_REF] Li | Collision of continental corner from 3-D numerical modeling[END_REF]. Details on the parametric study can be found in the Supplementary Material.

3D numerical model, setup and boundary conditions

All numerical experiments were performed with the code "I3ELVIS" [START_REF] Gerya | Introduction to Numerical Geodynamic Modelling[END_REF], solving the momentum, continuity and energy equations based on a staggered finite difference scheme combined with a marker-in-cell technique. Realistic and non-Newtonian visco-plastic rheologies [START_REF] Ranalli | Rheology of the Earth, Deformation and Flow Processes in Geophysics and Geodynamics[END_REF] are used in the model. The model domain (Figure 2) measures 2000x328x1000 km in the x, y and z dimensions, respectively. This domain is resolved by 501x165x197 grid points resulting in a resolution of 4, 2 and 5 km in the x, y and z dimensions, respectively. ~130 million randomly distributed markers are used for advecting the material properties and temperatures. The velocity boundary conditions are free slip at the top (y=0km) and at both the front and back boundaries (z=0km and z=1000km). The left and right boundaries (x=0km and x=2000km) use constant x-parallel velocities, which define the material influx.

Global mass conservation is ensured by material outflux through the lower permeable boundary (y=328km). The top surface of the lithosphere is calculated dynamically as an internal free surface through a 12km thick layer of "sticky air" [START_REF] Gerya | Introduction to Numerical Geodynamic Modelling[END_REF]. Surface processes are implemented using a highly simplified gross-scale erosion-sedimentation law according to which instantaneous erosion limits mountains height to 4 km above the reference water level (y=12km), whereas instantaneous sedimentation limits trench depth to 8 km below the water level [START_REF] Gerya | Robust characteristics method for modelling multiphase visco-elasto-plastic thermo-mechanical problems[END_REF]. The initial temperature gradient in the asthenospheric mantle is ~0.5 • C km -1 (adiabatic). The thermal boundary conditions are 0 • C for the upper boundary and nil horizontal heat flux across the vertical boundaries. An infinity-like external temperature condition [START_REF] Gerya | Introduction to Numerical Geodynamic Modelling[END_REF]) is imposed on the lower boundary. All the numerical experiments were performed using 20 cores on the BRUTUS cluster at the ETH-Zürich (the approximate run time of a simulation is two weeks).

A joint model of continental collision, oceanic subduction and slab tearing

Our numerical experiments were not designed to reproduce any particular case, but can be applied to the transition between a continental indenter flanked by active oceanic subduction. We stress that the ultimate goal of our numerical experiments is to investigate the effects of slab tearing on the mantle and surface strain, but the potential causes of a tear in the slab are beyond the scope of this work. As such, we impose mechanical weaknesses in the oceanic plate and spatial variations of the upper continental crust (as detailed below), but do not claim that, in natural contexts, slab tears are necessarily generated by these factors.

The reference model is set up as follows. Three continental plates (C1, C2 and C3) are included into the initial model domain (Figure 2a and Table 1). is the density, C p is the specific heat capacity, E a is the activation energy, V a is the activation volume, n is the stress exponent,  0 is the reference viscosity, Hr  eff is the effective internal friction angle. Cohesion is 1 MPa for each phase.

The crust of C2 and C3 are 50km and 45km thick, respectively. By varying the crustal thickness of the upper plate (C1) along the trench-parallel direction (z) one can control the slab descent rates [START_REF] Nikolaeva | Subduction initiation at passive margins: Numerical modeling[END_REF]) and therefore generate a tear in the slab. In particular, we imposed the crust of C1 as 35km thick where z≤490km and 45km thick where z>510km (linear interpolation in between). The initial thermal structure of C1, C2 and C3 is laterally uniform with 0 • C at the surface and 1300 • C at 90, 140 and 150 km depth respectively. An oceanic domain, characterized by a trench-parallel weak fracture zone to initiate subduction and a trenchperpendicular weak fratcure zone to allow for slab tearing, separates the three continental plates (Figure 2a,b). The thermal structure of the oceanic lithosphere is that of a half-space cooling age of 120Ma (e.g., [START_REF] Turcotte | [END_REF]. Uniform and constant in time x-parallel velocities equal to ~1.8cm yr -1 (convergence) are imposed to the x=2000km boundary. The evolution of the reference model is shown in Figure 3 and in Video S1, S2 and S3.

The details and timing of the model evolution depends upon the assumed thermal and mechanical parameters (see the Supplementary Materials), but the overall sequence of events (Figure 3) is robust and includes: 1) early subduction of the oceanic domain between C1 and C2 and 2) later subduction of the oceanic domain between C1 and C3. These major events are allowed by tearing of the slab and separated in time by the collision between C1 and C2. Trench retreat, back-arc extension and progressive approach/collision of the continental plates characterize the surface evolution during the entire model run.

Results

In this section, we outline the numerical results with particular focus on the relationships between mantle and crustal flow and the effects of slab tearing and continental collision on the surface strain.

Mantle vs. crustal flow

A toroidal motion (with a minor ascending poloidal component) induced by the downgoing slab dominates the deep dynamics during the first subduction event (Figure 4a). This toroidal flow occurs both below C1 and through the opening slab window. At the surface, extensional and compressional deformation occurs in the back-arc and collisional domain, respectively. Crustal and mantle flows are essentially decoupled and surface deformation seems to be substantially driven by incipient slab rollback and trench retreat within the oceanic domain between C1 and C3 and continental collision between C1 and C2. The rheological stratification of the continental plates is essentially unaffected by this early subduction event, implying that the traction due to mantle flow is imposed at the base of the lithosphere (Figure 4c). The model evolution during the second subduction event is highly conditioned by the precedent toroidal flow of hot asthenosphere. This flow warmed the incipiently subducting lithosphere close to the slab tear, enhancing here slab rollback and trench retreat, generating a non-cylindrical (i.e., not parallel to the initial upper plate margin) subduction (Figures 3,4b and 4e). This subduction geometry induces a lateral (i.e., with non-nil trench-parallel component of motion), horizontal to sub-horizontal return flow of hot asthenospheric material, which further enhances the asymmetry of subduction (Figure 4b). The surface deformation across C1 is characterized at this stage by a faster trench retreat in the proximity of the collisional domain and above the slab tear, while localized strike-slip deformation propagating toward the subducting plate as the model progresses accommodates the rotation of a crustal block (Figure 3b and4b).

The rheological stratification of the upper plate in the back-arc extensional domain is characterized by the absence of a rigid lithospheric mantle and the velocity-depth profile shows particularly fast mantle flow (Figure 4d), parallel in places to the surface velocity field (Figure 4b). This translates into shear stresses ( xy ) up to ~100MPa applied directly to the base of the crust by the asthenospheric flow (Figure 5).

Effects of slab tearing and continental collision on the surface strain

By comparing the reference model outputs to similar simulations, but not implying collision between C1 and C2 or slab tearing, one can decouple the first-order crustal and mantle contributions to the surface strain. It appears that strain localization along a major strike-slip fault zone and thus the onset of a distinctive block-like deformation and surface rotation only occurs when both collision between rigid continental plates and slab tearing take place (Figures 4b and6). In particular, the trench parallel component of motion within the upper plate is sensibly reduced in absence of continental collision or slab tearing. Besides, slab tearing allows/enhances rollback of the subducting plate, trench parallel/toroidal migration of the asthenosphere, trench retreat and extension within the upper plate. It is noteworthy that, when both continental collision and slab tearing allow for a blocklike deformation and surface rotation, the velocity fields within the crust and mantle are sensibly faster within the back-arc extensional region and close to the subduction trench than in the collisional domain (Figure 4b).

Discussion

Various combinations of processes such as extrusion tectonics, suction exerted by slab rollback, gravitational collapse of the collisional belt and mantle drag associated with regional plum-like upwelling were proposed to explain the observed surface kinematics across the Tethyan convergent domain (e.g., [START_REF] Tapponnier | Slip-line field theory and large scale continental tectonics[END_REF][START_REF] Royden | Surface Deformation and Lower Crustal Flow in Eastern Tibet[END_REF][START_REF] Armijo | Westward propagation of the north Anatolian into the northern Aegean: timing and kinematics[END_REF][START_REF] Allen | Late Cenozoic reorganization of the Arabia-Eurasia collision and the comparison of short-term and long-term deformation rates[END_REF][START_REF] Alvarez | Protracted continental collisions argue for continental plates driven by basal traction[END_REF][START_REF] Becker | Mantle conveyor beneath the Tethyan collisional belt[END_REF]. Most of the proposed models, however, neglected the role played by the dynamical asthenospheric flow originated by slab rollback and tearing. A major outcome of our numerical experiments is that the toroidal motion of the asthenosphere at the edge of a retreating tearing slab is able to modulate the thermal state of the lithosphere, in turn affecting the geometry and dynamics of later subduction events and associated mantle and surface strain. Tectonic inheritance is clearly not a new concept, but it most often refers to the mechanical properties of the lithosphere and mantle (e.g., Thomas 2006), while we show here that also the thermal and kinematic inheritances can significantly affect the surface evolution. Moreover, the mantle return flow following slab rollback and tearing is able to produce tectonically significant (i.e., up to ~100MPa; Bürgmann & Dresen 2008) shear stresses at the base of the lithosphere and crust. Thus, surface strain across subduction zones is likely to be jointly driven by slab rollback and mantle flow, especially across hot and thinned environments, characterized by the absence of a rigid mantle lithosphere such as the back-arc domains. It was demonstrated that the horizontal traction generated at the base of deep continental roots from regional thermal anomalies within the mantle dragged the continents together along the Tethyan axis and was able to protract convergence even after continental collision and brake-off of the Tethyan slabs [START_REF] Alvarez | Protracted continental collisions argue for continental plates driven by basal traction[END_REF]. The sub-horizontal mantle flow then seems to be able to actively drive the surface strain and kinematics across both thin and thick continental lithospheres. As such, the sub-horizontal couplings between the Earth's mantle and surface should benefit from the same attention that is paid to the vertical couplings and associated "dynamic topography" (e.g., Braun 2010).

Our numerical experiments show that the surface strain at the transition between a continental indenter and an oceanic subduction can be dominated by the suction exerted by slab rollback (enhanced by slab tearing) and the associated mantle flow, as demonstrated by faster crustal flow close to the subduction trench than to the collisional domain (Figures 4b). Although the continental indenter is necessary for setting a block-like deformation and surface rotation, its role is that of a passive bulwark if the overall convergence rates are slower than slab rollback.

Thus, the effectiveness of extrusion tectonics, crustal shortening and gravitational collapse at the margins of a collisional belt should be assessed in the light of the dynamics of the neighbouring oceanic subduction (Figure 7), since the "free boundary" is the most important factor controlling these mechanisms. We cannot resolve the timing of slab tearing across the Aegean-Anatolian or the eastern Indian-Eurasian margins nor can we address its relationships to slab rollback or continental collision through our modeling. However, trench migration rates are inversely correlated to the width of the subducting plate [START_REF] Schellart | Evolution and diversity of subduction zones controlled by slab width[END_REF]), which implies that the observed tears decoupling the actively subducting slabs from the mantle lithosphere below the collisional domains might have facilitated slab rollback and modulated trench retreat after continental collision, consistently with our numerical experiments. A common feature in the Mediterranean metamorphic core complexes exhumed in back-arc domains from the lower and middle crust is a highly asymmetric deformation with constant sense of shear -i.e., hangingwall motion relative to footwall - [START_REF] Jolivet | Subduction, convergence and the mode of backarc extension in the Mediterranean region[END_REF][START_REF] Jolivet | From mantle to crust: Stretching the Mediterranean[END_REF]. Our modeling results suggest that this asymmetric deformation is substantially controlled by the basal shear stresses imposed to the crust by the asthenospheric return flow.

The modeled surface velocity field shown in Figure 4b is also similar (regardless of the sense of spin) to the observed surface rotation across the Aegean-Anatolian and the eastern Indian-Eurasian margins (Figure 1). The majority of the observed surface rotation across the natural examples is accommodated by several fault zones amongst which the NAF, the Kun-Lun, Xianshuihe-Xiaojiang and the Altyn Tagh faults (e.g., [START_REF] Sengör | The North Anatolian transform fault: its age, offset and tectonic significance[END_REF][START_REF] Tapponnier | Propagating extrusion tectonics in Asia: New insights from simple experiments with plasticine[END_REF], whose kinematics is reasonably well reproduced in our numerical experiments by highly localised strike-slip deformation (Figures 3 and4b). A faster crustal flow close to the subduction trench than to the collisional domain in our reference model is consistent with GPS measurements across the Aegean-Anatolian region, and suggests that the Aegean extension and extrusion of Anatolia are principally driven by the subduction dynamics and associated mantle flow, instead of collision-related mechanisms (Figure 7). The arrival of the Australian plate to the Sumatra trench in middle-late Miocene times [START_REF] Ricou | Tethys Reconstructed-Plates, Continental Fragments and Their Boundaries since 260-Ma from Central-America to South-Eastern Asia[END_REF]) might have reduced the suction exerted by oceanic subduction. Nonetheless, the crust along portions of southeastern Tibet has been extending since at least ~4Myr (e.g., [START_REF] Royden | Surface Deformation and Lower Crustal Flow in Eastern Tibet[END_REF][START_REF] Wang | Present-Day Crustal Deformation in China Constrained by Global Positioning System Measurements[END_REF]. The many basins and low topographies across the Sichuan, Yunan, Indochina and Sunda provinces as well as the back-arc basins further testify widespread extension, implying that the Andaman-Sumatra and other western Pacific subduction zones substantially affected the surface strain across southeastern Asia during most of the Cenozoic [START_REF] Jolivet | Right-lateral shear along the Northwest Pacific Margin and the India-Eurasia Collision[END_REF]). In addition, tomographic models highlight the presence of hot and shallow asthenospheric material across southeastern Tibet (e.g., Li et al., 2008), which is analogous to the asthenosphere flowing at direct contact with the base of the crust and rotating around the edge of the collisional domain in our numerical experiments. Since a free boundary seems necessary to allow for continental extrusion (e.g., [START_REF] Tapponnier | Slip-line field theory and large scale continental tectonics[END_REF][START_REF] Tapponnier | Propagating extrusion tectonics in Asia: New insights from simple experiments with plasticine[END_REF][START_REF] Jolivet | Right-lateral shear along the Northwest Pacific Margin and the India-Eurasia Collision[END_REF][START_REF] Fournier | Backarc extension and collision: an experimental approach to the tectonics of Asia[END_REF]) and this free boundary is provided by oceanic subduction which we proved able to actively drive the surface strain more efficiently than the continental indentation, our experiment confirm that the rollback of the Andaman-Sumatra and the associated mantle return flow actively determined the surface strain east and southeast of the EHS (Figure 7).

Instantaneous global flow models calculated from mantle density anomalies, suggest that both the pull from the Tethyan slabs and mantle upwelling from underneath the Ethiopia-Afar and Reunion-Carlsberg regions are necessary to reproduce the crustal motion across the Tethyan collisional belts [START_REF] Becker | Mantle conveyor beneath the Tethyan collisional belt[END_REF]Faccenna et al. 2013). Our results suggest that the dynamical mantle flow associated solely with slab rollback and tearing can produce a distinct block-like surface rotation, without the intervention of an external dragging force such as the upwelling from the lower mantle. In our models, however, slab tearing is triggered along the preexisting weak oceanic fracture zones and modelled slabs are not tied to previously subducted plates as the Tethyan slabs probably were during most of the Cenozoic [START_REF] Ricou | Tethys Reconstructed-Plates, Continental Fragments and Their Boundaries since 260-Ma from Central-America to South-Eastern Asia[END_REF]. Mantle

upwellings, thus, have a major role in setting crustal motion at the scale of the whole Tethyan margins by dragging the northern and southern continents together along the Tethyan axis [START_REF] Alvarez | Protracted continental collisions argue for continental plates driven by basal traction[END_REF], warming and propelling the upper continental plates [START_REF] Becker | Mantle conveyor beneath the Tethyan collisional belt[END_REF]Faccenna et al. 2013) and, possibly, by thermally enhancing slab tearing, but the rollback of the Tethyan slabs and the associated mantle return flow are likely to be the primary drivers of the surface deformation at the edge of the Tethyan collisional belts.

Conclusions

3D Cartesian thermo-mechanical models of self-consistent oceanic subduction, slab tearing and continental collision calibrated on geologic, geophysical and geodetic observations and accounting for realistic crustal and mantle rheologies enabled us to conclude that:

 The dynamical asthenospheric flow originated by slab rollback and tearing is a primary factor controlling the surface evolution at the transition between collisional and subduction domains. In particular, the toroidal motion of the asthenosphere at the edge of a retreating tearing slab is able to modulate the thermal state of the lithosphere, in turn affecting the geometry and dynamics of later subduction events and associated mantle and surface strain. In addition, the sub-horizontal mantle return flow following slab rollback and tearing is able to produce tectonically significant shear stresses at the base of the lithosphere and crust. The surface strain across subduction zones is thus likely to be jointly driven by slab rollback and mantle flow, especially across hot and thinned environments such as the back-arc domains.

 The surface strain at the transition between a continental indenter and an oceanic (Gerya 2010) and experimentally derived deformation laws and parameters [START_REF] Ranalli | Rheology of the Earth, Deformation and Flow Processes in Geophysics and Geodynamics[END_REF]. The natural complexity, however, is strongly simplified and major uncertainties arise from the imposed boundary conditions, the input parameters and the inferred distribution of the continental plates within the initial model setup. We here present a few of the additional computations we performed to establish that the model predictions and conclusions discussed in the text are robust with respect to these choices. We performed numerical experiments with x-parallel velocities imposed to the x=2000km boundary ranging between 0 and 5.4cm yr -1 (three time faster convergence rates than those of the reference model). In Figure S2 we show the numerical experiment in which uniform and constant in time x-parallel velocities are equal to 5.4cm yr -1 . Although the timing and some details of the subduction dynamics differ from the reference model, the overall geodynamic evolution (see text)

is similar in all simulations. Conclusions are thus robust with respect to the imposed velocity boundary conditions.

Nikolaeva and co-workers presented an extensive investigation of the parameters that codetermine the initiation and evolution of an ocean-continent subduction on a 2D version of the numerical model used in this study [START_REF] Nikolaeva | Subduction initiation at passive margins: Numerical modeling[END_REF]. They found that the lithospheric and the crustal thickness of the upper continental plate is of primary importance. We ran simulations in which the initial 1300°C isotherm of the upper continental plate (C1) is at 50, 60, 70, 80, 100 and 110 km depth and the thickness of the crust varies between 25-55 km where z≤490km and 35-65km where z>510km (linear interpolation in between). Overall, increased thickness of the continental crust enhances subduction rates and leads to a faster geodynamic evolution. On the contrary, an increased lithospheric thickness leads to slower subduction rates and overall evolution. Modern GPS data were used to calibrate the numerical experiments with respect to the crustal and lithospheric thicknesses. The magnitude of the modeled surface velocities and GPS measurements is in a first approach similar when the lithospheric thickness is set to 90km and a the crustal thickness is equal to 35km where z≤490km and 45km where z>510km, with linear interpolation in between (i.e., the reference model). We show in Figure S3 a numerical experiment in which the initial 1300°C isotherm is at 70km depth and the crust is 30km where z≤490km and 40km where z>510km, with linear interpolation in between. In Figure S4, the initial 1300°C isotherm is at 100km depth and the crust is 35km where z≤490km and 45km where z>510km, with linear interpolation in between. Although the timing and some details of the subduction dynamics differ from the reference model, the overall sequence of events (see text) is similar in all simulations. Conclusions are thus robust with respect to these parameters. In Figure S5 we show a sensitivity test for the inferred distribution of the continental plates, with particular focus on the role of continental collision in setting the patterns of surface deformation. All thermal and mechanical boundary conditions as well as the rheological parameters of this numerical experiment are the same as those of the reference model, the only difference being that C1 is the only continental plate within the model domain. Similarly to the

Figure 1 :

 1 Figure 1: Maps of the Aegean-Anatolian and Indian-Eurasian systems. a,c) Topography and GPS velocities with 95% confidence ellipses[START_REF] Zhang | Continuous deformation of the Tibetan Plateau from global positioning system data[END_REF][START_REF] Gan | Present day crustal motion within the Tibetan Plateau inferred from GPS measurements[END_REF][START_REF] Reilinger | Geodetic constraints on the tectonic evolution of the Aegean region and strain accumulation along the Hellenic subduction zone[END_REF]). b,d) Tomographic models(Li et

Figure 2 :

 2 Figure 2: Model setup. a) The 3D model domain (2000x328x1000 km in the x, y and z directions, respectively) with colours showing different rock types: 1-sediments; 2-continental crust; 3-oceanic crust; 4-lithospheric mantle; 5asthenospheric mantle; 6-hydrated/serpentinized mantle (initially imposed "weak fracture zone"); 7-partially molten continental crust; 8-partially molten mantle. The top layer ("sticky air", y<~12 km) and other phases are cut off for clarity. Convergence is imposed by applying uniform and constant-in-time velocity equal to ~1.8 cm yr -1 on the

Figure 3 :

 3 Figure 3: Evolution of the reference model. a) Colours on the top panels show different rock types: 1-sediments; 2continental crust; 3-oceanic crust; 4-lithospheric mantle; 5-hydrated/serpentinized mantle (initially imposed "weak fracture zone"); 6-partially molten continental crust; 7-partially molten mantle. The top layer ("sticky-air", y<~12 km), the asthenosphere and other phases are cut off for clarity. The subducting plates shown on the bottom panels are visualised through an iso-viscosity contour equal to 10 22 Pa s. Red, yellow and black dashed lines on the lower panels show the plan view of C1, C2 and C3, respectively. b) Selected frames showing the 800°C isotherm (mostly corresponding to the crust-lithospheric mantle transition within the C1 domain), color-coded by the second invariant of the rate-of-strain tensor on this surface. A strike-slip fault zone (also visible at the surface, Figure 4b) is generated by the sub-horizontal mantle flow (red arrows) following slab rollback and tearing. Note that the timing of panels in (a) and (b) is different.

Figure 4 :

 4 Figure 4: a-b) Plan view of selected temporal frames. Black and red arrows represent the velocity vectors in the continental crust and asthenosphere, respectively. Colours on the left panel show the second invariant of the rate-ofstrain tensor. The subducting plates are shown on the right panel through an iso-viscosity contour equal to 10 22 Pa s. Red, yellow and black dashed lines on the right panels show the plan view of C1, C2 and C3, respectively. c-d) Depth-Velocity-Viscosity profiles (y-parallel). The profile location is shown by the yellow stars in (a) and (b). Solid lines, representing the magnitude of the velocity vector, follow the left-axis. Dashed lines, representing the viscosity, follow the right-axis. Shaded background colours show the phase at each depth: grey-continental crust, bluelithospheric mantle/slab, and purple-asthenospheric mantle. e) Temperature-depth profiles at 12.5 Ma (the green and blue star in (b) shows the locations of the dashed and solid profile, respectively.

Figure 5 :

 5 Figure 5: Shear stress (a) and strain rate (b) distribution at the base of the upper plate continental crust (C1) at 12.5Myr modeled time on the reference simulation. Red arrows show the velocity field within the asthenosphere. The area inside the black line is where the asthenosphere is at direct contact with the base of the continental crust.

Figure 6 :

 6 Figure 6: a) Plan view of a selected temporal frame of a simulation without slab tearing. b) Plan view of a selected temporal frame of a simulation without continental collision. Black and red arrows represent the velocity vectors in the continental crust and asthenosphere, respectively. Colours on the left panels show the second invariant of the rate-of-strain tensor. The subducting plates are shown on the right panel through an iso-viscosity contour equal to 10 22 Pa s. Red, yellow and black dashed lines on the right panels show the plan view of C1, C2 and C3, respectively.

Figure 7 :

 7 Figure 7: Illustrative representation of the dynamical interactions between continental collision, oceanic subduction, mantle flow and surface deformation.

  subduction can be dominated by the suction exerted by slab rollback (enhanced by slab tearing) and the associated mantle flow. The continental indenter acts as a passive bulwark if the overall convergence rates are slower than slab rollback. Thus, the effectiveness of extrusion tectonics, crustal shortening and gravitational collapse at the margins of a collisional belt should be assessed in the light of the dynamics of the neighbouring oceanic subduction. A faster crustal flow close to the Hellenic subduction trench than to the collisional domain (consistently reproduced by our numerical models) suggests that the suction exerted by the Hellenic subduction and associated mantle flow is the principal driver of the Aegean extension and extrusion of Anatolia, instead of currently proposed collision-related mechanisms. Our experiment also confirm that the rollback of the Andaman-Sumatra and the associated mantle return flow affected the surface strain east and southeast of the EHS as also suggested by the several basins and low topographies across the Sichuan, Yunan, Indochina and Sunda provinces, which testify widespread extension. Mantle upwelling drag and propel the northern and southern continents along the Tethyan axis and possibly thermally enhance slab tearing. The rollback of the Tethyan slabs and the associated mantle return flow, however, are likely to be primary drivers of the surface deformation at the edge of the Tethyan collisional belts.

Figure

  Figure S1: a) Same as the left panel of Figure 3b (reference model for comparison). Red, yellow and black dashed

Figure S2 :

 S2 Figure S2: Sensitivity test for the imposed velocity boundary conditions. For each temporal frame, the phase distribution (1-sediments; 2-continental crust; 3-oceanic crust; 4-lithospheric mantle; 5-hydrated/serpentinized

Figure S3 :

 S3 Figure S3: Sensitivity test for the geothermal gradient of the upper plate (C1). For each temporal frame, the phase distribution (1-sediments; 2-continental crust; 3-oceanic crust; 4-lithospheric mantle; 5-hydrated/serpentinized mantle (initially imposed "weak fracture zone"); 6-partially molten continental crust; 7-partially molten mantle) and the subducting plates are shown on the top and bottom panel, respectively. The top layer ("sticky air", y<~12 km), the asthenosphere and other phases are cut off for clarity. The subducting plates are visualised through an isoviscosity contour equal to 10 22 Pa s. Red, yellow and black dashed lines on the lower panels show the plan view of

Figure S4 :

 S4 Figure S4: Sensitivity test for the geothermal gradient of the upper plate (C1). For each temporal frame, the phase distribution (1-sediments; 2-continental crust; 3-oceanic crust; 4-lithospheric mantle; 5-hydrated/serpentinized mantle (initially imposed "weak fracture zone"); 6-partially molten continental crust; 7-partially molten mantle) and the subducting plates are shown on the top and bottom panel, respectively. The top layer ("sticky air", y<~12 km),the asthenosphere and other phases are cut off for clarity. The subducting plates are visualised through an isoviscosity contour equal to 10 22 Pa s. Red, yellow and black dashed lines on the lower panels show the plan view of C1, C2 and C3, respectively. The initial 1300°C isotherm is at 100km depth and the crust is 35km where z≤490km and 45km where z>510km, with linear interpolation in between.

Table 1 :

 1 Qz

. and Ol. correspond to the abbreviations of Quartzite and Olivine. k denotes the thermal conductivity,  0

Acknowledgements: This work has been financially supported by the Laboratoire d'Excellence (LABEX) VOLTAIRE of the University of Orléans and the European Research Council (ERC) under the 7 th Framework Programme of the European Union (ERC Advanced Grant, grant agreement No 290864, RHEOLITH). Numerical simulations were performed on the ETH-Zürich cluster BRUTUS. We are thankful to Claudio Faccenna and Thorsten Becker for stimulating discussions and to Kosuke Ueda for constructive comments on an early version of the manuscript and tips for the modelling. The paper benefited from insightful revision by Leigh Royden and two anonymous reviewers.

Supplementary Materials:

Videos S1-S3 show the evolution of the reference model.

Video S1-S3: Evolution of the reference model (video thumbnails)

Our models are based on state-of-the-art numerical solutions for geodynamic problems reference model, the overall evolution of subduction is characterized by slab rollback and tearing as well as ascending toroidal mantle flows. However, the deformation of the upper plate in Figure S5b is sensibly less localised than in Figure 5a (i.e., the reference model) due to the lack of continental collision occurring on one side of the model domain. This test suggests that while slab rollback, tearing and consequent mantle flow control the surface deformation, continental collision allows the localisation of this deformation along one major strike-slip lineament (Figures 3 and4), as discussed in the text. Red, yellow and black dashed lines on the right panels show the plan view of C1, C2 and C3, respectively.

In Figure S6 we show a sensitivity test for the process of slab tearing. All thermal and mechanical boundary conditions as well as the rheological parameters of this numerical experiment are the same as those of the reference model, the only difference being that we did not impose any weak fracture zone in the lithospheric mantle to allow for slab tearing. The deformation of the upper plate in Figure S6b is sensibly different than in Figure S6a (i.e., the reference model), which imply that slab tearing and associated mantle flow is a primary factor for the evolution of the surface strain and kinematics, as discussed in the text. Red, yellow and black dashed lines on the right panels show the plan view of C1, C2 and C3, respectively.