Jean Auriol

Federico Bribiesca-Argomedo

Delay-robust stabilization of an n + m hyperbolic PDE-ODE system

In this paper, we study the problem of stabilizing a linear ordinary differential equation through a system of an n + m (hetero-directional) coupled hyperbolic equations in the actuating path. The method relies on the use of a backstepping transform to construct a first feedback to tackle in-domain couplings present in the PDE system and then on a predictive tracking controller used to stabilize the ODE. The proposed control law is robust with respect to small delays in the control signal.

I. INTRODUCTION

In this paper we present a full-state feedback control design that delay-robustly stabilizes an interconnected system composed of a general system of coupled hetero-directional transport equations coupled through their unactuated boundary to a linear Ordinary Differential Equation (ODE). This delay-robust feedback law is obtained by partially leveraging the backstepping design in [START_REF] Di Meglio | Stabilization of coupled linear heterodirectional hyperbolic PDE-ODE systems[END_REF] and using the explicit mapping between hyperbolic systems and time-delay systems proposed in [START_REF] Auriol | An explicit mapping from linear first order hyperbolic PDEs to difference systems[END_REF]. Using such a time-delay formulation, we combine a full-state feedback of the Partial Differential Equation (PDE) states (to stabilize the PDE system) with a predictor-based tracking controller (that stabilizes the ODE).

The control of interconnected ODEs and hyperbolic PDEs is an extremely active research topic [START_REF] Bekiaris-Liberis | Compensation of wave actuator dynamics for nonlinear systems[END_REF], [START_REF] Bresch-Pietri | Adaptive output feedback for oil drilling stick-slip instability modeled by wave PDE with anti-damped dynamic boundary[END_REF], [START_REF] Di Meglio | Stabilization of coupled linear heterodirectional hyperbolic PDE-ODE systems[END_REF], [START_REF] Sagert | Backstepping and flatness approaches for stabilization of the stick-slip phenomenon for drilling[END_REF]. This class of systems naturally appears when modeling delays (that can be seen as first-order hyperbolic PDEs) in the actuating and sensing paths of ODEs [START_REF] Bekiaris-Liberis | Compensation of wave actuator dynamics for nonlinear systems[END_REF], [START_REF] Bresch-Pietri | Commande robuste de systèmes à retard variable: Contributions théoriques et applications au contrôle moteur[END_REF], [START_REF] Fridman | H∞-control of linear state-delay descriptor systems: an LMI approach[END_REF], [START_REF] Shivakumar | A generalized lmi formulation for input-output analysis of linear systems of odes coupled with pdes[END_REF], [START_REF] Sipahi | Stability and stabilization of systems with time delay[END_REF], [START_REF] Yue | Delayed feedback control of uncertain systems with time-varying input delay[END_REF]. For instance, one can consider the problem of the attenuation of mechanical vibrations in drilling applications (see [START_REF] Saldivar | A control oriented guided tour in oilwell drilling vibration modeling[END_REF] for a review of drilling vibrations models). For such a system, the hyperbolic PDEs model axial and torsional stress propagation (with finite speed of propagation) along the drill string. At the same time, the dynamics of the Bottom hole Assembly (BHA) can be represented by an ODE. The problem of delays acting on ODE systems has been tackled by the Smith predictor [START_REF] Smith | A controller to overcome dead time[END_REF]. Expressing these delayed terms as PDE states of transport equations [START_REF] Krstic | Backstepping boundary control for first-order hyperbolic PDEs and application to systems with actuator and sensor delays[END_REF], numerous related problems have been solved (non-constant and uncertain delays for instance) [START_REF] Bekiaris-Liberis | Nonlinear control under nonconstant delays[END_REF], [START_REF] Bresch-Pietri | Commande robuste de systèmes à retard variable: Contributions théoriques et applications au contrôle moteur[END_REF]. To solve the general problem of stabilizing an ODE with a system of first-order linear hyperbolic PDEs in the actuator path, a backstepping transformation has been introduced in [START_REF] Di Meglio | Stabilization of coupled linear heterodirectional hyperbolic PDE-ODE systems[END_REF] to map the original system into a cascade of exponentially stable subsystems. The major drawback of such a transformation is that it cancelled (among other terms) all the reflection terms at the actuated boundary. Although this is mathematically correct, such an approach completely neglects the impact on stability of small delays in the feedback loop (delay-robustness). Regarding the control of linear hyperbolic systems, it has been shown in recent contributions [START_REF] Auriol | Delay-robust control design for heterodirectional linear coupled hyperbolic PDEs[END_REF] inspired by [START_REF] Datko | An example on the effect of time delays in boundary feedback stabilization of wave equations[END_REF], [START_REF] Logemann | Conditions for robustness and nonrobustness of the stability of feedback systems with respect to small delays in the feedback loop[END_REF] (in which it has been observed that some feedback systems 1 Jean Auriol is with University of Calgary, Department of Petroleum Engineering, Calgary, Canada. jean.auriol@ucalgary.ca 2 Federico Bribiesca-Argomedo is with University of Lyon, INSA Lyon, Laboratoire Ampère CNRS UMR5005, F-69621, Villeurbanne, France possess a zero delay margin under linear state feedback) that some boundary reflections should not be compensated by the control action. Thus, a new approach, by leveraging the backstepping design given in [START_REF] Di Meglio | Stabilization of coupled linear heterodirectional hyperbolic PDE-ODE systems[END_REF] and complementing it with a predictor-type feedback has been proposed in [START_REF] Auriol | Delay-robust stabilization of a hyperbolic PDE-ODE system[END_REF] to guarantee the delay-robust stabilization of two linear PDEs coupled with an ODE. The control design and robustness analysis techniques developed in [START_REF] Auriol | Delay-robust control design for heterodirectional linear coupled hyperbolic PDEs[END_REF], [START_REF] Auriol | Delay-robust stabilization of a hyperbolic PDE-ODE system[END_REF] rely on a reformulation of the PDEs system as a difference equation. If only two coupled PDEs have been considered in [START_REF] Auriol | Delay-robust control design for heterodirectional linear coupled hyperbolic PDEs[END_REF], [START_REF] Auriol | Delay-robust stabilization of a hyperbolic PDE-ODE system[END_REF], such a mapping between hyperbolic PDEs and difference systems has been explicitly given in [START_REF] Auriol | An explicit mapping from linear first order hyperbolic PDEs to difference systems[END_REF] for an arbitrary number of PDEs, opening the perspective of adapting on hyperbolic PDE systems the stability analysis methods for time delay systems, such as those developed in [START_REF] Damak | Stability of linear continuoustime difference equations with distributed delay: Constructive exponential estimates[END_REF], [START_REF] Hale | Strong stabilization of neutral functional differential equations[END_REF], [START_REF] Niculescu | On delay-dependent stability under model transformations of some neutral linear systems[END_REF]. In recent contributions [START_REF] Aarsnes | Delay robust control design of under-actuated PDE-ODE-PDE systems[END_REF], an alternative approach has been considered to solve the problem of delay robust control design of an under-actuated PDE-ODE-PDE for which the PDEs are scalar (the last one being a transport equation). Finally, in [START_REF] Deutscher | Output feedback control of general linear heterodirectional hyperbolic ode-pde-ode systems[END_REF], the authors consider the stabilization of an ODE-PDE-ODE system (neglecting however the robustness properties of the feedback law).

The main contribution of this paper is the development of a new state-feedback control design (compared to [START_REF] Di Meglio | Stabilization of coupled linear heterodirectional hyperbolic PDE-ODE systems[END_REF]) for a PDE-ODE system, that ensures its delay-robust stabilization. Delay-robustness is guaranteed at the cost of preserving the proximal reflection terms in the target system. The approach of [START_REF] Auriol | Delay-robust stabilization of a hyperbolic PDE-ODE system[END_REF] cannot be directly extended due to the non-scalar structure of the problem. In this paper we proceed as follows: (i) Using two successive backstepping transformations (and associated feedback operators) adjusted from [START_REF] Auriol | An explicit mapping from linear first order hyperbolic PDEs to difference systems[END_REF], [START_REF] Coron | Finite-time boundary stabilization of general linear hyperbolic balance laws via Fredholm backstepping transformation[END_REF], we remove most of the in-domain couplings present in the PDEs. Without these in-domain couplings, we can prove that the stability properties of the original system are equivalent to the ones of an ODE coupled with a neutral system. (ii) The stabilization problem of the ODE-neutral system is reduced to the one of a linear ODE system with delayed input by solving a tracking problem. We then design a state-preictor feedback law. (iii) Finally, we study the delay-robustness properties of the system by algebraic means. The paper is organized as follows. In Section II we introduce the model under considerationand the notations. In Section III, we present the stabilization result: using successive backstepping transformations, we rewrite the original system as a distributed delay equation coupled with an ODE.After this reformulation, we design a stabilizing control law, shown to have a non-zero delay margin in Section IV. Finally, some simulation results are presented in Section V.

II. PROBLEM FORMULATION A. Definitions and notations

In this section we detail the notations used through this paper. For any integer p > 0, ||•|| R p is the classical euclidean norm on R p . We denote L 2 ([0, 1], R) the space of realvalued square-integrable functions defined on [0, 1] with the standard L 2 norm, i.e., for any

f ∈ L 2 ([0, 1], R), ||f || 2 L 2 = 1 0 f 2 (x)dx. The set L ∞ ([0, 1], R
) denotes the space of bounded real-valued functions defined on [0, 1] with the standard L ∞ norm, i.e., for any

f ∈ L ∞ ([0, 1], R), ||f || L ∞ = sup x∈[0,1] |f (x)|. In the following, for (u, v, X) ∈ (L 2 ([0, 1])) n+m × R p , we define the norm ||(u, v, X)|| = n i=1 ||u i || L 2 + m i=1 ||v i || L 2 + ||X|| R p . ( 1
)
The set C p ([0, 1]) (with p ∈ N ∪ {∞}) stands for the space of real-valued functions defined on [0, 1] that are p times differentiable and whose p-th derivative is continuous. The set T is defined as

T = {(x, ξ) ∈ [0, 1] 2 s.t. ξ ≤ x}. (2)
L ∞ (T ) stands for the space of real-valued L ∞ functions on T . For any (p, q) ∈ N, we denote M p×q (R) the set of matrices with p rows and q columns. The symbol I p (or I if no confusion arises) represents the p × p identity matrix. We use the notation f (s) for the Laplace transform of a function f (t), provided it is well defined. The set A stands for the convolution Banach algebra of BIBO-stable generalized functions in the sense of [START_REF] Vidyasagar | Input-output stability of a broad class of linear timeinvariant multivariable systems[END_REF]. A function g(•)

belongs to A if it can be expressed as

g(t) = g r (t) + ∞ i=0 g i δ(t -t i ), where g r ∈ L 1 (R + , R), i≥0 |g i | < ∞, 0 = t 0 < t 1 < . . . and δ(•) is the Dirac distribution. The associated norm is g A = g r L 1 + i≥0 |g i |.
The set  of Laplace transforms of elements in A is also a Banach algebra with associated norm ĝ  = g A .

B. System under consideration

We consider a class of systems consisting of an ODE coupled to a general linear heterodirectional first-order hyperbolic system in the actuation path.

More precisely, we consider systems of the form:

∂ t u(t, x) + Λ + ∂ x u(t, x) = Σ ++ (x)u + Σ +-(x)v, (3) ∂ t v(t, x) -Λ -∂ x v(t, x) = Σ +-(x)u + Σ --(x)v, (4) Ẋ(t) = AX(t) + Bv(t, 0), (5)
evolving in {(t, x) s.t. t > 0, x ∈ [0, 1]}, with the boundary conditions

u(t, 0) = Qv(t, 0) + CX(t), v(t, 1) = Ru(t, 1) + V (t), (6)
where X ∈ R p is the ODE state, and

u = (u 1 . . . u n ) T , v = (v 1 . . . v m ) T are the PDE states. The matrices Λ + and Λ - are diagonal (Λ + = diag(λ i ), Λ -= diag(µ i )) and their coefficients satisfy -µ m ≤ • • • ≤ -µ 1 < 0 < λ 1 ≤ • • • ≤ λ n .
The constant real boundary coupling matrices are defined by

Q = {q ij } 1≤i≤n,1≤j≤m and R = {ρ ij } 1≤i≤m,1≤j≤n .
The spatially-varying inside domain couplings matrices are defined as follows

Σ ++ (x) = {σ ++ ij (x)} 1≤i,j≤n Σ +-(x) = {σ +- ij (x)} 1≤i≤n,1≤j≤m Σ -+ (x) = {σ -+ ij (x)} 1≤i≤m,1≤j≤n Σ --(x) = {σ -- ij (x)} 1≤i,j≤m ,
their coefficients are assumed to belong in C 0 ([0, 1], R).

We assume that the diagonal terms of the matrices Σ ++ and Σ --are equal to zero (such coupling terms can be removed using a change of coordinates, see [START_REF] Coron | Local exponential H 2 stabilization of a 2× 2 quasilinear hyperbolic system using backstepping[END_REF] for details). The matrices A, B and C respectively belong to M p×p (R), M p×m (R) and M n×p (R) The initial conditions of the state (u, v) are denoted u 0 and v 0 and are assumed to respectively belong to

(L 2 ([0, 1], R)) n and (L 2 ([0, 1], R)) m .
The initial condition of the ODE ( 5) is denoted X 0 . The function V is an input function (control law) that has values in R m . The resulting system (3)-( 6) is well-posed [6, Theorem A.6, page 254]. Remark that this system naturally features several couplings that can be source of instabilities. Finally, we denote τ = 1 λ1 + 1 µ1 the sum of the largest transport times in each direction.

C. Control problem

Let us consider the two following assumptions Assumption 1: The pair (A, B) is stabilizable, i.e. there exists a matrix K such that A + BK is Hurwitz.

Assumption 2: The system defined for all i ∈ [1, m] by

z i (t) = 1≤l≤m 1≤k≤n R ik Q kl z l (t - 1 λ k - 1 µ l ). (7)
is exponentially stable In this paper, we show that if Assumption 1 and Assumption 2 are satisfied, then it is possible to explicitly design a feedback control law

V = K[(u, v, X)] (where K : (L 2 [0, 1]) 2 × R p → R
is a linear operator) that delayrobustly stabilizes (in the sense of [START_REF] Logemann | Conditions for robustness and nonrobustness of the stability of feedback systems with respect to small delays in the feedback loop[END_REF]) system ( 3)-( 6), i.e.:

• the state (u, v, X) of the resulting feedback system (3)-( 6) exponentially converges to its zero equilibrium (stabilization problem), i.e. there exist κ 0 ≥ 0 and ν > 0 such that for any initial condition

(u 0 , v 0 , X 0 ) ∈ (L 2 [0, 1]) 2 × R p ||(u, v, X)|| ≤ κ 0 e -νt ||(u 0 , v 0 , X 0 )||, t ≥ 0. (8)
• the resulting feedback system (3)-( 6) is robustly stable with respect to small delays in the loop (delayrobustness), i.e. there exists δ > 0 such that for any δ i (∈ [0, δ ]) m , the control law whose components are defined by (V (t -δ i )) i still stabilizes (3)-( 6).

The first assumption (stabilizability of the ODE subsystem) is necessary for the stabilizability of the whole system. It can be shown that violating the second assumption results in an open-loop transfer function with an infinite number of poles on the closed right half-plane.Consequently (see [START_REF] Logemann | Conditions for robustness and nonrobustness of the stability of feedback systems with respect to small delays in the feedback loop[END_REF]Theorem 1.2]), one cannot find any linear state feedback law V (•) that delay-robustly stabilizes (3)-( 6). Note that if the delays 1 λ k + 1 µ l are rationally independent a necessary and sufficient condition for this assumption to be satisfied is given by [START_REF] Hale | Introduction to functional differential equations[END_REF]Theorem 6.1] in terms of spectral radius.

III. DESIGN OF THE CONTROL LAW

In this section we derive a control law that guarantees the stabilization of ( 3)-( 6), following the methodology introduced above. Using two successive backstepping transformations, we map the original system to a simpler target system for which local coupling terms Σ •• have been moved to the boundary. The target system is recast as a neutral system for which a predictor-based control law is developed.

A. Backstepping transformations

We derive two backstepping transformations to remove the in-domain coupling terms of (3)- [START_REF] Bastin | Stability and boundary stabilization of 1-D hyperbolic systems[END_REF]. These transforms will introduce some non-local coupling terms yet keep pointwise coupling terms, ensuring a non-zero delay margin for the system.Let us consider the following Volterra transformation adjusted from the one defined in [START_REF] Auriol | An explicit mapping from linear first order hyperbolic PDEs to difference systems[END_REF], [START_REF] Coron | Finite-time boundary stabilization of general linear hyperbolic balance laws via Fredholm backstepping transformation[END_REF], [START_REF] Hu | Boundary exponential stabilization of 1-D inhomogeneous quasilinear hyperbolic systems[END_REF]

α = u- x 0 (K uu (x, ξ)u(ξ) + K uv (x, ξ)v(ξ)) dξ + γ 0 (x)X(t), (9)
β = v- x 0 (K vu (x, ξ)u(ξ) + K vv (x, ξ)v(ξ)) dξ, + γ 1 (x)X(t), (10)
where the kernels K uu , K uv , K vu , K vv belong to L ∞ (T ) (where T is introduced in (2)), while γ 0 and γ 1 are differentiable matrices defined on ([0, 1]). They satisfy the following set of PDEs

Λ + ∂xK uu + ∂ ξ K uu Λ + = -K uu Σ ++ (ξ) -K uv Σ -+ (ξ), (11)
Λ + ∂xK uv -∂ ξ K uv Λ -= -K uu Σ +-(ξ) -K uv Σ --(ξ), (12)
Λ -∂xK vu -∂ ξ K vu Λ + = K vu Σ ++ (ξ) + K vv Σ -+ (ξ), (13)
Λ -∂xK vv + ∂ ξ K vv Λ -= K vu Σ +-(ξ) + K vv Σ --(ξ), (14)
and ODEs

Λ + γ 0 (x) = -γ 0 (x)A + K uu (x, 0)Λ + C, (15)
Λ -γ 1 (x) = γ 1 (x)A -K vu (x, 0)Λ + C, (16)
with the boundary conditions

Σ ++ (x) -Λ + K uu (x, x) + K uu (x, x)Λ + = 0, ( 17
) Σ +-(x) -Λ + K uv (x, x) -K uv (x, x)Λ -= 0, (18)
Σ -+ (x) + Λ -K vu (x, x) + K vu (x, x)Λ + = 0, (19) Σ --(x) + Λ -K vv (x, x) -K vv (x, x)Λ -= 0, ( 20
) (K vv (x, 0)Λ -) ij = (K vu (x, 0)Λ + Q -γ 1 (x)B) ij (21) γ 0 (0) = 0, γ 1 (0) = 0, (22)
where in equation ( 21), we have 1 ≤ j ≤ i ≤ m. To ensure well-posedness of the kernel equations, we add the following artificial boundary conditions for

K uu ij (n ≥ i > j ≥ 1) on ξ = 0: K uu ij (x, 0) = k ij (x) ∈ C 1 ([0, 1]).
We have the following lemma Lemma 1: Consider system (11)- [START_REF] Krstic | Boundary control of PDEs: A course on backstepping designs[END_REF]. There exists a unique solution K uu , K uv , K vu and K vv in L ∞ (T ) and differentiable matrices γ 0 , γ 1 .

Proof: This result follows, with some minor adaptations, from [START_REF] Hu | Boundary exponential stabilization of 1-D inhomogeneous quasilinear hyperbolic systems[END_REF] and [START_REF] Di Meglio | Stabilization of coupled linear heterodirectional hyperbolic PDE-ODE systems[END_REF]Theorem 3.2]. The main idea consists in reinterpreting the ODEs as PDEs evolving in the triangular domain T with horizontal characteristic lines (since there is only an evolution along x) and then solving all the PDEs together. In the case n = m = 1, the proof can be found in [START_REF] Auriol | Delay-robust stabilization of a hyperbolic PDE-ODE system[END_REF]. Due to page limitation, we do not give the complete proof here. Since the transformation ( 9)-( 10) is a Volterra transformation, it is invertible (see e.g. [START_REF] Krstic | Boundary control of PDEs: A course on backstepping designs[END_REF]). The corresponding inverse kernels are denoted L αα , L αβ , L βα , L ββ , γ0 and γ1 . They have the same regularity.

u = α- x 0 L αα (x, ξ)α(ξ) + L αβ (x, ξ)β(ξ) dξ + γ0 (x)X(t), (23)
v = β- x 0 L βα (x, ξ)α(ξ) + L ββ (x, ξ)β(ξ) dξ, + γ1 (x)X(t). (24)
Defining the

L ∞ matrix G 1 (x) = γ 0 (x)B -K uu (x, 0)Λ + Q + K uv (x, 0)Λ -,
and the upper-triangular matrix G 2 (x) = -K vu (x, 0)Λ + Q+K vv (x, 0)Λ -+γ 1 (x)B, the invertible transformation ( 9)-( 10) maps the original system (3)-( 6) to the target system

∂ t α(t, x) + Λ + ∂ x α(t, x) = G 1 (x)β(t, 0), ( 25
) ∂ t β(t, x) -Λ -∂ x β(t, x) = G 2 (x)β(t, 0), (26)
Ẋ(t) = AX(t) + Bβ(t, 0), (27)
along with the boundary conditions

α(t, 0) = Qβ(t, 0) + CX(t), (28)
β(t, 1) = Rα(t, 1) + V (t) + (Rγ 0 (1) -γ1 (1))X(t) + 1 0 N α (ξ)α(t, ξ) + N β (ξ)β(t, ξ) dξ, ( 29
)
where

N α (ξ) = L βα (1, ξ) -RL αα (1, ξ) and N β (ξ) = L ββ (1, ξ) -RL αβ (1, ξ).
The associated initial condition, denoted (α 0 , β 0 , X 0 ), is related to the initial condition (u 0 , v 0 , X 0 ) by the transformation ( 9)- [START_REF] Bresch-Pietri | Adaptive output feedback for oil drilling stick-slip instability modeled by wave PDE with anti-damped dynamic boundary[END_REF]. Let us now consider the following transformation which is invertible due to its cascade structure

α(t, x) = w(t, x), (30)
β(t, x) = z(t, x) - 1 0 F (x, ξ)z(t, ξ)dξ, (31)
where the kernel F is a strict upper triangular matrix (i.e

F ij (x, ξ) = 0 if i ≥ j) defined on T 1 = {(x, ξ) ∈ [0, 1] 2 }
by the following set of PDEs

Λ -F x (x, ξ) + F ξ (x, ξ)Λ -= 0, (32)
along with the boundary conditions

F (x, 0) = G 2 (x)(Λ -) -1 , F (0, ξ) = 0. (33)
It is straightforward to show that this set of PDEs admits a unique solution

F ∈ L ∞ ([0, 1] × [0, 1]
) m×m whose components can explicitly be obtained by the method of characteristics [START_REF] Coron | Finite-time boundary stabilization of general linear hyperbolic balance laws via Fredholm backstepping transformation[END_REF]. Differentiating ( 30)-( 31) with respect to time and space, we obtain that the system ( 25)-( 29) (and consequently the system (3)-( 6)) is equivalent to

∂ t w(t, x) + Λ + ∂ x w(t, x) = G 1 (x)z(t, 0), ( 34
) ∂ t z(t, x) -Λ -∂ x z(t, x) = G 3 (x)z(t, 1), (35) Ẋ(t) = AX(t) + Bz(t, 0), (36)
with the boundary conditions w(t, 0) =Qz(t, 0) + CX(t), (37) z(t, 1) =Rw(t, 1) + V (t) + (Rγ 0 (1) -γ1 (1))X(t)

+ 1 0 N α (ξ)w(t, ξ) + N z (ξ)z(t, ξ)dξ, ( 38
)
where G 3 is the unique L ∞ solution of the equation G 3 (x) = F (x, 1)Λ -+ 1 0 F (x, ξ)G 3 (ξ)dξ (which can be directly solved due to the triangular structure of F ). and where N z (ξ) = F (1, ξ)+N β (ξ)-1 0 N β (ν)F (ν, ξ)dν. Note that G 3 is upper-triangular. For the control design of the target system (34)-(36) with the boundary conditions (37)-(38), we decompose the control input V (t) as V (t) = V 1 (t) + V 0 (t), where V 1 (•) will be designed in the next sectionsand where V 0 (t) is defined by,

V 0 (t) = -(Rγ 0 (1) -γ1 (1))X(t) - 1 0 (N α (ξ)w(t, ξ) + N z (ξ)z(t, ξ))dξ, (39)
so that (38) rewrites z(t, 1) = Rw(t, 1)+V 0 (t). Remark that, due to the invertibility of the backstepping transformations ( 9)-( 10) and ( 30)-( 31), V 0 (t) can be expressed in terms of u, v and X.

B. A neutral system satisfied by z(t, 1)

Using the specific structure of the target system ( 34)-( 38), it is possible to express the state z(t, 1) as the solution of a neutral system. We have the following theorem.

Theorem 1: There exists G 4 (•) an L ∞ ([0, τ ]) function such that for all 1 ≤ j ≤ m and all 1 ≤ i ≤ n and all t ≥ τ

z j (t, x) = z j (t - 1 -x µ j , 1) + m k=1 1-x µ j 0 (G 3 ) jk (x + µ j s) • z k (t -s, 1)ds, (40)
w i (t, x) = m k=1 (Q) ik z k (t - x λ i - 1 µ k , 1) + p k=1 C ik X k (t- x λ i ) + m l=1 τ 0 (G 4 ) il (x, η)z l (t -η, 1)dη. (41)
Consequently, combining this equation with (38), there exists m × m L ∞ ([0, τ ])-functions denoted G l j (•) (1 ≤ j, l ≤ m) such that for all t ≥ τ , for all 1 ≤ i ≤ m we have

z i (t, 1) = V 1 i (t) + n k=1 m l=1 R ik Q kl z l (t - 1 λ k - 1 µ l , 1)+ m l=1 τ 0 G l i (s)z l (t -s)ds + n,p k,l=1 R ik C kl X l (t - 1 λ k ) (42)
Proof: The proof follows the same steps as the ones of [START_REF] Auriol | An explicit mapping from linear first order hyperbolic PDEs to difference systems[END_REF]Theorem 4]. It relies on the method of characteristics. It is not given here, due to the page limitation. In the following we decompose the control law

V 1 (t) as V 1 (t) = V BS (t)+V ODE (t) where V BS is defined for all 1 ≤ i ≤ m by (V BS ) i (t) = - n k=1 p l=1 R ik C kl X l (t - 1 λ k ) - τ 0 G l i (s)z l (t -s)ds, (43)
while V ODE (•) has to be designed for the stabilization of the ODE dynamics. Remark that V BS can be expressed as delayed values of u, v and X. With this control law, equation (42) rewrites

zi(t, 1) = n k=1 m l=1 R ik Q kl z l (t - 1 λ k - 1 µ l , 1) + VODE(t). ( 44
)
As the principal part of system (44) generates an exponentially stable semigroup (due to Assumption 2), the control law V BS would guarantee the stabilization of z(t, 1) in the absence of the ODE. The objective now is to design V ODE such that the ODE state also converges to zero. If such a stabilizing control law only depends on X(•), then once the ODE has been stabilized, we are brought to the previous situation and have the convergence of z(t, 1) to zero, and consequently of the system (w, z), due to the cascade structure of system ( 34)-(38) (and using the invertibility of the backstepping transformations, the stability of the system (u, v)).

Remark 1: Contrary to [START_REF] Di Meglio | Stabilization of coupled linear heterodirectional hyperbolic PDE-ODE systems[END_REF] ,we do not cancel the principal part

n k=1 m l=1 R ik Q kl z l (t -1 λ k -1 µ l , 1
) in the control law to guarantee the existence of some delay-margins (see [START_REF] Auriol | Delay-robust control design for heterodirectional linear coupled hyperbolic PDEs[END_REF] for details). However, as it is done in [5, Theorem 5], we could cancel a small part of these reflection terms to improve the convergence rate.

C. A tracking problem solved by a predictor

To stabilize the ODE (36), we would like to directly control z(t, 0), and ideally we would like z(t, 0) = KX(t) (where K is defined in Assumption 1 and is such that A + BK is Hurwitz). This corresponds to a tracking problem. Using equation (35) (and the fact that G 3 is upper triangular), we have for all

1 ≤ i ≤ m z i (t, 0) = z i (t - 1 µ i , 1) + m k=i+1 1 µ i 0 (G 3 ) ik (µ i s)z k (t -s, 1)ds.
This yields

z i (t, 1) = z i (t + 1 µ i , 0) - m k=i+1 1 µ i 0 (G 3 ) ik (µ i s)z k (t + 1 µ i -s, 1)ds, (45)
and using the cascade structure of the integral part, we can express z(t, 1) as a function of future values of z(•, 0) which may lead to a non causal problem. Let us denote

1 μ = m k=1 1
µ k and z(t, 0) = z(t + 1 μ , 0). If instead of tracking z(•, 0), we choose to track z(•, 0), the problem becomes causal as equation ( 45) rewrites z i (t, 1) = zi (t-

1 μ + 1 µi , 0)- m k=i+1 1 µ i 0 (G 3 ) ik (µ i s)z k (t + 1 µi -1 μ -s, 1)
ds. Thus, if we want z(•, 0) to be equal to an arbitrary function ζ(•), it immediately implies

z m (t, 1) =ζ m (t - 1 μ + 1 µ m ), (46)
. . .

z 1 (t, 1) =ζ 1 (t - 1 μ + 1 µ 1 ) - m k=2 1 µ 1 0 (G 3 ) 1k (µ 1 s)z k (t + 1 µ 1 -s, 1)ds. ( 47
)
In what follows, the function ζ will be chosen to stabilize the ODE system. Using equations ( 46)-( 47), it becomes possible to express the state z(t, 1) as a function of ζ(•). More precisely, we have the following lemma Lemma 2: There exist L ∞ ([0, 1 μ ])-functions denoted H l j (•) (1 ≤ j, l ≤ m) that only depend on G 3 and on ζ, such that, if for all 1 ≤ i ≤ m we have

z i (t, 1) = χ(t) = ζ i (t - 1 μ + 1 µ i ) + m k=i+1 1 μ 1 μ -m l=k 1 µ l H ik (s)ζ k (t -s)ds, (48)
then zi (t, 0) = ζ i (t) and consequently z i (t, 0) = ζ i (t -1 μ ). Proof: The proof is straightforward by recursion, the initialization (i = m + 1) corresponding to (46). Let us assume that we have z(t, 0) = ζ(t -1 μ ). The ODE system (36) rewrites for

t ≥ τ + 1 μ Ẋ(t) = AX(t) + Bζ(t - 1 μ ), (49)
which is a finite-dimensional system with a delayed input ζ(•). Different methods [START_REF] Zhong | Robust control of time-delay systems[END_REF] can be used to design a control law that stabilizes equation (49). A classical result from [START_REF] Mirkin | Every stabilizing dead-time controller has an observer-predictor-based structure[END_REF] states that any control law that stabilizes such an equation is equivalent to a predictor. More precisely, let us define the function P by

P (t) = K e A μ X(t) + t t-1 μ e A(t-ν) BP (ν)dν , (50) It can be verified that if ζ(t) = P (t), then the system Ẋ = AX + Bζ(t -1 μ
) is exponentially stable (as in this case we have ζ(t) = P (t) = KX(t + 1 μ )). Consequently, we want z(t, 0) to converge to ζ(t -1 μ ) = P (t -1 μ ). This gives us the corresponding value of z(t, 1), using Lemma 2. Let us define for all 1 ≤ i ≤ m the control law (V ODE ) i (•) by

(V ODE ) i (t) =χ i (t) - n k=1 m l=1 R ik Q kl χ l (t - 1 λ k - 1 µ l ), (51)
where

χ i (t) = P i (t -1 μ + 1 µi ) + m k=i+1 1 μ 1 μ -m l=i 1 µ l H ik (s)P k (t -s)ds
, where the function H ik are defined in Lemma 2. We can now conclude this section with the following theorem. Theorem 2: The feedback law

V (t) = V ODE (t) + V BS (t) + V 0 (t), (52)
where V ODE V BS and V 0 are respectively defined by (51), ( 43) and (39), exponentially stabilizes (in the sense of equation ( 8)) the system (3)-( 6) to its zero-equilibrium.

Proof: Using the computations done above (and in particular equation (44)), and denoting e(t) = z(t, 1) -χ(t), we have for all

1 ≤ i ≤ m e i (t) = n k=1 m l=1 R ik Q kl (e l (t - 1 λ k - 1 µ l , 1)).
Thus, due to Assumption 2, this implies the convergence of the function e to zero and consequently, the convergence of z(t, 1) to χ(t). Using Lemma 2, we obtain the convergence of z(t, 0) to P (t -1 μ ). This implies the stabilization of X(t). Furthermore, the state-predictor P (t) converges to KX(t + μ), which implies that P (t) (and consequently χ(t) and V ODE (t)) exponentially converges to zero. Consequently, due to (44), z(t, 1) exponentially vanishes. Using (35), this implies that z(t, •) converges L 2exponentially to zero. This implies from [START_REF] Zhong | Robust control of time-delay systems[END_REF] and the boundary condition (37), that w(t, ) converges also L 2exponentially to zero. This yields the existence of 0 < κ 0 such that ||(w, z, X)|| ≤ κ 0 e -νt ||(w 0 , z 0 , X 0 )||. Thus, the control law V (t) ensures the exponential stabilization of ( 34)-(38). Due to the invertibility of the backstepping transformations ( 9)-( 10) and ( 30)- [START_REF] Smith | A controller to overcome dead time[END_REF], it is straightforward to prove the stabilization of ( 3)-( 6).

IV. DELAY-ROBUST STABILIZATION

In this section we prove the delay-robustness of the control law designed in the previous section. To do so, we show that the associated characteristic equation has all its zeros on the left-half plane. Let us consider a vector of positive delays δ = (δ

1 δ 2 • • • δ m )
T acting on the actuation input V (•) defined in (52) (each δ i acting on the corresponding component of V ). Using equations ( 40) and (41) and doing some change of variables (see [START_REF] Auriol | An explicit mapping from linear first order hyperbolic PDEs to difference systems[END_REF] for details) we have (

1 0 N α w(t, ξ) + N z z(t, ξ)dξ) i = m k=1 τ 0 (G 5 ) ik (s)z k (t - s, 1)ds + p l=1 τ 0 (G X ) il (s)X l (t -s)ds, where for all 1 ≤ i ≤ m, 1 ≤ k ≤ m and 1 ≤ l ≤ p, the functions (G 5 ) ik and (G X ) il are L ∞ functions defined in [0, τ ].
Thus, in presence of the delayed control law V (•), we now get for all 1 ≤ i ≤ m and all t ≥ τ + max{δ i }, denoting y(t) = z(t, 1), we have

yi(t) = n,m k,l=1 R ik Q kl y l (t - 1 λ k - 1 µ l ) + m l=1 τ 0 G l i y l (t -s)ds + n k=1 p l=1 R ik C kl X l (t - 1 λ k ) +
G l i y l (t -s -δi)ds - n,p k,l=1 R ik C kl X l (t - 1 λ k -δi) - p l=1 (Rγ0(1) -γ1(1))ilXl(t -δi) - m k=1 τ 0 (G5) ik (s)• y k (t -s -δi)ds - p l=1 τ 0 (GX ) il (s)X l (t -s -δi, 1)ds + χi(t -δi) - n k=1 m l=1 R ik Q kl χ l (t - 1 λ k - 1 µ l -δi) (53)
Taking the Laplace transfom of (53), we obtain

F (s, δ)(ŷ(s) -χ(s)) = (I -∆)(H 1 (s) X(s) + H 2 (s) χ(s)),
where for all 1 ≤ i, j ≤ m, for all 1 ≤ l ≤ p

(F (s, δ)) ij = η j i - n k=1 R ik Q kj e -( 1 λ k + 1 µ j )s -(1 -e -δis ) ( τ 0 (G 5 ) ij (s)e -νs dν + τ 0 G j i (s)e -νs dν), (54)
(H 1 (s)) il = ( n k=1 R ik C kl e -s λ k + (Rγ 0 (1) -γ1 (1)) il + τ 0 (G X ) il (s)e -sν dν), (55)
(H 2 (s)) ij = -(η j i - n k=1 R ik Q kj e -( 1 λ k + 1 µ j )s ) + τ 0 G j i (s)e -νs dν + τ 0 (G 5 ) ij (s)e -νs dν, (56)
with ∆ i,j = η j i e -δis , where η j i stands for the Kronecker symbol which is equal to one if i = j, and is equal to zero otherwise. To obtain the characteristic equation from (54), we need first to prove F is invertible on the Right Half Plane (RHP). From [32, Theorem 1], we know that F (•, δ) ∈ A has a unique inverse in A if and only if inf Re(s)≥0 | det(F (s, δ))| > 0. We have the following lemma on invertibility of F 1 (s, δ) in  (where the Banach algebra  is defined in section II-A).

Lemma 3: There exists δ ∈ (0, τ ] such that

inf δ∈([0,δ ]) m inf Re(s)≥0 | det(F (s, δ))| > 0. ( 57
)
Proof: The proof is analogous to the one of [Theorem 9][5] and is a consequence of Assumption 2 and of the fact that the δ i is arbitrarily small. We now obtain the characteristic equation associated to (54). We get from (35) thatẑ(t, 0) = M (s)ẑ(t, 1), where M is the upper-triangular matrix defined for all

1 ≤ i ≤ j ≤ m by (M (s)) ij = e -1 µ i s η j i + 1 µ i 0 (G 3 ) ij (µ i s)e -sν dν
. By construction (see equation ( 46)-( 47)), we also have χ(s) = M (s)e -s μ P (s) where P is the Laplace transform of the predictor state feedback given in (50), namely P (s) = K 0 (s) X(s) where

K 0 (s) = I -K(sI -A) -1 (I -e -(sI-A) 1 μ )B -1 Ke A μ .
Taking the Laplace transform of (36), we obtain for s ∈ C such that Re(s) ≥ 0

(sI -A -BK 0 (s)e -s μ ) X(s) = B ẑ(s, 0) -BK 0 (s)e -s μ X(s) = BM (s)(ẑ(s, 1) -χ(s)) = BM (s)(F (s, δ)) -1 (I -∆)(H 1 (s) + H 2 (s)M (s) K 0 (s)e -s μ ) X(s), (58)
as F is invertible. We are now finally able to prove that the control law V (t) as defined in (52) delay-robustly stabilizes the system (3)-( 6) by proving that the characteristic equation (58) does not have any root on the RHP. Theorem 3: The control law V (t) defined in Theorem 2 delay-robustly stabilizes the system (3)- [START_REF] Bastin | Stability and boundary stabilization of 1-D hyperbolic systems[END_REF]. That is, there exists δ > 0 such that, for all δ i ∈ ([0, δ ]) m , the control law defined for all 1 ≤ i ≤ m by V (t) = (V ODE (tδ i ) + V BS (t -δ i ) + V 0 (t -δ i )) i exponentially stabilizes the system (3)- [START_REF] Bastin | Stability and boundary stabilization of 1-D hyperbolic systems[END_REF].

Proof: The closed-loop characteristic equation associated to (58) can be written

p(s) = det(F 0 (s) -E 0 (s)(I -∆)H 0 (s)) = 0, (59)
where

F 0 (s) =sI -A -BK 0 (s)e -s μ , E 0 (s) = BM (s)(F (s, δ)) -1 , H 0 (s) =(H 1 (s) + H 2 (s)M (s)K 0 (s)e -s μ )
The holomorphic function F 0 (s) has all its roots in the left-half complex plane as the system Ẋ = AX + BP (t -1 μ ) is exponentially stable. As the holomorphic functions M, H 1 , H 2 and F -1 are bounded in the right-half complex plane (see Lemma 3 for the boundedness of F -1 ), so are E 0 (s) and H 0 (s) and the function det(E 0 (s)(I -∆)H 0 (s)). As the holomorphic function F 0 (s) goes to infinity for |s| large enough, there exists

M 2 > 0 such that ∀s ∈ Ω 1 = {s ∈ C, (s) ≥ 0 and |s| ≥ M 2 } | det(F 0 (s) -E 0 (s)(I -∆)H 0 (s))| > 0 Let us now consider s ∈ C ∈ Ω 2 = {s ∈ C, (s) ≥ 0 and |s| ≤ M 2 }.
By contradiction, assume that there exists s ∈ C, s = 0 and Re(s) ≥ 0, such that p(s) = 0. There exists η = 0 such that F 0 (s)η = E 0 (s)(1-∆)H 0 (s)η. This yields

η * F * 0 (s)F 0 (s)η = η * H * 0 (s)(I-∆) * E * 0 (s)E 0 (I-∆)H 0 (s)η, where * denotes the conjugate transpose. Since F (s) is non singular in C + , there exists M 3 > 0 such that M 3 < η * F * (s)F (s)η.
The function p(s) can only have a finite number of zeros in Ω 2 (isolated zeros theorem). As the quantity δ = min s δ r (s) is strictly positive, equation (60) holds for any δ i ≤ δ . This is a contradiction with the previous inequality. Consequently, there does not exist any s ∈ C + such that p(s) = 0. This implies delay-robust stability as the asymptotic vertical chain of zeros of p(s) can not be the imaginary axis (the principal term of p(s) being stable).

V. SIMULATION RESULTS In this section we illustrate our results with simulations Let us consider the unstable system (3)- [START_REF] Bastin | Stability and boundary stabilization of 1-D hyperbolic systems[END_REF] for which the coefficients are defined by

Λ + = 1, Λ -= 1 0 0 2 , Σ ++ = 0, Σ +-= ( 1
The parameters values are chosen such that • the ODE and the PDE open-loop system are unstable,

• the reflection terms Q and R satisfy Assumption 2. We consider the norm || • || defined by [START_REF] Aarsnes | Delay robust control design of under-actuated PDE-ODE-PDE systems[END_REF]. The initial condition is chosen as a C 1 function. The algorithm we use is adapted from the one proposed in [START_REF] Auriol | Minimum time control of heterodirectional linear coupled hyperbolic PDEs[END_REF]. Using the method of characteristics, we write the integral equations associated to the PDE-system (11)- [START_REF] Krstic | Boundary control of PDEs: A course on backstepping designs[END_REF]. Finally, the original system (3)-( 6) is simulated using a Godunov's discretization scheme. The predictor is adjusted from the one presented in [START_REF] Mondié | Finite spectrum assignment of unstable time-delay systems with a safe implementation[END_REF]. Figure 1 pictures the ||•||-norm of the state (u, v, X) in open-loop and using the control law (52) in presence of a 0.05-s delay.

VI. CONCLUDING REMARKS

In this paper, we have developed a delay-robust stabilizing feedback control law for a system composed of a PDE coupled with an ODE through its boundary. Our approach consists of a first backstepping-based feedback and a second prediction based feedback. This second feedback control is obtained after solving a tracking problem after reformulation of the PDE-ODE system as a neutral system. This mixed strategy is proved to robust to small delays as we preserve the proximal reflection terms in the target PDE systems used for the backstepping design. We will consider in future works the delay-robustness properties of the output-feedback controller.

0 (

 0 ik y k (t -s)ds + p l=1 τ GX ) il (s)X l (t -s, 1)ds -

  Similarly, H 0 (s) and E 0 (s) are bounded in C + , so that there exists M 4 > 0 such thatM 3 ≤ η * H * 0 (I -∆) * E * 0 E 0 (I -∆)H 0 η ≤ max i m|1 -e -δis | 2 M 4 .Construct δ r (s) = δ |s| , for some δ > 0 such that e δ < 1 + M3 mM4 . It follows that for any δ i ≤ δ r (s),|1 -e -δs | ≤ e δ -1 < M 3 mM 4 .

C

  = (0.1 0.1) .

Fig. 1 .

 1 Fig. 1.Evolution of the || • || norm of the system (3)-(6) for the parameters (61) in presence of a 0.05s delay.

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.


Siavash Nejadi

Nasser Kazemi

Jordan A Curkan

Jean Auriol

Paul R Durkin

Stephen M Hubbard

Kristopher A Innanen

Roman J Shor

Ian D Gates

The oil and gas industry, operating and service companies, and academia are actively looking for ways to see ahead of the drillbit while drilling to reduce the risks and costs of the operation and improve the well-placement process. Optimal drilling in the challenging and highly heterogeneous reservoirs, where geological interpretations overlook the high-frequency variations in the rock properties, requires reliable subsurface information from around and ahead of the drillbit. To provide this, we have developed a seismic-while-drilling imaging algorithm based on signal processing, drillstring modeling, and pre-stack wave-equation migration.

To extend the visibility ahead-of-the-bit, we use the drillbit as a seismic source and image the changes in acoustic properties of rocks both around and ahead of the drillbit. The common practice is to build a reverse vertical seismic profile (R-VSP) gathers. Here, we use a blind deconvolution algorithm to estimate the drillbit source signature from the data directly. Alternatively, we can estimate such a signature through drillstring modeling and top-drive measurements (i.e., force and velocity). The drillstring dynamics is modeled by using Riemann's invariants and a backstepping approach. Next, we input the estimated source signature to the pre-stack wave-equation depth imaging workflow. Our simulations show that providing drillbit source signature to the pre-stack wave equation depth migration consistently delivers reliable subsurface images around and ahead of the drillbit. The output of our workflow is a high-resolution subsurface image that provides vital information in oil sands reservoirs for placement of steam assisted gravity drainage (SAGD) well pairs. Compared to conventional practices, the proposed methodology images around and ahead of the drillbit enabling interactive decision making and optimal well-placement. The key feature of the presented methodology is that instead of cross-correlating the seismic-while-drilling data with the pilot trace and building R-VSP gathers, we use the estimated drillbit source signature and deliver high-resolution pre-stack depth migrated images. Through numerical modeling, we tested the potential impacts, validity, and challenges of the proposed methodology in drilling horizontal wells in SAGD settings with an emphasis on the McMurray Formation. We further compared the results with the conventional drilling practice. In contrast to existing tools that have limited depth of penetration, interpreting seismic-while-drilling data in real-time confidently maps key target features ahead of the drillbit. This imaging workflow provides sufficient time to precisely control the borehole trajectory and stay within the desired reservoir zone. Accordingly, it mitigates the risk of intersecting mud-filled channels and lean zones.

Introduction

The Lower Cretaceous McMurray Formation hosts the majority of bitumen in the Athabasca oils sandsthe largest known resource of bitumen. This formation is composed of large-scale point bars and other components of meander-belts that are highly heterogeneous [START_REF] Musial | Sub-surface and outcrop characterization of large tidally influenced point bars of the Cre-taceous McMurray Formation (Alberta, Canada)[END_REF]. The repeated erosional cut and fill events within the McMurray Formation have led to multi-stacked channel-belt deposits (Hein and Cotteril, 2013;[START_REF] Horner | The impact of Aptian glacio-eustacy on the stratigraphic architecture of the Athabasca Oil Sands[END_REF]. Associated point bar deposits are up to 50 m thick, and often exhibit thick packages of massive to cross-stratified sandstones near their base and interbedded sandstones and siltstones upwards, which are referred to as inclined heterolithic strata (IHS) [START_REF] Thomas | Inclined heterolithic stratificationterminology, description, interpretation and significance[END_REF]. More broadly, channel-belt deposits consist of an amalgam of depositional elements, resulting in complex sedimentary facies relationships, where rock types change both laterally and vertically over short distances (e.g., [START_REF] Hubbard | Seismic geomorphology and sedimentology of a tidally influenced river deposit, Lower Cretaceous Athabasca Oil Sands[END_REF]. Steam-Assisted Gravity Drainage (SAGD) is the most promising in-situ recovery process to produce bitumen resources from Athabasca oil-sands deposits (Butler, 1991;[START_REF] Strobl | Geological modelling of McMurray Formation reservoirs based on outcrop and subsurface analogues[END_REF]. SAGD well configurations typically consist of two parallel horizontal wells located at a short distance one above the other. The length of horizontal sections range from 500 to 1000m, and the vertical separation is 5-7m. As the name implies, the key production mechanism in SAGD is gravity. Hence, placing the producing well at high -non-optimal elevations reduces the ultimate recovery factor. Bitumen below the production well cannot be recovered at any significant distance from the well and therefore placement at or near the base of porosity is optimal. Another significant design limitation worth noting is controlling the borehole undulation, with less than 0.5 m unwanted horizontal deviations over its entire length targeted (i.e. 0.5/1000 m). Geological interpretations of the subsurface, together with engineering design, determine the optimal placement of the horizontal wells. In SAGD systems, the optimal field development requires not only an estimation of gross reservoir thickness but also the distribution of non-reservoir facies and connectivity. Fluvial processes related to point bar formation, intra-point bar erosion, counter-point bar formation, and channel abandonment contribute to the development of heterogeneities within the McMurray Formation [START_REF] Thomas | Inclined heterolithic stratificationterminology, description, interpretation and significance[END_REF][START_REF] Hubbard | Seismic geomorphology and sedimentology of a tidally influenced river deposit, Lower Cretaceous Athabasca Oil Sands[END_REF][START_REF] Martinius | Reservoir characterization and multiscale heterogeneity modeling of inclined heterolithic strata for bitumen-production forecasting[END_REF][START_REF] Yan | A 3D forward stratigraphic model of fluvial meander-bend evolution for prediction of point bar lithofacies architecture[END_REF][START_REF] Durkin | Three-dimensional reconstruction of meander-belt evolution, Cretaceous McMurray Formation, Alberta Foreland Basin, Canada[END_REF]. Bedding architecture and lithologic heterogeneities within different geobodies influence the connected hydrocarbon volume and recovery performance. Mapping inclined heterolithic strata (IHS) of sandwiched sand-shale/siltstone sequences and muddy abandoned channels, where possible, improve the understanding of flow barriers and baffles in the reservoir; that aids in both well placement and advanced well completion design for improved recovery. The use of three-dimensional (3D) seismic data for detailed subsurface imaging and interpretation is a common practice (e.g., [START_REF] Hubbard | Seismic geomorphology and sedimentology of a tidally influenced river deposit, Lower Cretaceous Athabasca Oil Sands[END_REF][START_REF] Durkin | Three-dimensional reconstruction of meander-belt evolution, Cretaceous McMurray Formation, Alberta Foreland Basin, Canada[END_REF]. 3D seismic data, a high density of wells and extensive cores are used to infer characteristics of sedimentary architecture, including the geometry of point bar elements. The workflow relies on the combination of data types to construct highresolution geological models of channelized reservoirs, which are the basis of development strategies. The elements of meander-belt deposits at different stratigraphic levels of the reservoir are deterministically modeled and combined with methodologies to infer reservoir rock parameters, such as facies and porosity trends, using seismic attributes [START_REF] Journel | Combining Knowledge from Diverse Sources: An Alternative to Traditional Data Independence Hypotheses[END_REF]. In general, several factors affect the resolution of seismic data. A 360-degree seismic ray path coverage is required to have an ideal and reliable seismic image of the subsurface. The 360-degree coverage guarantees uniform sampling and illumination of the subsurface structure. However, in the field, we only have access to seismic recordings at the surface. The complex stratigraphic architecture, high heterogeneity, low impedance contrast, and sub-seismic features introduce non-uniform illumination and uncertainties in seismic images. In the point-bar deposits of the McMurray Formation, the presence of sub-seismic mudstone drapes potentially impede steam chamber growth [START_REF] Richardson | The Effect of Small, Discontinuous Shales on Oil Recovery[END_REF][START_REF] Hartkamp-Bakker | Permeability patterns in point bar deposits: Tertiary Loranca Basin, central Spain, in the geological modelling of hydrocarbon reservoirs and outcrop analogs[END_REF]. These mudstone drapes cannot be imaged using the surface seismic method. Moreover, the deposits towards the base of the reservoir overlying the sub-Cretaceous unconformity surface are convoluted and show low impedance contrast, which makes detailed imaging of stratigraphic architecture difficult. The small-scale heterogeneities, convoluted layers, and irregular surface topography of the unconformity contribute to stratigraphic complexity and result in a poor seismic response in the region. To reduce the uncertainties inherent in surface seismic recordings, we need to use other sources of information. Conventional SAGD well placement practice includes acquiring data while drilling for optimal placement of the laterals. Measurement-while-drilling (MWD) and logging-while-drilling (LWD) have been adopted to monitor drilling parameters and acquire formation properties, in real-time, and transfer data to the surface. Applications of MWD systems include monitoring and controlling the directionaldrilling operation, assistance in detection of abnormal pressure zones, correlation logging, preliminary evaluation of some possible producing zones, monitoring of weight on bit (WOB), and drilling torque at the bit [START_REF] Dowell | Chapter 15 -Drilling-Data Acquisition[END_REF]. Accurate measurements of depth, inclination, and azimuth are essential for the path determination of horizontal SAGD wells. Among various LWD tools such as Gamma-radiation (GR), nuclear, acoustic, or NMR technologies, LWD resistivity measurements are the principal measurement for geosteering and optimal well placement; the depth of investigation of the LWD resistivity tools is significantly larger than that of any other LWD tools. LWD tools have sensors located immediately above or near the drillbit. Putting these tools at some distance behind the bit still provides an opportunity to first detect, and then avoid drilling into non-reservoir, which is a pseudo-look-ahead capability.

To place the wells in the best possible reservoir section, operators exploit look-ahead monitoring techniques to guide real-time well trajectory adjustment. Industry and academia have researched ways to use noise generated by the drillbit-rock interaction to monitor the rock properties ahead of the drillbit since the 1960s [START_REF] Poletto | Seismic while drilling: Fundamentals of drill-bit seismic for exploration[END_REF]. The drillbit-rock interaction generates significant P and S waves that are strong enough to bounce several times from the subsurface layer boundaries and be recorded at the surface or in the nearby wells. The methodology that uses the seismic waves generated by the drillbit-rock interaction to image the subsurface structure is called seismic-while-drilling (SWD) imaging. Authors also refer to SWD imaging as an inverse vertical seismic profiling (inverse VSP) [START_REF] Hardage | Vertical seismic profiling[END_REF][START_REF] Poletto | Seismic while drilling: Fundamentals of drill-bit seismic for exploration[END_REF]. This is because in VSP recordings, the seismic source is at the surface, and the recording devices, i.e., geophones, are placed in the wells. However, in SWD data, the drillbit-rock interaction in the well acts as a seismic source, and the geophones record drillbit noise at the surface. SWD has some advantages over the LWD and VSP methods. LWD provides images behind the bit; whereas SWD-based algorithms image around and ahead of the bit. Imaging ahead of the bit provides an opportunity to adjust the drilling parameters for improving the rate of penetration (ROP) and updating the well trajectory, which, in turn, ensures staying in the sweet spot, and avoiding risky and nonproducible zones. Also, compared to VSP, the SWD method does not interfere with drilling. Accordingly, the drilling crew can drill wells and record drillbit noise simultaneously [START_REF] Poletto | Seismic while drilling: Fundamentals of drill-bit seismic for exploration[END_REF]. Besides, in SWD, we have the option to record drillbit noise in nearby wells, allowing imaging between wells (i.e., crosswell imaging). Crosswell imaging using SWD gives vital information about the nearby subsurface structures through the development process. In other words, crosswell images improve reservoir characterization and are used to optimize subsequent well pair/pad placement. We use a new SWD imaging algorithm that provides high-resolution subsurface images. [START_REF] Kazemi | Illumination compensation with seismic-while-drilling plus surface seismic imaging[END_REF] show that if we understand drillbit-rock interactions and can model their seismic radiation patterns, prestack depth imaging of SWD data is possible. To estimate the drillbit signature, Kazemi et al. ( 2018) used a blind deconvolution method called SMBD [START_REF] Kazemi | Sparse multichannel blind deconvolution[END_REF]. Recently, [START_REF] Auriol | A Sensing and Computational Framework for Estimating the Seismic Velocities of Rocks Interacting with the Drill bit[END_REF] introduce a methodology that uses drill string modeling and top-drive measurements to model drillbit-rock interaction. These two techniques are used to estimate the drillbit source signature, which is an input to the SWD imaging algorithm developed by [START_REF] Kazemi | Illumination compensation with seismic-while-drilling plus surface seismic imaging[END_REF]. Here, we use this algorithm to image around and ahead of the drillbit. Subsequently, we use these images to update reservoir models and consider applications for improving well placement. The main focus of this study is to present a SWD methodology supported by several case studies inspired from the Surmont Phase 1 SAGD project. We investigate potential applications of look-aheadof the-bit information for an improved well pad drilling approach. The approach helps to reduce drilling risks and improves 1) the geosteering decisions for individual production well drilling including their vertical placement and 2) the knowledge of geological heterogeneity and reservoir boundaries in well pad design. The computational requirements of the methodology would not constrain its implementation for real-time operational decision making. The methodology provides high-resolution ahead-of-the-bit images of the subsurface that identify reservoir boundaries and highlight geological heterogeneities to be used for real-time decision making.

Methodology

Drillbit-rock interaction generates significant elastic energy in the subsurface. The drillbit energy propagates through or reflects from subsurface layers, and can commonly be recorded at the surface or nearby wells. These wave mods follow the wave equation [START_REF] Kazemi | Illumination compensation with seismic-while-drilling plus surface seismic imaging[END_REF] where P is pressure wavefield, s is slowness, ω is temporal frequency, f is the drillbit source signature, x s is source location, is Laplacian operator, and δ is the Dirac delta function. Solving the wave equation represented in Equation-1 provides the information necessary to generate a subsurface image. To solve Equation-1, we decompose the subsurface slowness and recorded pressure wavefield as follows

(2) where s 0 is the background slowness, m is the subsurface image, P 0 is the pressure wavefield satisfying the background medium, and ∆P is the perturbation in the pressure wavefield due to m. Then, by using the Born approximation (Born 1926), we get the subsurface image as

(3)
where d is the recorded SWD data at the surface, G 0 is the Green's function that satisfies the wave equation in the background medium, x ' s is source location, x ' r is receiver location, and P 0 is the sourceside wavefield. To generate the source-side wavefield, the drillbit source signature is estimated through the exploration of two methodologies. The first method uses a signal processing algorithm and directly estimates the drillbit source signature from the recorded SWD data. In this approach, we use the blind deconvolution method proposed by [START_REF] Kazemi | Sparse multichannel blind deconvolution[END_REF] to estimate the source signature directly from the data (i.e., SMBD algorithm). In the SMBD algorithm, data is modeled as the convolution of drillbit source signature with a reflectivity series (4)

where F is the convolution matrix built from the drillbit source signature f, r j and n j are the reflectivity series and noise term in channel j. [START_REF] Kazemi | Sparse multichannel blind deconvolution[END_REF] found the smallest common factor between the channels by solving an inverse problem where the solution fits the heterogeneous system of equations and provides a sparse reflectivity series. Then, the byproduct of the algorithm is the drillbit source signature. The SMBD algorithm estimates the smallest common factor between all of the channels in the SWD data by solving a heterogeneous system of equations. Let's convolve the data of channel p with the reflectivity of channel q and the data of channel q with the reflectivity of channel p. After some algebraic manipulations, we get (5) where D p and N p are the convolutional matrices built from d p and n p , respectively. After using all of the combinations similar to Equation-5 between all of the channels, we end up having a heterogeneous system of equations ( 6) where ( 7) and ( 8)

The term e is the concatenation of vectors similar to the right-hand side of Equation-6. To find the reflectivity series that fit the heterogeneous system of equations represented in Equation-6, SMBD solves the following cost function (9) where λ is a regularization parameter, and ε is a small number. The regularization term in the cost function is the Huber norm that promotes sparsity in the final solution. After finding the sparse reflectivity series, the drillbit source signature is estimated by a least-squares matching filter in the frequency domain that matches the data to the reflectivity series. Interested readers are referred to [START_REF] Kazemi | Illumination compensation with seismic-while-drilling plus surface seismic imaging[END_REF] for the details of the application of the SMBD algorithm in estimating the drillbit source signature and SWD imaging. The second approach for estimating the drillbit source signature relies on drill string modeling and topdrive measurements [START_REF] Auriol | A Sensing and Computational Framework for Estimating the Seismic Velocities of Rocks Interacting with the Drill bit[END_REF]. These authors showed that through modeling the drillstring dynamics and measuring the top-drive force and velocity time series, it is possible to estimate the drillbit source signature. They also use the estimated drillbit source signature in combination with the SWD data to estimate the velocity of rocks that are interacting with the drill. In this paper, we assume that the drillbit generates significant energy that can reflect several times from the layer boundaries and can be recorded at the surface or nearby wells. However, in situations where the drillbit does not generate significant noise, we should consider using a controlled source in the well. In both cases, i.e., using dillbit as a source or a controlled source, if the elastic energy reaches the surface or nearby wells, our proposed workflow is valid and provides high-resolution subsurface images.

Case Studies

This section presents different case studies, where the look-ahead-while-drilling could potentially improve SAGD well placement. The main objective is to assess value creation as a result of using SWD to sequential geosteering decisions. We illustrate the application of the proposed methodology in various geological settings inspired by Surmont Phase 1 developments. The optimal SAGD development requires not only an estimation of the gross oil-sand reservoir thickness, but also internal imaging of the reservoir, including the distribution of non-reservoir facies. Current best practice for SAGD well pad placement is to place the production well near the base of the pay zone. Development areas are selected according to various criteria including: net continuous bitumen; top gas and water, lean zones; cap rock thickness; and surface constraints including proximity to the central processing facility. In our study area, the base of the pay zone is the base of high porosity bitumen saturated sandstone facies, which is below the seismically-constrained youngest meander deposit. Typically, petrophysical well logs show that porosity and bitumen saturation is highest at these depths; in the absence of vertical flow barriers and baffles, it is a sweet spot for well pair placement. However, the lack of lithologic contrast and resolution of seismic data limits the interpretation of depositional elements at these depths. The amalgamation of sandy meander belt deposits makes correlations difficult and adds to the complexity of reservoir characterization. Mapping mudstone-filled abandoned channels and muddy IHS areas of the lower channel belt deposit(s) are critical to horizontal well pair placement. These impermeable elements block vertical steam chamber growth and limit depletion from the thick, overlying bitumen-saturated point bar deposits. Similarly, drilling and placing the wells inside these poor quality elements affect productivity or injectivity, leading to a partial loss of the drilled lateral length. Wireline well logs, core descriptions, and stratigraphic dip analysis are used to infer the distribution of channel elements in the underlying deposits [START_REF] Fustic | Dipmeter Application in Canadian Oil Sands Development -Evolution from Predicting Potential Block Slides to 3D Geological Interpretation and Modelling[END_REF]. This limited source of data incurs uncertainty in interpreting depositional elements and spatial distribution of reservoir vs. nonreservoir facies. Below the base of the meander-belt strata and above the sub-Cretaceous unconformity, the McMurray Formation is characterized by smaller upward fining fluvial channel deposits that are encased in finegrained deposits including floodplain, paleosol and lacustrine units. These smaller channels -less than 10 m thick and 100-150m wide -are related to relatively localized paleodrainages that incised Devonian strata. Imaging this amalgam of non-reservoir units -between the base of porosity and top of unconformity -is equally crucial to optimal well placement and field development planning.

Case A: Producer partially drilled in non-reservoir facies, reservoir boundary

This case examines the downsides of placing a horizontal producer in a suboptimal location. The producer is placed near the base of the pay zone to maximize recovery, but is partially drilled into a mudstone-filled channel associated with the lower meander-belt deposit. The proximity of the mudstone-filled channel to the highly reflective Devonian Unconformity challenges the precise interpretation and estimation of the three-dimensional distribution of this non-reservoir rock. Figure-3 shows a cross-section that encounters a horizontal well pair. Different colors highlight the inner architectural elements of the reservoir. The log example alongside both injector and producer show typical Gamma-radiation log responses for common sedimentary rock types; cool colors correspond to the high GR values. The presented model is constructed for the most likely interpretation of the reservoir, as this is thought to be the most likely outcome and most useful and realistic for field development purposes. For optimal field development planning, we recommend performing sensitivity studies and evaluating the uncertainty in both deterministic and stochastic modeling, although such studies are beyond the scope of this paper. The length of the producer is 860 m. Approximately the last 140 m of the lateral section shows a high GR value, highlighting suboptimal well placement and drilling of non-reservoir facies (mudstone). The corresponding injector is drilled according to the producible length of the producer and is 140 m shorter than the designed length and pad average. Figure-4A displays the production data for this well pair. The average monthly oil production is roughly 1,400 m 3 with a cumulative steam oil ratio (cSOR) of 3.6 m 3 /m 3 . This well pair has a lower oil production rate and a higher cSOR compared to other well pairs in the same pad, which is attributable to the shorter producible length and the smaller drainage volume as a result of proximity to the reservoir boundary (mudstone-filled channel or edge of the belts). Table-1 compares the production data with other case studies; the data for the pad average is also included as the baseline to highlight the difference between superior and low performing cases. Figure -5 shows the overlay of subsurface velocity profile and high-resolution SWD image corresponding to the cross-section represented in Figure -3. As clear, the mudstone boundary is remarkably imaged with the proposed SWD imaging algorithm. To generate the SWD image, we numerically simulated 15 equally spaced drillbit sources in the vertical section of the producer's well. Moreover, we considered that there is a vertical well on the right-hand side of the cross-section. Hence, we recorded the numerically simulated pressure wavefields at the right-hand side vertical well. By doing so, we show that the subsurface image of cross well SWD data could efficiently reveal the mudstone and reservoir boundary. Accordingly, the high-resolution SWD image is a vital source of information for optimal well placement and avoiding risky zones. In this case, where the horizontal producer reaches and penetrates the reservoir boundary, the optimal geosteering solution, in practice, is not steering the horizontal well above the non-reservoir zone. Whereas, SWD interpretations provide improved images of the boundary well in advance (Figure -5), with sufficient time to make a decision and steer the well to the sides or as the best choice, not to drill beyond a certain measured depth -in the mud. In practice, this application saves rig time and minimizes drilling costs; incurred expenses for inefficient 140 m dilling cover the cumulative costs of seismic-while-drilling.

Case B: Optimal drilling and well placement

This case presents the optimal well placement with little or no loss of producible length of the horizontal wells (Figure -6). The producer is placed near the base of porosity and both injector and producer are drilled in reservoir facies with good porosity and high bitumen saturation. The production performance of this well pair is well above the pad average and is one of the best well pairs in the pad (Figure -4B). The average monthly oil production is 3,600 m 3 with cSOR of 2.8 m 3 /m 3 (Table -1). It is interesting to highlight that the very last few meters (20 m) of the injector penetrates an element of the upper meander-belt. This element is a separate, older unit not related to the modeled point bar deposit, with high IHS and mudstone facies proportions. As a result, the GR response near the end of this injector is characterized by high values, which correspond to non-reservoir facies. Again, this demonstrates the need for detailed reservoir characterization and precise interpretation of point bar elements including the three-dimensional geometry and fill of channel bodies prior to field development planning in SAGD settings.

Case C: Producer is partially drilled in Paleosols/Basal Shales In this case, we discuss a scenario where portions of the horizontal producer is placed in non-reservoir rocks (Figure -7). Consider non-uniform 3D seismic resolution and a pad size of 1000 x 1000 m 2 ; Velocity (m/s)

Horizontal location [~950m]

Depth location [~85m] typically 10 vertical cored or logged stratigraphic wells are drilled in such pads. This leads to an average vertical well spacing (sampling) of 10 hectares/well. That is sufficient to map flat structures like the base meander-belt, but introduces significant uncertainty in 1) interpreting depositional elements and 2) estimating reservoir boundaries, as defined by the edge of the belts. As mentioned, optimal horizontal well placement in a SAGD setting necessitates placing the producer well as close as possible above the base of bitumen. A detailed map (<1 m uncertainty) of the highs and lows on the undulatory unconformity surface, facies variations, as well as basal water distribution is needed. In this example, because of poor vertical well sampling and low-resolution 3D seismic data near the base of the reservoir, the vertical position of the producer is suboptimal with approximately the first half of the producer is drilled in mudstone and IHS strata. This has significantly affected the production performance of the well pair. The average monthly oil production of this well pair is 2,100 m 3 and cSOR is 4.1 m 3 /m 3 , which is quite poor compared to the pad average (Figure -4C). The real-time SWD data potentially images the surface on top of non-reservoir basal shales, and it provides valuable information to support decision making, steering the well, and adjusting the landing point of the horizontal section.

Added value There have been limited reports on the application of SWD in SAGD developments. The excessive costs of acquisition and interpretation need to be clearly justified in order to convince operators of its benefits. Several reports (e.g. [START_REF] Chemali | Advanced Geosteering for Optimal Exploitation of Hydrocarbon Reserves[END_REF] clarify the benefits of SWD and its successful application in deep water and subsalt drilling environments. A discussion of missed or wasted resources during SAGD developments, considering both the cost of data acquisition and time spent on design and data interpretation, is not available at this time. In our case studies, we presented certain environments where look-ahead-of-the-bit information could drastically improve SAGD well placement, and ultimately, well deliverability. The production performance of adjacent well pairs shows the added value of optimal well placement and the potentials of SWD data.

A key design criterion in SAGD well placement is the target depth of the producer's horizontal track. Generally, horizontal wells are drilled using sparse vertical well sampling and surface seismic data. These sources of data have quite large uncertainty near the base of the pay zone, and as a result, laterals may not be optimally placed with regards to target petrophysical properties (Case C), and in some cases, even miss the desired target altogether (Case A). SWD can help to reduce this uncertainty by providing reflective information in real-time. Operators can update subsurface images and lessen the risks while the well is being drilled into the target zone.

Ahead-of-the-bit information helps geosteering engineers to mitigate the risks of drilling in unwanted zones without disturbing drilling operations. It provides indications of reflecting horizons/rocks and helps to detect non-reservoir facies. This supports both drilling and well completion design decisions. Ahead-of-the-bit information is highly valuable to optimize casing landing point for horizontal SAGD wells placement. Besides, the detailed and high-quality rock type information provides invaluable constraints to enhance engineering designs and optimize completion equipment (e.g. [START_REF] Nejadi | Optimization of Placement of Flow Control Devices under Geological Uncertainty in Steam Assisted Gravity Drainage[END_REF]. In the case of non-challenging SAGD development areas with flat structures near the base, seismicwhile-drilling information will provide detailed images of the inner elements of the deposit, resulting in reduced uncertainties in reservoir characterization step. The operator may use detailed geological models to fine-tune the well pair orientations in the future pads and place the pads perpendicular to the point bar packages -accretion surfaces and maximize productivity [START_REF] Su | SAGD Pad performance in a point bar deposit with a thick sandy base[END_REF].

Multi-Well Pad Drilling

Multi-well pad drilling reduces both time and costs by allowing for shared construction, equipment, and facilities, reducing the surface footprint while minimizing rig moving expenses. It further provides an opportunity for continuous process improvement. In this paper, we propose improvements to well placement strategies that can enhance production efficiency in the context of multi-well pad drilling. This approach fine-tunes the placement of horizontal wells with the help of drilling and subsurface data obtained in real-time while drilling previous laterals. A result may be that the sequence of lateral drilling may be modified as data is collected through the batch drilling process. The sequence starts with drilling the first lateral in the best-characterized, less uncertain part of the pad that inherently implements the best knowledge of the subsurface. Subsurface images are updated while drilling the horizontal wells. The real-time information reduces uncertainties and updates the current understanding of the presence, scale, and continuity of low-quality facies. This approach follows the idea of just-in-time well placement and design, which is to specify the next horizontal well, in the drilling sequence, according to the risks and uncertainties, instead of a pre-defined order using a priori knowledge. The selection criteria and placement of the next laterals are modified in real-time using the lessons learned from the previous drilling. Our proposed workflow challenges the practice of merely skidding the rig over a few meters and drilling the next lateral in a pre-designed manner that solely relies on sparse vertical well sampling and three-dimensional seismic with non-uniform resolution in the depths of lateral placement.

Conclusions

This paper discussed the value of look-ahead-of-the-bit information to multi-well pad drilling in SAGD settings. We use a novel SWD imaging algorithm that provides high-resolution subsurface images. Imaging ahead of the bit provides an opportunity to adjust the drilling parameters, updating the well trajectory and improves 1) well placement and 2) drilling performance. The approach consists of estimating the drillbit source signature and using it as an input to the SWD imaging algorithm. Later on, we use these images to update reservoir models and improve well placement. Implementation of SWD information in the context of SAGD well pad drilling in the presence of poor three-dimensional seismic quality and geological uncertainty demonstrates reasonable improvement in optimal well placement. Application of SWD information potentially addresses the problem inherent in horizontal producer placement, and minimizes risks of drilling in non-reservoir facies; it implements an improved image of the subsurface and reduces the uncertainties in reservoir size, distances to the boundary, and formation properties. In addition to improved well placement, the near-real-time updated subsurface images provide sufficient time to fine-tune downhole well completion designs.

Nomencluture

Figure 1 -

 1 Figure 1 -A) Paleogeographic reconstruction of the Early Cretaceous Western Canada Sedimentary Basin (WCSB), including the location of the ancient meander belt deposit of this study. B) Stratigraphic setting of the McMurray Formation in the lower Mannville Group. C) Simplified stratigraphic column for the McMurray Formation, which unconformably overlies Devonian carbonate strata in the study area and Gamma-radiation (GR) log through the entire McMurray Formation. (Modified from Hubbard et al., 2011, and Durkin et al., 2017)

Figure 2 -

 2 Figure 2 -A) 3D seismic-amplitude stratal slice and B) zone model of the seismically constrained youngest meander-belt deposit in the focused study area. (Modified from Durkin et al., 2017)

Figure 3 -

 3 Figure 3 -Vertical reservoir cross section along Case A well pair. Different colors show the zone model and and the log show Gamma-radiation along the laterals. The horizontal scale is 1:2800, and the vertical exaggeration is 4.5.

Figure 4 -

 4 Figure 4 -Oil production and cumulative steam oil ratio for different case studies.

Figure 5 -

 5 Figure 5 -Overlay of subsurface velocity profile and high-resolution SWD image corresponding to the cross-section represented in Figure-3. Black arrows show the mudstone boundary.

Figure 6 -

 6 Figure 6 -Vertical reservoir cross section along case B well pair. Different colors show the zone model and the logs show Gamma-radiation along the laterals. The horizontal scale is 1:2800, and the vertical exaggeration is 4.5.

Figure 7 -

 7 Figure 7 -Vertical reservoir cross section along case C well pair. Different colors show the zone model and and the log show Gamma-radiation along the laterals. The horizontal scale is 1:2800, and the vertical exaggeration is 4.5.

  Figure-8 shows a seismic depth section of the meander belt in the study area. As discussed and illustrated in different case studies, in the area used for this study, the horizontal well pairs are often placed between the base of the upper meander-belt deposit and top Devonian strata. Seismic profiles at these depths show non-uniform resolution rendering geological interpretations, yielding a challenging task. Due to the strong reflectance at the contact between poorly consolidated sediments of the McMurray Formation and competent underlying Devonian limestones, the surface seismic resolution is relatively poor up to 10 m above the unconformity.

Figure 8 -

 8 Figure 8 -A) Seismic-time cross section of the McMurray Formation in the study area. B) Line trace showing the elements of meander-belt deposits at different stratigraphic levels. 3D seismic resolution is relatively poor up to 10 m above the unconformity. (Modified from Durkin et al., 2017)

  term in SMBD f: Seismic source function m: subsurface model n: Noise in the data r: Subsurface reflectivity s: Slowness A: Channels combination matrix for SMBD

Acknowledgements

The authors would like to acknowledge the University of Calgary's Canada First Research Excellence Fund (CFREF) program on unconventional resources for funding support and ConocoPhillips Canada for allowing use of the Surmont geomodel. PETREL software is provided by Schlumberger.