Jérôme Franc

Radu State

Thomas Engel

Olivier Festor
email: olivier.festor@inria.fr

Enforcing Security with Behavioral Fingerprinting

Keywords:

Although fingerprinting techniques are helpful for security assessment, they have limited support to advanced security related applications. We have developed a new security framework focusing especially on the authentication reinforcement and the automatic generation of stateful firewall rules based on behavioral fingerprinting. Such fingerprinting is highly effective in capturing sequential patterns in the behavior of a device. A new machine learning technique is also adapted to monitor high speed networks by evaluating both computational complexity and experimented performances.

I. INTRODUCTION

Fingerprinting automatically infers specific information from monitored devices: operating system [START_REF] Zalewski | P0f[END_REF], [START_REF] Lyon | Nmap Network Scanning: The Official Nmap Project Guide to Network Discovery and Security Scanning[END_REF], protocol used [START_REF] Ma | Unexpected means of protocol inference[END_REF] or device type [START_REF] Franc ¸ois | Automated behavioral fingerprinting[END_REF], [START_REF] Franc ¸ois | Machine Learning Techniques for Passive Network Inventory[END_REF], [START_REF] Franc ¸ois | PTF: Passive Temporal Fingerprinting[END_REF] (or protocol stack). The latter infers automatically the types of devices: brand name, series of hardware or name and version of software.

Such tools are helpful for security assessment by revealing devices having security flaws in a network like for instance operating systems or software having known vulnerabilities. Moreover, an attacker can hide itself by spoofing normal device type. To detect that, our approach takes benefit from our previous fingerprinting method [START_REF] Franc ¸ois | Automated behavioral fingerprinting[END_REF], [START_REF] Franc ¸ois | PTF: Passive Temporal Fingerprinting[END_REF] which characterizes a device by its behavior i.e., the way it interacts with others. This kind of fingerprinting can also catch attack behaviors to create automatically attack patterns which might be inputted as signatures for an Intrusion Detection System (IDS).

However, a main issue is the real-time applicability of such techniques for safeguarding efficiently users, hosts and network infrastructure. Thus, the contribution of this paper is threefold. Firstly, a new method is proposed for managing behavioral fingerprints. Secondly, we propose the Behavioral Security Framework using fingerprinting for security purposes like vulnerable device discovery, spoofing detection or automated attack signature creation. Finally, complexity issues are evaluated analytically and experimentally. Regarding our previous work [START_REF] Franc ¸ois | Automated behavioral fingerprinting[END_REF], [START_REF] Franc ¸ois | PTF: Passive Temporal Fingerprinting[END_REF], we have reused the construction of the behavioral fingerprints as well as the SVM classification.

This paper is structured as follows. The framework is described in section II. Section III details the fingerprinting technique. The security applications are assessed in section IV as well as the complexity in section V. Section VI highlights main related work. Conclusion is given in section VII.

II. FINGERPRINTING BASED SECURITY FRAMEWORK

The Behavioral Security Framework depicted in figure 1 is based on behavioral fingerprinting. This paper introduces two main applications: authentication reinforcement (authentication module) and attack/intrusion detection (attack module), and voluntary omits applications related to security vulnerabilities assessment based on a general inventory of connected devices, which may be directly derived from fingerprinting [START_REF] Franc ¸ois | Machine Learning Techniques for Passive Network Inventory[END_REF]. Figure 1 shows the architecture including a first level of protection (firewall and/or an IDS) screening the traffic as a first step towards a more rigorous authentication.

A. Authentication reinforcement

The first module of our security framework reinforces the user authentication. Obviously, the latter cannot be exclusively based on fingerprinting but should use an authentication mechanism such as password protection or certificate based solutions but they are still faced with security issues. For example, stealing passwords using a dictionary, social engineering or phishing is easily possible [START_REF] Van Der Merwe | Characteristics and responsibilities involved in a phishing attack[END_REF]. Thus, device fingerprinting is helpful in detecting an anomalous device type even if the user succeeds to authenticate. Anomalies could be related to attacker activities but can also reflect a normal behavior change. Therefore, the authentication module within the behavioral security framework (figure figure 1) only provides results to the original authentication module responsible to trigger counter-measures. For example, if the fingerprinting detects an anomaly, the user could be required to do an additional check (second password, email confirmation, usage of a TLS certificate, alarm reported to the administrator...). Moreover, if the second check succeeds, the user profile has to be updated.

In fact, anomalies might be detected in parallel in two manners. On one hand, a difference between the announced type, since most of protocols define a user agent field [START_REF] Rosenberg | SIP: Session Initiation Protocol[END_REF], [START_REF] Fielding | Hypertext Transfer Protocol -HTTP/1.1[END_REF], and the type automatically inferred can be detected. The attackers, who execute specific tools, prefer to announce faked user agents for avoiding detection (type spoofing) since they often specific tools. Thus, the first authentication reinforcement has to compare the inferred device type by fingerprinting with the announced device type. On the other hand, the user behaviors may be monitored over time to profile the types they use and track deviations from them.

B. Attack module

The attack module (figure figure 1) might be trained with attack generators such as KiF [START_REF] Abdelnur | KiF: A stateful SIP Fuzzer[END_REF], which produces and discovers automatically new attacks for SIP [START_REF] Rosenberg | SIP: Session Initiation Protocol[END_REF] (Session Initiation Protocol). The fingerprint of the attacks can then be automatically generated and expressed with stateful firewall rules. If the firewall is able to filter the traffic regarding such rules, its configuration is updated, otherwise it has to delegate the process to the attack module itself which will then only output the action to trigger.

Otherwise, the attack module can be trained with normal behaviors to act as an anomaly intrusion detection system by looking for behaviors too much divergent from a profile.

III. BEHAVIORAL FINGERPRINTING

Behavioral fingerprinting was introduced in [START_REF] Franc ¸ois | PTF: Passive Temporal Fingerprinting[END_REF]. The key idea is to represent the devices by their behaviors: the way they interact with other ones. Since the original work includes all theoretical details, this section gives an overview and provides some extensions: a new technique based on K-TRACE and the definition of the result confidence.

A. TR-FSM model

The behavior of a device type is represented as a Temporal Random Parameterized Tree Extended Finite State Machine (TR-FSM) and corresponds to the manner a device sends the requests and how it replies to the requests from other ones. Moreover, the behavioral fingerprinting takes in account the message delays. A TR-FSM is a tree based representation of the device specific protocol state machine including delays as shown in figure 2. In fact each node represents the type of the message which is emitted (prefixed by !) or received (prefixed by ?). An edge between two nodes indicates that they have been successively monitored in the same session. The delay between them is put as an attribute on the edge. Basically, for each observed sequence of messages, a branch is created. Then, the branches are aggregated using a generic root node and the common prefix are merged such as the node !INVITE in figure 2(a).

The protocol knowledge is reduced to the message types, which can be automatically inferred as highlighted in [START_REF] Franc ¸ois | PTF: Passive Temporal Fingerprinting[END_REF] or based on a message parser. It eases the applicability of our behavioral security framework to any protocol. Furthermore, the time measurements are dependent of the network locations. Hence, the round-trip time has to be measured and filtered out. Otherwise, if the measurements are always done at the same location, this step can be avoided.

B. Classification

Once the TR-FSMs are constructed, they have to be classified. Support Vector Machines (SVM) are good trade-offs between accuracy and computational complexity [START_REF] Hsu | A comparison of methods for multiclass support vector machines[END_REF] and have proved their efficiency [START_REF] Franc ¸ois | PTF: Passive Temporal Fingerprinting[END_REF]. To extend the evaluation, this paper introduces another method, TRACE/K-TRACE [START_REF] Cifarelli | A kernel k-means procedure for classification[END_REF], which is, like SVM, dedicated to separate data points which are not linearly disjoint and we propose a very low computational version hereafter. The limited overhead is also highlighted in our evaluation while the accuracy is not degraded. A common classification problem with non linearily separable data is depicted in figure 3 where traditional clustering methods fails to separate the different clusters.

1) SVM: Supposing P different kinds of devices which form the set D = {d 1 , . . . , d P }, L = {l 1 , . . . , l M } is the training set of M TR-FSMs. Multi-class SVM resolve optimization problems to determine separating hyperplanes in a higher dimensional space where points have been projected using a function φ (figure 4), and with maximum margins for each pair of classes (types). Hence, assuming two types, d l and d k , the hyperplane h lk is characterized by a subset of points of L, namely support vectors, SV lk , belonging to type d l or d k . The decision function is directly derived from these supports vectors. Thus, f lk (t x ) returns the predicted type for t regarding h lk . This corresponds of the side-position of a point to assign, $ in figure 4. Since there are several hyperplanes (one per pair of types), the predicted type corresponds to mostly chosen one by all the decision functions. As defining the projection function φ is hard, a kernel function is applied to a pair of points, K(t i , t j ), and has to follow strong properties [START_REF] Cristianini | An introduction to support Vector Machines: and other kernel-based learning methods[END_REF]. Details may be found in our prior work [START_REF] Franc ¸ois | PTF: Passive Temporal Fingerprinting[END_REF] proving the correctness of a TR-FSM based kernel function which mainly relies on counting the number of common paths (sequence of nodes from the root to any other node) weighted by the delays.

2) K-TRACE: K-TRACE [START_REF] Cifarelli | A kernel k-means procedure for classification[END_REF] is a kernel-based extension of the Total Recognition by Adaptive Classification Experiments (TRACE). TRACE algorithm is executed in a supervised manner and is well suited for piece-case linearly separable data points for which traditional clustering methods fail because they look for continuous clusters. To discover the multiple clusters corresponding to a single class (device type), TRACE finds multiple barycenters (figure 4) thanks to a learning process (also called training). The initialization creates one cluster per class containing all points from L of this class. Hence, the set C = {c 1 1 , . . . , c 1 P } is built from the trees t i of the training set:

c 1 j = {t i |real(t i ) = d j } (1)
where real(t i ) returns the real identity of a tree in the learning set L. Actually, c j i represents the j th clusters of the type i. Then, the algorithm computes the barycenters of the clusters and divides a cluster in case of incoherency, i.e. the closest barycenter from a point is not the one of its own cluster. More precisely, it works as follows:

1) compute each barycenter, B = {b j i }, for each cluster c j i of C 2) compute all Euclidean distances d(t i , b l k ) for all t i belonging to the training set and all b l k in B

3) selects all points whose the closest barycenter is not the one of the corresponding class:

S = {t i |∃b l k , b s r , t i ∈ c l k , d(t i , b s r ) < d(t i , b l k )} (2)
4) stop if S is empty 5) select the tree of S, b + , which is the farthest from its barycenter: 6) get all the clusters, C ′ , containing points of the type real(b + ) 7) apply k-means on C ′ by using the original barycenters and b + as seeds. As K-means is usually employed with unlabeled data, the seeds are randomly selected unlike our approach. In this case, our approach and its evaluation is not dependent on a specific initial configuration. 8) since a new barycenter is added, a new cluster of the same type (real(b + )) is created. The set C is updated with the constructed clusters by k-means 9) go to step 1 The testing stage is based on the final set of barycenters B. In fact, each new point to identify is affected to the type corresponding to the closest barycenter (figure 4).

Applying a projection function, like SVM, helps to increase the separability of the points and so to reduce the number of TRACE iterations and the set of final barycenters, which speeds up the classification during the testing stage. This algorithm, K-TRACE, is a kernel-based method and is depicted in figure 4. Therefore, the algorithm is exactly the same than TRACE but it is executed on the projected data points φ(t i ). Similarly to SVM, the definition of this function is not straightforward and justify the use of the kernel trick: K(t i , t j ) = φ(t i ).φ ( t j ) . However, the Euclidian distance normally based on φ(t i ) has to be computed as follows between t i and the barycenter b l k :

d(ti, b l k ) = K(ti, ti)- 2 |β l k | x∈β l k K(ti, x)+ 2 |β l k | 2 x∈β l k y∈β l k K(x, y) (3)
with:

β l k = {t i |t i ∈ c l k } (4)
Identifying a new device representation needs to compute its distance from all barycenters and so from all other data points from the training set as shown in equation ( 3). Thus, the complexity is very high whereas the goal is to speed up the classification. Consequently, we propose a simple adaptation to keep the original interest of TRACE by representing each subcluster by an unique point. In fact, the closest point of each computer barycenters during the training stage is kept. They form the set B ′ .Then, the distances between a new tree, t x , to classify and each of them, t j , is computed in the higher dimensional using the kernel function similar to [START_REF] Cifarelli | A kernel k-means procedure for classification[END_REF]:

d(t x , t j ) = K(t x , t i ) + K(t j , t j ) -2K(t x , t j ) (5)
Unlike equation ( 3), the calculation does not require to compute the kernel functions with all other points of the cluster.

C. Confidence level

Security applications usually rely on confidence levels in order to trigger correct counter-measures without blocking too many normal connections (false positives). Thus, we propose in this section a way to determine the confidence level of the predicted type.

Assuming the TR-FSM t x , it is assigned to the type d = assigned(t x ) based on several decision functions using SVM. In fact, the sign of f lk (t x ) indicates if the device t x belongs to the type l or k. Then the class which is mainly chosen is considered as the good one. Considering the distribution

V (t x = d p ) = |{f dp , f dp (t x ) > 0}| |{f lk }| , the identification result d corresponds to arg max k V (t x = d p ).
In fact, the distribution represents the proportion of voters for each possible device type. Hence, the value V (t x = d) helps to determine the confidence level. The result should have more impact if the vote is not floating. A floating vote result is equivalent to a low standard deviation. Considering σ V as the standard deviation of the distribution V , the confidence level is defined as follows:

c(t x ) = α × V (t x = d) + β × σ V (6)
α and β can be tuned. In this paper, they are fixed to 0.5 for computing the usual mean. Regarding K-TRACE, the higher the difference between the distance from the closest barycenter to the point to assign and the distance from others, the higher the confidence should be. So, to keep the same definition as prior, the distances are normalized between 0 and 1 and the following distribution is defined:

V (t x = d p ) = 1 -min l d(t x , b l p ) b l k 1 -d(t x , b l k )
This refined distribution definition is then used in [START_REF] Franc ¸ois | PTF: Passive Temporal Fingerprinting[END_REF].

IV. EVALUATION

The experiments are based on SIP [START_REF] Rosenberg | SIP: Session Initiation Protocol[END_REF] which is well used by VoIP operators since it is designed for managing multimedia sessions. Due to its popularity, a broad range of VoIP threats appeared last years [START_REF] Keromytis | Voice over IP: Risks, threats, and vulnerabilities[END_REF] whereas it becomes one major service in the Internet. Therefore, designing efficient solutions for VoIP service is necessary but the approach described in this paper remains generic.

A. Accuracy assessment

Our evaluation is based on the same datasets than in [START_REF] Franc ¸ois | PTF: Passive Temporal Fingerprinting[END_REF] (table I)) with network traces from our own testbed and also from a real VoIP operator (T1, T2) that we assume originally free of attacks. 40% of the data were used for training the system with the testbed dataset because of its relative small size, although we have evaluated our system in a critical way when using real data as only 10% have been used during the training. In addition, each experiment is run ten times where message sequences are shuffled prior.

General statistics about the datasets and fingerprinting accuracy are given in table I. Regarding our context, the accuracy is the proportion of devices which are assigned to the correct types. Assuming a set P of N predictions represented as pred ij meaning that a prediction of type j was made for a device of real type i, the accuracy is:

predij ∈P,i==j 1 N Performances are better with the testbed dataset because the learning is executed on a bigger proportion of data. The operator dataset includes some device types under represented, which leads to degrade the results. Moreover, K-TRACE is less accurate but is still acceptable since the accuracy varies between 0.73 and 0.84 depending on the dataset. The standard deviation highlights that results are quite stable. Hence, as it is claimed in this paper, this proves that fingerprinting can enforce security applications but have to be carefully taken in account due to the errors entailed by classification techniques. Additional results about SVM and the tuning of parameters (especially the number of message sequences to use) can be found in the original paper [START_REF] Franc ¸ois | PTF: Passive Temporal Fingerprinting[END_REF].

B. Calibration

The following experiment evaluates the confidence level described in section III-C in two ways: when the identification is right or when it is wrong. The goal is to calibrate the system and to define a threshold τ evaluating the correctness of the identification: a confidence level lower than τ indicates that the result needs to be considered carefully. Due to space constraint, this section presents SVM-based results but conclusions are identical for K-TRACE.

After running multiple experiments and calculating the average confidence levels for right and wrong results, the minimal and maximal values can be computed as shown in table II. It shows that the minimal average confidence for a correct identification is always greater than the maximal average confidence of wrong identifications. The confidence level of the right and wrong identifications have close average confidence levels (around 0.56665) for the testbed dataset and with a standard deviation of 0.00005. This experiment helps to determine the threshold τ which can be set for example to the average between the two values represented. For example, τ = 0.5456 is suitable for the dataset T1. In this case, only right identifications are kept. Hence, the proposed confidence level definition clearly improves the viability of fingerprinting for supporting security applications. Choosing always a high threshold, close to the maximal observed values for the correct identifications, is another possible choice to be very cautious.

C. Authentication reinforcement

Attackers use specific tools not used by legitimate users. Hence, they masquerade as legitimate device types. Figure 5 summarizes the functioning of the detection. In fact, each device announces its type (user-agent) which is compared to the fingerprinting result. However, due to fingerprinting errors, the confidence level is taken in account. Hence, when the confidence level is high, a similar announced and fingerprinted type implies no alert (level 0) whereas a divergence is synonym of an intruder pretending to be another kind of device (maximal alert level). When the confidence level is not high, the decision strength is lower and then alert level is one when the identity does not diverge and two otherwise.

Attack traffic was generated on the testbed using the bot presented in [START_REF] Nassar | VoIP Malware: Attack Tool & Attack Scenarios[END_REF] generating phone calls. The bot pretends to belong a known device type from the training set. Table III highlights the statistics from the attack traffic generated. The traffic from the bots is mixed with 60% of the normal traffic.

τ is set to 0.5667 as it is an intermediary value highlighted on table III. The distribution of benign and attack traffic regarding the alerts is represented on figure 6. The benign traffic entails many alerts of level zero and one whereas the attack traffic is split into level two and three. Therefore, the attacks are clearly distinguishable even if ten percents of the benign traffic is about equally divided into alerts level two and three. However, no attack has an alert level lower than two.

The ROC (Receiver Operating Characteristic) curve is plotted in figure 7 by varying τ . Only the identifications having a confidence higher than τ are plotted (true or false positives). The scale highlights a huge increase of the true positives rates comparing with the false positives since detecting 100% of the attacks implies ten percents of false positives. Thus, if the system triggers an additional check for a priori malicious users, ten percents of normal users will be affected. Using K-TRACE, the results are affected in a similar way than in section IV-A.

Training the system with attack tools implies that 99.1% of them are directly detected afterwards with a high confidence value (higher than 0.5667 for SVM, which is the previously determined value) whereas the normal traffic is never considered as an attack.

D. Attack detection

Automated learning of stateful firewall rules is one of the main application of the attack module of our behavioral security framework. In fact, stateful firewall rules can be directly generated from the TR-FSM of an attacking device. The attack is run multiple times in order to collect several possible sessions which are aggregated into a TR-FSM. Figure 8(a) represents the results from a denial of service attack reported in CVE (Common Vulnerabilities and Exposures) 2007-4459 against a Cisco 7940. The attack uses only ten messages and the device responses are often errors (400, 404) due to the invalidity of the requests. Each individual message is legitimate and thus a simple signature detection will not discard it. Therefore, only the global sequence of sent messages can achieve the attack. Then, the TR-FSM can be considered as a stateful rule for a firewall but it can be restricted to messages from the attacker (prefixed by !) for discarding the normal device behavior. For example, if an attacker tries this attack on many devices without knowing their real types, the responses sent can differ but not the attack messages. Moreover, building several TR-FSM with different devices allows to construct a more complete stateful rule. Thus, by removing responses from figure 8(a), all branches are identical and composed of INVITE and OPTIONS in the same order leading to aggregate them into a single sequence based signature depicted in figure 8(b). This rule is injected into SecSip, a stateful firewall [START_REF] Lahmadi | SecSip: A Stateful Firewall for SIP-based Networks[END_REF] based on a new language, VeTo [START_REF] Lahmadi | Veto: An exploit prevention language from known vulnerabilities in sip services[END_REF], for tracking SIP device behavior along its lifetime. The first step defines the events ev_Invite and ev_Options corresponding to the different kinds of emitted messages (INVITE and OPTIONS):

definition SIPMessages begin when SIP:request.method @ match "ˆINVITE" -> let: event ev_Invite; when SIP:request.method @ match "ˆOPTIONS" -> let: event ev_Options; definition end Then, the attack is described as a sequence of events (regular expressions are authorized) equivalent to the messages in figure 8 Because responses are not defined as events, they are ignored by SecSip. Hence, when the previously defined sequence of events is detected, the message corresponding to the last event is dropped to stop the attack. For improving the detection, the delays can also be included in VeTo rules. Including delays in stateful rules is a user choice but can be clearly an added value for a certain kind of attack like flooding attacks as SpIT (Spam over Internet Telephony) can be.

Obviously, this attack is a simple example to illustrate the process and the rule could be created easily by hand but more complex attacks can include many messages sequences (alternative paths in the tree) for which an automatic tool such as proposed in this section is required.

V. REAL-TIME APPLICABILITY

For assessing the real-time applicability, this section includes analytical facts as well as experimental results based on a single computer (Intel Core 2 Duo 3GHz).

A. Message parsing

The first task is to determine the types of the messages when building a TR-FSM. Thus, the complexity is highly dependent on the protocol grammar and a straight-forward manner is to parse entirely the message to construct its syntactic tree as also used for device fingerprinting in [START_REF] Franc ¸ois | Machine Learning Techniques for Passive Network Inventory[END_REF]. However, the type is generally a field which is identified by a certain keyword or a certain position and so it can be easily found without constructing the complete syntactic tree: direct access to the position in the message or one of the first parsed fields as for instance in SIP. This implies a constant complexity since the type is given by few bytes encapsulating in a small range of the first bytes. In contrast, a syntactic parser able to deal with huge grammars presenting potential ambiguities [START_REF] Tomita | Lr parsers for natural languages[END_REF], [START_REF] Tomita | Generalized LR parsing[END_REF] is known to have a linear complexity in most cases.

We present results of evaluating the time to parse the messages and to retrieve the necessary information in figure 9. For most cases, the time function follows the theoretical complexity mentioned previously but highlights an incompressible time equivalent to 13,1 milliseconds (intercept) if the complete syntactic tree is built, although that the time for a simple parsing is well concentrated around the median value (1.097 × 10 -5 seconds) as shown in table IV. Hence, syntactic fingerprinting as proposed in [START_REF] Franc ¸ois | Machine Learning Techniques for Passive Network Inventory[END_REF] is not suited to the context of high speed monitoring.

B. TR-FSM construction

Constructing a TR-FSM is divided in two times: (1) message sequences must be monitored and (2) have to be aggregated to build a TR-FSM. The running time of the first step cannot be avoided but also dependent on the network throughput and the device activity. So, associating a complexity to this one is meaningless. However, the second stage has to agglomerate such messages into a TR-FSM which implies a linear complexity regarding the total number of messages composing the tree [START_REF] Franc ¸ois | PTF: Passive Temporal Fingerprinting[END_REF]. This is checked in table IV. Fortunately, the entire TR-FSM construction is very fast with a median value equals 0.024 ms (table IV).

C. Machine learning

By design, the kernel function computation requires to compute the intersection between paths of the trees to compare. From an analytical complexity point of view, this is equivalent to compare syntactic trees [START_REF] Franc ¸ois | Machine Learning Techniques for Passive Network Inventory[END_REF] except that timing differences are also computed, i.e. the complexity is quadratic with respect to the number of nodes. This is confirmed in figure 10 and in table IV.

K-TRACE and SVM both iterate multiple times to find a way to separate data points. Hence, the computation time of the training stage is highly dependent on the data (separability of data points) and the initialization of the SVM related optimization problem. Furthermore, evaluating precisely the complexity of this stage is not the most important since it can be done offline and only the identification itself (testing stage) has to be done online. However, our experiments show that the learning time is not the bottleneck (maximum one hour for a daily update).

The main advantage of the learning techniques is to extract a subset of training data points which will be then compared

reduction f actor =      N train N SV if SV M N train N bary if K -T RACE
The less efficient the training is, the lower this factor is (minimum is one). This factor is plotted in figure 11 and increases while the learning percentage also. Actually, when trees are added to the training samples, the probability to have redundant information, i.e., similar trees, increases also. Thus, more trees will be discarded and the reduction factor is better. Figure 11 highlights the great efficiency of the K-TRACE method comparing with SVM. Thus, K-TRACE is more capable to discard redundant information when the training set is expanded although SVM kept most of points as shown by a reduction factor close to 1.

D. Identification time

Assuming t parsing , the time for parsing a message, t build , the time for constructing a TR-FSM and t comp the time for comparing two TR-FSMs, the total execution time of one online identification is:

t = t parsing × card + t build + t comp × N samples
where card is the average number of nodes in a TR-FSM and is varying depending on monitored sessions. Hence, the global time corresponds to parse the different messages, to construct the TR-FSM and finally to compare this tree with Worst case Fig. 12: Online identification time previously selected training trees. In order to keep pace with the network throughput, a device has to be identified before the next identification starts which means that t ≤ inter packet time × card where inter packet time is the inter arrival time between two packets. Hence, we obtain:

t parsing + t build card + t comp card × N samples ≤ inter packet time
The different terms of this equation depend on the protocol and the user activities. The operator and the testbed network are similar in terms of TR-FSM sizes and message lengths (table I) which are the main factors impacting the computation time. Hence, regression functions in table IV can be considered for extrapolating the results to the operator dataset. Therefore, figure 12 displays the global computation time depending on the number of selected TR-FSM (N SV or N bary ) and the average value for card (table IV) together with the worst case when card = 1. The operator dataset (T1) has an average throughput of 110 messages per second equivalent to 9 milliseconds between two messages and so excludes the SVM method because the training process selects more than 10,000 N SV (support vectors) whereas K-TRACE computes about 110 distinct N bary (barycenters) which is acceptable (figure 12). In this case, fingerprinting requires around 0.2ms per packet which corresponds to a maximal network througput of 5000 packets per second. Considering the statistics in table I, this is also equivalent to 2.2 GB/second. Assuming our configuration and the same reasoning, SVM, which provides the best accuracy, is only viable if the througput for SIP is lower than 21 MB/second.

VI. RELATED WORK

There are different kinds of fingerprinting methods. Passive techniques only monitor the devices like p0F [START_REF] Zalewski | P0f[END_REF] which aims at finding the operating system. NMAP [START_REF] Lyon | Nmap Network Scanning: The Official Nmap Project Guide to Network Discovery and Security Scanning[END_REF] uses active techniques i.e., probing the devices with specific requests. Although there are different levels of granularity (operating systems: [START_REF] Zalewski | P0f[END_REF], [START_REF] Lyon | Nmap Network Scanning: The Official Nmap Project Guide to Network Discovery and Security Scanning[END_REF]; protocol: [START_REF] Ma | Unexpected means of protocol inference[END_REF]; device type: [START_REF] Scholz | SIP Stack Fingerprinting and Stack Difference Attacks[END_REF], [START_REF] Franc ¸ois | Automated behavioral fingerprinting[END_REF]), none of them are leveraged for advanced security purposes especially device type fingerprinting. The behavioral fingerprinting was introduced in [START_REF] Franc ¸ois | Automated behavioral fingerprinting[END_REF], [START_REF] Franc ¸ois | PTF: Passive Temporal Fingerprinting[END_REF] and does not require a full syntactic [START_REF] Franc ¸ois | Machine Learning Techniques for Passive Network Inventory[END_REF] or semantic [START_REF] Franc ¸ois | Digital Forensics in VoIP networks[END_REF] knowledge about the protocol.

As argued in this paper behavioral fingerprinting can be easily extended to security purposes. IDS research has skyrocketed in parallel of the Internet. Different surveys, as for example [START_REF] Patcha | An overview of anomaly detection techniques: Existing solutions and latest technological trends[END_REF], highlight the two main classes of IDS. Signature based IDSs such as Snort [START_REF] Roesch | Snort -lightweight intrusion detection for networks[END_REF] identify an attack based on a signature. Anomaly based IDSs detect an attack based on its deviation from a normal profile. Such systems have a higher ability to detect unknown attacks but the number of false positives is higher. Hence, the trend is to use more and more sophisticated techniques as for example [START_REF] Mahoney | Learning rules for anomaly detection of hostile network traffic[END_REF].

Obviously, the idea of modeling behavior was explored in the past. User behavior was generally modeled as its activity in a system i.e., the executed commands. For example, the authentication reinforcement presented in this paper is close to masquerade detection as for example by analyzing UNIX command sequences [START_REF] Seo | Masquerade detection based on svm and sequencebased user commands profile[END_REF]. System calls could be also analyzed to detect anomalies [START_REF] Cha | Anomaly intrusion detection for system call using the soundex algorithm and neural networks[END_REF].

The authors in [START_REF] Sekar | Specification-based anomaly detection: a new approach for detecting network intrusions[END_REF] extend this concept to network protocol state machine by combining manual building of a subpart of the protocol state machine and by learning statistical properties about transitions. However, such an approach needs to know and manually analyze each protocol specification. NetSTAT [START_REF] Vigna | Netstat: a network-based intrusion detection system[END_REF] also employed state machine for detecting intrusion but needs to manually define prior scenario. A close work is [START_REF] Sengar | Voip intrusion detection through interacting protocol state machines[END_REF] which did a similar analysis for VoIP intrusion detection but relies on fully described attacks unlike ours which automatically builds the signatures of the attacks for any protocol. Finally, the use of Markov chain was also explored in order to model the transitions, between the messages, labeled with probabilities as for example in [START_REF] Estevez-Tapiador | Stochastic protocol modeling for anomaly based network intrusion detection[END_REF].

Therefore, our framework combines state machine based IDS and machine learning based IDS (no manual analysis). Actually, the TR-FSM is a specific state machine representation built by monitoring the network with a limited computational overhead [START_REF] Franc ¸ois | PTF: Passive Temporal Fingerprinting[END_REF]. This is also due to a limited parsing of messages unlike deep packet inspection. Moreover, retrieving such representation does not need a full specification knowledge [START_REF] Franc ¸ois | Automated behavioral fingerprinting[END_REF]. Thus, this paper demonstrates the capacity of fingerprinting regarding different security use-cases and proves its online applicability.

VII. CONCLUSION

We extend the original model of behavioral fingerprinting to being able to evaluate the confidence level and experimental results show its interest in the security context. We also shown that stateful firewall rules can be automatically inferred. Furthermore, the device fingerprinting in real time was evaluated and experimented practically with VoIP networks. By proposing a new classification algorithm, we have shown that monitoring high speed network is feasible. Our future work will focus on improving the firewall rule generation by merging and reducing the TR-FSM sizes.

Fig. 1 :

 1 Fig. 1: Behavioral security framework

Fig. 2 :Fig. 3 :

 23 Fig. 2: TR-FSM examples. Two shared paths are grey colored

Fig. 5 :

 5 Fig. 5: Intrusion detection based on fingerprinting

Fig. 6 :Fig. 7 :

 67 Fig. 6: Traffic distribution regarding the alert levels

Fig. 9 :Fig. 10 :Fig. 8 :

 9108 Fig. 9: Computation time for recovering the message typ regarding the length of the message

  (b): veto attack_cisco uses SIPMessages begin (ev_Invite,ev_Options[ * 2],ev_Invite,ev_Options[ * 3], ev_Invite,ev_Options[ * 2]) -> drop; veto end

Fig. 11 :

 11 Fig. 11: Reduction factor with different training set size

TABLE I :

 I Experimental datasets fingerprinting results

		Testbed	T1	T2
	Min. confidence of right identification	0.56677	0.54706	0.54580
	Max. confidence of wrong identification	0.56658	0.54410	0.54153

TABLE II :

 II Maximal and minimal confidence levels

TABLE III :

 III Attack traffic statistics

  = 9.911 × 10 -7 x + 1.478 × 10 -06 Tree comparison 1.097 × 10 -5 2.694 × 10 -5 5.198 × 10 -5 O(#nodes × #nodes) y = 5.348 × 10 -7 x 2 + 2.033 × 10 -05

	Stage	Execution time (seconds) 25% percentile median 75% percentile	complexity	regression
	Parsing	1.001 × 10 -5	1.097 × 10 -5	1.192 × 10 -5	Constant	y = 1.087 × 10 -5
	Tree construction	1.597 × 10 -5	2.384 × 10 -5	3.409 × 10 -5	O(#messages)	y

TABLE IV :

 IV Performance evaluation

ACKNOWLEDGMENT

This work is partly funded by OUTSMART, a European FP7 project under the Future Internet Private Public Partnership programme. It is also supported by MOVE, a CORE project funded by FNR in Luxembourg.

Note: figures include wikimedia contents (visit http://commons.wikimedia.org for license information).


G Morra

D A Yuen

L Boschi

P Chatelain

P Koumoutsakos

P J Tackley

Morra et al. simulated the fate of subducting plates interacting with the mid-lower mantle with variable viscosity

Keywords: Subduction, slab-mantle interaction, spin transition, mid-mantle viscosity

A c c e p t e d M a n u s c r i p t 2

Abstract

In the last two decades it has been proposed several times that a non-monotonic profile might fit the average lower mantle radial viscosity. Most proposed profiles consist in a more or less broad viscosity hill in the middle of the mantle, at a depth roughly between 1,200 km and 2,000 km. Also many tomographic models display strong signals of the presence of "fast" material lying at mid mantle depths and a recent spectral analysis of seismic tomography shows a very clear transition for degree up to around 16 at a less than 1,500 km depth. Finally latest works, both theoretical and experimental, on the high-to-low spin transition for periclase, have suggested that the high-spin to low-spin transition of Fe++ might lie at the heart of all these observations. To verify the dynamical compatibility between possible mantle profile and observed tomographic images and compare them with possible mineral physics scenarios, such as the spin transition, we employ here a recently developed Fast Multipole-accelerated Boundary Element Method (FMM-BEM), a numerical approach for solving the viscous momentum equation in a global spherical setting, for simulating the interaction of an individual slab with a mid-mantle smooth discontinuity in density and viscosity. We have focused on the complexities induced to the behaviour of average and very large plates O (2,000 km -10,000 km), characteristic of the Farallon, Tethys and Pacific plate subducting during the Cenozoic, demonstrating that the a mid mantle density and/or viscosity discontinuity produces a strong alteration of the sinking velocity and an intricate set of slab morphologies. We also employ the Kula-Farallon plate system subducting at 60 Ma as a paradigmatic case, which reveals the best high-resolution tomography models and clearly suggests an interaction with a strong and/or denser layer in the mantle. Our 38 models show that a plate might or might not penetrate into the lowest mantle and might stall in the mid lower mantle for long periods, depending on the radial profiles of density and viscosity, within a realistic range (viscosity 1, 10 or 100 times more viscous of the rest of the mantle, and a change of differential density in the range -2% to 2%), of a transitional layer of 200 km or 500 km. We conclude that a layer with high viscosity or negative density would naturally trigger the observed geodynamic snapshot. We finally propose a scenario in which the long time accumulation of depleted slabs in the mid-mantle would give rise to a partially A c c e p t e d M a n u s c r i p t

Introduction

The radial profile of the lower mantle viscosity is still largely unknown due to incomplete knowledge of composition, grain size, volatiles content, and to the lack of laboratory experimental data at such high pressure and temperatures and the difficulty of performing exhaustive ab-initio models of diffusion and dislocation creep for polycrystalline materials. However, geophysical inversion data based on geoid and post-glacial rebound [START_REF] Forte | Deep-mantle high-viscosity flow and thermochemical structure inferred from seismic and geodynamic data[END_REF][START_REF] Mitrovica | Radial profile of mantle viscosity: Results from the joint inversion of convection and postglacial rebound observables[END_REF]) has suggested at the first order a peak of high viscosity in the middle of the lower mantle at around 2000 km and in a range between 10 and 100 times than the one in the upper mantle.

Few viscosity radial profiles for the lower mantle have been put forward. Among them [START_REF] Wuming | Inferring the mantle Viscosity and Its Three Dimensional Structure From Geoid ,Topography and Plate Velocity[END_REF] proposed that a peak of viscosity in the mid mantle might explain a set of geoid and topographical data. From geophysical considerations [START_REF] Wen | Slabs, hotspots, cratons and mantle convection revealed from residual seismic tomography in the upper mantle[END_REF] imposed a chemical barrier at 1200 km. [START_REF] Kellogg | Compositional Stratification in the Deep Mantle[END_REF] introduced the concept of a barrier at 1800 km depth based on geochemical and radiogenic heat production arguments, combined with tomographical evidences. Later a number of works came to similar conclusions, for example [START_REF] Forte | Deep-mantle high-viscosity flow and thermochemical structure inferred from seismic and geodynamic data[END_REF] who combined satellite and ground data, [START_REF] Ito | Pressure effect of self-diffusion in periclase (MgO) by molecular dynamics[END_REF] have found a peak of activation energy for vacancy diffusion at mid mantle pressures (both experimentally and numerically from large-scale molecular dynamics simulations) and Peltier (personal comm., 2008) looked at the GRACE dataset and found that the hypothesis of the viscosity hill in the mid mantle is still compatible with the most up-to-data gravity data, although this presents a controversy, as other studies, e.g., [START_REF] Soldati | Inferring radial models of mantle viscosity from gravity (GRACE) data and an evolutionary algorithm[END_REF] indicated that lower mantle viscosity structure can hardly be constrained from gravity data.

Tomographical Evidences

A clearer picture of how anomalous the mid-lower mantle appears from mantle tomography. One of the better seismically resolved areas of the Earth's mantle is under North America, where the Farallon slab is arguably the best-known sudbuction system characterized by likely penetration into the lower mantle, e.g. [START_REF] Grand | Mantle shear structure beneath the Americas and the surrounding oceans[END_REF]. The set of A c c e p t e d M a n u s c r i p t panels on the left of Figure 1 illustrates the paleo plate reconstructions for the period from 80 Ma to 40 Ma [START_REF] Müller | Age, spreading rates and spreading asymmetry of the world's ocean crust[END_REF] of the splitting of the wide oceanic plate with over 8,000 km width subducting under North-America into two subplates called Farallon and Kula. A recent high resolution tomography [START_REF] Sigloch | Two-stage subduction history under North America inferred from multiple-frequency tomography[END_REF], reproduced in figure 1, right panel, offers a very detailed image of the Farallon slab at mid mantle depth. The most prominent features are the thickening of the slab at depth between 1,200 km and 1,500 km and the segmentation into two plates, most probably the signature of the division of the paleo Farallon-Kula plate, also confirmed by the combination of paleo-tectonic reconstruction and global tomography, as in [START_REF] Ren | Understanding seismic heterogeneities in the lower mantle beneath the Americas from seismic tomography and plate tectonic history[END_REF].

Such anomalous mid-lower mantle behaviour is confirmed by the visual inspection of a global tomography as the the one shown in Figure 2 where the isosurface corresponding 0.8% heterogeneity, associated with the shear-velocity model of [START_REF] Grand | Mantle shear structure beneath the Americas and the surrounding oceans[END_REF]) is originally displayed with the Paraview software [START_REF] Henderson | ParaView Guide, A Parallel Visualization Application[END_REF]. While also here the intriguing pattern of the Farallon slab is visible, with its shape that largely broadens at 1,500-2,000 km depths, other similar patterns appear in other regions of the lower mantle. Finally the picture shown in figure 1 and 2 for the Farallon slab is shown to be robust and independent on the tomographic models assumption through the comparative analysis of five independent global shear velocity models in Figure 3: tx2007 [START_REF] Simmons | Constraining mantle flow with seismic and geodynamic data: A joint approach[END_REF], pri-s05 [START_REF] Montelli | A catalog of deep mantle plumes: New results from finite-frequency tomography[END_REF], rmsl-s06 [START_REF] Li | A new global model for $P$ wavespeed variations in Earth's mantle[END_REF], saw642an [START_REF] Panning | A three dimensional radially anisotropic model of shear velocity in the whole mantle[END_REF] and smean [START_REF] Becker | A comparison of tomographic and geodynamic mantle models[END_REF]. All such models convey very clearly the picture of a Farallon slab that does not cross straightaway in the lower mantle, but instead flattens at depth of 1500-2000 km, as it would encounter an obstacle.

The 3-D seismic structure of the lower mantle can be also analyzed globally, looking at the spectral signal of all the slabs together. Detailed analysis of spherical harmonics results as the one shown in the Figure 4 (modified from figure 1 of [START_REF] Boschi | On the statistical significance of correlations between synthetic mantle plumes and tomographic models[END_REF] illustrates an unexpected non-monotonic radial pattern. The logarithm of the ratio of positive-to-negative shear-velocity spectra from model SMEAN [START_REF] Boschi | On the statistical significance of correlations between synthetic mantle plumes and tomographic models[END_REF] recently also observed in [START_REF] Houser | The relative wavelengths of fast and slow velocity anomalies in the lower mantle: Contrary to the expectations of dynamics?[END_REF] elucidates, in fact, how the positive (blue) fast anomalies dominate until the depth of ~1500km up to spherical harmonic degree 15, while the situation is reversed below where negative (red) values A c c e p t e d M a n u s c r i p t 6 dominate at lower depths. This unquestionable tomographical observation is the premise for searching the causes and consequences of the phenomenon.

Mineral physics insights

The dichotomy between the deepest 1000km and the rest of the lower mantle has been already proposed to be of compositional nature [START_REF] Van Der Hilst | Compositional heterogeneity in the bottom 1000 km of Earth's mantle: towards a hybrid convection model[END_REF].

However this picture has been enriched by recent works that combine the interpretation of seismic data with the knowledge of the elastic properties of silicates whose spin varies at mid mantle conditions [START_REF] Da Silva | The composition and geotherm of the lower mantle: constraints from the elasticity of silicate perovskite[END_REF]. Such novel works depart from the most recent experiments on ferropericlase during its spin transition [START_REF] Marquardt | Elastic Shear Anisotropy of Ferropericlase in Earth's Lower Mantle[END_REF] and exploit the surprisingly transient and anomalous character change of the bulk modulus during the transition, that largely softens. Because the observed pattern of vertical gradient of bulk Velocity (dVp/dz) is instead smooth and apparently adiabatic [START_REF] Cammarano | Inferring upper-mantle temperatures from seismic velocities[END_REF][START_REF] Matas | On the bulk composition of the lower mantle: predictions and limitations from generalized inversion of radial seismic profiles[END_REF], the only explanation of this phenomena is a mid-lower mantle thermo-chemical transition at the same depth of the spin transition in order to compensate for the bulk-modulus weakening [START_REF] Cammarano | Role of iron-spin transition in ferropericlase on seismic interpretation: A broad thermochemical transition in the mid mantle?[END_REF].

If this hypothesis proves to be correct, this groundbreaking interpretation would then demand a thorough investigation of its geodynamical consequences and above all, of the conditions that might have produced it. Such a mid-mantle transition is clearly not as sharp as the upper lower mantle one but it is believed to be smooth. The inversion of the laboratory data indicates a thickness of several hundred kms [START_REF] Marquardt | Elastic Shear Anisotropy of Ferropericlase in Earth's Lower Mantle[END_REF], which is less than the prior spin transition thickness i.e. up to 1000 km [START_REF] Da Silva | The composition and geotherm of the lower mantle: constraints from the elasticity of silicate perovskite[END_REF]. The combination of seismic and mineralogical investigations strongly suggests that the thermo-chemical should happen at depths around 1200-1600 km [START_REF] Cammarano | Role of iron-spin transition in ferropericlase on seismic interpretation: A broad thermochemical transition in the mid mantle?[END_REF][START_REF] Marquardt | Elastic Shear Anisotropy of Ferropericlase in Earth's Lower Mantle[END_REF].

All these indications are converging toward a scenario in which the viscosity profiles in the lower mantle might be non-monotonic, which would be entirely different from the canonical picture one gets from an Arrhenius type of activation energy, while it would be better described by compositional differentiation. A fundamental point to clear is the origin of the hill shape of the mantle viscosity arising from some inversion.

If the rising of viscosity can be associated with the differentiated composition, i.e. with more viscous rocks at higher depth, the following lowering viscosity would be simply With the exception of the work by [START_REF] Wen | Slabs, hotspots, cratons and mantle convection revealed from residual seismic tomography in the upper mantle[END_REF], who postulated a chemical layer at around 1200 km depth based on geoid inversion, this is probably the first attempt to fit the observed mantle observable with an explicit compositional well defined radial transition. The details of the generation of such a three-layered mantle require careful study and are still speculative, however the strength of the tomographic indications listed above are so neat that require models able to explain them. Yet there is only one conclusive evidence of the existence of such dramatic alteration, i.e. [START_REF] Cammarano | Role of iron-spin transition in ferropericlase on seismic interpretation: A broad thermochemical transition in the mid mantle?[END_REF] and more independent studies are required to confirm this observation. Finally definitive estimates of this new radial profile have not been presented.

Geodynamic scenario and our parameterization

If the above mineral-physics and seismological interpretations are both correct, it remains to be explained the geodynamic origin of the lower mantle transition and what is its impact for rising plumes and sinking slabs. This work aims at understanding just the fate of a slab in the lower mantle. The novel geodynamic scenario here suggested is the one in which a sinking slab must cross three distinct mantle layers in its way to the core: the upper mantle, the external-lower mantle (i.e. above 1500-2000 km) and the interior-lower mantle (i.e. below 1500-2000 km). In particular, we propose that the viscosity hill often suggested in the literature was indicating such a mid-lower mantle viscosity transition.

We consider 4 fundamental parameters influencing the fate of the slab. Other factors taken fixed in all models are the thickness of the plate, set to 200 km, which should account for the broadening that the plate undergoes due to diffusion mechanisms and maybe thickening in the early stages of subduction. The plate starts each model from an initial straight shape, inclined of 20 degrees from the vertical, just above the hill.

The first examined parameter is the differential density between slabs and mantle.

Because the slab composition cannot change in its way to the core (except some devolatilization phenomena, mostly in the upper mantle), and the spin transition A c c e p t e d M a n u s c r i p t 8 happens in the slab at shifted depths compared to the rest of the mantle, due to lithosphere's cooler temperature, we consider the entire range of possible differential densities, 2%, 0% and -2%, from the lower mantle during the transition. Before and after the mid-mantle transition we assume instead unvaried 2% negative buoyancy for the slab (Figure 5). Real values can be different for each real case, but our end-member scenarios should encompass all possible events.

The second parameter is mantle viscosity. We investigate a non-homogeneous mantle with a mid mantle peak viscosity 1, 10 and 100 times of the far field mantle, as they are end members and are of lower or equal order of magnitude to the proposed mantle viscosity hills (Figure 5). We do not vary the mantle-lithosphere viscosity ratio.

Such ratio has been amply studied in most geodynamic simulations of subduction that give a value of about 100 [START_REF] Capitanio | Dynamic models of downgoing platebuoyancy driven subduction: Subduction motions and energy dissipation[END_REF] in the upper mantle. In the lower mantle such a value is typically estimated as being much lower. We assume therefore a slab with a viscosity ratio of 10 between slab and mantle. Interesting tests of how a homogenously viscous slab might evolve in the lower mantle have , PhD Thesis) showing a complexity of shapes, also seismically detected and interpreted [START_REF] Ren | Understanding seismic heterogeneities in the lower mantle beneath the Americas from seismic tomography and plate tectonic history[END_REF].

In the technical implementation of the viscosity hill, as explained in the appendix, only the high shear stresses due to the highly viscous mantle have been implemented.

Several reasons justify this assumption: a) the spin transition should take place at different depths inside and outside the slab, due to the differential temperature between slab and mantle; b) we follow the conceptual model that the viscosity rises in the mantle due to a compositional variation, not directly due to the spin change (however, this is not clear yet); c) we are interested in the fate of the slab in the particular context of their ability to cross a mid mantle viscosity peak and reach the core-mantle boundary contributing to the D'' composition or to a more or less thick differentiated mid mantle layer, and this is due more to the mantle shear stresses applied to the slab then to the slab strength itself. In other words, we are interested only in the case in which the lithosphere remains a coherent entity, while we can only detect, and not directly model, the case in which it becomes unstable, thus producing fragments that drop in the underlying mantle. The last two parameters are geometrical, i.e. the thickness of the transition, in which two end-members cases are tested: a thin 200 km and a thick 500 km. and two different slab width cases: very wide, 10,000 km and relatively narrow, 2,000 km, corresponding respectively to the largest plates on Earth and to a typical subduction zone.

The total of 36 possible combinations gives origin to a very large set of possible outcomes that we summarize in three categories: (1) Stokes flow, for slabs that directly penetrate through the transition and do not feel the presence of mantle peak effects, (2) transient, for slab that sensibly slow down at the mid mantle transition and then go through it arriving at the lowest mantle and (3) stalling for slabs that remain entrapped in the mid-mantle and never cross the transition.

Numerical Simulations

The simulation of the flow of creeping systems with large and sharp viscosity and density variations presents formidable challenges, due to the extremely complex evolving heterogeneities. We employ an indirect boundary integral formulation representing the boundary between a strong slab a less viscous mantle. The resulting formulation, based on bounded Stokes flow laws, is a set of Fredholm integral equations of the second kind. While the methods based on boundary integral equations [START_REF] Brebbia | The Boundary Element Method for Engineers[END_REF] reduce the effort necessary for representing complex geometries, their standard formulation is computationally very expensive for large problem due to the construction of dense matrices. In order to accelerate the solution, we replace the calculation of the dense matrix with the implementation of an algorithmic "matrix multiplier operator", based on the multipole expansion of the integral terms [START_REF] Greengard | A Fast Algorithm for Particle Simulations[END_REF] as detailed in the appendix.

In our numerical scheme, the slab surface is discretised into linear triangular elements. It has been shown that the linear system arising from its related equation set is well-conditioned and dense [START_REF] Zhu | A 3-D indirect boundary element method for bounded creeping flow of drops[END_REF]. The system is then solved employing the GMRES Algorithm which exploits the potentialities of the FMM fast multipole approach [START_REF] Greengard | A Fast Algorithm for Particle Simulations[END_REF]. Once the single and double layer integrals are calculated, the position of each node is updated with a simple explicit forward step in time. The heat or chemical diffusion equation is not solved in this work because it would unnecessarily complicate our procedures, but it can be straightforwardly Standard implementations of the boundary element method allow only the study of homogeneous flow as discussed in [START_REF] Morra | Large Scale Three-Dimensional Boundary Element Simulation of Subduction[END_REF] [START_REF] Morra | Earth curvature effects on subduction morphology: Modeling subduction in a spherical setting[END_REF] and illustrated in the first part of the appendix. In practical terms this means that for each element, the density difference and the viscosity ratio between slab and mantle are constant. In this work, for the first time, we employ a non homogenous viscosity and density boundary formulation, which allows us to model the interaction of the slab with a non-homogeneous mantle, characterized by non-monotonous radial density and viscosity profiles. Both requirements are exploited employing different approaches, described in more detail in the technical Appendix: a) the viscosity ratio is straightforwardly mapped on each boundary element separately, depending on the local position of the boundary, but the implementation requiring the calculation of the viscosity ratio between inner and outer slab requires a careful implementation. The approximation here employed is explained in detail in the appendix, formula (9), and has been tested with excellent results, giving only local errors in proximity to the mantle-hill zone. In practice, it corresponds to a mean-field approximation, neglecting the effects of the far field viscosity transition; b) the non-homogeneous differential density requires the calculation of a one dimensional radial finite difference implementation of the pressure, applied to the boundary element formulation, replacing the component of the stress in the Stokeslet component, as explained in recalculating the formula (10). This is done at each step, with minimum computational time consumption, by integrating over the entire mantle density profile.

Model Results

In this section we describe in detail the results of our experiments in numerical modelling; we provide in the Discussion section a geophysical explanation of our theoretical insights. Let's consider first the interaction of a small plate (2000 km x 2000 km) that encounters a purely 500 km thick viscosity hill, without considering any density variation, i.e. only the case in which Δρ hill is 2% higher in the slab compared to the surrounding mantle. The transition from the 10x to the 100x viscosity hill model is The slab in both cases crosses the viscosity hill, but the time required for going through it is one order of magnitude higher for the 100x model. Furthermore the impact with the high viscosity layer changes dramatically the morphology of the down welling, which converts from a slab like structure into a cylindrical drop, growing as Rayleigh-Taylor (RT) instability. We call the case shown above figure 6 as being Stokes sinking, while the one on the lower mantle is defined as being transient.

In order to show the effect of a density variation superposed to a viscosity hill, the following figure 7 illustrates the time sequences of two 10x viscosity models with a small 2,000 km x 2,000 km plate testing 2 different smooth density variations Δρ hill = ρ litho -ρ mantle at the peak either 0% or -2% compared to surrounding mantle. The 0% model, whose time evolution is illustrated in figure 7 (top panels), displays an intermediate outcome between the two models shown in figure 6, with the final morphology (after the hill is crossed) resembling partly such as slab and partly such as a RT instability, falling in a ring shape. The lower panel of figure 7 shows that in this simple configuration, the opposing buoyancy model (Δρ hill =-2%) is able to halt the plate at the viscosity transition. The two models here displayed are paradigmatic for a transition from a transient to a stalling case.

The figure 8 illustrates the final stage of 6 different modeled slabs, focusing on the role of differential density, summarizing their morphology. Only the mild viscosity hill (10x) is considered as they show a larger variation compared to the 100x cases. Here upper and lower panels represent the comparison between a small (2,000 km x 2,000 km, top panels) and a very wide plate (10,000 km x 2,000 km, bottom panels). Although the timing of the deformation is very different for small and large plates (see figure 9), the final configurations are similar. The most surprising result is that, while a variation from 2% to 0% triggers little effect, the negative-buoyancy (Δρ hill =-2%) model is characterized by a rounded long-term shape of the plate, with down welling at the sides of the plate itself. As it is shown in Figure 9, such a down welling arises a long time after the plate has met the mid-mantle hill. It is therefore foreseeable that for a morphological complex hill this instability might grow from segments inducing several mesoscale down welling [START_REF] Cizkova | Layered convection with an interface at a depth of 1000 km: stability and generation of slab-like downwellings[END_REF], but possibly also being an A c c e p t e d M a n u s c r i p t origin for upwelling plumes. In general, the models with inverse density hill (able to invert the sign of the slab buoyancy) display an unstable upside-down effect for which the slab splits and then fall laterally in form of drops.

All models are finally summarised in figure 9 where it is displayed the radial position of the center of mass of each plate versus a common non-dimensional time (a 2,000 km x 2,000 km plate in top panels and a 10,000 km x 2,000 km plate in the bottom panels).

The trajectories of the centre of mass for the small and the 5x wider plate display a similar pattern but the mean-velocity is about double for the wider plate (see x-axis of the upper and lower panel). This can be explained by the prediction of [START_REF] Capitanio | Dynamic models of downgoing platebuoyancy driven subduction: Subduction motions and energy dissipation[END_REF] based on the classical work of [START_REF] Happel | Low Reynolds Number Hydrodynamics[END_REF], that the drag of a plate in the mantle is proportional to [1+ln(L/S)] where L is the maximum length of the plate and S is the shortest. In our case D L =1+ln(10,000 km / 2,000 km)=1+ln(5)=2.6 while D S =1+ln(2,000 km / 2,000 km)=1+ln(1)=1. The sinking velocity is given by the ratio between plate total buoyancy (volume*Δρ) and drag, v L ~(5x1)/2.6~1.9 and v S ~(1x1)/1=1, therefore the ratio v L /v S =1.9 explains the observed sinking velocity.

The non-dimensional time employed in figure 9 generates a range of possible sinking rates that depend on the physical parameters involved in the system. In order to detect the equivalent order of magnitude for the modelled plates, we can convert the nondimensional time using the velocity renormalization factor 2/9*Δρga 2 /η and the time renormalization 9/2*η/Δρga, where η is the mantle viscosity, Δρ is the differential density, g is the gravity and a is the effective length of the slab a=S*[1+ln(L/T)], where S and L are as above, 2,000 km and 2,000-10,000 km. Assuming a lower mantle background viscosity in the range of 10 22 -3x10 22 Pa s and an average differential density in the range 1-2% (~50 kg/m 3 ) for the sinking lithosphere, one finds that a time unit is equivalent to 2 to 6 Myrs, which can be linearly renormalized for a different lower mantle viscosity. This implies that we predict a crossing time of the lower mantle of the order of 20-60 Myrs for very wide plates, as Farallon, Tethys or Pacific, and of 40-120 Myrs for small plates as single pacific slabs, both for the models that directly cross the mid mantle transition. Comparing this time with the diffusive timescale of a plate in the mantle, this could explain why only the largest plates have been detected in the deepest mantle.

A c c e p t e d M a n u s c r i p t 13 For both, wide and narrow plates, the trajectories can be grouped in three families (Figure 9 and table 1), defined roughly by the same parameter subset. The first family is specified by the cluster of trajectories that cross the hill in an almost unperturbed form and is called Stokes. For a hill of only 200 km, such models are defined by the ones that cross a small or no viscosity transition hill and have a lithosphere-mantle density difference positive or null, while for a hill thickness of 500 km already when the density difference is null the trajectory is perturbed. The second family is composed by the plates that remain permanently or for extremely long time stalled at the hill. We call this dynamics Stalling. This involves all the models with a hill-thickness of 500km and negative Δρ and high viscosity and negative Δρ (see table 2 for better visualizing the transition). Collected into a third family there is then a large set of models that are characterized by moderate anomalous viscosity and moderate density anomaly, displaying a transient behaviour. They sensibly slow around the viscosity/density hill for periods of the same order of magnitude of the ones that penetrate down to the bottom of the lower mantle and then sink down, often displaying a new morphology.

Those models are called Transient in the tables. Although their definition is here slightly arbitrary, the direct analysis of the models shows their actual existence. Their importance, as discussed in the next section, resides in the possibility that the long transient time spent around the mid mantle transition might act as a reservoir that in the life of the Earth evolution would originally cause the compositional anomaly observed in the lower mantle, as better discussed in the next session.

Discussion and new scenarios for the Earth's mantle

Convergent results taken from gravity data, and more recently mineral physics and seismology, indicate that a non-monotonic profile of viscosity/density mighty better describe the mid lower-mantle dynamics. In particular, it has been proposed that a high viscosity layer several hundred km thick might exist at depths 1,200-1,600 km. More recently a study that compares mineral physics and seismological data introduces evidences for a compositionally distinct layering below ~1,500 km depth [START_REF] Cammarano | Role of iron-spin transition in ferropericlase on seismic interpretation: A broad thermochemical transition in the mid mantle?[END_REF].

We find little difference between the dynamics of a very wide plate and a small plate, except for the sinking-velocity, which respects a logarithmic law [START_REF] Capitanio | Dynamic models of downgoing platebuoyancy driven subduction: Subduction motions and energy dissipation[END_REF] A c c e p t e d M a n u s c r i p t 14 [START_REF] Morra | Earth curvature effects on subduction morphology: Modeling subduction in a spherical setting[END_REF]. This terminal velocity diversity explains why only very large-scale plates (Farallon, Tethys and Pacific), while smaller slabs would tend cross, although being slowed down, by the non-monotonic radial pattern of the lower mantle viscosity. This filtering effect could be enhanced by a mid mantle layering that induces a transient dynamics (Figure 9), which induces small slabs to stall in the mid-mantle after a period of several tens Myrs, diffusing heat and thereby dissipating away their differential buoyancy.

We modelled numerically the effects of a wide but realistic range of thicknesses, densities and viscosities for a non-homogeneous layer in the lower mantle and on the lower-mantle subduction of a relatively small and a very wide plate. Depending on such parameters, a variety of different behaviours is possible. The change in hill thickness induces alone stalling if a sufficiently strong negative buoyancy (to the mantle) is applied, while the transient behaviour is common to both 200 km and 500 km models (table 1), changing only the timing of the transition. We conclude that even a mild viscosity and density transition, if broad enough, can dramatically change the fate of the lower mantle slabs.

An unexpected outcome of the model is the fate of slabs that during their sinking encounter a combination of a high viscosity zone with a negative density anomaly. We predict that they might stall a long time in the mid-mantle transition zone, furthermore they might tear and/or fragment because of the opposite buoyancy forces encountered and then amplified by the local high viscosity of the mantle. This would reduce the wavelength of the seismic anomalies observed below in the deep mantle, which can explain the observation of Figure 4.

Plotting the hill-mantle viscosity ratio vs. slab-hill differential density in table 2 shows the most distinct transition between the three dynamic domains that we show in this work, with a banded configuration that illustrates how both the negative buoyancy and the high viscosity ratio are able to create a transient and a stalling dynamics. This somehow predictable result is here directly quantified, illustrating how a) if the viscosity peak proposed before for the mid mantle hill is realistic (100x) then the slab must either stall or have a long residential time in the mid mantle (see bottom line of the table 2); b) if the compositional heterogeneity in the deep lower mantle is more dense A c c e p t e d M a n u s c r i p t then the above one (as it is predictable) the dynamics of the lower portion of the lower mantle will be strongly inhibited and the slabs will always display a transient or stalling behaviour (see right column of the table 2); c) if the two effects act in a combined way, with close to neutral buoyancy and a mild hill-mantle viscosity ratio, the thickness of the transitional layer will play an important role controlling the slab behaviour and also the slab width will be relevant because inducing a larger mantle flow.

If a density and viscosity non-homogenous layer of thickness of 200 km or 500 km exists in the lower mantle, this must be necessary due to a compositional heterogeneity, either within locally in the layer, or between a sandwich consisting of a top and a bottom layer in the lower-mantle, because the spin transition in the ferropericlase alone cannot induce it [START_REF] Wentzcovitch | Anomalous compressibility of ferropericlase throughout the iron spin cross-over[END_REF]. Furthermore, the bulk modulus softening in ferropericlase during the spin transition is not observed in the dVp/dz data, which also implies that a compositional transition might compensate for the lack of dVp/dz jump. Such a configuration, if confirmed, requires a geodynamic system that at steady state is able to create such compensation. We envisage here this new mechanism, proposing that many sinking slabs would have transient time or completely stall in the mid lower mantle. Such slabs represent a compositional differentiated reservoir, as it has been for example already proposed for an hazburgite component in the D'' layer [START_REF] Hirose | The fate of subducted basaltic crust in the Earth's lower mantle[END_REF]. In the long term, a differentiation in the mid lower mantle would dynamically develop by a feedback due to the already existing compositional lowest mantle segregation, in a self-sustained form. This mechanism could induce over a long term the viscosity-density hill itself, causing a radial separation in two layers of the lower mantle, therefore preventing global mantle convection and favouring a more layered form of mantle dynamics.

Appendix. Fast Multipole Boundary Element Method for Stokes Flow

In order to capture the most fundamental modes that characterize how a slab can or not cross a density and/or a viscosity transition we simplified the system to its fundamentals. The effective viscosity (T,C) [START_REF] Kellogg | Compositional Stratification in the Deep Mantle[END_REF] (1)

We solve the generalized Stokes equations that comprise the momentum conservation and incompressibility condition:

(2)

For our problem, this amounts to subdividing the mantle into several closed regions, each characterized by a homogenous density and viscosity. When such an approximation is acceptable, it is possible rewrite (2) as a boundary integral equation for each sub-domain. This constitutes a decrease in the problem dimensionality and therefore potential computational gains. Moreover, this approach offers inherent multiscale capabilities as the surface mesh resolution can vary dynamically and track the physics of interest.

We show later how a system characterized by a perturbed viscosity  and/or density  in function of space:

      = '(z)       = '(z)
as is our case, can be solved using a perturbative approach to our boundary equations.

A c c e p t e d M a n u s c r i p t

Boundary Equations

We consider Stokes flow as described by the equations ( 1) and ( 2). Assuming a constant  in each domain, the velocity of each point inside the domain can be written as a sum of two surface integrals [START_REF] Pozrikidis | Boundary Integral and Singularity Methods for Linearized Viscous Flow[END_REF], called single and double layer integrals respectively, representing the effects of the forcing  ik (x)n k and velocity u i (x)

(5)

where G ij T ijk are the steady Green's functions for velocity and stress respectively, also known as the Stokeslet and the Stresslet:

(6) (7)
In turn, the above equation is implemented both for the inner and the outer fluid, in this way the boundary equations are cast into a form more appropriate for a quasi-steady multiphase flows in the presence of a gravity field. Hence for a point x on the surface S that separates different fluids, we obtain: (8)

where PV denotes the principal value of the integral,  0 is the viscosity of the external fluid, taken as a reference and  int  0 is the viscosity ration between inner and outer fluid and f is a normal stress jump that in our case accounts for gravity.

The equations have been widely used and tested for homogenous media. In this work we considered two types of inhomogeneities parameterized into a spatial variation of the coefficients 1) viscosity ratio  int  ext ; 2) density differential  int - ext We assumed that a zone of thickness between 200 km and 500 km, within a domain of several 1000's km, could have viscosities of 1-2 orders of magnitude higher then the surrounding and  0 =0 or even negative. In practice for each panel point x 0 at the center of the Panel(x 0 ),  has been rescaled as:

(9)

The implementation of non-homogenesous density is slightly more complex.

Exploiting a first-order hydrostatic approximation the non-uniform density is implemented as a space dependent shift of the pressure term in the single layer integral (in our case, the space variable is the radius and the pressure is defined by a profile of the pressure of the mantle to the lithosphere). We recalculated the entire radial pressure profile at each time step and the external forces applied in the single layer integral to the system. Hence the pressure P at the radius r=(x ii 2 ) 1/2 becomes:

(10) that is applied in the single layer of the boundary equation, replacing f with Pn.

Multipole Approach

The sinking slab surface S in Figure 5 and the supported quantities velocity u, and stress tensor at the boundary  are discretized with boundary elements (also called panels). The boundary integral equation thus becomes a linear system (11)

Many approaches carry out the construction of the matrix; which scales as N pan 2 both memory-and computation time-wise though, making it impractical for large systems.

We use a Fast Multipole Method (FMM) [START_REF] Barnes | A hierarchical O(N log N) force-calculation algorithm[END_REF][START_REF] Greengard | A Fast Algorithm for Particle Simulations[END_REF] in eq (8). The FMM algorithm, illustrated in figure 10, dramatically reduces the complexity of matrix-vector multiplication involving a certain type of dense matrix, which can arise out of many physical systems.

The FMM scales as N log(N), which is far more tractable and still allows the use of a Generalized Minimized Residual method GMRES or any Krylov space based method that does not rely on the storage of the full matrix. By treating the interactions between asymptotic basis functions using the FMM, the corresponding matrix elements do not need to be explicitly stored, resulting in a significant reduction in required memory. moments of the source distributions (quadrupoles). The reader is referred to [START_REF] Barnes | A hierarchical O(N log N) force-calculation algorithm[END_REF][START_REF] Greengard | A Fast Algorithm for Particle Simulations[END_REF][START_REF] Warren | A parallel hashed Oct-Tree N-body algorithm[END_REF] for general information on multipole methods and the work of [START_REF] Tornberg | A fast multipole method for the threedimensional Stokes equations[END_REF] for the transposition of harmonic multipoles to the evaluation of the Stokeslet and Stresslet.

The FMM algorithm thus sorts the sources in a tree structure whose cells contain the moment integrals--or multipoles--and carries out a field evaluation through a tree traversal (figure 11). The refinement of the interactions is determined by a tree traversal stopping criterion based on a prescribed tolerance.

The FMM-BEM drastically improves the computational cost of the method. For the coarse resolutions, the method displays the nominal N 2 scaling of a direct interaction [START_REF] Sigloch | Two-stage subduction history under North America inferred from multiple-frequency tomography[END_REF]. The [SG] line traces the Farallon-Kula paleo-plate boundary. At 1500km depth the two plates superpose, due to less space at depth [START_REF] Morra | Large Scale Three-Dimensional Boundary Element Simulation of Subduction[END_REF]. The broadening and flattening of the slab in the mid lower mantle confirms the the lower resolution outlook of the global tomography models displayed in figures 2 and 3.

Fig. 2. 3D isosurface of the 0.8% anomaly to the radial profile for the global tomographic model of [START_REF] Grand | Global seismic tomography: A snapshot of convection in the Earth[END_REF]. The three dimensional morphology of the Farallon slab is well distinguishable in the top portion of the figure. It is complex, flattening laterally and apparently reaching the core only in a very confined location.

There are also other prominent features of large and broad flat anomalies at mid mantle depth. [START_REF] Simmons | Constraining mantle flow with seismic and geodynamic data: A joint approach[END_REF], pri-s05 [START_REF] Montelli | A catalog of deep mantle plumes: New results from finite-frequency tomography[END_REF], rmsl-s06 [START_REF] Li | A new global model for $P$ wavespeed variations in Earth's mantle[END_REF], and saw642an [START_REF] Panning | A three dimensional radially anisotropic model of shear velocity in the whole mantle[END_REF]. The smean [START_REF] Becker | A comparison of tomographic and geodynamic mantle models[END_REF] model is an average the main proposed models and summarizes them. The cusp of the Farallon slab is a common feature of most tomographic models. Also the broadening of the image (indicateing slab flattaening) at depth between 1,500km and 2,000km is independent on the inversion approach employed. [START_REF] Boschi | On the statistical significance of correlations between synthetic mantle plumes and tomographic models[END_REF]. It is reproposed here as the main evidence of a dramatic transition at mid-lower mantle depths. It shows the logarithm of the ratio of positive-to-negative shear-velocity spectra from model SMEAN of [START_REF] Becker | A comparison of tomographic and geodynamic mantle models[END_REF], as a function of harmonic degree l (horizontal axis) and depth from Moho to core-mantle boundary (vertical axis). (Boschi Soo also the discussion proposed by [START_REF] Houser | The relative wavelengths of fast and slow velocity anomalies in the lower mantle: Contrary to the expectations of dynamics?[END_REF]. Positive (blue) values correspond to dominance of fast anomalies at a given depth and harmonic degree; vice-versa for negative (red) values: at relatively low harmonic degrees, fast anomalies dominate the pattern of seismic heterogeneity in the mid mantle, but the situation is reversed in the lower mantle, with the transition at ~1,500 km depth. Reyleigh Taylor instability from the center of the stalling plate, with a ring shape. The time evolution shown in the lower panels illustrates how negative buoyancies (i.e. hill density higher then lithosphere density) inhibits the crossing inducing slab stalling. In the long term, if the differential buoyancy remains respect to the lowest mantle, how small scale or even a subduction-like unstability triggered, as illustrated in Figure 8.

  , starting from the less prominent high-spin to low-spin contribution on the basis of mantle density and rheology.

  rising temperature, i.e. to the Arrhenius law again. In this sense a selection of stronger rocks at higher depth is the natural outcome of a very long-term dynamics in which the deep earth is naturally less mobile.

  A c c e p t e d M a n u s c r i p t 10 implemented either with a MP-BEM solver or coupling the code with FD finitedifference or FEM finite-elements methodologies, as it is in the implementation stage.

  A c c e p t e d M a n u s c r i p t 11 paradigmatic of the dramatic change due to the pure mantle shear stresses to the slab sinking. Three steps whose time evolution for the two models are illustrated in Figure 6.

  that characterizes rheological behaviour of subducted plates has been replaced by two mean values:  int inside the plate and  ext for the external mantle. Using such rheological definition the gradient of velocity and  is indicated by the viscous stress tensor:

  exploits the decay of the kernel to convolve and makes a controlled approximation. It does this by expanding the system Green's function using a multipole expansion, which allows one to group sources that lie close together and treat them as if they are a single source. More explicitly, let us compute (12) We consider the contribution from D i , a part of D that is far enough from our evaluation point x 0 and proceed with a Taylor expansion of the kernel G about x c D i We note that the equation involves successive moments of the  distribution in D i . The present FMM code can handle convolutions with the Green's functions for the Poisson equation, the Stokeslet or the Stresslet. It employs up to the second order

FiguresFig. 1 .

 1 Figures

Fig. 3 .

 3 Fig.3. A comparison of 5 global S-vawe tomographic models around the Farallon slab at 4 different depths. Models include tx2007[START_REF] Simmons | Constraining mantle flow with seismic and geodynamic data: A joint approach[END_REF], pri-s05(Montelli et al.,

Fig. 4 .

 4 Fig.4. This figure is mostly reproduced from figure 1 of[START_REF] Boschi | On the statistical significance of correlations between synthetic mantle plumes and tomographic models[END_REF]. It is

  , 2008) computed independent harmonic expansions of negative and positive velocity anomalies from SMEAN, and took the logarithm of the ratio of the resulting spectra.

Fig. 5 .

 5 Fig.5. Viscosity and density radial profiles for the lower mantle. Above are shown two models proposed in the past by (Ricard and Wuming, 1991) and (Forte and Mitrovica, 2001), both displaying a non-monotonous viscosity profile in the lower mantle, with a maximum at mid mantle depth. The lower panel illustrates instead the viscosity and density profiles that we test in this work. Out lower mantle model exhibit a less prominent viscosity peak respect to the one proposed above, as we observe in our model that 100x has already a strong influence on the geodynamic evolution. The density profile shown in bottom right figure illustrates both the hill-mantle and litho-mantle density differential. Those values are end-member, corresponding to typical seismic fluctutions. We model all the nine combinations.

Fig. 6 .

 6 Fig.6. Model results for the subduction of a small plate (2,000km x 2,000km) that encounters a mid lower mantle discontinuity. The differential density between slab and mantle is kept constant, 2%, while the viscosity peaks are tested, 10x and 100x, for a peak width of 500km. The viscosity transition is able to stop the slabs for a certain time, and retard largely the penetration into the lowest mantle, but the penetration in the deep mantle is very different for the viscosity peak of 10x or 100x, being the former a slab avalanche where the slab shape of the plate is maintained, while the last one is a Rayleigh Taylor instability type. The 100x is a much slower process, producing much longer residential times in the mid lower mantle.

Fig. 7 .

 7 Fig.7. Evolution of a sinking slab through a viscosity peak with two different lithospheric-mantle differential densities  for a not very large plate (2,000km) and very broad peak zone (500km). The viscosity hill is chosen 10x the background lower

Fig. 8 .

 8 Fig.8. The final crossing morphology of 6 different models, 3 characterized by a small 2,000km wide plate and 3 by a wide one (10,000km), all crossing a 200km thick transition zone and viscosity hill of 10x the background mantle. The models from left to right change through differential lithosphere-mantle density. The positive to null lithosphere-hill differential density comparison shows little effect, while the negative differential density (Hill Density > Slab Viscosity) instead triggers in the long term long residency times in the mid mantle and an unstable behaviour, coherently found for small or very wide plates, at different scales.

Fig. 9 .

 9 Fig.9. Summary of all models tested, through the radial trajectory of center of mass of the sinking slab. Small plates (2,000km wide) are displayed in the panel above, wihle very wide plate (10,000km) in the lower panel. The three families of models also recalled in the tables are (a) Stokes-like, which sink down to the core almost not interacting with the viscosity hill; (b) Stalling, with plates entrapped close to the viscosity-density hill for a period of the order of lower mantle Stokes crossing time or longer; (c) Transient, spend a moderate time around the viscosity hill, typically no longer of the lower mantle Stokes crossing time. Meaning of the legend: the line style indicates the thickness of the hill: dotted 500km (T=5) and continuous 200km (T=2).Colors instead display the hill peak viscosity ratio: blue 100 (V=100), green 10 (V=10) and red 1 (V=1). The dimensionless lithosphere-hill differential density is shown by symbols, which translated to standard differntial density values becomes: circle +2% (D=1), triangle 0% (D=0) and rhombic -2% (D=-1). Externally to the transition the litho-mantle density difference (D=1) and the viscosity ratio (V=1) are kept constant.

Fig. 10 .

 10 Fig.10. A general sketch of the multipole method employed. The summation on each surface is done using local (Point) or Multipole terms, depending on the relative position of different centers. Right Panel: The component of the summation are the integrals on each panel (=boundary element). Summing the integral contribution within a certain domain, the contribution of a set of panels is compressed into the one of one pole. Left Panel: There are three possible summation methods. Local(=Point)-Multipole, Local-Local, Multipole-Multipole. They depend on the relative contribution of the green function.

Fig. 11 .

 11 Fig.11. Sketch of a oct-tree and the hierarchical organization of its poles. Each internal node has up to eight children in order to partition the space by recursively subdividing it into eight octants until the required resolution is reached.

  Figure(s)

Table 1 :

 1 R = Hill  Mantle ) Summary of the results of the 36 experiments.Color version of the table 2: Litho  Hill  R=  Hill  Mantle

	Tables
	Hill thickness (km)	Viscosity Ratio ( Differential density ( Litho  Hill )	Slab width (km)	Fate of the slab (Stokes, Transient or Stalling)
		1	-2%	2,000	Transient
		1	-2%	10,000	Transient
		1	0%	2,000	Stokes
		1	0%	10,000	Stokes
		1	2%	2,000	Stokes
		1	2%	10,000	Stokes
		10	-2%	2,000	Transient
		10	-2%	10,000	Transient
		10	0%	2,000	Stokes
		10	0%	10,000	Stokes
		10	2%	2,000	Stokes
		10	2%	10,000	Stokes
		100	-2%	2,000	Stalling
		100	-2%	10,000	Stalling
		100	0%	2,000	Transient
		100	0%	10,000	Stalling
		100	2%	2,000	Transient
		100	2%	10,000	Transient
		1	-2%	2,000	Stalling
		1	-2%	10,000	Stalling
		1	0%	2,000	Stokes
		1	0%	10,000	Stokes
		1	2%	2,000	Stokes
		1	2%	10,000	Stokes
		10	-2%	2,000	Stalling
		10	-2%	10,000	Stalling
		10	0%	2,000	Transient
		10	0%	10,000	Transient
		10	2%	2,000	Stokes
		10	2%	10,000	Stokes
		100	-2%	2,000	Stalling
		100	-2%	10,000	Stalling
		100	0%	2,000	Stalling
		100	0%	10,000	Stalling
		100	2%	2,000	Transient

Table 2 :

 2 Synthetic summary of the resultsBW version of the table 2

:

   Litho  Hill  R=  Hill  Mantle

				
				
				Transient (T=200km)
	1	Stokes	Stokes	Stalling (T=500km)
			Stokes (T=200km)	Transient (T=200km)
	10	Stokes	Transient (T=500km)	Stalling (T=500km)
			Transient (T=200km)
	100	Transient	Stalling (T=500km)	Stalling

Table 2 :

 2 Synthetic summary of the results

Acknowledgements

Comments by two anonymous reviewers greatly helped to improve the manuscript.

We thank discussions with Fabio Cammarano, Anne Hofmeister, Matt Knepley, Marc Monnereau, Fabio Capitanio. This work, as part of the Eurohorcs/ESF European Young Investigators Awards Scheme, was supported by funds from the National Research Council of Italy and other National Funding Agencies participating in the 3rd Memorandum of Understanding, as well as from the EC Sixth Framework Programme.

The FMM-BEM has been parallelized using MPI. The parallel efficiency has been tested on a Opteron cluster with Quadrics connections. The scaling is very good up to 64 CPUS, still keeping 90% of efficiency. In its current implementation the FMM-BEM uses a shared (not distributed) tree, thus reducing the communication load at the expense of memory requirements.