Russel Nzekwa
email: russel.nzekwa@inria.fr

Romain Rouvoy
email: romain.rouvoy@inria.fr

Lionel Seinturier
email: lionel.seinturier@inria.fr

Towards a Stable Decision-Making Middleware for Very-Large-Scale Self-Adaptive Systems

Keywords: context-awareness, self-adaptation, stabilization, very-large-scale systems

With the development of the communication infrastructures, the number of applications collaborating at large scale increases. To maintain, and continue to deliver services of good quality to the end-users, very-large-scale applications continuously adapt themselves, depending on the changes in their surrounding. Stabilization of the system thus becomes a keystone issue in the adaptation process, in order to reduce the system reconfiguration cost. Existing approaches, for the stabilization of very-large-scale systems, provide solutions that are partially efficient. For example, learning-based stabilization algorithms give good results in predicting application behaviors, but still suffer from their weak reactivity. In this paper, we propose an approach of combining different goals-oriented stabilization algorithms, in order to provide sustainable and efficient stabilization for large-scale systems.

Introduction

Because of the huge amount of data that they are intended to process, very-large-scale systems need to retrieve only relevant information from the data flow, in order to perform the system adaptation. To optimize the system adaption process, stabilization mechanisms are required. The concept of context region defined in Quality of Object (QuO) [START_REF] Atighetchi | Building auto-adaptive distributed applications: The quo-apod experience[END_REF], is an example of a stabilization mechanism that advocates a strict partition of context space, to reduce the adaptation side effects. However, context regions do not handle application behavior during the transition of the system from one state to another. On the contrary, some stabilization algorithms based on learning strategies, such as Bayesian Network (BN) are much more efficient in handling that issue. But, despite some good results in predicting application behavior, even during the state transition, they still suffer from weak reactivity. The stabilization process can be more efficient, if the stabilization mechanisms can be adapted to the processed data. We call stabilization mechanisms, the set of algorithms or techniques aiming at regulating the responsiveness of adaptive system. This paper describes an adaptive stabilization approach for composing stabilization algorithms in order to increase their efficiency and accuracy.

The rest of the paper is organized as follows. Section 2 presents the motivations and challenges of the stabilization for very-large-scale systems. Section 3 discusses some related works. Section 4 describes our stabilization approach. We show some simulation results in Section 5. Finally, we conclude in Section 6.

Motivations & Challenges

Motivations With the growing complexity of software nowadays, the cost (time, money) of maintaining applications in a changing environment increases exponentially. Therefore, there is a high demand for self-managing system achieving desired quality requirements with a reasonable cost. Self-adaptive software answers to this demand. Such systems are able to monitor, to detect significant changes, to decide how to deal with, and to act consequently. When changes in the application environment frequently occur, it is important to limit system reconfigurations in order to keep a good tradeoff between the adaptation cost and the system performance. This is specially true, when dealing with large-scale self-adaptive applications where the communication cost can be expansive. On this specific issue, stabilization mechanisms play a key role allowing to reduce adaptation cost. As an example of a motivation scenario, and to easily point out challenges of the stabilization issue for large-scale self-adaptive application, we consider the "Trucks Tracking" scenario. This scenario is the subject of SALTY ANR project 1 .

Trucks Tracking scenario: We consider a self-adaptive application for the management of the truck fleet (800.000 to 1.000.000 trucks), of a company specialized in the transport of fragile products. All the trucks do not have the same characteristics. Some are equipped with good air conditioning system, while other provide robust slip system. In the same way, transported goods do not have the same requirements. Some are very sensitive to the temperature variation, when other need high security system. The overall objective of the system is to make sure the trucks reach their destination on time. The application must be able to detect stop times, temperature variation inside containers, security violation access and truck position. The application is also connected to remote services like the weather service or the city traffic service. Self-adaptive application, must also notify destination logistic platforms about truck arrival. If one truck gets into trouble, the system can remap the route of other trucks to help it. All aforementioned information is sent to the central location platform which processes and decides which adaptation process can be triggered. Figure 1, provides an illustration of the scenario.

From the analysis of data exchanged between the nodes of the application, several adaptations can result, such as, the frequency of the positioning requests, or the level of allocated resources. Given the large amount of data involved, and the distributed nature of the application, a naive technique consisting of sending all data for analysis cannot be effective. Therefore, stabilization mechanisms must be implemented.

Challenges In our opinion, two important challenges are to be met, when dealing with stabilization of self-adaptive applications in a large-scale distributed environment. The first challenge is to find the right level of stabilization, so that the application can be reactive enough to reflect important variations in the surroundings. Secondly, when selfadaptive applications are running, new features (data, services) can be added to fill specific needs. In this particular case, to be efficient, the stabilization mechanism need to be adapted in order to take into account the new added features. We believe that an adaptive strategy for the composition of stabilization algorithms is the beginning of the solution to move towards more accurate and efficient stabilization of self-adaptive system. In the next section, we present our proposition on this issue.

Related Work

Stabilizations techniques are used in the project TEA [START_REF] Albrecht | Multi-sensor context-awareness in mobile devices and smart artifacts[END_REF]. The principal limitation of this project is that implemented stabilization techniques are tightly coupled to application architecture and cannot easily be replaced or modified. [START_REF] Folliot | Phoenix: a self-adaptable monitoring platform for cluster management[END_REF], presents the PHOENIX project which targets the management and monitoring of distributed system (clusters) processing many heterogeneous data. To improve the performance of the system monitoring, and detect emergency situation just when they happen, PHOENIX introduces delta operator (DO) with other first order logic operators. DO allow evaluation of signal amplitude variation, so the data flow is relayed only when it reach the required threshold value. Despite good results recorded by the method, one of the main drawbacks is that the PHOENIX platform does not support extension of new operators, and that the stabilization method is targeted for quantitative numerical data. In [START_REF] Padovitz | Maintaining continous dependability in sensor-based context-aware pervasive computing systems[END_REF], Padovitz suggested an architecture to handle stabilization issue in self-adaptive systems on the base of Kalman's filter learning algorithm. The main limitation of that approach was the weak reactivity of the system compared to other approaches.

From the "Truck Tracking" scenario presented in the previous Section 2, we can notice the following points concerning the system reconfiguration.The amount of information to process for the geographical localization of trucks, is huge and variable: (i)-the number of trucks to track , (ii)-variable frequency of the requests , (iii)-variable precision of the positioning information. In order to provide a good service, the system must implement dynamic adaptation mechanisms. Implementing learning-based stabilization algorithm like Dempster-Shaffer (DST) [START_REF] Wu | Sensor data fusion for context-aware computing using dempster-shafer theory[END_REF] or Bayesian Network (BN), will improve the knowledge of the system about the environment, but at the same time will introduce latency in detecting new changes. The method proposed by Padovitz et al. [START_REF] Padovitz | Maintaining continous dependability in sensor-based context-aware pervasive computing systems[END_REF], in the case of the scenario is not suitable, because of the complexity of the mathematical expression determining the instability of the system, which increases with the number of arguments.

To the best of our knowledge, a unique stabilization algorithm that could meet all requirements of self-adaptive applications does not exist. In fact, in the context of very-large-scale systems, the stabilization strategies require to be flexible enough in order to adapt to the evolution of the environment. We believe that, in the perspective of improving the management of stabilization process for self-adaptive systems, a solution can be found in the flexible combination of several existing approaches.

Flexible Context-Aware Architecture

This section presents our stabilization approach for self-adaptive systems. We consider that in order to boost flexibility and efficiency, stabilization issue should be handled separately from decision-making basics concerns. For the general case of an autonomic control loop [START_REF] Kephart | The vision of autonomic computing[END_REF], we suggest to integrate stabilization blocks at two levels of the chain processing information data: Firstly, after acquisition of raw data from different probes or sensors, and secondly, after the decision-making block. Figure 2 depicts this architecture. The purpose of the first stabilization block, located between data collectors and decision-making blocks, is to filter data coming from collectors and to forward to decision-making only significants values that could lead to system reconfigurations. This considerably reduces access to decision-making module, resulting in a gain of performance for the system, since running decision-making processes is quite often a costly procedure for the application. The purpose of the second stabilization block, located in our architecture between decision-making block and actuator, is to optimize the work of actuators by finding the suitable moment when reconfigurations of the system can be executed. This issue has a significant impact when dealing with distributed applications, where the decision-making structure and the application are not hosted on the same machine.

In order to meet flexibility requirements for an efficient stabilization of self-adaptive applications, we propose to define a model dedicated to the composition of stabilization algorithms.

Composition Model

Our composition model defines how to combine different stabilization mechanisms in order to meet flexibility in the application. The model consists of two modalities of composition: Horizontal and Vertical composition.

Horizontal Composition

Learning-based stabilization algorithms, like Dempster-Shaffer (DST) [START_REF] Wu | Sensor data fusion for context-aware computing using dempster-shafer theory[END_REF] or Bayesian Network (BN), while providing a good efficiency in predicting contextual changes, introduce a latency in detecting variations in the application environment. This weakness can be filled in, by associating to learning-based algorithms other algorithms with a higher reactivity speed.

Horizontal composition consists in executing concurrently several stabilization algorithms. The idea behind this concept is to benefit from passive and reactive stabilization techniques by combining them adequately. The detection of irregular application's behavior can be done by combining a reactive algorithm like Delta Operator (DO) , and less reactive algorithm like Dempster-Shaffer (DST) algorithm, through a composition rule (CR). A CR can be some simple rule like " max(v n , v n+1 ) ⇒ proceed value", where v n , v n+1 are context values, or a more complex rule involving Quality of Context (QoC) of the processed data. Hence, using this composition model can help to improve accuracy of the stabilization process while keeping a reasonable level for the reactivity of the system. We will further detail to this example in section 5.

Vertical Composition

Vertical composition in our approach consists in applying sequentially two or more stabilization algorithms on the same data set. Some works [START_REF] Sekkas | Context fusion through imprecise reasoning[END_REF], suggest that it can be interesting in terms of performance for stabilization of context information to apply successively several algorithms of stabilization on the same sample of data. A good illustration of that idea is the work of Sekkas et al. [START_REF] Sekkas | Context fusion through imprecise reasoning[END_REF] where the authors compare the efficiency of using Bayesian Network (BN), Dynamic Bayesian Network (DBN), and Fuzzy Logic (FL) alone or in a combined way. In our approach we believe that, the combination of algorithms using vertical composition can increase efficiency of the stabilization. In order to limit the overhead introduced by the use of several algorithms, "cheap" (low execution cost) algorithms are found at the beginning (bottom) of the stabilization chain while "expensive" (high execution cost) algorithms are on the top of the architecture. This association rule is mostly justified by the fact that, the amount of processed data decreases from the bottom to the top, thus more costly algorithms at the top of the architecture would have to process less data, decreasing by the same way the overall cost of the stabilization process, which is tightly bound to the amount of context information processed.

Figure 3 gives an illustration of the composition model. From left to right, we have horizontal composition, then vertical and finally the combination of the both types of composition. Combination of both primitive composi- tion models in our approach can be beneficial for the stabilization process in order to meet accuracy and efficiency properties.

Evaluation and Simulations

In the previous sections, we presented a novel approach in order to meet flexibility when dealing with stabilizations issues in pervasive applications. We suggested to use composition of algorithms based on the presented composition model. Algorithms composition offers the possibility to maximize the efficiency of the stabilization process at a relative low cost for the application. While in most of the existing systems, the decision-making module applies a generic algorithm to processed data, our approach offer the possibility to adapt the stabilization mechanism to the data being processed. To evaluate our proposition, we have made some experimental simulations that illustrate the feasibility of our approach. Our simulations are based on COSMOS [START_REF] Rouvoy | Software architecture patterns for a context-processing middleware framework[END_REF](Context Entities Composition and Sharing). We also use the simulation engine SIAFU [START_REF] Miquel | A generic large scale simulator for ubiquitous computing[END_REF] to generate contextual informations.

The scenario for this simulation is similar to the one presented in Section 2, because the application performs some adaptations depending on the changes in the environment. In the current scenario, we have temperature sensors simulated by the simulator engine that send data concerning the temperature of the environment. Time elapsed between two consecutive readings from sensors is 0.3 second. Collected values have an accuracy close to 0.1. The overall objective of the application process will be to detect all the relevant variations in the application environment, around the range of 24 • C to 25.5 • C. We consider that a relevant variation in the environment, is a variation that persist along three consecutive readings from the sensors.

The first implemented algorithm, is "Delta operator" (DO) [START_REF] Folliot | Phoenix: a self-adaptable monitoring platform for cluster management[END_REF]. The choice of the threshold value is motivated by the data sample's variance, in order to choose the most appropriate value. The second implemented algorithm is Kalman filter. The central point for Kalman's filter algorithm is the determination of the matrix of transition A, in this case we defined the transition matrix as following: Given A the transition matrix, given x k-1 , z k-1 , R the prediction, the value of measured variable at k -1 step and the variance respectively. We have defined ε, ∆ then A as follows: 4 shows the curves of the application behavior without stabilization, stabilized with KF, and stabilized with DO. As expected, application stabilized with DO is very reactive, we can notice that around t = 2s for example, the curves of the application without stabilization and stabilized with DO are almost indistinguishable. Figure 5 shows application's behavior curves when combining the both algorithms using our model, in this case an horizontal composition. On Figure 5, the curve of the behavior of the application stabilized with both KF and DO algorithms is represented. On that graph, we can notice two things. Firstly, the detection of the first event for the application occurs earlier (zone B) than for an application stabilized by KF only. That is due to the fact that the application benefit from the reactivity property of DO algorithm. Secondly, we can notice that when an event is detected the behavior of the system does not change suddenly (zone A or zone B) even when the next raw value processed by the system does not belong to the target area. That last point is a property that the system inherited from KF algorithm. Obviously, the composition model offers a good compromise between a reactive and less reactive system. Finally, to characterize our system we have measured the filtration rate, the results of the measurement are depicted in 6. We can notice that resulting application (DO+KF) has the same filtrating rate as the application application stabilized only with KF.

∆ = x k-1 -z k-1 , ε = R * z k-1 A k =            A k-1 f or |∆| ≤ ε z k -ε x k-1 f or |∆| > ε and ∆ < 0 z k +ε x k-1 f or |∆| > ε and ∆ > 0 Figure

Conclusion

In this paper, we have proposed an adaptive approach for context stabilization in a distributed environment. For that purpose, we suggest two modalities of composing stabilization algorithms in order to increase accuracy and efficiency of the stabilization. We showed a simulation example to illustrate the feasibility of our approach. In the future, we plan to probe this approach on a real very-large-scale platform like an experimental or a production grid.

Figure 1 :

 1 Figure 1: Self-Adaptive Distributed Application

Figure 2 :

 2 Figure 2: Flexible context-aware architecture

Figure 3 :

 3 Figure 3: Composition Models

Figure 4 :

 4 Figure 4: Stabilized application behavior -Kalman filter and Delta operator

Figure 5 :

 5 Figure 5: Comparatives curves of stabilized application behavior

Figure 6 :

 6 Figure 6: Filtration rate measurements

SALTY ANR project: a research project funded by the french research agency. http://salty.unice.fr


Kota Mukumoto

Yann Capdeville

Sneha Singh

Takeshi Tsuji

Accounting for sub-wavelength heterogeneities in full waveform inversion based on wavefield gradient measurements

teaching and research institutions in France or abroad, or from public or private research centers.

INTRODUCTION

Recent advances in seismic instrumentation have enabled the use of wavefield gradient measurements to more comprehensively characterise ground motions and further constrain various seismic phenomena, including earthquake mechanisms, the directionality of seismic waves and anisotropic structures. However, it is important to carefully examine the feasibility of directly applying seismic methods that have been used for translational measurements to wavefield gradient measurements.

The ability to acquire seismic wavefield gradient measurements has existed for decades. Tiltmeters, which measure the tilt of the ground surface (spatial derivatives of displacement), have traditionally been used to investigate tidal tilt [START_REF] Nishimura | On earth tides[END_REF][START_REF] Melchior | The earth tides[END_REF]. Seismic observations have also been made using tilt meters [START_REF] Sassa | On phenomena forerunning earthquakes[END_REF]; however, measurements of seismic rotational motion have been a long-standing challenge due to their much smaller amplitudes compared with the translational motions that are excited by earthquakes. Recent advances in sensor technologies, such as ring laser gyroscopes [START_REF] Mcleod | Comparison of standard and ring laser rotational seismograms[END_REF], have enabled rotational motion measurements, and subsequent studies have shown the benefits of using rotational sensors in various geological settings. For example, collocated measurements of rotation and translational motions have been analysed to retrieve the local phase velocity [START_REF] Mikumo | Determination of local phase velocity by intercomparison of seismograms from strain and pendulum instruments[END_REF][START_REF] Igel | Rotational motions induced by the m8. 1 tokachi-oki earthquake, september 25[END_REF][START_REF] Fichtner | Sensitivity densities for rotational ground-motion measurements[END_REF]. Recently, the theory to measure the dispersion relation using both translational and wavefield gradient measurements at single observation point in weakly anisotropic media has been introduced [START_REF] Tang | Single-point dispersion measurement of surface waves combining translation, rotation and strain in weakly anisotropic media: theory[END_REF]. Detailed strain measurements from distributed acoustic sensing (DAS) arrays have led to a significant increase in the number of geophysical studies that implement dynamic subsurface strains [START_REF] Willis | Quantitative quality of distributed acoustic sensing vertical seismic profile data[END_REF][START_REF] Zhan | Distributed acoustic sensing turns fiber-optic cables into sensitive seismic antennas[END_REF]). These recent advances in sensor technology and the increased use of rotational sensors and DAS arrays have enabled a transition towards actively observing and utilising measurements of seismic wavefield gradients.

Many seismic analyses based on rotational sensor and DAS observations have been conducted.

Phase-velocity and polarisation analyses using rotational sensors [START_REF] Wassermann | Toward a single-station approach for microzonation: Using vertical rotation rate to estimate love-wave dispersion curves and direction finding[END_REF][START_REF] Keil | Single-station seismic microzonation using 6c measurements[END_REF], and structural imaging [START_REF] Parker | Active-source seismic tomography at the brady geothermal field, nevada, with dense nodal and fiber-optic seismic arrays[END_REF], seismic-wave monitoring [START_REF] Dou | Distributed acoustic sensing for seismic monitoring of the near surface: A traffic-noise interferometry case study[END_REF] and earthquake studies (Ajo-Franklin et al. 2019) using DAS cables have made important advances in seismology. These studies have demonstrated the potential of using wavefield gradient measurements to address fundamental seismological questions; however, it is important to pay attention to the distinct nature of the wavefield gradients compared with translational motions when analysing wavefield gradient measurements.

Wavefield gradients can be significantly influenced by local small-scale heterogeneities in the subsurface that distort both the phase and amplitudes of the wavefield gradients [START_REF] King | Tidal tilt measurement in europe[END_REF][START_REF] Harrison | Cavity and topographic effects in tilt and strain measurement[END_REF][START_REF] Van Driel | Strain rotation coupling and its implications on the measurement of rotational ground motions[END_REF][START_REF] Singh | Correcting wavefield gradients for the effects of local small-scale heterogeneities[END_REF][START_REF] Muir | Wavefield-based evaluation of das instrument response and array design[END_REF]. This issue has been acknowledged as a critical problem in tilt measurements [START_REF] Melchior | The earth tides[END_REF][START_REF] Lennon | The earth tide signal and its coherency[END_REF][START_REF] Goulty | Strainmeters and tiltmeters in geophysics[END_REF]. A previous study reported sudden changes in tilt values due to local subsurface heterogeneities at the sensor, such as cavities, discontinuities and elastic contrasts, and also the local topography, and noted discrepancies between the theoretical and observed tilt values [START_REF] Meertens | Tilt observations using borehole tiltmeters: 2. analysis of data from yellowstone national park[END_REF]). Similar effects have been observed in rotational sensor and DAS measurements [START_REF] Singh | Correcting wavefield gradients for the effects of local small-scale heterogeneities[END_REF][START_REF] Muir | Wavefield-based evaluation of das instrument response and array design[END_REF]. [START_REF] Singh | Correcting wavefield gradients for the effects of local small-scale heterogeneities[END_REF] obtained a strong fit between the observed and synthetic displacement waveforms using a smooth Earth model whereas the observed rotational waveforms from the same locations did not fit the synthetic waveforms well; they attributed these findings to the effects of small-scale heterogeneities on measurements of the wavefield gradients. Two-scale homogenisation theory, which has recently been introduced in seismology, can mathematically explain the effects of small-scale heterogeneities on wavefield gradient measurements (Capdeville et al. 2010a;Capdeville et al. 2010b;[START_REF] Cupillard | Non-periodic homogenization of 3-d elastic media for the seismic wave equation[END_REF]. This theory shows that coupling of the effective strain to gradient measurements can result in the effects of small-scale heterogeneities. [START_REF] Singh | Correcting wavefield gradients for the effects of local small-scale heterogeneities[END_REF] proposed a method for inverting the coupling-term coefficients to address the discrepancies between the observed and simulated rotational data. These coefficients can then be used to correct the simulated data and explain the observed wavefield gradients [START_REF] Singh | Correcting wavefield gradients for the effects of local small-scale heterogeneities[END_REF][START_REF] Muir | Wavefield-based evaluation of das instrument response and array design[END_REF]. Therefore, both the coupling-term coefficients (hereafter "correctors") and a smooth tomographic model are required to accurately simulate wavefield gradient observations, as the ability to obtain high-resolution tomographic models with small-scale structures is generally challenging.

Methods have been developed to account for the effects of small-scale heterogeneities when forward modelling wavefield gradients. However, few studies have investigated the effects of small-scale heterogeneities on wavefield gradients in seismic inverse problems, such as full waveform inversion (FWI). Although it may be straightforward to assume that a FWI using wavefield gradients can recover the small-scale structures around the receivers, thereby generating higherresolution images than those obtained using translational motions, the impacts of small-scale heterogeneities on wavefield gradient-based FWI have not been thoroughly investigated. This study aims to investigate the effects of wavefield gradient measurements on FWI results and address the impact of small-scale heterogeneities. Our findings demonstrate that small-scale heterogeneities can strongly affect the FWI inversion results. Given the sensitivity of FWI frameworks to these small-scale effects, we propose a new FWI framework that can account for small-scale heterogeneities and improve the accuracy of FWI results.

FULL WAVEFORM INVERSION AND HOMOGENIZATION THEORY

Context

We consider a two-dimensional (2D) elastic domain Ω with absorbing boundary conditions on its boundary ∂Ω. We use a Cartesian coordinate system in Ω, where the x-axis is the horisontal axis and the z-axis is the vertical axis. The particle displacement vector u(x, t) in Ω is driven by the elastic wave equation:

ρ∂ tt u -∇ • σ = f s , (1)
σ = c : ϵ(u) , (2)
where

ϵ kl (u) = 1 2 (∂ k u l + ∂ l u k ) is the strain operator, (k, l) ∈ {x, z} 2 , ρ(x) is the density, c(x) is
the elastic tensor, σ(x, t) is the stress tensor, f is the source term and x = (x, z) is the position vector in Ω. We use point sources, such that:

f s (x, t) = F s δ(x -x s )g(t) (3)
for a point force, or:

f s (x, t) = -M s • ∇δ(x -x s )g(t) (4)
for a point moment tensor, where F s is the force vector, M s is the moment tensor, x s is the source location and g(t) is the source time wavelet.

The mechanical properties (ρ(x), c(x)) are considered heterogeneous at all spatial scales, as expected in geological media. The two-scale homogenisation mathematical framework is therefore appropriate to deal with such heterogeneities.

Two-scale homogenisation

Two-scale homogenisation is a mathematical framework that can deal with multiscale problems.

It was originally developed for periodic and stochastic elastostatic and elastodynamic problems [START_REF] Sanchez-Palencia | Non homogeneous media and vibration theory[END_REF][START_REF] Bensoussan | Asymptotic Analysis of Periodic Structures[END_REF]) before being extended to deterministic multiscale heterogeneous media with no scale separation for wave propagation (Capdeville et al. 2010a;Capdeville et al. 2010b;[START_REF] Cupillard | Non-periodic homogenization of 3-d elastic media for the seismic wave equation[END_REF].

Homogenisation theory considers scales of heterogeneities as either small-scale (microscopic scale) or large-scale (macroscopic scale) features. Small-scale heterogeneities are meant to be homogenised and replaced by effective properties. For periodic heterogeneities, the small scale is the periodic structure, and the large scale is the constant background. In seismology, there is no periodicity in the media and we rely on an arbitrary scale λ 0 to separate small for large scales.

Homogenisation is an asymptotic method based on a small parameter ε 0 , which is defined as follows:

ε 0 = λ 0 λ min . ( 5
)
ε 0 measures the position of the separation between the small and large scales with respect to the minimum wavelength λ min . Note that any value can be chosen for ε 0 . However, ε 0 = 0.5 (meaning that heterogeneities half the size of the minimum wavelength and smaller are considered smallscale heterogeneities) is a good approximation that has generally been employed for most cases.

The two-scale homogenisation method relies on two space variables: the regular space variable x (also called the large-or macroscopic-scale variable) and a new space variable y = x ε 0 (also called the small-or microscopic-scale variable). A new mathematical problem is derived, whereby all of the physical quantities depend on both the macro-scale spatial variable x and the micro-scale spatial variable y. The details of the derivation of this homogenisation problem can be found in [START_REF] Capdeville | Machine Learning in Geosciences[END_REF].

The displacement u, which is a solution of the wave equation (equations ( 1)-( 2)), is a powerseries approximation in ε 0 :

u ε 0 (x, y, t) = u * (x, y, t) + ε 0 u 1 (x, y, t) + ε 2 0 u 2 (x, y, t) + ... , (6)
where the homogenised coefficients u i need to be determined. Similarly, the strain can also be expanded as a power-series approximation in ε:

ϵ ε 0 (x, y, t) = ϵ 0 (x, y, t) + ε 0 ϵ 1 (x, y, t) + ε 2 0 ϵ 2 (x, y, t) + ... (7)
where ϵ i are the strain homogenised expansion coefficients. Note that the exponent i is the power for ε 0 and an index for u i and ϵ i in equations ( 6) and ( 7). We can also define a homogenised expansion for the rotational field ω = 1 2 ∇ × u:

ω ε 0 (x, y, t) = ω 0 (x, y, t) + ε 0 ω 1 (x, y, t) + ε 2 0 ω 2 (x, y, t) + ... . (8)
The true u, ϵ and ω values are then found via:

u(x, t) = u ε 0 (x, y = x ε 0 , t) , (9) ϵ
(x, t) = ϵ ε 0 (x, y = x ε 0 , t) and (10) ω(x, t) = ω ε 0 (x, y = x ε 0 , t) . (11)
The main results of the homogenisation development are the following:

• The displacement of the leading-order approximation (in ε 0 ) is independent of the small-scale variable y:

u * (x, y, t) = u * (x, t) . ( 12
)
• u * is the solution of the effective wave equation:

ρ * ∂ tt u * -∇ • σ * = f * s , (13)
σ * = c * : ϵ * , (14)
where ϵ * = ϵ(u * ) and (ρ * , c * ) are the effective mechanical properties (see below), and f * s is the effective source [START_REF] Burgos | Homogenized moment tensor and the effect of near-field heterogeneities on nonisotropic radiation in nuclear explosion[END_REF][START_REF] Capdeville | Homogenization of seismic point and extended sources[END_REF]).

• The first-order homogenised coefficient can be written as follows:

u 1 (x, y, t) = χ(x, y) : ϵ * (x, t) + u 1 (x, t) , ( 15
)
where χ is the first-order corrector (a third-order tensor) and ⟨u 1 ⟩ can generally be neglected.

• The strain of the leading-order approximation depends on the small-scale variable y:

ϵ 0 (x, y, t) = (I + ϵ y (χ)(x, y)) : ϵ * (x, t) , ( 16
)
where I is the fourth-order identity tensor and

[ϵ y ] ij (v) = 1 2 (∂ y i v j + ∂ y j v i ) for any v.
• Similarly, the leading-order rotational coefficient can be written as:

ω 0 (x, y, t) = ω * (x, t) + 1 2 (∇ y × χ(x, y)) : ϵ * (x, t) , (17)
where ω * = 1 2 ∇ × u * . Note that ω * and ω 0 are generally different. • The homogenised quantities of ρ * (x) and c * (x), as well as the corrector χ(x, y), can be computed through the homogenised operator H:

(ρ * , c * , χ) = H(ρ, c) . ( 18
)
H implies solving a set of equations known as the "cell problem," which is a set of static elasticity equations with a set of source terms (Capdeville et al. 2010b;[START_REF] Capdeville | Fast fourier homogenization for elastic wave propagation in complex media[END_REF]. The solution only depends on ε 0 and the elastic model; neither the time nor the seismic source incluence H.

The presence of normal stress-free boundary conditions requires the use of effective boundary conditions [START_REF] Capdeville | Shallow layer correction for spectral element like methods[END_REF][START_REF] Capdeville | Residual homogenization for seismic forward and inverse problems in layered media[END_REF] to homogenise the small-scale variations in both the local topography and structures just beneath the subsurface. However, we omit the effective boundary condition in this study for simplicity. H is a non-linear operator that often yields non-intuitive homogenisation results. For example, the effective properties are almost always anisotropic, even if the true model is isotropic. Solving the cell problem and then finding the effective mechanical properties generally involve a numerical solver, with the layered-medium case being the only case that leads to an analytical solution. Homogenisation theory yields the well-known results of [START_REF] Backus | Long-wave elastic anisotropy produced by horizontal layering[END_REF] for the layered-medium case. We do not expand further on this aspect because it is not required for our discussion; see [START_REF] Capdeville | Machine Learning in Geosciences[END_REF] for a complete discussion. ε 0 ≤ 0.5 is considered a reasonable choice to obtain homogenised quantities and preserve the solution accuracy for most geological media (Capdeville et al. 2010b).

Finally, we can also find the relationship between the effective fields, and the true rotational and gradient of the wavefield. The leading-order approximation of the wavefield gradient is:

∇u(x, t) = ∇u * (x, t) + G(x) : ϵ * (x, t) + O(ε 0 ) , where: G(x) = ∇ y χ(x, y)| y= x ε 0 , (19)
and ∇ y is the gradient operator that is applied to the y variable.

The leading-order approximation of the rotational field is:

ω(x, t) = ω * (x, t) + J(x) : ϵ * (x, t) + O(ε 0 ) , where: J(x) = 1 2 ∇ y × χ(x, y)| y= x ε 0 . (20)
In practice, once the effective medium (ρ * , c * ) is obtained, the effective wave equation (equations ( 13)-( 14)), which is the standard wave equation, can be then solved to obtain the effective displacement u * and its derivative ϵ * , ∇u * and ω * . Nevertheless, if using u * we obtain directly the true displacement with a good approximation via equation (12), this is not the case for the strain and rotational wavefields respectively. Equations ( 19) and (20) show that an extra term is needed to obtain the correct leading-order approximations for the effective strain and rotational fields. This extra term depends on y, thereby suggesting that the waveforms of the wavefield gradients are affected by the small-scale structures of the medium. It also implies that the gradient and rotational field are sensitive to small-scale structures, whereas the displacement is not. In practice, this small-scale heterogeneity causes the strain components to couple to the expected gradient of the rotational components. This effect has already been observed in both synthetic and real data [START_REF] Van Driel | Strain rotation coupling and its implications on the measurement of rotational ground motions[END_REF][START_REF] Singh | Correcting wavefield gradients for the effects of local small-scale heterogeneities[END_REF][START_REF] Muir | Wavefield-based evaluation of das instrument response and array design[END_REF].

Earth models and associated spaces

We first define the space of possible Earth models prior to introducing the inverse problem. An

Earth model m is defined by its mechanical properties (ρ(x), c(x)) at any location x in the domain Ω. We define the admissible model space M, which contains all of the physically admissible Earth models. If we assume that Ω is heterogeneous at all scales (similar to the real Earth), then M is infinite-dimensional.

We can define the homogenised counterpart of each model m ∈ M at a given frequency band as p * = H(m). Equation (18) shows that p * is more than just the effective mechanical properties m * = (ρ * , c * ), as it generally contains the corrector χ. m * is the restriction of p * to the mechanical properties: m * = p * | m . Note that the actual definition of the effective model generally depends on the receiver geometry, type of observation and homogenisation order. If there are only displacement observations and the homogenisation order is 0, then there is no need for correctors. However, if the observations are rational, even for zeroth-order homogenisation, then the corrector curl is needed at the receiver locations (see equation ( 20)). We therefore define the effective model space M * , which contains all of the possible effective models p * . Symbolically, we have M * = H(M), where M * is finite-dimensional.

Full waveform inversion based on displacement data

We consider an idealised FWI problem whose objective is to recover the true mechanical properties (ρ t , c t ) of the elastic domain Ω using seismic waveforms. We assume that the true model m t = (ρ t , c t ) belongs to M. We consider Ω with N r receivers located at x r , r ∈ {1, ..., N r } and N s sources located at x s , s ∈ {1, ..., N s }. The signal from each source is recorded for a time duration T at each receiver. We assume that the displacement data set is accurately modelled by solving the wave equation (equations ( 1)-( 2)). We define the least-squares misfit function and associated minimisation problem as:

E(m) = r,s T 0 (d s (x r , t) -u(x r , t; x s , m)) 2 dt , m = argmin m∈M E(m) , (21)
where d s (x r , t) represents the displacement data generated by source number s and recorded by receiver number r. Under perfect conditions, the inverse problem has a unique solution and m = m t [START_REF] Nachman | Reconstructions from boundary measurements[END_REF][START_REF] Nakamura | Global uniqueness for an inverse boundary problem arising in elasticity[END_REF]. However, equation ( 21) is generally an unsolvable problem because M is an infinite-dimensional space and the numerical cost of solving the wave equation is not bounded.

The classical strategy to obtain a solvable inverse problem is to limit the data frequency band.

Therefore, a low-pass filter F fm is introduced, such that the filtered data d fm s = F fm (d s ) have no signal beyond the maximum frequency f m in the frequency domain. A new misfit is then defined from this filtered data set as follows:

E fm (m) = r,s T 0 (d fm s (x r , t) -u fm (x r , t; x s , m)) 2 dt , (22)
where u fm are the synthetic data, which are computed with a low-pass-filtered source time function

g fm s = F fm (g s ).
The maximum frequency f m ensures that the wavefield has a minimum wavelength λ min for most media (see Zhao et al. (2016) and [START_REF] Capdeville | Machine Learning in Geosciences[END_REF] for exceptions). Based on the common assumption that displacement wavefields are insensitive to spatial heterogeneities that are smaller than λ min , it is assumed that it is possible to introduce a finite-dimensional model space M h , whereby at least one model that minimises (22) can be found. Furthermore, M h can be designed in such a way that the wave equation (equations (1)-( 2)) can be solved in a bounded time for every one of its models. This assumption and model design yield the following inverse problem:

mh = argmin m∈M h E fm (m) (23)
that can be solved in practice, where h is a 'resolution' parameter, which is a number that characterises the discretisation made in the finite-dimensional approximation. For example, if the elastic model is spatially represented by constant velocity blocks, then h can be the size of these blocks.

h is generally directly related to λ min ; however, other information, such as the illumination angles, offset ranges or data coverage, can influence the choice of h.

One serious drawback arises when the introduction of frequency-limited data and M h make the inverse problem solvable: the relation between the true model m t and the inverse problem solution mh is unknown. Furthermore, there are many different ways to design M h for a given experiment geometry and λ min , with a different mh for each of them.

Capdeville & Métivier (2018) outlined an approach to solve this problem as follows. For a given signal maximum frequency f m , we search for a solution in M * , as defined in the previous section. We therefore define a new inverse problem for a fixed f m :

E * (m * ) = r,s T 0 (d fm s (x r , t) -u * (x r , t; x s , m * )) 2 dt , m * = argmin m * ∈M * E * (m * ) , (24)
where u * is the leading-order homogenised displacement, which is defined in equation ( 6). Here, we only use the displacement observations and zeroth-order homogenisation, which implies that no corrector is needed in the definition of M * . This simplifies u * to the displacement obtained in the effective media m * , such that u * = u fm (m * ). As shown in equation ( 12), u * is insensitive to small-scale heterogeneities.

For any model m and its homogenised version m * = H(m), it can be shown that:

E fm (m) = E * (m * ) + O(ε 0 ) . (25)
The true m t minimises E and also

E fm ; therefore, m * t = H(m t )
m * = H(m t ) . (26)
This important result shows we can access the effective true model by solving (24), which is a solvable inverse problem. In practice, we can only set up an explicit parametrisation of M * in the layered model case. We need to rely on an approximate space model M * h for the general case and use the homogenisation operator H to project the solution into M * . In practice, the construction of M * h is very similar to a classical construction of M h . However, we need to have M * ⊂ H(M * h ), which often implies that M * h is fully anisotropic. The numerical tests shown in [START_REF] Capdeville | Elastic full waveform inversion based on the homogenization method: theoretical framework and 2-d numerical illustrations[END_REF] show that the solution of the HFWI problem is indeed the homogenised true model.

Full waveform inversion based on rotational waveforms

We now consider the case where the observations d are rotational waveforms. The FWI inverse problem is very similar to the displacement case, and it also needs to be based on limited frequency band data to be solvable:

E fm (m) = r,s T 0 (d fm s (x r , t) -ω fm (x r , t; x s , m)) 2 dt , mh = argmin m∈M h E fm (m) . (27)
Although this problem is solvable, the relationship between the solution mh and true model m t is lost.

Following what was done for the displacement case in the previous section, we can define a homogenised inverse problem based on rotational data as follows:

E * (m * ) = r,s T 0 (d fm s (x r , t) -ω * (x r , t; x s , m * )) 2 dt , m * = argmin m * ∈M * E * (m * ) . (28)
However, ω fm ̸ = ω * for the leading-order approximation, which means that

E fm (m) ̸ = E * (m * ).
This implies that we cannot relate the solution m * of the inverse problem to the true model m t .

In practice, a solution to the above problem can be found, but it contains errors and/or differences compared with the homogenised true model m * t that cannot be controlled. We can correct this problem by relying on ω 0 instead of ω * :

E * (p * ) = r,s T 0 (d fm s (x r , t) -ω 0 (x r , t; x s , p * )) 2 dt , p * = argmin p * ∈M * E * (p * ) , (29)
where ω 0 (x r , t;

x s , p * ) = ω * (x r , t; x s , p * | m ) + J r : ϵ * (u * ), with J r = J(x r ), from equation (20).
Here M * contains all of the possible models p * = ((ρ * , c * ), J r , r ∈ {1, ..., N r }). and each inverted mechanical property is expanded within each element. Here we name such a parametrisation P N n (I). We use a fully anisotropic parametrisation for M * h , thereby implying that each c ij and ρ is inverted. Therefore, we have 7×n 2 ×(N +1) 2 free parameters to characterise M * h if no correctors are inverted.

Gauss-Newton iterative inversion

We iteratively update the model parameters using a Gauss-Newton inversion scheme to minimise the misfit function. The i-th iteration model m i for the rotational wavefields is updated as follows:

m i+1 = m i + ((F i ) T • F i + λ i ) -1 [(F i ) T • (d -ω(m i ))] , (30)
where F i is the partial derivative matrix and λ i is the damping parameter, which is used to stabilise the inversion. The damping parameter λ i is optimised in each numerical test. The partial derivative matrix F i with respect to the model parameters were evaluated using the adjoint technique [START_REF] Tarantola | Inversion of seismic reflection data in the acoustic approximation[END_REF][START_REF] Tarantola | Theoritical background for the inversion of seismic waveforms, including elasticity and attenuation[END_REF][START_REF] Pratt | Gauss-newton and full newton methods in frequency domain seismic waveform inversion[END_REF]. The detailed calculations of the partial derivatives for both the displacement and rotational wavefields are provided in appendix A.

We set up two criteria for convergence:

• A 98% decrease in the misfit between the observed and synthetic data, where the initial model in the inversion is a homogeneous model.

• A <1% decrease in the misfit from the previous iteration.

FWI with correctors

The model space for a rotational-based FWI includes both the effective mechanical properties and correctors. Adding the corrector parameters to the inverted parameters can be done via the same framework that was used for the elastic parameters alone. The Gauss-Newton scheme remains the same for this case; the main modifications are an increase in the size of the partial Hessian matrix to be inverted and slightly more complex programming to keep track of the inverted parameters. The adjoint sources for the rotational data are moment tensors, as explained in appendix A. Therefore, moment tensors of the adjoint sources for the rotational data should be corrected with inverted correctors [START_REF] Burgos | Homogenized moment tensor and the effect of near-field heterogeneities on nonisotropic radiation in nuclear explosion[END_REF][START_REF] Capdeville | Homogenization of seismic point and extended sources[END_REF]).

Here we rely on a two-stage approach to simplify the computations. During each Gauss-Newton iteration, the mechanical parameters are inverted to obtain the current model m i in the first stage, and the correctors are then independently inverted following [START_REF] Singh | Correcting wavefield gradients for the effects of local small-scale heterogeneities[END_REF] in the second stage. We minimise the squared error by summing over all of the sources for a given rotational receiver location x r as follows:

E(J r ) = s T 0 (d s (x r , t) -ω 0 (x r , t; x s , p * )) 2 dt ,
where:

ω 0 (x r , t; x s , p * ) = ω * (x r , t; x s , p * | m ) + J r : ϵ * (u * ) , (31)
and ω * and ϵ * are the synthetic data and strain computed in the current mechanical model m i , respectively. The misfit is simple quadratic with respect to J r and can be inverted via a simple least-squares approach:

J r = (F t r F r ) -1 F t r b, (32)
where F r is the partial derivative matrix that is limited to specific receiver, and b is the residual between the target and synthetic data. This calculation is a small computation (three and six free parameters for 2D and three-dimensional (3D) problems, respectively) and under-parameterised, which implies that damping is not needed to obtain a stable result. Fig. 2 shows the overview of the FWI with correctors employed in this study.

This two-step inversion process is expected to slow down the Gauss-Newton inversion convergence compared with a single-step inversion scheme. However, we generally do not observe a significant difference in the number of iterations compared with the cases where correctors are not needed.

Experiment design

Target model

The target model that the FWI tries to estimate is a 2D elastic medium, such as that shown in figure 1. The Ω domain is a 40×40 km region that consists of a homogeneous background with a 16 × 16 km heterogeneous area embedded at the centre. This heterogeneous area is made of 40 × 40 square anomalies, with each square possessing constant mechanical properties that are randomly generated. The background model is isotropic with shear and compressional velocities (V S and V P , respectively), and a density of 2.8 km/s, 5.0 km/s and 3.0 t.m -3 , respectively. The random mechanical properties of the heterogeneous region are generated using a flat probability distribution that is within 60% of the elastic tensor and density of the homogeneous background values.

Spectral element simulation setup

We use the spectral element method [START_REF] Komatitsch | The spectral element method: an effective tool to simulate the seismic response of 2D and 3D geological structures[END_REF] to solve the wave equations (equations ( 1)-( 2)) for both generating the data to be inverted and solving the inverse problem.

The absorbing boundaries for each of our experimental cases are based on perfectly matched layers [START_REF] Festa | The newmark scheme as velocity-stress time-staggering: an efficient implementation for spectral element simulations of elastodynamics[END_REF] that are used all edges of the domain. Two different spectral element meshes are used:

• We rely on a square 100 × 100 element mesh to generate the data to be inverted. We use a polynomial degree of 10 in each tensorial direction of each element. The mesh is designed such that all of the material discontinuities are honoured by an element boundary, which is important to obtain good accuracy.

• A different spectral element mesh than the previous one is used when computing the full wavefields from the sources and adjoint sources to compute the partial derivatives. This mesh consists of 50 × 50 elements with a polynomial degree of 4.

The different spatial meshes for generating the target and synthetic data sets during the inversion process allows us to avoid an "inverse crime" because the numerical errors due to spatial discretisation are different. Furthermore, we do not introduce any a priori knowledge of the heterogeneity geometry in the inversion modelling mesh.

FWI acquisition geometry

The source and receiver geometries used in our numerical experiments are shown in Fig. 3. The sources are located all around the inverted area, with each consisting of a single vector force that produces a Ricker wavelet. This wavelet possesses a 0.25-Hz central frequency and 0.7-Hz maximum source frequency. The associated minimum wavelength in the background medium is 4 km, which is 10 times larger than the lengths of the small anomalies in the target model.

We use two different geometries for the receivers:

• Configuration A (Fig. 3): The receivers are located around the inverted area and are far from the heterogeneities.

• Configuration B (Fig. 3): The receivers are located within the inverted area and are in contact with the heterogeneities.

This setup is designed to investigate the impact of small-scale heterogeneities on FWI using rotational wavefields recorded using two receiver configurations: closely located receivers (configuration B) or not (configuration A).

Representation and evaluation of the FWI results

We first need to define what will be presented and how the results will be assessed prior to conducting the numerical experiments.

The FWI results are always a full elastic tensor c after homogenization. We present the results by first projecting the anisotropic elastic tensor c * to the nearest isotropic tensor c *,iso following [START_REF] Browaeys | Decomposition of the elastic tensor and geophysical applications[END_REF]. The compressional-and shear-wave velocities are defined as follows:

V * P (x) = (c *,iso 1111 (x)/ρ(x)) , (33)
and:

V * S (x) = (c *,iso 1212 (x)/ρ(x)) , (34)
respectively, and the total anisotropy is defined as:

aniso(x) = ijkl (c iso ijkl (x) -c ijkl (x)) 2 ijkl (c iso ijkl (x)) 2 . ( 35
)
One objective of our work is to compare the rotational-based FWI results to the displacementbased FWI results, particularly to assess the accuracy and resolution of the two approaches. A classical way to perform this analysis is to analyse the diagonal of the partial Hessian and compute the point spread functions [START_REF] Fichtner | Resolution analysis in full waveform inversion[END_REF]. However, our use of a homogenisation step in the FWI, which is a non-linear step, makes it impossible for us to work directly on the partial Hessian. We overcome this difficulty by performing a brute force statistical analysis. We generate 200 different data sets for both the displacement and rotational receiver configuration cases by adding random noise to the data. This added noise possesses a standard deviation of 10% of the original data set and spans the same frequency band as the original data set. Two hundred inversions are then performed, one for each data set. We define the inversion error as follows:

E 2 m (ε 0 , x) = n=N n=1 ij (c * ,n ij (x) -c * ,t ij (x)) 2 N ij (c * ,t ij (x)) 2 , ( 36
)
where c * ,n and c * ,t are the homogenised FWI model results for noise number n and the true homogenised elastic tensor, respectively, for a given ε 0 . We finally define the inverted area-averaged inversion error as:

Ē2 m (ε 0 ) = I n=N n=1 ij (c * ,n ij (x) -c * ,t ij (x)) 2 dx N I ij (c * ,t ij (x)) 2 dx , (37)
Analysing Ē2 m (ε 0 ) as a function of ε 0 allows us to find an estimate of the smallest ε 0 that can be used. This smallest ε 0 can be seen as the overall resolution limit of the inversion. Constraining ε 0 in this manner highlights the compromise between the resolution and error in the inversion: using a smaller ε 0 would generate higher-resolution images with larger error, whereas using a larger ε 0 instils good confidence in the accuracy of the imaged structures at a lower spatial resolution.

We also assess the inversion results using three statistical measures: the standard deviations, mean values and resolutions (correlations) of V * S and V * P (obtained from equations ( 33) and (34), respectively), and density. The mean values are the same as the noise-free inversion results and are used to compare the inversion results and target model for a given ε 0 . We note that the mean values and noise-free inversion results are not exactly the same, but it does not change the point of the paper and we ignore this differences. For a given ε 0 , the covariance of mechanical property p * (V * P , V * S or the effective density) is written as:

Cov ε 0 (x, y) = N n=1 (p * n (x) -µ p (x))(p * n (y) -µ p (y)) N , (38)
where µ p is the mean value of p, and N is the number of samples. The deviation σ ε 0 can be derived in terms of the covariance as follows:

σ ε 0 (x) = Cov ε 0 (x, x) , (39)
such that the resolutions R ε 0 are then written as follows:

R ε 0 (x, y) = Cov ε 0 (x, y) σ ε 0 (x)σ ε 0 (y) . ( 40
)
Here we provide the V * S results for our presented inversion and statistical analyses. Comprehensive results, which include the V * P and effective density results, are provided in the appendix. Note that both σ ε 0 and R ε 0 depend on ε 0 . This co-dependence introduces a trade-off between σ ε 0 and R ε 0 . Therefore, the resolution limit and overall error cannot be measured independently for this inversion process.

A simple forward modelling example of a small-scale heterogeneity effect on the data

We first illustrate the different impacts of small-scale heterogeneities on displacement and rotational measurements prior to introducing the FWI results. Example of rotational and displacement waveforms for the two source-receiver configurations are shown in Fig. 4. Here, four receivers are located in close proximity to each other at two distant locations, A and B, whereby they are spaced at 100-m horizontal and vertical intervals from each other. We acquired both the displacement and rotational waveforms using a source that generates a Ricker wavelet with a 0.7-Hz maximum frequency and 4-km minimum wavelength. It should be noted that station intervals and lengths of the anomalies are much smaller than the 4-km minimum wavelength. The displacement waveforms from the four receivers at location A, which is within the heterogeneous area, are the same, as expected, whereas there are observable differences among the rotational waveforms. These differences in the rotational waveforms are due to the corrector in equation ( 20). Conversely, the rotational waveforms from the four receivers at location B, which is ∼2 km from the heterogeneous area, are the same. This is because the effect of the corrector is limited to the heterogeneous area, with no small-scale effects propagating outside of the heterogeneous area.

This example clearly highlights that a smooth tomographic model would have a hard time explaining data that are contaminated by small-scale heterogeneities. It also underscores why it is necessary introduce correctors to account for such effects.

FWI results without correctors

Receivers outside of the heterogeneous area

Here we present displacement and rotational inversion results for receivers positioned outside the heterogeneous area ("configuration A" in Fig. 3).

The displacement-based inversion results, which are presented in Fig. 5, show that the misfit decreases steadily throughout all of the Gauss-Newtown iterations and satisfies the convergence criteria at iteration 19. The strong fit between the target and synthetic data indicates that the final model can adequately explain the target data. We have used ε 0 = 1.0 to compute the homogenised models, which is a good compromise between resolution and accuracy, as we will discuss later.

The residual is calculated as the deviation from a background V S value of 2.8 km/s. A comparison of the modelled V * S structures shows that the final model exhibits a broadly similar pattern to that in the target model, with the exception of some high-velocity anomalies in the target model that are not well-recovered in our final inverted model.

The rotational inversion results are shown in Fig. 6. The chosen receiver configuration makes it possible to conduct a rotational-based FWI without any expected small-scale heterogeneity effects because the rotational waveforms are only influenced by small-scale heterogeneities that are in close proximity to the receiver locations, as explained in equation ( 20) and demonstrated in Fig. 4. The presented waveform comparisons show a strong fit between the waveforms for the target and final models, which indicates that the rotational-and displacement-based FWI models are comparable in their abilities to explain the target data. The V S residuals between the target and final models for the rotational waveform case are also comparable to those for the displacement case.

Our assessment of the optimal ε 0 value, which is based on equations ( 37) and ( 37), is shown in Fig. 7. We used the mean model from 200 samples, as opposed to a single solution, for the model error Ēm , which was computed for models that used ε 0 in the 0.3-2.9 range. The model error Ēm increases exponentially as ε 0 decreases, with the transition to a rapid increase in Ēm occurring around ε 0 = 1.0 for both the displacement-and rotation-based FWI. This result suggests that the displacement-and rotation-based FWI have almost the same resolution limit. We therefore select ε 0 = 1.0 as the optimal choice that provides the best overall resolution and accuracy of the final model results.

Our statistical comparisons between the displacement and rotational cases, which consist of the error of mean values, standard deviation, and resolution between the target and final models, are based on ε 0 = 1.0 and shown in Fig. 8. The resolutions were computed by selecting a single grid point y i and computing the correlations against all of the x points based on equation (40). Strong correlations (approaching unity) indicate that these points cannot be determined separately; therefore, the areas with strong correlations correspond to the spatial resolution of the final model. All of the statistical quantities exhibit comparable V * S trends between the displacement and rotational cases.

Receivers within the heterogeneous area

Here we present the displacement-and rotational-based inversion results using receivers positioned within the heterogeneous area ("configuration B" in Fig. 3). The forward modelling test presented in section 4.2 indicates that it is necessary to invert for the corrector to obtain good results. Nevertheless, we ignore this fact here and perform the displacement-and rotational-based FWI without correctors. The displacement-based inversion results, which are presented in Fig. 9, show that the misfit steadily decreased throughout all of the iterations and satisfied the convergence criteria at iteration 11. A strong fit between the target and synthetic data is confirmed, thereby indicating that the final model can explain the target data. Our model comparisons reveal that the errors between the target and final models for the receivers within the heterogeneous area are slightly better than the model comparisons when all of the receivers are positioned outside of the heterogeneous area. The nearby heterogeneities have a small effect on the displacement measurements, and subsequently the FWI process, as theoretically expected (equation 12).

The rotational-based inversion results are presented in Fig. 10. The misfit reductions satisfied the convergence criteria at iteration 46, despite a temporary increase in the misfit at iteration 13.

The comparison between the target and synthetic models suggests that the final model is unable to fully explain the target model, which indicates that the inversion did not converge to the best solution in this case. Furthermore, the V S residuals in the rotational case are larger than those in the displacement case, despite having the same source and receiver configurations. We did not observe significant misfits in the waveforms and model parameters when the receivers were positioned outside of the heterogeneous area. The larger residuals between the target and final models and larger misfits between the waveforms are likely associated with the presence of smallscale heterogeneities in close proximity to the receiver locations, as theoretically expected.

The trade-off between the error of the mean model Ēm and ε 0 are shown in Fig. 11. The overall trend of Ēm is similar to that when the receivers are placed within the heterogeneous area.

However, the Ēm results show that errors from the rotational-based FWI are consistently larger than those from the displacement-based FWI, regardless of the ε 0 value. The effects of small-scale heterogeneities can be clearly seen, even if the scale is limited to very large scales.

Fig. 12 presents our statistical comparison of the displacement-and rotational-based FWI results for ε 0 = 1.0. The error of the mean value of V S for the rotational case indicates that the inversion could not determine the model parameters as accurately as in the other cases. The deviations and resolutions also show peculiar results for the rotational case. Unlike the displacement case, the deviations for the rotational case are not uniform, and the largest anomalies are dis-tributed. Furthermore, the spatial resolutions for the rotational case are no longer valid, as there are strong correlations at distant locations.

FWI results with correctors

Here we present the rotational-based FWI results, with correctors, for receivers located within the heterogeneous area ("configuration B" in Figure 3). The correctors are obtained during the inversion based on equation ( 32). Incorporating correctors in the rotational-based FWI procedure is theoretically expected to enhance the results compared with the FWI procedure without correctors.

The inversion results, which are presented in Fig. 13, show that the misfit decreased in all of the iterations and met the convergence criteria at iteration 12. The strong fit between the waveforms for the target and synthetic models indicates that the final model with correctors can adequately explain the target data, whereas the final model without correctors is unable to achieve such a fit.

Our model comparisons yield smaller residuals between the models with correctors than those without, thereby suggesting that including correctors led to improvements in the final model. Ēm than that for the displacement case was observed when correctors were not included in the rotational-based FWI (Fig. 11). However, the rotational-based FWI with correctors possesses an error that is almost the same as that for the displacement case. The resolution limits for both the rotational and displacement cases can be chosen to be 1.0 based on the Ēm results, which indicates that there is no difference in the resolution limit between the rotational and displacement cases.

Fig. 15 presents the statistical results based on ε 0 = 1. As discussed in the previous section, the FWI without correctors resulted in large errors in the mean values and peculiar outcomes for the deviations and resolutions. However, we obtained similar statistical results to those in the case for receivers positioned outside of the heterogeneous area when correctors were included in the inversion process. This indicates that the inclusion of correctors in the FWI procedure mitigates the effects of small-scale heterogeneities.

Table 1 provides a summary of the misfits, which are defined by equation ( 37) for a single model, in each case. Here we use single solution, as opposed to the mean model, to compute when rotational data are used and all of the receivers are located within the heterogeneous area.

However, incorporating correctors in the inversion process significantly reduces the misfit, with a reduction from 0.095 (without correctors) to 0.049 (with correctors) observed. The misfit of 4.86 for the rotational-based FWI with correctors is comparable to the FWI results when receivers are placed outside of the heterogeneous area.

We present comprehensive inversion and statistical results, including the V * P , density * and anisotropy * results for three cases (FWI using displacement receivers positioned within the heterogeneous area; FWI, without correctors, using rotational receivers positioned within the heterogeneous area; and FWI, with correctors, using rotational receivers positioned within the heterogeneous area), in the appendix.

A short investigation of additional FWI results using strain receivers is also presented in appendix. The conclusion of the strain-based FWI results is same as that for the rotational case:

the FWI does not converge to the best model when the effects of small-scale heterogeneities are present. It is therefore necessary to include correctors in strain-based FWI to effectively mitigate any small-scale effects.

JOINT DISPLACEMENT-AND ROTATIONAL-BASED FWI

Gradient sensors are often collocated with seismometers for more comprehensive field studies.

Here we perform a numerical experiment where both displacement and rotational sensors are used. We used the same receiver and source configurations as in the previous section, whereby the receivers are positioned within the heterogeneous area, and assume that both displacement and rotational data are recorded at each receiver location. We therefore incorporated data from 52 sensors (26 rotational and 26 displacement sensors) for the inversion. We utilised a least-squares misfit function by combining equations ( 22) and ( 29) in this numerical experiment. A data covariance matrix that is different from the identity matrix is used to balance the significant amplitude differences between the displacement and rotational data to ensure that both data types are equally weighted in the inversion [START_REF] Tarantola | Inverse problems= quest for information[END_REF]. The results are evaluated using the model misfit (equation ( 37)), which are based on ε 0 = 1, and presented in table 2. The results show that the misfits in the model parameters are minimised when correctors are incorporated. These misfits are even smaller than those in the case when only displacement data are used; this is likely related to the inclusion of additional data in the joint displacement and rotational data case. This numerical experiment indicates that the inclusion of correctors in the inversion process is necessary, even when translational motions are used with wavefield gradients.

USING CORRECTORS TO OBSERVE SUB-WAVELENGTH HETEROGENEITIES

We have demonstrated that the rotational-based FWI is hampered by the sensitivity of the rotational wavefield to small-scale structures rather than improving the spatial resolution. Although incorporating correctors in the FWI process can help mitigate the effects of small-scale heterogeneities, it can only produce tomographic images that are at the same resolution as displacementbased FWI results, as shown in the statistical evaluations. Therefore, rotational-based FWI can only estimate the effective mechanical properties at best, and are therefore unable to capture finescale heterogeneities in the tomographic images, despite the sensitivity of the rotational wavefield to small-scale heterogeneities. However, the inverted correctors in our new FWI framework may contain information about any small-scale structures that may be present. Here we demonstrate how small-scale structures can be observed in the inverted correctors.

We randomly place 500 rotational receivers within a 1.8 × 1.8 km area (Fig. 16) for this numerical test; the absolute values of the inverted correctors at the final iteration for each component are illustrated in Fig. 16. The diagonal components of the correctors appear to only possess high values in the vicinity of the vertices of the small-scale anomalies, whereas the off-diagonal components of the correctors more clearly show the defined boundaries of the small-scale anomalies. Therefore, the inverted correctors in our new FWI framework can provide indications of the presence of strong small-scale discontinuities in the subsurface that are difficult to image using conventional seismic tomography methods.

DISCUSSIONS AND CONCLUSIONS

We have first shown the necessary difference in model parametrisation between displacement-and wavefield gradient-based FWI. The leading-order approximation in the two-scale homogenisation mathematical framework highlights that the displacement-based FWI only requires effective elastic mechanical properties to explain the data, whereas the wavefield gradient-based FWI also requires the corrector in the model space to explain the data. Current FWI applications implicitly rely on the FWI solution space being within the M * subspace, such that the solution space is only represented by the effective elastic mechanical properties. However, FWI applications that are based on wavefield gradients, such as rotational and strain wavefields, require the solution space to extend to the entire M * , thereby including the effective elastic mechanical properties and correctors, as theoretically expected. Our numerical experiments have confirmed that this extension of the solution space is necessary, as the effective true model was only retrieved in the wavefield gradient-based FWI when correctors were included.

Another important result is that the resolution limit is close to ε 0 = 1 for both the displacementand rotational-based FWI, such that there is no difference in terms of resolution. The model error exponentially increases as ε 0 decreases, with a rapid increase observed when ε 0 < 1. The resolution images that were obtained for different ε 0 values are also helpful for assessing the resolution limit. Fig. 17 shows the resolution images for all of the cases and ε 0 = 0.5 , 1 , 1.5. The size of the area with a strong correlation is almost same between ε 0 = 0.5 and ε 0 = 1 for all of the cases, and then increases when ε 0 = 1.5. We therefore suggest that the resolution limit is around ε 0 = 1 for both the displacement-and rotational-based FWI. Finally, we note that ε 0 = 1 is generally not a valid value. Other configurations, such as either more receivers or more scattering, may allow this value to be lowered to around ε 0 = 0.5. ε 0 = 0.5 is sufficient to model data in most geological media, such that it is unnecessary to lower ε 0 below this limit.

We have found no advantage of using wavefield gradients instead of displacement data during the FWI process in terms of accuracy and resolution. We note that this conclusion is only true for excellent data coverage. We have not investigated cases with partial data coverage; therefore, it is still possible that gradient data could make a significant difference in such contexts.

We have shown that the inverted correctors can serve as indicators of small-scale heterogeneities, which can lead to high-resolution mapping of complex geological structures, as long as such features are in the vicinity of the receivers. Therefore, we may be able to use arrays of many receivers to locate abrupt subsurface discontinuities, such as faults, if they are located near the sensors. DAS cables may be a viable option since they can easily obtain a significant number of very short interval strain data, which can assist in imaging small, abrupt, subsurface discontinuities via the corrector.

It is important to note that we have only inverted part of the correctors in this study. Indeed, the parametrisations of M * h that we have been using contain small-scale heterogeneities (continuity between elements is not enforced). Therefore, the correctors shown in Fig. 16 cannot be directly compared with the corrector that would be obtained from the direct homogenisation problem. Such a comparison requires the correctors that are derived from the inversion process to be combined with those obtained from applying the homogenisation operator H to the inverse solution.

Although the FWI with correctors can mitigate the effects of small-scale heterogeneities, it cannot completely eliminate these effects. Our approach assumes that evaluating up to the firstorder term in equation ( 6) is sufficient, such that the effect of the higher-order terms are negligible and can be ignored. However, it should be noted that the effects of these higher-order terms become more pronounced as the elastic contrasts increase. Therefore, our method is still considered to be affected by small-scale heterogeneities to some extent. Nevertheless, all of the synthetic tests in this study show that a better solution is obtained using correctors.

This study demonstrates that there is no intrinsic gain in using gradient sensors instead of conventional seismometers for FWI studies, as rotational-based FWI results can only recover information of the same quality as translation-based FWI results if receiver-coupling correctors are included in the inversion process. Nevertheless, these correctors contain sub-wavelength heterogeneities that can potentially be exploited to obtain a resolution at the sensor array scale.

C2. Receivers within the inversion area for the FWI without correctors

Figures A.8 and A.9 illustrate the inversion results when the receivers are positioned within the heterogeneous area and correctors are not included. The misfit decreased until iteration 2, but then increased in subsequent iterations. Furthermore, there is a large degree of misfit between the synthetic and target waveforms, as well as large residuals between the target and final models.

Similar to the rotational-based FWI case, the strain-based FWI without correctors fails to converge to the best model and is therefore incapable of explaining the target data when the effects of smallscale heterogeneities are present in the target data.

C3. Receivers within inversion area for the FWI with correctors

We then performed the FWI with correctors using receivers that were positioned within the heterogeneous area. Table 3 provides the misfits in the elastic tensors based on equation ( 37) for a single model in each case, with the incorporation of correctors leading to an improvement in the accuracy of the inverted elastic parameters. Therefore, it is necessary to incorporate correctors in the FWI process to obtain accurate results when the strain observations are used. using receivers positioned within the heterogeneous area ("Displacement inside"). Rotational-based FWI results, without correctors, using receivers positioned within the heterogeneous area ("Rotation inside").

Rotational-based FWI results, with correctors, using receivers positioned within the heterogeneous area ("Rotation inside (corrector)").

  is a solution of the homogenisation FWI (HFWI) problem (24). If we assume this solution is unique (this a strong assumption; see the discussion in[START_REF] Capdeville | Elastic full waveform inversion based on the homogenization method: theoretical framework and 2-d numerical illustrations[END_REF]), then the solution of the HFWI problem is related to the true model as follows:

  The last inverse problem definition leads to E fm (m) = E * (p * ) + O(ε 0 ), which makes sure that we can relate the inverse problem solution m * = p * | m to the effective true model m * t . We use the same inversion framework as in Capdeville & Métivier (2018), whereby we spatially parametrise M * h by dividing the squared inversion domain I ∈ Ω into non-overlapping n × n square elements: I = ∪ n 2 e=1 I n e . An example of inversion domain I with 22 × 22 elements (I 22 e ) is shown in Fig. 1. A N × N degree polynomial basis is then used for each element I n e ,

Fig. 14

 14 Fig.14 shows relationship between the error of the mean model Ēm and ε 0 . A larger error

  Zhao, M., Capdeville, Y. & Zhang, H., 2016. Direct numerical modeling of time-reversal acoustic subwavelength focusing. Wave Motion, 67, 102-115. Appendix the corresponding residuals between the target and final models are shown in Fig. A.7, with the residual values being comparable to those obtained via the displacement-based FWI.

  Fig. A.10 shows the misfit reduction and waveform comparisons. The misfit steadily decreased throughout the iterative process, and a strong fit between the synthetic and target waveforms was observed, thereby indicating that FWI with correctors can adequately explain the target data. Fig.A.11 presents the model comparisons and residuals for each parameter, whereby the residuals show improvements in all of the parameters compared to the results without correctors.

Figure 1 .Figure 2 .Figure 4 .Figure 5 .Figure 6 .Figure 7 .Figure 9 .Figure 10 .Figure 11 .

 12456791011 Figure 1. (A) Schematic diagram of the I 22 e mesh structure. The shaded grey area shows the target area to be inverted. Black lines define the element boundaries. (B) V S (left), V P (centre) and density (right) structures of the target model.

Figure 12 .Figure 13 .

 1213 Figure 12. Statistics of the displacement-and rotational-based FWI results, without correctors, using receivers positioned within the heterogeneous area. (A) Error of the mean models of V * S . (B) Standard deviation of V * S . (C) Resolution of V * S . The error and standard deviation values are presented as fluctuations from a background value of 2.8 km/s.

Figure 16 .Figure 17 .

 1617 Figure 16. Numerical demonstration to show indications of small-scale structures via the inverted correctors. (A) Source (black stars) and receiver (black diamonds) configuration for the numerical experiment. Five-hundred rotational receivers are positioned at location A. (B) V S and (C) V P heterogeneities around location A. Black dots represent the receiver positions. Inverted diagonal components, (D) J r,11 and (E) J r,22 . (F) Inverted off-diagonal component, J r,12 .

Figure A. 1 .

 1 Figure A.1. Comparison of the FWI results. The V * S , V * P , density * and anisotropy * structures, which are homogenised using ε 0 = 1.0, are shown. Target model ("Target model"). Displacement-based FWI results

Figure

  Figure A.2. Comparison of the error of mean models of the FWI results: The V * S , V * P and density * structures, which are homogenised using ε 0 = 1.0, are shown. Mean model for displacement receivers positioned within the heterogeneous area ("Displacement inside"). Mean model, without correctors, for rotational receivers positioned within the heterogeneous area ("Rotation inside"). Mean model, with correctors, for rotational receivers positioned within the heterogeneous area ("Rotation inside (corrector)").

Figure

  Figure A.3. Comparison of the standard deviations in the FWI results from the target model. The V * S , V * P and density * structures, which are homogenised using ε 0 = 1.0, are shown. Deviation for displacement receivers positioned within the heterogeneous area ("Displacement inside"). Deviation, without correctors,for rotational receivers positioned within the heterogeneous area ("Rotation inside"). Deviation, with correctors, for rotational receivers positioned within the heterogeneous area ("Rotation inside (corrector)").

  Figure A.4. Comparison of the resolutions of the FWI results. The V * S , V * P and density * structures, which are homogenised using ε 0 = 1.0, are shown. Resolution for displacement receivers positioned within the heterogeneous area ("Displacement inside"). Resolution, without correctors, for rotational receivers positioned within the heterogeneous area ("Rotation inside"). Resolution, with correctors, for rotational receivers positioned within the heterogeneous area ("Rotation inside (corrector)").

Figure

  Figure A.5. Source (red stars) and receiver (black diamonds) configurations for the synthetic tests of the strain-based FWI. (A) Configuration with 8 sources and 26 receivers positioned outside of the heterogeneous area. (B) Configuration with 8 sources positioned outside of the heterogeneous area and 150 receivers positioned within the heterogeneous area.

Figure A. 6 .

 6 Figure A.6. Strain-based FWI results using receivers positioned outside of the heterogeneous area. (A) Misfit reduction. (B) Waveform comparisons between the target and synthetic data. Source-receiver positions are shown in the left panel, with the green source (star) and receiver (diamonds) positions indicating the source-receiver pairs used to generate the presented target (black) and synthetic (red) waveforms in the right panel. The waveforms are normalised by the maximum amplitude of the selected receivers.

Figure A. 8 .

 8 Figure A.8. Strain-based FWI results, without correctors, using receivers positioned within the heterogeneous area. (A) Misfit reduction. (B) Waveform comparisons between the target and synthetic data. Sourcereceiver positions are shown in the left panel, with the green source (star) and receiver (diamonds) positions indicating the source-receiver pairs used to generate the presented target (black) and synthetic (red) waveforms in the right panel. The waveforms are normalised by the maximum amplitude of the selected receivers.

  Figure A.10. Strain-based FWI results, with correctors, using receivers positioned within the heterogeneous area. (A) Misfit reduction. (B) Waveform comparisons between the target and synthetic data. Sourcereceiver positions are shown in the left panel, with the green source (star) and receiver (diamonds) positions indicating the source-receiver pairs that were used to generate the presented target (black) and synthetic (red) waveforms in the right panel. The waveforms are normalised by the maximum amplitude of the selected receivers.

Table 1 .

 1 Model misfits, which are based on equation (37), for the numerical FWI results using displacement and rotational data.

	Case

the model error Ēm . The misfits for the FWI without correctors indicate that the misfit is largest

Table 2 .

 2 Model misfit based on equation (37) for the FWI results of the joint displacement and rotational data experiments.

	Case	Ēm (ε 0 = 1)
	FWI without correctors
	Joint displacement and rotational data	0.066
	FWI with correctors
	Joint displacement and rotational data	0.038

Table 3 .

 3 Model misfit based on equation (37) for results of numerical experiment using strain data.

	Case

ACKNOWLEDGEMENTS

This work was supported by the Japan Society for the Promotion of Science (JSPS) through a Grant-in-Aid for JSPS Research Fellowship (22KJ2397) and JSPS Overseas Challenge Program for Young Researchers. This work was also supported by JSPS KAKENHI (Grant Numbers JP20H01997). This work was completed in part with resources provided by the GLiCID (www.glicid.fr) Computing Center. Yann Capdeville's work was partly supported by the ANR ASIRISplus-SDS grant ANR-19-CE22-0015. This work has benefited from discussions within the SPIN-ITN European project.

DATA AVAILABILITY

This study used only synthetic data.

Appendix A: Adjoint partial derivative calculations

A.1 Preliminary remarks

Here we use elastic wave equation in the frequency domain for simplicity:

-ω 2 ρ(x)u i (x, ω) -∂ j c ijkl (x)∂ k u l (x, ω) = f i (x, ω) .

(A.1)

The Green's function for a single point force at x s is described by f p i (x, ω) = δ(x -x s )δ ip , which satisfies the equation:

The Green's function for a moment component point source f pq i (x, ω) = -δ ip δ jq ∂ j δ(x -x s )

satisfies the equation:

-ω 2 ρ(x)G pq i (x, ω; x s ) -∂ j c ijkl (x)∂ k G pq l (x, ω; x s ) = -δ ip δ jq ∂ j δ(x -x s ) . (A.3)

Finally, we note that G M i for the moment tensor force described in equation ( 4) can be written as:

We then omit ω to denote the reciprocity theorem as follows: .5) and:

We also have:

A.2 Partial derivatives for a single scattering point in x

The Born approximation reads:

Equation (A.8) can be reformulated via integration by parts and the reciprocity principle as follows:

We can then deduce the equation for a gradient component:

We substitute equation (A.7) into equation (A.11) to obtain:

The coupling tensor J, with [K qi ] pn = J qipn , can be written as:

We can utilise the coupling tensor J and Levi-Civita symbol θ to derive the adjoint source of the moment tensor for rotational wavefields as:

(A.13) Furthermore, we can derive the adjoint source of the moment tensor for the strain wavefields as follows:

We stored the calculated full wavefields for both sources and adjoint sources on disk and computed the partial derivatives by convolving each source-receiver pair.

Appendix B: Supplementary inversion results

We present comprehensive inversion and statistical results for three different cases: displacementbased FWI using receivers positioned within the heterogeneous area; rotational-based FWI, without correctors, using receivers positioned within the heterogeneous area; and rotational-based FWI, with correctors, using receivers positioned within the heterogeneous area. We show the target Appendix C: Strain-based FWI

Here we present numerical FWI experiments using strain data. We followed the inversion scheme that is explained in section 3 for this numerical experiment. We performed numerical tests with two different receiver configurations. One configuration involved placing all of the receivers outside ). We only used the vertical component along the vertical axis of the strain data at each receiver location in these numerical experiments. This setup was designed to simulate a DAS-based FWI; however, we did not consider the gauge length and strain rate of DAS measurements in our numerical experiments.

C1. Receivers outside of the inversion area