Jakob F Ulf

Aarsnes

Jean Auriol

Florent Di Meglio

Roman J Shor

Estimating friction factors while drilling

We design and test an algorithm for estimating the bit angular velocity and side forces on a drill string using only topside measurements fed into an adaptive observer. To derive the algorithm gains the drill string is modeled as a distributed wave equation coupled with Ordinary Differential Equations (ODEs) at the boundaries which represent the lumped dynamics of the top-drive and BHA, respectively. Here, the algorithm uses only measurements from the top-drive, while the non-linear side-forces are lumped at the BHA. The obtained observer gains are then implemented on a previously validated torsional drill string model with a distributed BHA model and distributed Coulomb side forces implemented as an inclusion. This design approach combines a high fidelity model with recent theoretical developments on estimation of PDEs to find appropriate feedback gains to ensure a fast and robust tuning. The feasibility of the approach is illustrated by showing convergence in estimated friction factors to the true values when testing on simulated data. When testing on full scale field data, friction factors converge to the same value for different initial values which indicates both good estimates and a robustness of the approach. The estimates produced by the algorithm can be displayed to a driller in real time in an advisory system, and the result can be built on to help optimize the drilling operation, detect faults and unwanted incidents, aid on-site decision making, and improve control of directional drilling.

Introduction

Exploration and production of oil and gas reserves in the deep subsurface necessitates the precise placement of a narrow wellbore, between 10 and 50 cm diameter, of lengths up to 10,000 m through a varied subsurface geological environment. Modern wellbores are rarely straight and typically must follow preplanned wellplans, ranging from simplier horizontal or deviated wells to complex three dimensional paths. In these complex wellbores, torque and drag due to friction between the drillstring and borehole wall become significant and increase the complexity of steering operations. Steering is typically achieved through the use of rotary steerable systems -which use actuated downhole tools to physically point or push the bit in a predefined direction -or bent motor assemblies -which use a permanent bend in a downhole mud motor to impart a deviation in hole angle if the motor bend is held in a fixed orientation. The latter system is often used for land wells due to its low cost, but presents challenges for the directional drillers due the complexity of the drillstring dynamic behavior and uncertainty in borehole torque and drag modelling.

Drillstring dynamic behavior often manifests itself as drill string vibrations, and has been extensively studied.

In particular, the torsional oscillations known as stick slip are considered to be one of the most prevalent vibrations. These stick-slip oscillations are characterized by a series of stopping -"sticking" -and releasing -"slipping"events of the bit. These models, if properly formulated, lend themselves to estimation problems and can be used to estimate friction parameters, and in turn, directional abilities of the drilling system.

A distributed model was used to characterize stick slip oscillation as caused by Coulomb friction induced side forces in [START_REF] Aarsnes | Torsional vibrations with bit off bottom: Modeling, characterization and field data validation[END_REF]. Such stick slip can also occur off-bottom, and does not require a velocity weakening in the bit-rock interaction. Off-bottom stick slip is a well known phenomenon from the field, and when mentioned in literature is hypothesized to be caused by a negative difference between static and kinetic along-string Coulomb-type friction [START_REF] Brett | Uses and Limitations of Drillstring Tension and Torque Models for Monitoring Hole Conditions[END_REF][START_REF] Halsey | Drillstring Vibrations: Comparison Between Theory and Experiments on a Full-Scale Research Drilling Rig[END_REF][START_REF] Zhao | Abnormal Down Hole Pressure Variation by Axial Stick-Slip of Drillstring[END_REF]. This cause is distinct [START_REF] Aarsnes | Stick-slip and torsional friction factors in inclined wellbores[END_REF] from stick-slip caused by the regenerative effect in the bit-rock interaction [START_REF] Aarsnes | Dynamics of a distributed drill string system: Characteristic parameters and stability maps[END_REF][START_REF] Detournay | A phenomenological model for the drilling action of drag bits[END_REF][START_REF] Richard | Self-excited stickslip oscillations of drill bits[END_REF].

The distributed effects of the drill string becomes an increasingly prominent feature of torsional dynamics in long wells, and perhaps especially wellbores with highinclination laterals. The nonlinear nature of the Coulomb friction can excite a wide range of frequencies where higher order modes become essential for representing the dynamics of the system, in particular for long wells. Hence, lumped approximations of the drill string easily fall short, and it is desirable to employ distributed model representations of the torsional dynamics, along the lines of [START_REF] Aarsnes | Linear stability analysis of self-excited vibrations in drilling using an infinite dimensional model[END_REF][START_REF] Aarsnes | Axial and torsional self-excited vibrations of a distributed drill-string[END_REF][START_REF] Germay | Multiple mode analysis of the self-excited vibrations of rotary drilling systems[END_REF][START_REF] Aarsnes | Torsional vibrations with bit off bottom: Modeling, characterization and field data validation[END_REF].

The approach taken in the present work is to use a distributed model of the drill string with the lumped representation of the BHA used in [START_REF] Aarsnes | Torsional vibrations with bit off bottom: Modeling, characterization and field data validation[END_REF] where the Coulomb friction is given as a source term implemented as an inclusion, with the states propagated according to the Filippov solution [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF]. The distributed nature of the model enables it to capture a large part of the possible dynamics for a significantly extended range of well lengths and surveys compared to the more standard lumped approximations. At the same time, the relative mathematical simplicity of the formulation, a 1-D wave equation with a source term with ODEs as boundary conditions, enables us to employ new results from control theory [START_REF] Di Meglio | Robust output feedback stabilization of an ODE-PDE-ODE interconnection[END_REF] to effectively derive a soft sensor as an Luenberger Observer. Using only topside measurements, this soft sensor can in turn estimate kinetic and static friction factors as well as the distributed states of the system including drill string angular velocity and angular strain. This have several applications 1. Effectively estimating static and kinetic friction factors while drilling can help optimize the drilling operation, detect faults and unwanted incidents and aid on-site decision making. 2. An adaptive on-line observer opens up several promising new venues for avoiding torsional oscillations. 3. Improved estimates of drill string orientation have important applications for directional drilling [START_REF] Skyles | Converting static friction to kinetic friction to drill further and faster in directional holes[END_REF][START_REF] Cheng | Theoretical model for calculating pulling loads for pipes in horizontal directional drilling[END_REF].

This paper begins with a derivation of the torsional drillstring model in Section 2, including validation with two sets of field data. The semi-lumped approximation needed for observer design is then given in Section 3, and the observer design itself in Section 4. Finally, in Section 5, convergence behavior of the observer is demonstrated on both synthetic and field data.

Model

In this section we present the model under consideration. The validity of the proposed model is illustrated by considering the open-loop fit of the model to full scale field data at the end of the Section.

Torsional dynamics of drill string

The distributed model we use in this paper is similar to the one proposed in [START_REF] Aarsnes | Dynamics of a distributed drill string system: Characteristic parameters and stability maps[END_REF][START_REF] Aarsnes | Linear stability analysis of self-excited vibrations in drilling using an infinite dimensional model[END_REF][START_REF] Germay | Multiple mode analysis of the self-excited vibrations of rotary drilling systems[END_REF], the main difference being that we only consider here the torsional dynamics. That is, for the angular motion, we denote the angular velocity and torque as ω(t, x), τ (t, x) which are functions of (t, x) evolving in {(t, x)| 0 < t < T, x ∈ [0, L]} (where L is the total length of the drill string and T a positive time), see Fig. 1. The torque is found from shear strain, given as twist per unit length. Letting φ denote the angular displacement (see Fig. 2) in the string s.t. ∂φ(t,x) ∂t = ω(t, x), we have τ (t, x) = JG φ(t, x) -φ(t, x+dx) /dx, where J is the polar moment of inertia and G is the shear modulus. Hence the equations for the angular motion are given by

∂τ (t, x) ∂t + JG ∂ω(t, x) ∂x = 0, (1)
Jρ ∂ω(t, x) ∂t + ∂τ (t, x) ∂x = S(t, x), (2)
where the source term is due to frictional contact with the borehole and is modeled as

S(t, x) = -k t ρJω(t, x) -F(ω, τ, x), (3)
where k t is a constant damping term representing the viscous shear stresses between the pipe and drilling mud, and F(ω, τ, x) is a differential inclusion, that still has to be described, representing the Coulomb friction between the drill string and the borehole.

J p J c τ + ,ω + τ -,ω - L p L c

Discontinuities of a multiple sectioned drill string

The lowermost section of the drill string is typically made up of drill collars which may have a great impact on the drill string dynamic due to their added inertia. In particular, the transition from the pipes to collars in the drill string will cause reflections in the travelling waves due to the change in characteristic line impedance [START_REF] Aarsnes | Linear stability analysis of self-excited vibrations in drilling using an infinite dimensional model[END_REF].

We split the drill string into a pipe section with polar moment of inertia and lengths J p , L p and a collar section with the same parameters given as J c , L c . We use τ + , ω + to denote the strain and velocity at the top of the drill collar and τ -, ω -at the bottom of the pipe, see Fig. 3. The boundary conditions at the transition are given by the following continuity constraints

ω + = ω -, τ + = τ -. (4)

Coulomb friction as an inclusion

The Coulomb friction is modeled as an inclusion

     F(ω, x) = µ k F N (x), ω > ω c , F(ω, x) ∈ [-µ s F N (x), µ s F N (x)], |ω| < ω c , F(ω, x) = -µ k F N (x), ω < -ω c , (5)
where F N is the normal force acting between the drill string and the borehole wall, ω c is the threshold on the angular velocity where the Coulomb friction switches from static to kinetic, µ s , µ k are the static and kinetic friction factors, respectively. Here, we distinguish between the kinetic and static friction factors: Static friction is friction between two or more solid objects that are not moving relative to each other. Kinetic friction, also known as dynamic friction or sliding friction, occurs when two objects are moving relative to each other and rub together. The function [START_REF] Aarsnes | Dynamics of a distributed drill string system: Characteristic parameters and stability maps[END_REF] and take the boundary values ±µ s F N (x) if this relation does not hold. The normal force profile of the drill string should be obtained from an appropriate torque and drag model, see e.g. [START_REF] Menand | AADE-13-FTCE-21 Borehole Tortuosity Effect on Maximum Horizontal Drilling Length Based on Advanced Buckling Modeling[END_REF][START_REF] Weijermans | Drilling with Rotary Steerable System Reduces Wellbore Tortuosity[END_REF][START_REF] Aarsnes | Torsional vibrations with bit off bottom: Modeling, characterization and field data validation[END_REF].

F(ω, τ, x) ∈ [-µ s F N (x), µ s F N (x)] denotes the in- clusion where
F(ω, τ, x) = - ∂τ (t, x) ∂x -k t ρJω(t, x) ∈ [-µ s F N (x), µ s F N (x)],

Boundary condition

At the topside boundary, the top drive is actuated by a motor delivering a torque τ m which we assume to be the control input. The topdrive has the inertia I T D and hence satisfies the dynamics

∂ω 0 ∂t = 1 I T D τ m -τ (t, x=0) , (7)
and finally, the angular velocity at the top of the drill string is equal the top drive velocity ω(t, x=0) = ω 0 .

Derivation of Riemann invariants

The Riemann invariants of a Hyperbolic PDE are the states corresponding to a transformation of the system which has a diagonalized transport matrix, i.e. the system can be written as a series of transport equations only coupled through the source terms [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]. On the set {(t, x)| 0 < t < T, x ∈ [0, L]}, we define the Riemann invariants as

α = ω + c t JG τ, β = ω - c t JG τ, (8)
where c t = G ρ is the velocity of the torsional wave. This transformation enables us to rewrite (1)-( 2) in variables α, β as the diagonalized PDE system

∂α ∂t (t, x) + c t ∂α ∂x (t, x) = S(t, x), (9)
∂β ∂t (t, x) -c t ∂β ∂x (t, x) = S(t, x). ( 10
)
where the source term S is defined by

S(t, x) = S(t, x) Jρ = k t (α(t, x) + β(t, x)) + 1 Jρ F(t, x), (11)
where F is the inclusion given in (5) (expressed as a function of α and β). In the Riemann coordinates, the boundary conditions (4) rewrite as follows

β + = 1 1 + Z α + (1 -Z) + 2 Zβ -, (12)
α -= 1 1 + Z 2α + + ( Z -1)β -, (13)
where we have denoted the relative magnitude of the impedance as

Z = c t JG + c t JG - . (14)
The boundary condition [START_REF] Auriol | Two-sided boundary stabilization of heterodirectional linear coupled hyperbolic pdes[END_REF] rewrites

∂ω 0 ∂t = 1 I T D τ m + GJ c t (β(t, 0) -ω 0 (t)) . (15)
In the case of the same material being used at both sides of the discontinuity, the only change is in the polar moment of inertia. That is, for a pipe-collar sections of e.g. steel, we have, following Fig.

3 Z = J c J p . ( 16
)
Note the meaning of ( 12)-( 13) as reflections of incoming waves from both sides, as they are split into an upward and a downward travelling wave.

Model validity

The effectiveness of this modeling approach is explored in some detail in [START_REF] Aarsnes | Torsional vibrations with bit off bottom: Modeling, characterization and field data validation[END_REF]. Here we illustrate this by briefly considering the open-loop fit of the model to full scale field data shown in Fig 4 and5. Details of the numerical implementation is given in Appendix A.

In both the cases of Fig 4 and5, the model accurately replicates the stick slip oscillations of the field data. The cases covers two qualitatively different behaviour of the stick slip, demonstrating a certain degree of versatility of the model. In particular we note that, for Fig. 4 where down-hole data is available, a good replication of the angular BHA velocity is also achieved. Down-hole data for the well in Fig. 5 is not available.

These fits were achieved by tuning the friction factors of the model to get the simulation results to match the data. The results in [START_REF] Aarsnes | Torsional vibrations with bit off bottom: Modeling, characterization and field data validation[END_REF] verify that the model is capable to replicate the full scale field dynamics in most cases when this tuning is good. The goal of the present paper is to obtain an automated algorithm, the adaptive observer (a.k.a. soft sensor), which performs this tuning online. The design and testing of this adaptive observer is done in the following.

Semi-lumped approximation

In this section we derive a model for observer design by using a lumped approximation of the drill collar section following the approach of [START_REF] Aarsnes | Avoiding stick slip vibrations in drilling through startup trajectory design[END_REF]. This amenable model approximation will be used for the observer design.

Lumped BHA

The approximation entails lumping the effect of the source term [START_REF] Detournay | A phenomenological model for the drilling action of drag bits[END_REF] into the lumped dynamics of the BHA. This is a reasonable approximation for many drill-strings as much of the torque acting on the drill string will come from stabilizers located in, or close to, the BHA [START_REF] Aarsnes | Torsional vibrations with bit off bottom: Modeling, characterization and field data validation[END_REF]. The inertia of the lumped BHA is

I BHA = ρL c J c , (17)
and hence the angular velocity of the BHA is governed by where d(t) accounts for the now lumped effect of the distributed source term, i.e.:

∂ω Lp ∂t = 1 I BHA τ (t, x=L p ) -d(t) , (18)
d(t) ≈ L 0 S(ω, x). (19)
Here, ( 19) is meant for illustration and is not used directly. When we later employ the flat formulation, facilitated by this approximation, for estimation d(t) will be treated as an uncertain disturbance.

Using this lumped approximation of BHA, we obtain what we will refer to as the semi-lumped formulation, given by the distributed wave-equation for the pipe section

∂τ p (t, x) ∂t + JG ∂ω p (t, x) ∂x = 0, (20)
J p ρ ∂ω p (t, x) ∂t + ∂τ p (t, x) ∂x = 0, ( 21
)
defined on {(t, x)| 0 < t < T, x ∈ [0, L p ]} with the boundary conditions ω p (t, x=0) = ω 0 , ω p (t, x=L p ) = ω Lp governed by ∂ω 0 ∂t = 1 I T D τ m -τ p (t, x=0) , ( 22
)
∂ω Lp ∂t = 1 I BHA τ p (t, x=L p ) -d(t) . (23)
The fit of this approximation and how to quantify any resulting error is discussed in [START_REF] Aarsnes | Avoiding stick slip vibrations in drilling through startup trajectory design[END_REF].

In Riemann invariants

Using the Riemann Variables (α, β) as states, the semilumped system (20)-(23) defined on {(t, x)| 0 < t < T, x ∈ [0, L p ]}, rewrites:

∂α p ∂t (t, x) + c t ∂α p ∂x (t, x) = 0, ( 24
)
∂β p ∂t (t, x) -c t ∂β p ∂x (t, x) = 0. ( 25
)
with the boundary conditions

β p (t, L p ) = 2ω Lp (t) -α p (t, L p ), (26)
α p (t, 0) = 2ω 0 (t) -β p (t, 0), ( 27
) ∂ω 0 ∂t (t) = a 0 (-ω 0 (t) + β p (t, 0)) + 1 I T D τ m (t) (28) ∂ω Lp ∂t (t) = a Lp (-ω Lp (t) + α p (t, L p )) - 1 I BHA d, (29)
where the frequency constants given as a 0 =

α p (t, x=L p ) = α p (t -t D , x=0) (30)
β p (t, x=0) = β p (t -t D , x=L p ), (31)
where t D = L p /c t . Thus, we note that this system is characterized by three time constants 1/a 0 : Top drive time constant.

1/a Lp : BHA time constant.

t D : Drill string travel time.

Soft sensor (Observer) design

The core of the approach is a 'soft-sensor' combining measurements from physical sensors with a model of the system dynamics to provide estimates of states and sideforces. The soft-sensor is based on an observer algorithm. Observers are dynamical systems whose states are the estimates of the aforementioned variables. The observer state dynamics are a combination of • 'open-loop' terms based on the equations of physics describing the behaviour of the system.

• 'closed-loop' terms that correct the dynamics based on the value of the measurements. Typically, the value from the sensors are compared with their estimates from the observer, and a function of the discrepancy between the two are added as a correction term.

The task of design an observer is to engineer the combination of 'open-loop' and closed-loop correction terms.

Often, the correction terms are linear function of the estimation error of the measured states. These correction terms are then used to update the estimated states as new measurement data become available. In our case, they appear in the observer dynamics as spatially-varying and lumped gains called observer gains. We will use the specific observer gains derived in [START_REF] Di Meglio | Robust output feedback stabilization of an ODE-PDE-ODE interconnection[END_REF] using a backstepping method which guarantees fast and robust convergence of the estimates assuming the model is correct 1 . The soft sensor is adaptive in that it estimates the side forces and then updates the observer model kinetic or static friction factor itteratively for each time step depending on if the drill string BHA is stuck or sliding. The algorithm presented here is not able to simultaneously update the transition velocity ω c in (5), hence this value has to be fixed. We remark that it was shown in [START_REF] Aarsnes | Torsional vibrations with bit off bottom: Modeling, characterization and field data validation[END_REF] that this transition velocity can have a significant impact on model bahvior in some cases, which might limit the performance of the proposed algorithm.

Observer system equations

The observer equations are given as a copy of the plant equations plus the correction terms. The measured output of the system corresponds to the top drive angular velocity ω 0 .

We denote an estimated variable with the • superscript. We define the measured estimation error as: e = ω0 -ω 0 . In the following p α (•), p β (•), p 0 , p 1 , P 0 , P 1 , l s and l k denote the observer gains, which still have to be defined. The estimation of the top drive angular velocity is given by ω0 = a 0 βp (t, 0) -ω0 + 1

I T D τ m -p 0 e, (32)
For the pipe section, the estimation of the Riemann invariant is given by

∂ αp ∂t (t, x) + c t ∂ αp ∂x (t, x) = Ŝp (t, x) -p α (x)e, (33)
∂ βp ∂t (t, x) -c t ∂ βp ∂x (t, x) = Ŝp (t, x) -p β (x)e. ( 34
)
For the collar section, the estimation of the Riemann invariant is given by

∂ αc ∂t (t, x) + c t ∂ αc ∂x (t, x) = Ŝc (t, x) -p 1 e, (35)
∂ βc ∂t (t, x) -c t ∂ βc ∂x (t, x) = Ŝc (t, x) -p 1 e. ( 36
)
Finally, the boundary conditions are

αp (t, 0) = 2ω 0 (t) -βp (t, 0) -P 0 e, (37) βp
(t, L p ) = αp (t, L p )(1-Z) + 2 Z βc (t, L p ) 1 + Z -P 1 e (38) αc (t, L p ) = 2α p (t, L p ) -βc (t, L p )(1-Z) 1 + Z (39) βc (t, L) = -α c (t, L) (40)
The source term in each section are computed from the estimated states and friction factors Ŝi (t, x) = Ŝ(t, x)

J i ρ = k t (α i (t, x) + βi (t, x)) + 1 J i ρ F(t, x), (41)
where F(t, x) is the inclusion given in ( 5), and the estimates of the friction factor is updated according to

μs (t) = -l s e, |ω Lp | ≤ ω c , 0, |ω Lp | > ω c , (42)
μk (t) = 0, |ω Lp | ≤ ω c , l k e, |ω Lp | > ω c , (43)
Finally, the following saturation is used to improve robustness of the method:

μs = max(μ s , μk ). ( 44
)
The initial condition of ( 32)-( 40) can be arbitrarily chosen.

Semi-lumped approximation and error system

We want to design the observer gains using the approach proposed in [START_REF] Di Meglio | Robust output feedback stabilization of an ODE-PDE-ODE interconnection[END_REF] to ensure the convergence of the observer state (solution of (32)-( 40)) to the real state (solution of ( 9)-( 15)). To do so, we need to rewrite the observer system in a suitable form (i.e. without the inclusion). This is done using the lumped approximation of the drill collar section introduced in Section 3.1. More precisely, the observer dynamics (32)-( 40) can be rewritten

ω0 (t) = a 0 βp (t, 0) -ω0 (t) + 1 I T D τ m -p 0 e, (45)
∂ αp ∂t (t, x) + c t ∂ αp ∂x (t, x) = -p α (x)e, ( 46
)
∂ βp ∂t (t, x) -c t ∂ βp ∂x (t, x) = -p β (x)e, ( 47
)
ωLp (t) = a Lp αp (t, L p ) -ωLp (t) - 1 I BHA d -p 1 e, ( 48
)
with the boundary conditions αp (t, 0) = 2ω 0 (t) -βp (t, 0)

-P 0 e, ( 49
) βp (t, L p ) = 2ω Lp (t) -αp (t, L p ) -P 1 e, (50)
Similarly to what has been done in Section 3.1, the term d accounts for the now lumped effect of the distributed source term (obviously depending on the expression given in (41)). For the design of observer gains, we consider the term d as a constant lumped disturbance, normalized by

I BHA L 0 F N (x)dx
(thus, d is also assumed to be constant). Subtracting system (24)-(29) from system (45)-( 50) and denoting the error variables with the • superscript (i.e αp = αp -α p for instance), we get the following error system ω0 (t) = a 0 βp (t, 0) -e -p 0 e, (51)

∂ αp ∂t (t, x) + c t ∂ αp ∂x (t, x) = -p α (x)e, (52)
∂ βp ∂t (t, x) -c t ∂ βp ∂x (t, x) = -p β (x)e, ( 53
) ωLp (t) = a Lp αp (t, L p ) -ωLp (t) -p 1 e - 1 I BHA d, ( 54
) ḋ = 0, ( 55
)
with the boundary conditions αp (t, 0) = 2ω 0 (t) -βp (t, 0) -P 0 e, (56) βp (t, L p ) = 2ω Lp (t) -αp (t, L p ) -P 1 e, (57)

Control dual problem

The backstepping approach proposed in [START_REF] Di Meglio | Robust output feedback stabilization of an ODE-PDE-ODE interconnection[END_REF] provides an explicit method to design a robust output feedback boundary control law for an ODE-PDE-ODE interconnection. However, it has been proved in [START_REF] Auriol | Two-sided boundary stabilization of heterodirectional linear coupled hyperbolic pdes[END_REF] that the gains of such a control law correspond to the gain of the observer of the dual problem (this has actually been proved in the case of a system of n + m PDEs, but the proof can be easily adjusted to an ODE-PDE-ODE interconnection). More precisely, adjusting the methods proposed in [START_REF] Auriol | Two-sided boundary stabilization of heterodirectional linear coupled hyperbolic pdes[END_REF] we can prove that the system (51)-(57) has the same stability properties as those of the system defined on

{(t, x)| 0 < t < T, x ∈ [0, 1]} by ∂ψ ∂t (t, x) - c t L p ∂ψ ∂x (t, x) = 0, (58)
∂φ ∂t (t, x) + c t L p ∂φ ∂x (t, x) = 0, (59) χ0
(t) = 2 c t L p ψ(t, 0) -a 0 χ 0 -p 0 χ 0 (t) -p 1 χ 1 (t) -P 0 c t ψ(t, 0) -P 1 c t φ(t, L p ) - Lp 0 p α (ξ)ψ(t, ξ) + p β (ξ)φ(t, ξ)dξ, (60)
χ1 (t) = 2 c t L p φ(t, 1) -a Lp χ 1 (t), (61) χ2
(t) = - 1 I BHA χ 1 (t), (62)
with the boundary conditions

ψ(t, 1) = -φ(t, 1) + L p c t a Lp χ 1 (t), (63) φ(t
, 0) = -φ(t, 0) + L p c t a 0 χ 0 (t) (64)
This new system (which is the adjoint of ( 51)-( 57)) has exactly the structure which is considered in [START_REF] Di Meglio | Robust output feedback stabilization of an ODE-PDE-ODE interconnection[END_REF]. This immediately gives us the expression of the observer gains (we choose to not rewrite them here for sake of clarity). More precisely, these observer gains only depends on three tuning parameters (namely, κ 1 , z and p) that can be tuned to shape the observer top-side reflection and place the poles of the down-hole "error system". Considering the adaptive update law, the gain l k , is derived from the observer assuming a constant lumped disturbance, and then normalized by

I BHA L 0 F N (x)dx .
The gain for the update law of the static friction factor is then chosen to be two times this l s = 2l k . With these gains we can ensure the convergence of the observer states to the real states in presence of a constant disturbance. One must be aware that this assumption is only an approximation and does not hold in the real case. However, due to the robustness properties of the observer (see [START_REF] Di Meglio | Robust output feedback stabilization of an ODE-PDE-ODE interconnection[END_REF] for details) we can still have the convergence of the different states to their real values when considering the real inclusion term. Considering the estimations of the friction terms, we do not have any guarantee of convergence, but, due to the robustness properties of the observer, they should converge close to the real friction terms. This is experimentally validated in the next section.

Field Application

The observer detailed in the previous section may be used to provide online parameter estimation for the friction parameters of the drillstring model presented in Section 2 but can also be used as a method to estimate BHA rotation and torque. This can be of particular usefulness in directional drilling scenarios where real-time estimation of tool face angle -BHA angular orientation -is essential, and in feedforward stick-slip mitigation systems.

The envisioned industrial implementation of the estimatino algorithm is briefly described in this section. The flow of data is show in Figure 6 and a flowchart describing the process is presented in Figure 7.

Output of the estimation algorithm may include an estimate of the drillstring state as a function of measure depth, and a time series of the friction coefficient, estimated drillstring twist and the soft sensor (observer) gains.

We now consider a test of the estimation algorithm (32)-(43). We want to evaluate the real time estimation of the friction coefficients and of the BHA and top drive rotation in two different situations: [START_REF] Aarsnes | Torsional vibrations with bit off bottom: Modeling, characterization and field data validation[END_REF] comparison against a simulation model (2) comparison against field data. In each situation, we run the observer using high frequency measurements (100 Hz) but also downsampled 5 Hz and 1 Hz measurements. These examples are of particular note since a majority of supervisory control systems currently deployed on drilling rigs in the field operate at 5 to 10 Hz.

Test against simulation model

We test our observer against the simulation model described in [START_REF] Aarsnes | Axial and torsional self-excited vibrations of a distributed drill-string[END_REF] with the wellbore survey shown in Figure 8, using the numerical implementation described in Appendix A. The kinetic friction is chosen to be equal to 0.187, while the static friction is chosen to be equal to 0.6 which is similar to values reported using traditional friction tests in the field. The well represents a simple build

Step 6: Output friction coefficient and drillstring state

Step 5: Compute estimated drillstring state including angular displacement and velocity

Step 4: Compute estimate of static and dynamic coefficients of friction

Step 3: Compute or update observer gains (if applicable) and propagate model estimate and update based on data

Step 2: Read RPM and Torque from sensors on drilling rig

Step 1: Obtain well plan (wellbore survey) and drillstring configuration and hold well used throughout the world. More discussion of this synthetic example may be found in [START_REF] Aarsnes | Torsional vibrations with bit off bottom: Modeling, characterization and field data validation[END_REF].

We have pictured in Fig. 9 the real value of the BHA angular velocity and of the top drive angular velocity as well as the estimations obtained using the observer in three different situations. In the first plot, the measurements are available every 0.01 s (100 Hz), while in the two next plots they are only available every 0.2s (5 Hz), resp. 1s (1 Hz)). This means that for the cases with lower sampling rates, the observer correction term is less frequently updated, which we would expect to potentially reduce convergence rate or accuracy. As is seen in Fig. 9, the 100 and 5 Hz cases show a rapid convergence of the estimated states to the correct trajectories, while at 1 Hz sampling rate the performance has started to degrade.

We have pictured in Fig. 10, the estimation given by the observer of the friction parameters for different starting points in the three different situations (100 Hz, 5 Hz and 1 Hz measurements). In the two first cases, we have a rapid convergence of the estimations regardless of the initial guess. However, the convergence value (in particular for the static friction term) may not exactly correspond to the exact one (even if it remains close). In the last case (1 Hz measurements), the estimations of the friction terms oscillate around the exact value. The estimations remain good however.

Test on field data

We now test our observer against field data obtained from an onshore well with the wellpath shown in Fig. 4. Surface rpm and torque data was recorded at 100Hz using a high frequency data recorder and downhole data in the BHA was recorded on a memory tool as a 10 second windowed minimum, maximum and average value. RPM data was obtained using an encoder at surface and via magnetometer/accelerometer data downhole. Surface torque data was obtained from current integration in the inverter controller the topdrive and pipe torque was inferred by removing inertial acceration torque. Top drive inertia was obtained using a chirp test. We have pictured in Fig. 11 the real mean value of the BHA angular velocity and of the top drive angular velocity as well as the estimations obtained using the observer, again for the three cases of 100 Hz, 5 Hz and 1 Hz data sampling rates. Once again, the observer is run at a frequency of 100 Hz, possibly using the same value of the output for successive iterations in the two latter cases. In the three cases, the estimations provided by the observer tends to convergence to the real values. Note that the field data can be unsynchronized and that we don't take here the possible offset into account. Then we cannot really compare the phases of the field data/estimated angular velocities but only their respective magnitudes.

We have pictured in Fig. 12 the estimation given by the observer of the friction parameters for different starting points in the three different situations (100 Hz, 5 Hz and 1 Hz measurements). In the three cases the estimations converge to the same values (0.6 for the static friction and 0.35 for the kinetic friction) regardless of the initial guess.

Conclusion

We have proposed a model and algorithm to estimate both drill bit velocity and friction factors while drilling, using only surface measurements. The approach relies on an observer designed for the linearized system, combined with an update law of the static and kinemtic friction factors used in a nonlinear coulomb-friction model. The estimation is robust to poor initial estimates of the states and parameter, as have been illustrated on both synthetic and field data. The adaptive estimation of the friction factors enables the states estimates of the observer to converge (close to) actual values in real time. This means that we are able to continously monitor the downhole behavior of the drill string, specifically the dynamic angular velocity and torques, even when drilling through stick slip. Furthermore, we envision that these state and parameter estimates could enable measures to counteract unwanted drill string behaviour, e.g. using a techniques such as proposed in [START_REF] Aarsnes | Avoiding stick slip vibrations in drilling through startup trajectory design[END_REF], and represent a step towards automated closed loop geosteering. and at cell # j and time step # k

Figure 1 :

 1 Figure 1: Schematic indicating the distributed drill string lying in deviate borehole.

Figure 2 :

 2 Figure 2: Infinitesmal drill string element.

Figure 3 :

 3 Figure 3: Collar-Pipe transition.

Figure 4 :

 4 Figure 4: Recorded and simulated a) RPM and b torque of a drill-string at a bit depth of 1,733 m in a well with the survey shown in c), using the friction parameters: µs = 0.34, µ k = 0.187, ωc = 19 (RPM).

Figure 5 :

 5 Figure 5: Recorded and simulated a) RPM and b) torque of a drill-string at a bit depth of 2,506 m in a well with the survey shown in c), using the friction parameters: µS = 0.43, µ k = 0.185, ωc = 17 (RPM).

  GpJp ctI T D and a Lp = GpJp ctI BHA . They are expressed in seconds -1 and they repectively represent the inertia of the top-drive and BHA, relative to the line impedance of the pipe section of the drill string ζ p = GpJp ct . The solution to (24)-(25) can be written as the delay equations:

Figure 6 :

 6 Figure 6: Flow of data from the surface sensors, into the estimation algorithm and the resulting estimated output.

Figure 7 :Figure 8 :

 78 Figure 7: Flow chart describing the flow of information and steps used in the potential application.

Figure 9 :

 9 Figure 9: Time evolution of the real and estimated angular velocity in the test against simulation model for, a) 100 Hz, b) 5 Hz and c) 1 Hz measurements.

1 )

 1 -α j-1 )+ (A.1)β j + c t ∆t ∆x (β j+1 -β j ) + ∆tk t (α k j + β k j ) , (A.2) j ) min(| Fk j |, µ s F N ), j ) min(| Fk j |, µ k F N ), -∆tk t (α k j + β k j ) j -β k j+1 ) -∆tk t (α k j + β k j ) -F k j .

Figure 10 :

 10 Figure 10: Estimation of static and kinetic friction factors for different initial conditions in the test against simulation model for, a) 100 Hz, b) 5 Hz and c) 1 Hz measurements.

Figure 11 :

 11 Figure 11: Real and estimated angular velocity of the BHA and the top drive in the test against field data for, respectively, a) 100 Hz, b) 5 Hz and c) 1 Hz measurements.

Figure 12 :

 12 Figure 12: Estimation of static and kinetic friction factors for different initial conditions in the test against field data for, a) 100 Hz, b) 5 Hz and c) 1 Hz measurements.

Table 1 :

 1 Nomenclature (•), p β (•), p 0 , p 1 , P 0 , P 1 , l s and l k are observer gains.

	Parameters
	a 0 , a Lp top drive, BHA frequency coefficient
	c t	torsional wave velocity
	F	Coulomb component of side force
	G	shear modulus
	I T D	top drive inertia
	I BHA	BHA inertia
	J p , J c	pipe, collar polar moment of inertia
	k t	viscous component of side force
	L	drill string length
	L p , L c pipe, collar drill string length Z collar-pipe relative impedance
	µ k , µ s kinetic, static friction coefficient
	ρ	drill string density
	p α Dependent variables
	d	lumped side force
	F N normal force on drill string
	S	torque source term (side force)
	S	Normalized source term
	α, β drill string Riemann invariants
	τ	drill string torque
	τ m	motor torque
	φ	drill string orientation
	ω	drill string angular velocity
	ω 0	Top drive angular velocity
	Theˆsuperscript denotes estimasted quantities.
	The˜superscript denotes estimation residual.
	Independent variables
	t time
	x position relative to top of drill string

More precisely, a proof of asymptotic convergence of the estimates to the true states is provided in[START_REF] Di Meglio | Robust output feedback stabilization of an ODE-PDE-ODE interconnection[END_REF] for the linearised dynamics, i.e. close to the constant homogeneous velocity profile.

Acknowledgments

This work was in part supported by the Research Council of Norway, ConocoPhillips, AkerBP, Statoil and Wintershall through the research center DrillWell (203525/O30) at IRIS, and by the Marie Sk lodowska-Curie (MSCA) -Individual Fellowships through the Robust Estimation and Control of Infinite Dimensional Systems (RECIDS) Project. Additional support came from the University of Calgary's Canada First Research Excellence Fund Program, the Global Research Initiative in Sustainable Low Carbon Unconventional Resources.

Appendix A. Numerical implementation

In the numerical implementation of the model the wave equation ( 1),( 2) is transformed into transport equations discretized using a first order upwind scheme. This choice is made to ensure numerical robustness and to avoid spurious oscillations, as higher order schemes perform poorly due to the temporal discontinuities introduced by the distributed differential inclusions which are used to represent the Coulomb friction. Numerical accuracy can then be ensured by having a sufficiently fine spatial grid, and this is an amenable approach due to the fact that simulation speed is not of critical importance for the present study.

In all simulations a spatial grid of 500 cells is used for the drill string and the time-step is chosen to enforce the Courant-Friedrichs-Lewy (CFL) condition [START_REF] Courant | On the partial difference equation of mathematical physics[END_REF].

In the numerical treatment of the model, F is implemented as follows. For cell size ∆x and and time step ∆t,


Jean Auriol

Nasser Kazemi

Roman J Shor

Kristopher A Innanen

Ian D Gates

A sensing and computational framework for estimating the seismic velocities of rocks interacting with the drill-bit

We have developed a sensing and computational framework to estimate seismic velocities of rocks interacting with the drill-bit during the drilling process. The performance of drilling depends on our knowledge of the subsurface. The interaction between the drill-bit and rock can introduce severe vibrations in the drill-string and result in safety and performance issues. However, we can use seismic waves radiated from drill-bit-rock interactions to determine seismic velocities of the rocks interacting with the drill-bit. Our approach consists of a distributed (wave equation) representation of the dynamics of the drill-string for which we show (using Riemann's invariants and a backstepping approach) that it is possible to express the force-on-bit as a function of the top-drive force and the topdrive velocity, without requiring explicit information about the subsurface properties. We also show that seismic waves generated by drill-bit-rock interaction can be modelled as functions of the force-on-bit and of rock velocities. The rock velocity independent formulation of the force-on-bit, along with modelling of the seismic waves generated by drill-bit-rock interaction as a function of force-on-bit and rock velocities allow us to estimate seismic velocities of rocks interacting with the drill-bit. We use the alternating minimization algorithm to estimate the velocities. Numerical examples on simulated data are indicators of the validity of the approach. The proposed methodology is the first step towards a subsurface-aware drilling system.

I. INTRODUCTION

Optimizing the drilling process requires that subsurface material properties be accurately known. Drill-bit-rock interaction can cause severe vibrations in the drilling system and result in an inefficient Rate Of Penetration (ROP). It may also raise safety issues. Control mechanisms may be used to reduce these vibrations in the system [START_REF] Aarsnes | Avoiding stick slip vibrations in drilling through startup trajectory design[END_REF], [START_REF] Aarsnes | Benchmarking of industrial stick-slip mitigation controllers[END_REF]. However, such control mechanisms often depend on the rock properties that the system is drilling into. Hence, having access to the accurate subsurface information in near real-time will allow us to design an efficient drilling system. Here, we only focus on estimating the seismic velocities of rocks interacting with the drill-bit in the drilling process. To do so, the first step is to understand the drilling mechanism.

The drilling of an oil well consists of creating a borehole up to several thousand meters deep into the ground until an oil reservoir is reached. The mechanical part of a drilling system is made of three components: the rotating mechanism Corresponding author: J. Auriol (email: jean.auriol@ucalgary.ca).

(usually a rotary table or a top drive), the drill-string and the Bottom Hole Assembly (BHA). The dynamical behaviour of drill-strings is complex and involves many dynamic phenomena (vibrations, bending and twisting quasi-static motion, bit-rock interactions [START_REF] Kapitaniak | Unveiling complexity of drill-string vibrations: Experiments and modelling[END_REF], [START_REF] Spanos | Oil and gas well drilling: a vibrations perspective[END_REF]). To describe the dynamics of interest, distributed parameter models have been proposed in the literature [START_REF] Meglio | A distributed parameter systems view of control problems in drilling[END_REF], [START_REF] Saldivar | A control oriented guided tour in oilwell drilling vibration modeling[END_REF]. In this class of models, the axial and torsional efforts can be modeled by a set of hyperbolic Partial Differential Equations (PDEs). For this class of representation, special attention must be paid to the boundary conditions and in particular to the bit-rock interaction law at the downhole boundary [START_REF] Boussaada | Analysis of drilling vibrations: A time-delay system approach[END_REF]- [START_REF] Richard | Self-excited stickslip oscillations of drill bits[END_REF] that depends on the drilled rock properties. The drill-string interaction with the borehole gives rise to a wide variety of non-desired oscillations [START_REF] Dunayevsky | Application of Stability Approach to Bit Dynamics[END_REF]- [START_REF] Saldivar | Stick-slip oscillations in oillwell drilstrings: Distributed parameter and neutral type retarded model approaches[END_REF] which have been extensively studied. Among them, torsional oscillations known as stick-slip are considered to be one of the most prevalent vibrations. These stick-slip oscillations are characterized by a series of stopping -"sticking" -and releasing -"slipping" -events of the bit. These vibrations may lead to a reduction of the ROP, cause fatigue on the equipment and wellbore instability. They can also lead to premature failure of the bit and may cause catastrophic failures, or, at least, wear to expensive components of the drill-string [START_REF] Kriesels | Cost savings through an integrated approach to drillstring vibration control[END_REF]. Such oscillations can be the consequence of specific downhole conditions (rock composition or small borehole diameter). Numerous models assume that stick-slip is a consequence of the non-linear frictional force actuating at the bit by contact with the rock [START_REF] Leine | Stick-slip Whirl Interaction in Drillstring Dynamics[END_REF], [START_REF] Nandakumar | Stability analysis of a state dependent delayed, coupled two DOF model of drill-string vibration[END_REF]. More precisely, in these models, the stick-slip phenomenon is assumed to be related to the velocity-weakening of the frictional force at the bit (Stribecklike effect) associated with typical dry friction profiles (static friction and dynamic friction) [START_REF] Kapitaniak | Unveiling complexity of drill-string vibrations: Experiments and modelling[END_REF], [START_REF] Brett | The Genesis of Bit-Induced Torsional Drillstring Vibrations[END_REF]. However, the bit-rock interaction is not the only cause of stick-slip; for instance, it would not explain the occurrence of such oscillations in the case of bit off bottom [START_REF] Brett | Uses and Limitations of Drillstring Tension and Torque Models for Monitoring Hole Conditions[END_REF]- [START_REF] Zhao | Abnormal Down Hole Pressure Variation by Axial Stick-Slip of Drillstring[END_REF].

A clear understanding of drill-string dynamics appears to be crucial to control these vibrations and to obtain an improvement of the performance of drilling systems (ROP), preventing damage and reducing safety risks. Moreover, as the cost of drilling a well is in large part determined by the time it takes, strategies for reducing non-productive time and increasing ROP are sought. Different models of friction at the contact point between the cutting device and the drilling surface have been proposed in the literature [START_REF] Saldivar | A control oriented guided tour in oilwell drilling vibration modeling[END_REF]. These models usually depend on the nature of the drilled rock (dependence of the drill-string dynamic response on bit-rock interaction and rock properties were verified in [START_REF] Shor | For better or worse: applications of the transfer matrix approach for analyzing axial and torsional vibration[END_REF]). They can later be used to express the downhole boundary condition of the distributed parameter models previously mentioned, such a boundary condition being for instance found by expressing the weight-on-bit, which are the result of cutting forces and friction forces. Estimating the seismic velocities of rocks while drilling is not an easy task as downhole sensors are expensive and may raise potential technical risks. A viable option is to remotely record the seismic energy generated by the drillbit rock interactions and use the recorded measurements to estimate the seismic velocities of the rocks interacting with the drill-bit.

Since the 1960s, recording and interpreting the seismic energy generated by the drill-bit-rock interaction has been an active area of research [START_REF] Meehan | Drill bit seismic: A drilling optimization tool[END_REF]- [START_REF] Anchliya | A review of seismic while drilling (swd) techniques: A journey from 1986 to 2005[END_REF]. This methodology is called seismic while drilling (SWD). The SWD data are rich in information, and, provided that we understand the nonimpulsive and correlative source signature of the drill-bit rock interaction, they can enhance our knowledge of the subsurface materials. In recent years, due to the advances in modelling the radiation patterns of the drill-bit rock interaction and increased computational power in the field, new methodologies have been developed [START_REF] Meehan | Rekindling interest in seismic while drilling[END_REF]- [START_REF] Bharadwaj | Focused blind deconvolution[END_REF]. Here, we make use of such data to estimate the seismic velocities of rocks interacting with the drill-bit.

In this paper, we present a new algorithm for estimating the seismic velocities of rocks interacting with the drillbit. The proposed algorithm uses far-field radial and angular displacements (obtained by using seismic sensors) combined with the measurements of the drill-string top-drive force and top-drive velocity. Our methodology only requires knowledge of the physical properties of the drilling system such as its Young's modulus, pipe density, etc, but does not require any knowledge about the subsurface properties. Employing a distributed axial dynamics of the drill-string, we show that the force exerted on the bit can be expressed as a function of the top-drive force and top-drive velocity (which are measured) without depending on the downhole boundary condition or on the subsurface. Combining this relationship with a seismicwhile-drilling approach, we can express the far-field radial and angular displacement as a function of the force-on-bit and rock velocity, the only unknown parameters now being the seismic velocities of the rock. We implement alternating minimization approach to estimate the seismic velocities of the rocks interacting with the drill-bit in near real-time. This is the first step towards a subsurface-aware drilling system.

The paper is organized as follows. In Section II, we derive a distributed model to describe the axial dynamics of the drill-string. Seismic-while-drilling techniques are described in Section III. In Section IV, we derive the Riemann invariants associated with the distributed model of the drill-string. Then, using a backstepping approach, we express the drill-bit axial force as a function of the top-drive force and velocity. This relationship is used in Section V to develop a new algorithm that estimates the seismic velocities of the rocks ahead of In what follows, given a set Ω ⊆ R, its characteristic function will be denoted by

1 Ω (θ) = 1 if θ ∈ Ω 0 otherwise.

II. DRILL-STRING MODEL

In this section we present the model under consideration to describe the mechanical dynamics of the drilling system.

A. Distributed Axial dynamics of the drill-string

Let us consider a vertical well, as pictured in Figure 1. The drilling system consists of three parts: the rotary table, the drill-string and the drill-bit. We consider the case of a rollercone drill-bit.

The rotary table sets the drill-string into a rotary motion around its main axis and the dynamics of interest can be derived by assuming elastic deformations and using equations of continuity and state. More precisely, a distributed model has been proposed in [START_REF] Meglio | A distributed parameter systems view of control problems in drilling[END_REF], [START_REF] Germay | Nonlinear Drillstring Dynamics Analysis[END_REF] to describe the evolution of the axial displacement ξ(t, x) of the drill-string which is a function of (t, x) evolving in {(t, x) | 0 < t < T, x ∈ [0, L]} (where L is the total length of the drill-string and T a positive time). The axial motion satisfies the following wave Partial Differential Equation

∂ 2 ξ ∂t 2 (t, x) -c 2 ξ ∂ 2 ξ ∂x 2 (t, x) = -k ξ ∂ξ ∂t (t, x), (1)
where c ξ = E ρ , ρ being the pipe mass density E its Young's modulus and k ξ is a damping coefficient representing the Fig. 1. Schematic representation of a drilling system (modified from [START_REF] Leine | Stick-slip Whirl Interaction in Drillstring Dynamics[END_REF]) viscous shear stresses acting on the pipe. The axial force associated to ξ can be found from the strain, given as the local relative compression:

w(t, x) = AE (ξ(t, x) -ξ(t, x + dx)) dx , (2)
A being the cross-sectional area of the drill-string and dx → 0 the infinitesimal axial position increment. The axial velocity satisfies v(t, x) = ∂ξ(t, x) ∂t .

It can easily be shown from (1) that the axial force and velocity satisfy the following set of PDEs

∂w(t, x) ∂t + AE ∂v(t, x) ∂x = 0, (3)
∂v(t, x) ∂t + 1 Aρ ∂w(t, x) ∂x = -k ξ v(t, x). (4)
The topside weight on the drill-string, w(t, 0), corresponds to the system actuation. The downhole boundary condition at x = L is obtained from a force balance on the lumped Bottom-Hole Assembly (BHA). These two boundary conditions will be discussed in the following.

B. Topside boundary condition

In the considered drilling system, the drill-string is connected at the top to the top-drive suspended over the drill floor by the traveling block. This block is connected by several drill lines with one attached to the deadline anchor and the other being spooled on a drum controlled by AC induction motors [START_REF] Cayeux | On the Importance of Boundary Conditions for Real-Time Transient Drill-String Mechanical Estimations[END_REF]. Thus, we can assume that the operator control the weight on the drill-string w 0 (t) = w(t, 0). This yields

-EA ∂ξ(t, x) ∂x = w 0 (t), (5)
where w 0 is considered as an input.

C. Downhole boundary condition: bit-rock interaction

The lowermost section of the drill-string is typically made up of drill collars which may have a great impact on the drill-string dynamic due to their added inertia. In particular, the transition from pipes to collars in the drill-string will cause reflections in the traveling waves due to the change in characteristic line impedance [START_REF] Aarsnes | Linear stability analysis of self-excited vibrations in drilling using an infinite dimensional model[END_REF]. Rigorously speaking, we should write a new wave PDE for the Bottom-Hole Assembly (BHA). This PDE would be analogous to (1) but using the BHA density and Young's modulus. However, as the length of the BHA (≈ 200m) is much smaller than the one of the drill-string (≈ 2000m), its effect can be lumped into an ODE coupled with the drill-string [START_REF] Meglio | A distributed parameter systems view of control problems in drilling[END_REF]. Thus, the downhole boundary condition at x = L can be obtained from a force balance on the lumped BHA. This yields

M b ∂v ∂t (t, L) = -w b (v(t, L), w(t, L)) + w(t, L), (6)
where M b is the mass of the lumped BHA and w b (•, •) the force acting from the rock on the BHA through the drilling bit, known as the weight on bit. From [START_REF] Germay | Nonlinear Drillstring Dynamics Analysis[END_REF], [START_REF] Cayeux | On the Importance of Boundary Conditions for Real-Time Transient Drill-String Mechanical Estimations[END_REF], the weight on bit can be related to the bit velocity by considering the combined depth of cut per revolution [START_REF] Richard | A simplified model to explore the root cause of stick-slip vibrations in drilling systems with drag bits[END_REF]. More precisely, the combined depth of cut per revolution d(t) satisfy the following equation,

d(t) = v(t, L) ω bit ,
where the bit angular velocity ω bit is assumed constant here. The cutting force is finally expressed as

w b (v(t, L), w(t, L)) = w f + Kd(t) = w f + aζ d(t), (7)
where w f is a friction weight independent of the bit velocity (and therefore constant) while K = aζ , with a being the bit radius, ζ a characterization of the cutting angle and the intrinsic specific energy of the rock [START_REF] Richard | A simplified model to explore the root cause of stick-slip vibrations in drilling systems with drag bits[END_REF], [START_REF] Detournay | A phenomenological model for the drilling action of drag bits[END_REF]. Thus, the boundary condition ( 6) rewrites as

M b ∂ 2 ξ ∂t 2 (t, L) = - aζ ω bit ∂ξ(t, L) ∂t -w f -EA ∂ξ(t, L) ∂x . (8)
Due to the complexity of this boundary condition, controlling the weight on the drill-string to achieve optimal drilling performance is not an easy task. It has to be updated in realtime to adapt to the changing operating conditions (different type of rocks for instance). To do so, classical control procedures rely on topside drilling data and either assume known the nature of the drilled rock or use simple PID controllers. This explains why the knowledge of the nature of the subsurface could lead to an improvement in the performance of the control mechanism and thus optimize ROP while drilling.

III. SEISMIC WHILE DRILLING

Each layer in the subsurface is considered to be a homogenous and isotropic medium. In this model, each layer is recognized by its seismic velocities. The seismic velocities depend on Lamé parameters. In the drilling process, the impact of drill-bit on the rock generates vibrations. These vibrations vary depending on the drilling configurations and the Lamé parameters of the drilled rock. For example, the pure axial motion of the drill-bit in vertical direction generates pressure wavefield in the direction of motion and SV (Shear wave vertical component) wavefield in the direction perpendicular to the drill-bit motion. The radiation patterns of such motion are depicted in Figure 2. The vertical only motion of drill-bit does not generate SH (Shear wave horizontal component) wavefield. However, if the drill-bit rotates in the drilling process the impact of rotating motion on the layer can generate SH waves. These radiated seismic wave energies can be recorded by the receivers at the surface. The recorded seismic while drilling data provides vital information about the subsurface materials provided that we understand the behavior of the seismic radiation patterns of the drill-bit rock interaction.

A. SWD radiation patterns

While drilling, the drill-bit rock interaction radiates significant elastic, P-and S-wave, energy. In isotropic and homogeneous media, these radiations are functions of the drill-bit point force and the seismic velocities of rocks. It has been shown in [START_REF] Rector | Radiation pattern and seismic waves generated by a working roller-cone drill bit[END_REF], that for a roller-cone drill-bit, the seismic radiation pattern proceeding from the axial component drill-bit impacts can be modeled as a transient, monopolar point force acting along the axis of the borehole. Let us denote U r the Pwave radiation and U φ the S-wave motion. These radiations can be measured at the surface using the drill-bit seismic methods described above. More precisely, the far-field radial displacement resulting from the point force w(t, L) satisfies the following relation

U r (r, φ, t) = A 1 cos(φ) ρ f V 2 p r w(t - r V p , L), (9)
and the far-field angular displacement satisfies

U φ (r, φ, t) = A 1 sin(φ) ρ f V 2 s r w(t - r V s , L), ( 10
)
where r is the straight line distance from the source to the wavefront, ρ f is the rock density, V p is the rock compressional velocity, V s the rock shear velocity, A 1 is a constant, and the angle φ is measured relative to the direction of the point force (i.e. relative to the direction of axial drill-tooth impacts at the bottom of the borehole). The rock density ρ f can be expressed as a function of the rock compressional velocity V p . More precisely, we have (see [START_REF] Brocher | Empirical relations between elastic wavespeeds and density in the earth's crust[END_REF] for details)

ρ f = 1.74V 0.25 p , (11)
where V p is expressed in km.s -1 and ρ in g.cm -3 .

IV. EXPRESSION OF THE DRILL-BIT AXIAL FORCE

In this section, we express the axial force exerted at the drillbit as a (delayed) function of the topside velocity and hook load. This is done deriving the Riemann invariants associated to (1) and using a backstepping approach.

A. Derivation of Riemann invariants

The Riemann invariants of a hyperbolic PDE are the states corresponding to a transformation of the system which has a diagonalized transport matrix, i.e. the system can be written as a series of transport equations only coupled through the source terms [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]. On the set {(t, x)| 0 < t < T, x ∈ [0, L]}, we define the Riemann invariants as

u(t, x) = ∂ ∂t ξ(t, x) -c ξ ∂ ∂x ξ(t, x), (12)
z(t, x) = ∂ ∂t ξ(t, x) + c ξ ∂ ∂x ξ(t, x). (13)
This transformation enables us to rewrite (1) as the PDE system

∂ ∂t u(t, x) + c ξ ∂ ∂x u(t, x) = - k ξ 2 (u(t, ξ) + z(t, ξ)), (14)
∂ ∂t z(t, x) -c ξ ∂ ∂x z(t, x) = - k ξ 2 (u(t, ξ) + z(t, ξ)). ( 15
)
The boundary condition (5) rewrites

u(t, 0) = z(t, 0) + 2c ξ EA w 0 (t), (16)
while the boundary condition (8) rewrites

z(t, L) = -u(t, L) + 2 Ẋ(t), (17)
Ẍ(t) = - aζ M b ω bit Ẋ(t) - w f M b - EA 2c ξ M b (z(t, L) -u(t, L)). (18)

B. A backstepping transformation

In this section, we use a backstepping transformation [START_REF] Krstic | Boundary Control of PDEs: A Course on Backstepping Designs[END_REF] to map the system ( 14)-( 15) into a diagonal target system. This method consists in performing an integral change of variables (usually using a Volterra transformation) that maps the original system to a so-called "target system" for which the analysis is easier. Provided that the transformation is invertible (which is always the case for a Volterra transformation [START_REF] Krstic | Boundary Control of PDEs: A Course on Backstepping Designs[END_REF]), the original system and the corresponding target system have equivalent properties [START_REF] Krstic | Boundary Control of PDEs: A Course on Backstepping Designs[END_REF]. Let us first consider the following change of coordinates

û(t, x) = u(t, x)e k ξ 2c ξ x , ẑ(t, x) = z(t, x)e - k ξ 2c ξ
x , (19) so that the system ( 14)-( 15) rewrites

∂ ∂t û(t, x) + c ξ ∂ ∂x û(t, x) = - k ξ 2 e k ξ c ξ x ẑ(t, x), (20)
∂ ∂t ẑ(t, x) -c ξ ∂ ∂x ẑ(t, x) = - k ξ 2 e - k ξ c ξ x û(t, x). (21)
We consider the following Volterra transformation which is adjusted from the one defined in [START_REF] Coron | Local Exponential H 2 Stabilization of a 2 x 2 Quasilinear Hyperbolic System Using Backstepping[END_REF], [START_REF] Vazquez | Backstepping boundary stabilization and state estimation of a 2 x 2 linear hyperbolic system[END_REF] α(t, x) = û(t, x) +

x 0 (K uu (x, y)û(t, y)

+ K uz (x, y)ẑ(t, y))dy, (22)
β(t, x) = ẑ(t, x) + x 0 (K zu (x, y)û(t, y) + K zz (x, y)ẑ(t, y))dy, (23)
where the kernels K uu , K zu , K uz and K zz are continuous functions defined on the domain

T = {(x, y) ∈ [0, L] 2 , y ≤ x}.
They are defined by the following set of PDEs

∂ ∂x K uu (x, y) + ∂ ∂y K uu (x, y) = k ξ 2c ξ e - k ξ c ξ y K uz , (24)
∂ ∂x K uz (x, y) - ∂ ∂y K uz (x, y) = k ξ 2c ξ e k ξ c ξ y K uu , (25)
∂ ∂x K zu (x, y) - ∂ ∂y K zu (x, y) = - k ξ 2c ξ e - k ξ c ξ y K zz , (26)
∂ ∂x K zz (x, y) + ∂ ∂y K zz (x, y) = - k ξ 2c ξ e k ξ c ξ y K zu , (27)
with the boundary conditions

K uu (x, 0) = 0, K uz (x, x) = k ξ 4c ξ e k ξ c ξ x , (28)
K zz (x, 0) = 0, K zu (x, x) = - k ξ 4c ξ e - k ξ c ξ x . ( 29
)
It is shown that the system ( 24)-( 29) admits a unique continuous solution on T [START_REF] Coron | Local Exponential H 2 Stabilization of a 2 x 2 Quasilinear Hyperbolic System Using Backstepping[END_REF]. The proof is quite classical (see [START_REF] Hu | Control of Homodirectional and General Heterodirectional Linear Coupled Hyperbolic PDEs[END_REF], [START_REF] John | Continuous dependence on data for solutions of partial differential equations with a prescribed bound[END_REF] and [START_REF] Whitham | Linear and nonlinear waves[END_REF]) and consists in writing the integral equations associated to equations ( 24)- [START_REF] Djikpesse | Reducing uncertainty with seismic measurements while drilling[END_REF]. These integral equations are then solved using the method of successive approximations [START_REF] Krstic | Boundary Control of PDEs: A Course on Backstepping Designs[END_REF]. It is also possible to explicitly solve equations ( 24)-( 29) (see [START_REF] Vazquez | Marcum -functions and explicit kernels for stabilization of linear hyperbolic systems with constant coefficients[END_REF] for details). For all (x, y) ∈ T , we have

K uu (x, y) = k ξ 2c ξ e k ξ (x-y) 2c ξ yI 1 ( x 2 -y 2 k ξ 2c ξ ) x 2 -y 2 , ( 30
)
where I 1 is the modified Bessel function of first kind. The kernel K uv can be obtained using equations [START_REF] Greenberg | Seismic while drilling keeps bit turning to right while acquiring key real-time data[END_REF]. Finally the remaining kernels can be found by changing k ξ in -k ξ . Differentiating ( 22)-( 23) with respect to time and space and integrating by parts, one can show that the states α and β satisfy the following set of transport PDEs

∂ ∂t α(t, x) + c ξ ∂ ∂x α(t, x) = -c ξ K uz (x, 0)β(t, 0), (31)
∂ ∂t β(t, x) -c ξ ∂ ∂x β(t, x) = c ξ K zu (x, 0)α(t, 0). ( 32
)
with the topside boundary condition

α(t, 0) = β(t, 0) + 2c ξ EA w 0 (t). (33)
We choose not to write the boundary condition at x = L as this one is not required in the following. As the transformation ( 22)-( 23) is a Volterra transformation, it is invertible [START_REF] Krstic | Boundary Control of PDEs: A Course on Backstepping Designs[END_REF]. The inverse transformation is defined by

û(t, x) = α(t, x) + x 0 (L αα (x, y)α(t, y) + L αβ (x, y)β(t, y))dy, (34)
ẑ(t, x) = β(t, x) + x 0 (L βα (x, y)α(t, y) + L ββ (x, y)β(t, y))dy, (35)
where the kernels L αα , L βα , L αβ and L ββ are continuous functions defined on the domain T = {(x, y) ∈ [0, 1] 2 , y ≤ x}. They satisfy the following Volterra equation [START_REF] Krstic | Boundary Control of PDEs: A Course on Backstepping Designs[END_REF] defined for 0 ≤ y ≤ x ≤ 1.

L(x, y) = -K(x, y) - x ξ K(x, ν)L(ν, y)dν, (36)
where we have defined the functions K(•, •) and L(•, •) as K(x, y) = K uu (x, y) K uv (x, y) K vu (x, y) K vv (x, y) , L(x, y) = L αα (x, y) L αβ (x, y) L βα (x, y) L ββ (x, y) .

These inverse kernels are required to express the force exerted on the bit as a function of the topdrive force and velocity.

C. A delayed equation

In this section, we express for all x ∈ [0, L] the Riemann invariants α(t, x) and β(t, x) as functions of α(•, 0) and β(•, 0). Then, using the inverse transformation ( 34)- [START_REF] Richard | A simplified model to explore the root cause of stick-slip vibrations in drilling systems with drag bits[END_REF], it becomes possible to express for all x ∈ [0, L] the functions u(t, x) and z(t, x) as functions of α(•, 0) and β(•, 0). Applying the method of characteristics on (31)- [START_REF] Germay | Nonlinear Drillstring Dynamics Analysis[END_REF], for any x ∈ [0, L] and for any t ≥ x c ξ , we get

α(t, x) = α t - x c ξ , 0 -c ξ x c ξ 0 K uz (x -c ξ s, 0)β(t -s, 0)ds, (37)
β(t, x) = β t + x c ξ , 0 -c ξ x c ξ 0 K zu (x -c ξ s, 0)α(t + s, 0)ds. ( 38
)
Combining these equations with [START_REF] Aarsnes | Linear stability analysis of self-excited vibrations in drilling using an infinite dimensional model[END_REF], we obtain for all t

≥ L c ξ û(t, L) = α t - L c ξ , 0 -c ξ L c ξ 0 K uz (L -c ξ s, 0)β(t -s, 0)ds + L 0 L αα (L, y)α(t - y c ξ , 0) + L αβ (L, y)β(t + y c ξ , 0)dy -c ξ L 0 y c ξ 0 L αα (L, y)K uz (y -c ξ s, 0)β(t -s, 0)dsdy -c ξ L 0 y c ξ 0 L αβ (L, y)K zu (y -c ξ s, 0)α(t + s, 0)dsdy.
This implies, using Fubini's theorem and a simple change of variable

û(t, L) = α t - L c ξ , 0 -c ξ L c ξ 0 K uz (L -c ξ s, 0)β(t -s, 0)ds + c ξ L c ξ 0 L αα (L, c ξ y)α(t -y, 0) + L αβ (L, c ξ y)β(t + y, 0)dy -c ξ L c ξ 0 L c ξ s (L αα (L, y)K uz (y -c ξ s, 0)dy) β(t -s, 0)ds -c ξ L c ξ 0 L c ξ s L αβ (L, y)K zu (y -c ξ s, 0)dy α(t + s, 0)ds.
Finally, using a new change of variable, we obtain

û(t, L) = α t - L c ξ , 0 + L c ξ -L c ξ c ξ f u 1 (s)α(t -s, 0)ds + L c ξ -L c ξ c ξ f u 2 (s)β(t -s, 0)ds, (39)
where f u 1 and f u 2 are L ∞ functions defined on the interval

[-L c ξ , L c ξ ] by f u 1 (s) = 1 [0, L c ξ ] (s)(L αα (L, c ξ s)) -1 [-L c ξ ,0] (s) ( L -c ξ s L αβ (L, y)K zu (y + c ξ s, 0)dy), (40)
f u 2 (s) = 1 [-L c ξ ,0] (s)(L αβ (L, -c ξ s)) -1 [0, L c ξ ] (s)(K uz (L- c ξ s, 0) + L c ξ s L αα (L, y)K uz (y -c ξ s, 0)dy). (41)
Similar computations can be done on the state ẑ to obtain for all

t ≥ L c ξ ẑ(t, L) = β t + L c ξ , 0 + L c ξ -L c ξ c ξ f z 1 (s)α(t -s, 0)ds + L c ξ -L c ξ c ξ f z 2 (s)β(t -s, 0)ds, (42)
where f z 1 and f z 2 are L ∞ functions defined on the interval

[-L c ξ , L c ξ ] by f z 1 (s) = 1 [0, L c ξ ] (s)(L βα (L, c ξ s)) -1 [-L c ξ ,0] (s)( L -c ξ s L ββ (L, y) • K zu (y + c ξ s, 0)dy + K zu (L + c ξ s, 0)), (43)
f z 2 (s) = 1 [-L c ξ ,0] (s)(L ββ (L, -c ξ s)) -1 [0, L c ξ ] (s) • L c ξ s L βα (L, y)K uz (y -c ξ s, 0)dy . (44)
Using these expressions we are now able to express the axial force on the drill-bit as a function of the top-drive force and velocity.

D. Expression of the weight-on-bit

Our objective is now to express the axial force at the drillbit as a function of the topside velocity and hook load. We have

w(t, L) = -EA ∂ξ(t, L) ∂x = EA u(t, L) -z(t, L) 2c ξ = EA e - k ξ L 2c ξ û(t, L) -e k ξ L 2c ξ ẑ(t, L) 2c ξ . ( 45
)
Using equations ( 39) and ( 42), and noticing that α(•, 0) = u(•, 0) and β(•, 0) = z(•, 0), we obtain

w(t, L) = EA e - k ξ L 2c ξ u(t -L c ξ , 0) -e k ξ L 2c ξ z(t + L c ξ , 0) 2c ξ + L c ξ -L c ξ (f u (s)u(t -s, 0) + f z (s)z(t -s, 0)dy) ( 46
)
where the functions f u and f z are respectively defined by

f u (s) = EA 2 (e - k ξ L 2c ξ f u 1 (s) -e k ξ L 2c ξ f z 1 (s)) (47)
f z (s) = EA 2 (e - k ξ L 2c ξ f u 2 (s) -e k ξ L 2c ξ f z 2 (s)). ( 48
)
Using the expressions of the Riemann invariants u(t, x) and z(t, x) given in ( 12)-( 13) and the definition of the axial force and velocity, we immediately obtain

w(t, L) = EA 2c ξ (e - k ξ L 2c ξ v(t - L c ξ , 0) -e k ξ L 2c ξ v(t + L c ξ , 0)) + 1 2 (e - k ξ L 2c ξ w(t - L c ξ , 0) + e k ξ L 2c ξ w(t + L c ξ , 0)) + L c ξ -L c ξ (f u (s) + f z (s))v(t -s, 0)ds + L c ξ -L c ξ c ξ EA (f u (s) -f z (s))w(t -s, 0)ds. ( 49
)
It is interesting to note that this expression does not directly depend on the downhole boundary condition. In the next section, we describe a procedure to use this relation to estimate the compressional and shear velocities.

V. ROCK VELOCITY ESTIMATION Considering a single channel sensor, the radial component of the far-field recorded seismic-while-drilling data d obs = U obs r (r, φ, t) can be modelled by using equations ( 9) and ( 11) as

d cal = U cal r (r, φ, t) = A 1 cos(φ) 1.74V 0.25 p V 2 p r w(t - r V p , L), ( 50
)
where w is the calculated force-on-bit with equation [START_REF] Kazemi | Imaging with a seismic-whiledrilling dataset[END_REF].

As clear, the recorded data is the delayed and amplitude modulated version of the force-on-bit. Hence, in noise free data, simple cross-correlation of the calculated force-on-bit with the recorded data can provide a rough estimate of the velocity. In other words, by measuring the similarity between the recorded data and the function w, the estimated lag, i.e., r Vp , between the two data sets gives a rough estimate of the rock velocity V p . In the case of noisy data, assuming that the data is contaminated with white Gaussian noise, we propose to solve

{ Â1 , Vp } = argmin A1,Vp d obs -d cal 2 2 . ( 51
)
The optimization problem of equation ( 51) can be rewritten as

{ Â1 , Vp } = argmin A1,Vp N i=1 (U obs r (r, φ, t i ) -U cal r (r, φ, t i )) 2 , ( 52
)
where N is the number of time samples in the data. It is worth mentioning that equation ( 52) is only valid for direct arrival of the primary wave in the case of single isotropic and homogeneous rock. However, the SWD data contains different wave modes. In dealing with such dataset, to estimate the velocity of rock interacting with the drill-bit, similar to the work of Rector and Hardage [START_REF] Rector | Radiation pattern and seismic waves generated by a working roller-cone drill bit[END_REF] we should match the amplitude of the direct arrival in the source compensated SWD dataset with the amplitude of the modeled radiation pattern of the direct arrival given by equation [START_REF] Kazemi | Least-squares rtm of a seismic-while-drilling dataset[END_REF]. Nonetheless, source estimation is a necessary step in both cases. In this paper, we consider a single isotropic and homogenous rock simulation so the only wavefield generated is the direct arrival. Accordingly, we use equation ( 52) for the rock velocity estimation process. The cost function of equation ( 52) can be solved with alternating minimization algorithm (Algorithm 1). As the alternating minimization method is a local minimization technique, we use the rough velocity estimated by the cross-correlation approach as an initial solution. Steps 1 and 2 in the algorithm is solved with least-squares parameter estimation technique. Note that the proposed approach can be adjusted to estimate the shear velocity V s using [START_REF] Dunayevsky | Application of Stability Approach to Bit Dynamics[END_REF].

VI. SIMULATION RESULTS

In this section, we evaluate through simulations the performance of the estimation of the weight-on-bit obtained with equation [START_REF] Kazemi | Imaging with a seismic-whiledrilling dataset[END_REF] in the presence of noisy measurements. Then, we analyze the performance of the rock velocity estimation algorithm designed in Section V (i.e., Algorithm 1). We test the efficiency of the algorithm on single channel synthetic data, and then, extend the approach to a multichannel recording.

Algorithm 1 Alternating minimization algorithm for rock velocity estimation Require: d obs , w, r, L, φ, k = 0 Initialize:

V 0 p = r lag(w d obs )
While not converged 1:

Update Amplitude

A 1 k+1 = argmin A1 N i=1 (U obs r (r, φ, t i ) - A1 cos(φ) 1.74(V k p ) 0.25 (V k p ) 2 r w(t i -r V k p , L)) 2:
Update Velocity

V k+1 p = argmin Vp N i=1 (U obs r (r, φ, t i ) - A k+1 1 cos(φ) 1.74Vp 0.25 Vp 2 r w(t i - r Vp , L)) 3: Update k ← k + 1 If converged Output Vp ← V k p Â1 ← A k 1 Param. Value Param. Value E 2 × 10 11 Pa ρ 8000 kg.m -3 c ξ 5000 m.s -1 k ξ 0.3s -1 A 3.5 × 10 -3 m 2 L 2000 m M b 12000 kg ω bit 1.05s -1 a 0.1m ζ 0.6 w f 71280 N φ 0.

rad TABLE II NUMERICAL VALUES OF THE PARAMETERS

A. Estimation of the drill-bit source

Consider a drill-string whose axial motion satisfies the PDE given by equation (1) subject to the boundary conditions given in equations ( 5)- [START_REF] Germay | Multiple mode analysis of the self-excited vibrations of rotary drilling systems[END_REF]. Table II shows the numerical values of parameters for the drill-string system. The drill-bit is interacting with a sedimentary rock ( = 57 Jcm -3 ). The PDE system is simulated using a classical finite volume method based on a Godunov scheme [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]. The control input w 0 is chosen as a sinusoidal function. The corresponding force-onbit is denoted F real bit . By using the top-drive velocity and force measurements and equation ( 49), we can estimate the force-on-bit F est bit . To do so, first, we need to compute the kernels K •• . A numerical solution to equations ( 24)-( 25) can be reached with the method of characteristics and a fixed point algorithm (see [START_REF] Auriol | Robust design of backstepping controllers for systems of linear hyperbolic pdes[END_REF] for details) or the explicit expression given by equation [START_REF] Kazemi | Illumination compensation with seismic-while-drilling plus surface seismic imaging[END_REF]. The different integrals are computed using a trapezoidal method with adjustable precision. We consider the case of noisy measurements, the noise is modeled by a white Gaussian noise that is characterized by its signal-to-noise ratio (SNR). More precisely, similar to what is done in [START_REF] Kazemi | Surface-consistent sparse multichannel blind deconvolution of seismic signals[END_REF], we define the SNR as

SN R = a 2 rms σ 2 n , (53)
where a rms is the root-mean-square amplitude of the noisefree signal, and σ 2 n is the variance of the noise. Choosing a value for the SNR gives the corresponding amplitude for the white Gaussian noise. Note that the value of the SNR may be different for the velocity and force measurements.

The estimated force-on-bit F est bit is compared with the ground truth F real bit in Figure 3. We use the sedimentary rock and noisy top-drive measurements with SN R = 5. The bit is not rotating, only bouncing axially. We notice that after a time L c ξ , the functions are comparable (they are not identical due to the presence of the noise). Note that the estimated force-on-bit is not estimated for t ≤ L c ξ or t ≥ T f -L c ξ (where T f is the final time of the simulation) due to the fact that, on these intervals, we cannot compute all terms present in equation [START_REF] Kazemi | Imaging with a seismic-whiledrilling dataset[END_REF].

As is done in [START_REF] Kazemi | Surface-consistent sparse multichannel blind deconvolution of seismic signals[END_REF], to examine the performance of the proposed method, we introduce the quality of the reconstruction metric for the estimated force-on-bit. Let us denote y 0 as true generic signal (in our case y 0 = F real bit ) and let us denote by y its estimation (in our case y = F est bit ). We define the quality of the reconstruction Q as follows

Q = 10 log y 0 2 2 y 0 -y 2 2 , ( 54
)
where • 2 is the 2 norm. Figure 4 shows the variation of Q for different values of the SNR of top-drive velocity and the top-drive force (the SNR being chosen between 0.5 and 10).

We have also plotted a dashed line that corresponds to the limit value of Q = 10. We have a symmetrical effect between the SNR of top-drive velocity and the SNR of top-drive force. Moreover, the critical value of Q = 10 is not reached for low SNR values (SN R < 2). In our simulations, we find that when the quality factor is larger than 10, the rock velocity estimation algorithm results in satisfactory performances, even in the presence of severe noise in the SWD data. It is worth noting that the noise level in top-drive measurements is different from the noise level in the SWD data. Next, we evaluate the performance of the rock velocity estimation algorithm on single and multichannel SWD datasets.

B. Single sensor

We now use Algorithm (1) to estimate the velocity of the rocks that are interacting with the drill-bit. Let us consider that only one seismic sensor is available with its location being characterized by the angle φ and the distance r = L cos(φ) (see Section III for details). The numerical values of the different parameters corresponding to the drill-string system are given in Table II. We analyse the efficiency of the proposed algorithm on three different types of rocks with different hardness, i.e., different compressional/shear velocities,

• Unconsolidated sands (i.e. water saturated) for which V p = 1750 m.s -1 , ρ f = 2 g.cm -3 and = 11 J.cm -3 . • Sedimentary rocks (e.g. sandstone) for which V p = 3500 m.s -1 , ρ f = 2.4 g.cm -3 and = 57 J.cm -3 . • Igneous or metamorphic rocks (e.g. granite) for which V p = 5750 m.s -1 , ρ f = 2.7 g.cm -3 and = 151 J.cm -3 .

Using the simulated weight-on-bit and the (known) value of V p , we compute the function U obs r to which we add white Gaussian noise with SNR equal to 10. Moreover, by using the top-drive force and velocity measurements, which are given by the simulation of the PDE, we estimate the weight-on-bit (see equation ( 49)).We have plotted in Figure 5, the time evolution of the far-field radial displacement provided by almost noisefree seismic sensors for different kind of rocks (unconsolidated sands, sedimentary rocks, metamorphic rocks). The angle φ is equal to 10 degrees.

Let us test our approach in the case of unconsolidated sands. We have plotted in Figure 6 the cross correlation between U r and w. The largest spike occurs at t a = 1.155s which gives us the corresponding estimate

V 0 p = r t a = L t a cos(φ) = 1758 ms -1 ,
which is extremely close to the real value (1750 ms -1 ). Finally, for each type of rock, we give in Table III the different estimates for V p obtained using our Algorithm. We compare these results with the ones obtained using only the crosscorrelation estimate. One can notice that even if they are extremely similar, the proposed algorithm slightly improves the quality of the estimations.

C. Multiple sensors

Let us now consider that several sensors are available. For each sensor, the procedure described in the previous section can still be applied. Then, we can compute the mean and the standard deviations of the obtained set of data. Table IV shows the mean and the standard deviation of the obtained estimations for different SNR for the three different kinds of rocks using 100 different sensors. One can check that even in the presence of strong noise (SNR=0.1), the mean remains close to the real value while the standard deviation remains small. This highlights the quality of the proposed approach.

Our results indicate that the drilling process with a pure axial motion can provide a better estimation of the velocities in the case of unconsolidated rocks, i.e., near surface, than the hard rocks. Moreover, as the SNR decreases the performance of the method is deteriorated.

VII. FUTURE WORK

In this paper, we focused on the axial motion of the drill-bit. However, the torsional motion of the drill-bit also contributes to the overall signature of the drill-bit source. To properly model and estimate the elastic source signature both axial and torsional motions should be incorporated in the estimation system. Moreover, the main focus of the proposed algorithm was to provide a good estimate of the drill-bit source signature and to invert the direct arrival for the velocity of the rocks that are interacting with the drill-bit. Rector and Hardage [START_REF] Rector | Radiation pattern and seismic waves generated by a working roller-cone drill bit[END_REF] show that the radiation patterns for the direct arrival in the field data, after removing the source signature and compensating for the quality factor of the earth, can be approximated by equation [START_REF] Richard | Self-excited stickslip oscillations of drill bits[END_REF]. However, the SWD data contains a variety of wave modes such as up-going and down-going reflected wave fields, converted waves, head waves, rig arrivals, and so on. To properly see ahead of the drill-bit and monitor the rock properties, the full waveform should be considered and inverted. Recently, the authors proposed an imaging workflow for SWD data [START_REF] Kazemi | Imaging with a seismic-whiledrilling dataset[END_REF], [START_REF] Kazemi | Least-squares rtm of a seismic-while-drilling dataset[END_REF]. The results show that if the drill-bit source signature is known, the pre-stack depth imaging of such dataset is achievable. Moreover, authors also combined the SWD data with the surface seismic dataset and jointly image the subsurface [START_REF] Kazemi | Illumination compensation with seismic-while-drilling plus surface seismic imaging[END_REF]. Accordingly, the approach presented in this paper for estimating the drill-bit source signature can be used as an input to the SWD imaging workflows. The fullwaveform inversion (FWI) is also possible. FWI is an approach that directly inverts the full recorded seismic wavefield for the rock velocities of the subsurface [START_REF] Virieux | An overview of full-waveform inversion in exploration geophysics[END_REF]. As the FWI algorithm is highly nonlinear, the initial estimate of the rock velocity is an important factor for preventing the solution to be trapped in local minima. The source signature is also one of the input parameters of the conventional FWI method. The background velocity information usually is borrowed from the kinematics analysis of the surface seismic data. Hence, FWI inversion of the SWD dataset, after estimating the drill-bit source signature and using the background velocity information, is possible. In the case of SWD data, the FWI inversion result of surface seismic data could also be used as an initial estimate. Note that the background velocity information or the velocity calculated within surface seismic FWI can also be used as an initial estimate of the rock velocity that is interacting with the drillbit (i.e., initial velocity in Algorithm 1). The full-waveform inversion of the SWD dataset, in theory, can provide the subsurface velocity around and ahead of the drill-bit. To date, we are not aware of any developed FWI algorithm for SWD measurements.

VIII. CONCLUDING REMARKS

We have derived a new algorithm that provides a near realtime estimation of the seismic velocities of rocks interacting with the drill-bit, using seismic-while-drilling and drillstring top-drive force and top-derive velocity measurements. The proposed algorithm combined a distributed representation of the drill-string dynamics, seismic-while-drilling measurements, and classical parameter estimations methods. The interesting feature of the presented methodology is that it does not require any explicit knowledge about the subsurface properties, in so far as only surface measurements are required. The approach provided an efficient and reliable estimation of the seismic velocities of rocks interacting with the drill-bit, enabling a more precise characterization of the subsurface. Integrating the velocity estimation with the existing control techniques can, in theory, lead to an improvement in the drilling performance by reducing non-productive time, increasing the ROP, and reducing the effect of undesired torsional oscillations. This is an important step towards a subsurfaceaware drilling system. Testing the algorithm against real-field data is the focus of our future contributions.

Fig. 2 .

 2 Fig. 2. Seismic-while-drilling radiation patterns.

Fig. 3 .

 3 Fig.3. Comparison of the simulated force-on-bit and of the estimated one (computed using (49)) using noisy top-drive measurements in the case of sedimentary rock (SNR=5).

Fig. 4 .

 4 Fig.4. Sensitivity analysis of the drill-bit source estimation to different levels of noise in the top-drive measurements. We use 75 realizations for each noise level combination and report the average value of the quality of the reconstruction Q. The dashed line corresponds to Q = 10 contour which we consider it to be the success limit for the drill-bit source estimation.

Fig. 5 .

 5 Fig. 5. Comparison of the far-field radial displacement Ur (φ = 10 degrees) with respect to time for different rocks (unconsolidated sands, sedimentary rocks, metamorphic rocks).

Fig. 6 .

 6 Fig.6. Cross correlation between d obs and w for unconsolidated sands, the largest spike occurs at t = 1.155s (SNR=10).

ACKNOWLEDGMENTS This work was supported by the University of Calgary's Canada First Research Excellence Fund Program, the Global Research Initiative in Sustainable Low Carbon Unconventional Resources.