Reduced basis method for managed pressure drilling based on a model with local nonlinearities

To circumvent restrictions of conventional drilling methods, such as slow control actions and inability to drill depleted reservoirs, a drilling method called managed pressure drilling (MPD) has been developed. In MPD, single‐phase flow processes can be modeled as a feedback interconnection of a high‐order linear system and a low‐order nonlinear system. These nonlinearities appear locally both inside and at the boundaries of the computational domain. To obtain a fast simulation platform for real‐time purposes (eg, online model‐based controller implementation), model order reduction is required for MPD. However, the local nonlinearities render applying model order reduction techniques challenging. In this study, a new approach is proposed to deal with such nonlinearities within the reduced basis (RB) context and it is successfully tested on a model for MPD. Contrary to the classical RB technique, the proposed approach not only does not generate nonphysical spikes at the locations of these local nonlinearities but also yields high speedup factors. The obtained reduced‐order model can be used for efficient online simulation and controller design for drilling systems with MPD.


INTRODUCTION
The design of pressure control systems for drilling with Managed Pressure Drilling (MPD) requires an accurate model of the hydraulics of the system. Moreover, training novice drillers for drilling operations, simulating well control incidents in a controlled environment and also monitoring drilling wells are other important reasons for hydraulics modeling in drilling. However, the desire for accurate modeling of drilling systems typically leads to highly complex models. In particular, this type of modeling for MPD automation gives rise to models involving parameterized, nonlinear, nonconservative hyperbolic Partial Differential Equations (PDEs) completed by nonlinear and implicit boundary conditions. [1][2][3] The governing equations can become even more complicated when the flow path cross-section is discontinuous along the well. All these features of the model render its numerical simulation computationally expensive and also make the controller design for the system cumbersome. Moreover, an accurate and fast assessment of the downhole pressure is crucial in drilling processes. A tool that facilitates this assessment can be a significant help/support for the drilling procedure. For instance, in case of emergencies, an easy-to-use and much faster than real-time simulation framework can provide insights into effects of a potential reaction to contingencies occurring during drilling. This reduces the probability of making hazardous decisions and hence increases the safety of the drilling rig. To provide such efficient model-based simulation tool, model order reduction (MOR) 4 can be very useful.
MOR techniques aim at the automatic construction of reduced-complexity models that combine high predictive capacity and low complexity. These techniques drastically reduce the size of the problem and, subsequently, its computational cost at a price of a minimal and quantifiable loss of model accuracy. The resulting low-complexity models can then be used for simulation, control, optimization, parameter estimation, and inverse modeling. 4 Among several existing techniques for reducing the dimension of expensive and complex computational models, the reduced basis (RB) 5 method is an efficient approach for dealing with parameterized systems. In particular, the MPD model described in Reference 2 can be considered as a parameterized system where the (varying) parameters can be the properties of the drilling mud or the wellbore geometry.
The use of the RB method for MPD models is not straightforward due to the time-dependent, nonlinear, and state-dependent boundary conditions. The underlying idea of the RB method is to find the solution of the system as a combination of some precomputed functions, called RB functions. Handling boundary conditions becomes challenging due to the global nature of the RB functions and their independency over time. 6 To resolve this issue, a lifting method 7 is implemented here and tested for an industry-relevant MPD application. Within this lifting method, the effect of the boundary conditions is segregated from the internal domain and is incorporated later as a state dependent and locally acting input.
Besides the nonlinearities at the boundaries, local nonlinearities can occur inside the computational domain. One realistic example is represented by the discontinuities in the cross-section of the flow path in a drilling system. 8 These local nonlinearities generate discontinuities in the solution of the state variables. While various reduction techniques for distributed nonlinearities, such as Empirical Interpolation Method (EIM) 9 and Discrete EIM, 10 have been developed, little attention, to the best of our knowledge, has been paid to the issue of local nonlinearities in the RB setting. To take into account the local nonlinearities inside the domain, the RB functions are enriched with some local basis functions with compact support (henceforth called the localized enrichment). This method shares similarities with the approaches introduced in the balanced truncation community. 1,11 In this paper, to the best of our knowledge, for the first time the RB method is applied to a drilling system with MPD and the ansatz introduced in Reference 7 is tested for dealing with implicit and nonlinear boundary conditions governing the hydraulics dynamics of the drilling system. Moreover, a new methodology is introduced to automatically capture the local nonlinearities inside the computational domain by enriching RB functions with specific local basis functions.
The outline of this paper is as follows. In Section 2, a model for MPD is briefly introduced. In Section 3, the new approach to use the RB method for capturing local nonlinearities in MPD is proposed. In Section 4, numerical results comparing the full-order model and the reduced-order model are presented. Finally, conclusions and future works are presented in Section 5.

PROBLEM STATEMENT
The industrial problem under investigation is a drilling system, with a special focus on MPD. The configuration of the system is shown in Figure 1. A drilling liquid known as mud is pumped into a pipe at high pressure; this pipe is called the drillstring. At the bottom of the drillstring, the mud leaves the drillstring through nozzles located inside the drill bit and enters the area between the drillstring and the wellbore known as the annulus. It then flows up through the annulus and carries the rock cuttings out of the well. In MPD, the annulus is sealed off from the surroundings at the top with a Rotating Control Device and the mud circulates out of the well through a choke valve. The circulation path of the mud can be observed by following the green arrows in Figure 1. Furthermore, the diameters of the drillstring and the wellbore experience sudden variations along the well; the diameter profile of the drillstring and the wellbore are shown in red (denoted by d) and blue lines (denoted by D out ), respectively, in Figure 1. The parameters L and denote the well length and the well inclination with respect to the horizontal line. In this paper, a single-phase flow model for MPD is studied. For a more comprehensive explanation of MPD systems and multiphase flow modeling, the reader is referred to Reference 2. The mathematical formulation employed to model this system is described below.
F I G U R E 1 A schematic configuration of an managed pressure drilling system

Mathematical modeling for MPD
To model the drilling hydraulics, the drillstring and annulus can be treated as two pipes connected at the bit and the hydraulics within MPD can be modeled by the so called U-tube modeling approach. 2 Different models have been introduced for the mathematical representation of the single-phase flow inside each pipe. 2,12 One of the most widely used models is obtained by simplifying the isothermal Euler equations and accounting for the area variation of the drillstring and the wellbore. This model is described by the following system of PDEs: where (t, x), v(t, x), and A(x) are density, velocity, and the cross-sectional area of each pipe, respectively. Here, t and x denote the temporal variable and the spatial variable, respectively. Parameters T and L, respectively, denote the temporal and spatial length of the computational domain. In addition, is the momentum of the flow, and superscript (⋅) T shows the transpose operator). The source term S(u, x) consists of a friction term F(u, x), a gravitation term G(u, x) and the effect of variable cross-sectional area. The flux function f(u, x) attributed to (1) is f (u, x) = [ Av, Ac 2 l ] T where c l is the speed of sound in the fluid medium. In case of laminar flow, we have where D h (x) is the hydraulic diameter of each pipe, 0s is the reference density and is the viscosity of the liquid. For gravity, we set where g is the gravitational acceleration and is the inclination of the pipe. To calculate pressure, the equation of state for the fluid can be used, which reads as follows: where p(t, x) is the pressure and p 0 is the reference pressure. Although system (1) is parameter-dependent (eg, the inclination and the properties of the fluid such as c l and 0s can vary), the parameter dependency is not shown for the sake of simplicity. The parameter dependency is made explicit later in Section 3.
Remark 1. System (1) is solved once for the drillstring and once for the annulus. These two solutions then interact at the bit. For the simulation of the hydraulics in the drillstring, we set D h (x) = d(x) where x = 0 and x = L correspond, respectively, to the pump and bit location in the drillstring. Additionally, for gravitational source term, we choose "+" sign in (3). For the annulus, we set D h (x) = D out (x) − d(x) with x = 0 indicating the bit location and x = L indicating the choke location with "−" sign in the gravity term (3).
Remark 2. Area discontinuities, such as those presented in Figure 1, leads to an impulse in the right-hand side of the momentum equation in (1), leading to a discontinuity in the density and velocity. At these locations, where v is discontinuous and v∕ x is an impulse, the simplified isothermal Euler equations as in (1) are not accurate. At these locations, the original isothermal Euler equations as in Reference 13 are exploited to solve the fluid flow.
The finite-volume (FV) method is usually employed to numerically solve this type of system of PDEs. 14 In the following section, the numerical framework to solve system (1) is discussed.

FV discretization
To numerically solve (1), a Gudonov-type scheme is employed as follows where the discrete field variable U n i = [ n i , n i v n i ] T , i ∈ {1, · · · N }, n ∈ {1, · · · N t }, is the average of field variables u = [ , v] T over the ith grid cell at the time instant t n ∶= nΔt. Moreover, N and N t are, respectively, the number of spatial and temporal grid cells of the discretized computational domain. Also,  (⋅, ⋅) is the scheme-specific numerical flux function (to be introduced in Section 2.4). Here, Δt and Δx refer to the temporal and spatial discretization step sizes, respectively. It should be noted that parameter dependency is not mentioned in (5) for the sake of notational simplicity.
The effects of area variation are incorporated by the starred variables (U n, * i+1 , U n, * i−1 ) appearing in (5). These variables can be computed by solving the following system of equations at each time step (refer to Reference 15 for details), After obtaining n, * i+1 and v n, * i+1 , U n, * i+1 = [ n, * i+1 , n, * i+1 v n, * i+1 ] T is constructed at each time step n (U n, * i−1 is formed in the same manner). Since this approach is developed only to consider the effects of diameter discontinuity, system (6) is solved only at the locations of area variations. At other locations in the spatial domain where the area does not change over the interfaces, we set U n, * i+1 = U n i+1 and U n, * i−1 = U n i−1 . Remark 3. The set of equations (6) are derived according to the isothermal Euler equations, not (1), see Remark 2.

Initial and boundary conditions
In order to be able to numerically solve system (1), initial conditions and boundary conditions should be imposed. Initial conditions are chosen to the steady-state solution of system (1) corresponding to a selected set of parameters. Then, by changing the inputs of the system (ie, the boundary inputs), the dynamics of the system are excited. System (1) consists of two first-order PDEs; hence, two physical boundary conditions for each pipe should be defined. For the drillstring, the governing equation of a pump comprises the left boundary condition while the drilling bit equation governs the right boundary condition. For the annulus, on the other hand, the drilling bit constitutes the left boundary, while the right boundary is specified by a choke valve, see Figure 1. The governing equations of the pump, bit and choke are summarized in Table 1. There, q p , A p , and v p represent the volumetric flow rate of the pump, the cross-sectional area at the pump, and the velocity of the liquid at the pump, respectively. Moreover, Δp b , b ,ṁ b , A N , and C D denote the pressure

Pump
Bit Choke TA B L E 1 Governing equations of the pump, bit, and choke drop over the bit, density at the drillstring side of the bit, the mass flow rate through the bit, the total area of the nozzle holes of the bit, and the nozzle coefficient of the bit, respectively. Finally, q c , K c , z c , c , p c , and p atm are the volumetric flow rate through the choke, the choke constant, the choke opening, the density at the choke inlet, the pressure at the choke inlet, and atmospheric pressure, respectively. Control inputs of the system are q p and z c . For a more comprehensive description of the boundary conditions, refer to Reference 2.
Remark 4. It should be noted that a nonreturn valve (NRV) is always installed above the drilling bit inside the annulus. This valve allows the drilling mud to flow only from the drillstring into annulus, not vice versa. The NRV model is introduced in Reference 2.

Upwind scheme
All numerical test cases in this paper are performed by an upwind scheme. 16 The corresponding numerical flux function in (5) holds where By inserting (7) into (5), the following equation is attained by compressing the notation (F n i ∶= F and similarly for G n i ): Using (9), one can find the vector of discrete field variables over the ith cell at each time step. Augmenting all U n i , i ∈ {1, · · · N }, n ∈ {1, · · · N t }, over all cells and taking into account the variation of the diameter (by solving (6)) yields the following algebraic equations where , m, D h ∈ R N are, respectively, the vector-valued averages of the density and momentum and the vector containing diameter values at all grid cells gathered in a vector of dimension N (number of grid cells). Matrix I is the square identity matrix of dimension N . The function "diag" in (11) generates a square matrix with the diagonal terms according to its argument. Furthermore, n 0 and m n 0 are the left boundary values of each pipe at time instant with index n while n N +1 and m n N +1 are the right boundary values of each pipe at time instant with index n. We recall that in the term 32 ∕( 0s D h 2 ) in (11), we set D h = d for the drillstring and D h = D out − d for the annulus (d is the discrete version of d(x) and similarly for D out ). Also, we have In (10) and (11), the variables with the superscript "+", (⋅) +,n ∈ R 2N d with N d number of diameter discontinuities along the spatial domain, are defined as: where the variable (⋅) n is subtracted from the starred variable (⋅) * ,n computed by (6). The operators L * 1 and L * 2 incorporate the nonlinearities due to discontinuities in the diameter. These operators are tridiagonal matrices with zero diagonal terms and nonzero entities at the off-diagonal entities corresponding to the locations of area discontinuities. These nonzero terms are, respectively, a copy of the off-diagonal terms of L 1 and L 2 at the neighboring cells of the interface where the area discontinuity occurs; the other off-diagonal entities are set to zero. The operator ∈ R N ×2N d in (10) and (11) is composed of the canonical vectors at the neighboring cells of the interface with the area discontinuity. For instance, if the area variation only occurs at the interface x i+1/2 of the discretized domain, we have: where e i ∈ R N is the ith canonical vector (a zero vector with only one nonzero element at the i-th entry).
To show the internal dynamics (10) and (11), boundary equations and the area discontinuity solver (6) in a compact form, a description of the discretized system dynamics is shown below in the form of a feedback interconnection: where h(⋅, ⋅, ⋅) is the solver of the dynamics occurring at the boundaries, u c contains the control inputs acting at the boundaries (in MPD, u c = [q p , z c ]), and (⋅) is the solver of (6). We have set The other matrices such as A 11 are defined according to (10) and (11) and y n ∈ R N y is a vector of the outputs of interest with the dimension N y . The other outputs of Σ lin , y n V and y n W (correspondingly the output matrices C V and C W ), are, respectively, required for the computations of the boundary conditions and the area discontinuities. This interconnection is shown in Figure 2.
It should be noted that the area variation changes the wave reflection pattern inside the system. It also significantly changes the velocity profile and therefore the pressure drop due to the friction. The formulation (15) allows to capture these effects even when the order of the system is reduced by the RB method. Remark 5. Before going to the MOR part, let us add the parameter dependency into (15), where is a vector in Here,  is a continuous space containing varying parameters of the system. To summarize, system (16) is linear except at the locations of area discontinuities and at the boundaries. The high-dimensional nature of the system dynamics Σ is only present in Σ lin and not in the nonlinear part Σ nl and this motivates the system complexity reduction by reducing only the linear part. Thus, developing a MOR method to reduce the linear part and simultaneously account for this type of local nonlinearities is essential, which is the topic of the next section. We in particular consider the case in which the computational effort to solve for the dynamics of Σ lin is significantly higher than that required for Σ nl . We care to stress that this is the case in many application scenarios in which Σ lin expresses the discretized, large-scale, dynamics of an internal domain and Σ nl expresses nonlinear, though low-dimensional, dynamics (or even static relationships) characterizing local nonlinearities. In such relevant scenarios, substantial computational efficiency gains can be expected after reduction. The reduction approach can still be pursued if this assumption is violated, but the computational gain may be less significant. In support of the MOR paradigm, we adopt the following assumptions. Assumption 1. To enable efficient reduction of (16), all operators in Σ lin ( ) should be affine with respect to the parameters , that is, where Θ j ( ) are scalar functions of and A j 11 are matrices independent of the parameter . Remark 6. If Assumption 1 is not satisfied, the operators in (16) can be approximated by an affine expression by using EIM. 9

REDUCED ORDER MODEL FORMULATION
The RB method targets parameterized problems requiring repeated evaluations (many-query analysis) or (faster than) real-time simulations (real-time analysis). It typically consists of an offline and an online stage. During the offline stage, an RB space Φ is generated and several parameter-independent operators related to system (16) are evaluated to be used in the online stage. In the online stage, the reduced solution can be obtained much faster than the one with a classical (eg, FV) numerical simulation with a reasonable accuracy.
In this paper, we propose an additional stage, called "intermediate" stage, necessary to sufficiently handle the internal nonlinearities of system (16) (ie, nonlinearities due to the discontinuous diameter of the drillstring and the annulus). In this section, we explain the three mentioned stages of the method: offline, intermediate, and online stage. The requirements of each stage are explained within each section. Finally, the computational decomposition of these stages is illustrated.
Remark 7. In this study, during the offline stage, we propose to consider a reference domain (constant cross-sectional area along the well) excluding nonlinearities in the internal domain and including the discontinuous cross-sectional area in the intermediate stage to be used in the online stage, which is further explained in the sequel.

Offline stage
During the offline stage, a set of N RB functions is computed. Although N is typically much smaller than the degrees of freedom of classical numerical scheme (eg, FV), this stage can be computationally expensive, since the RB functions are obtained via proper orthogonal decomposition (POD) on the solutions of (16) for a specific set of parameter values. The selection of these parameter values is performed by the Greedy algorithm. The entire procedure is called POD-Greedy 5 and is explained in Section 3.1.3.

Proper orthogonal decomposition
Let assume that the solution of (16) for an arbitrary parameter value has been computed. POD extracts the most energetic modes of the snapshots which represent the set of RB functions denoted by . Mathematically speaking, POD solves a minimization problem as follows 17 where U = [U 0 , U 1 , · · · , U N t ] with U n ∈ R N is the numerical solution at all spatial grid points at time instant n,  is the Hilbert space and ‖⋅‖ F is the Frobenius norm. The POD algorithm (to be used in a Greedy fashion later) is summarized in Algorithm 1, where is computed based on the discrete field variables U and the number of desired RB vectors. In Algorithm 1, the number of POD modes is either given by the user or determined by the decay of singular values. As we use this algorithm in a Greedy fashion, we show POD algorithm for the former option. Before discussing the POD-Greedy approach, we discuss the incorporation of the nonlinear boundary conditions.

Nonlinearities at the boundaries
To handle the nonlinearities at the boundaries, the following RB ansatz is used by lifting the reduced solution 7 where U n B is a vector that enables the RB solution Û n ( ) to satisfy the boundary conditions (Û can be either̂orm).
In other words, the vector U n B encodes the exact satisfaction of the boundary conditions. The RB vectors Φ = { i , i = 1, · · · , N} are the set of N RB functions evaluated at the discrete spatial points, which are obtained by the POD (Algorithm 1) applied to modified snapshots as below. To incorporate the boundary conditions correctly, i should vanish at the location of boundaries, see (19). To do so, we modify the snapshots as below (step 5 of Algorithm 2): Find U n ( k ) as the solution of (16), n = {1, · · · , N t },
Finally, a i , i = {1, 2, · · · , N}, which is the new state variable of the ROM, is the modal coordinate corresponding to i . The governing dynamics of the modal coordinates are explained in Section 3.3.
The vector U n B should be exact at the boundaries and can be an interpolation in the internal domain. For instance, for the density variable, we set where X = [1 · · · N ] T ∈ R N and 1 ∈ R N is a vector of ones. This expression at the boundaries leads to n B | x=0 = n 0 and n B | x=L = n N +1 , which satisfies the requirements for the new ansatz (19). A similar expression can be used for m n B . Remark 8. A similar lifting approach can be applied to two-(2D) and three-dimensional (3D) cases. The modified snapshots and correspondingly the RB functions then should vanish at the boundary lines or surfaces.

POD-Greedy approach
The POD-Greedy Algorithm is summarized in Algorithm 2. To enable the implementation of Algorithm 2, we need to introduce a finite parameter domain  h ⊂  which serves as a computational surrogate for  (which is the continuous parameter domain). We assume that the set of first i parameter values have been selected and the first set of RB vectors in Φ has been computed. The next parameter value i+1 is the one corresponding to the worst approximated RB solution of (22) among all the members of the discrete parameter domain  h (which is determined in step 10 of Algorithm 2). For this parameter value, we solve the full-order model (15) and gather the snapshots of density and momentum. The worst approximation is determined by evaluating the output error between the RB solution of system (22) and the full order solution of (16). An error bound or an error estimate can be useful to drastically alleviate the computational effort of the error exploration among the domain  h 7 (if the error estimate can be computed at low computational cost, ie, the computation of the error estimate does not scale with the dimension of the full-order model). Other empirical error estimates can also be useful for nonlinear problems, such as those in References 18,19. In this study, however, during the offline stage, the actual error in the output is considered as the measure to select the next parameter. Using the actual error during the Greedy algorithm is called the strong Greedy approach 20 and it generally is computationally expensive. It is known that the Greedy approach has a higher convergence rate compared to POD and nested POD algorithms. 21 Moreover, saving all snapshots for all members of the parameter domain for all time steps induces stringent memory requirements while for the strong Greedy algorithm, we only need to save the output of the system. The POD-strong Greedy algorithm terminates when the maximum error is smaller than a given tolerance or the number of RB vectors is equal to a given number of RB vectors.
Remark 9. The POD-strong Greedy is not computationally feasible in general, but the goal of the paper is to show that it is possible to use a reduced model to compute accurate solutions, not accounting for an effective offline stage (ie, using weak greedy approach).
By inserting the ansatz (19) into (15), applying a standard Galerkin projection and taking into account the orthonormality of Φ, a reduced-order model is obtained as below: ] T + C V C [ a n a n m ] T , where the hatted operators( ⋅) are obtained by Galerkin projection of the full-order operators on the space spanned by the orthogonal RB vectors Φ and Φ m , and .
Operators such as CC 0 ∈ R N y ×2 are of low dimensions and can be defined in the offline stage. For more information, refer to Reference 7. The initial conditions for the system are also obtained by a Galekin projection of the initial condition; for instance, a 0 = Φ T 0 . It should be recalled that in the offline stage, the diameter is considered constant along each pipe.
Note that in step 8 of Algorithm 2, the POD is performed on the interpolation error instead of the modified snapshots to consistently add new information at each step of the POD-strong Greedy procedure. To apply POD (step 8 of Algorithm 2), a singular value decomposition (SVD) is performed on the modified snapshots as demonstrated in Algorithm 1. The set of RB vectors is enriched with the first left singular vector. For more information, refer to Reference 5. At this point, if the maximum number of RB vectors is not yet reached or the expected maximum error among members of the discrete parameter domain  h is not achieved, we proceed in the same way to compute i+2 .

Intermediate stage
To incorporate the effect of the spatially dependent diameters, an extra stage is added to the standard RB method, which is an intermediate stage coming in between the offline and the online stages. Although the computation of this step scales with the dimension of the full-order mode, this stage has to be performed when the diameter profile is changed. If this stage is not included, whenever the diameter profile is changed, the entire offline stage should be repeated. In the intermediate stage, operators related to the diameter are defined and the RB vectors are enriched with some local basis vectors with compact support. Then, the enriched set of basis vectors is orthonormalized. This is explained in Algorithm 3 and in the following section.

Algorithm 3. Intermediate stage including local enrichment
Inputs: Continuous RB vectors Φ, diameter profile d, D out , FV operators as in (10)

Localized enrichment
For capturing the internal nonlinearities (discontinuous diameter inside the domain), a local enrichment of RB vectors is introduced here, which also shares similarities with the EIM approach for the distributed nonlinearities. The generic EIM procedure consists of an offline stage to generate the collateral basis functions and detect some interpolation points. 9 Then, the distributed nonlinearities are approximated by linear interpolation of these collateral basis functions. The coefficients of the interpolation are computed such that the interpolation is exact at the interpolation points. Generally, EIM is tailored for distributed nonlinearities, not for local nonlinearities. Here, by a similar formulation of the problem suited to EIM (see (10) and (11) where operator appears in the equations), the nonlinearity computation is decoupled from the linear subsystem.
In the current problem setting, the MOR method is supposed to be flexible and deal with any generic cross-sectional area of the domain. The nonlinearities do not enter the problem until the area discontinuity is introduced by the user in the intermediate stage. Here, unlike the procedure in the EIM approach, the interpolation points are already known and selected at the locations of area discontinuity. Moreover, the local basis vectors are the same as introduced in (10) and (11).
The locations of area discontinuities can be detected as Therefore, the interpolation points are defined as Then, the local basis vectors are defined as the collocation of canonical vectors corresponding to the locations in x m , that is, = {e k , e k+1 |D h (x k ) ≠ D h (x k+1 ), k = 1, · · · , N − 1}, where e k is the canonical vector with the nonzero entity at kth element.
If no area variation occurs, the local basis vector is a null vector. The local nonlinearities of area variation,  ∈ R N for instance for the density, can be approximated by where +,n d is calculated at the locations of x m as in (13). In other words,  is a vector with nonzero elements at the locations of area discontinuities and zero elements elsewhere. In this way, we approximate the local nonlinearities by an affine combination of precomputed basis vectors and local values of the field variables. This also enables the exact interpolation at the location of area variation.
At this stage, RB vectors are enriched with local basis vectors to take into account the effect of area variations and also the reduced operators are redefined. To account for the discontinuities in the area, the set of continuous RB vectors, obtained during the offline stage, are enriched by local discontinuous basis vectors with compact support. For instance, for an area variation at the discrete point x i+1/2 , one of the two basis vectors shown in Figure 4 is added to the set of RB vectors and then the orthogonalization is applied. These discontinuous basis vectors account for the discontinuities generated in the solution due to the discontinuous diameter and these should also vanish at the boundaries. This is explained in Algorithm 3. The motivation behind choosing the discontinuous basis vectors shown in Figure 4 is explained in the following remark. Remark 10. In the case of a discontinuous diameter, the snapshots in U, and therefore the RB vectors in , are also discontinuous. Since, in the offline stage, local discontinuities are ignored, we should account for these nonlinearities in the intermediate stage. Assume that we have the steady-state solution as the initial condition in the online stage. Then, we should enrich the current set of RB vectors such that the following cost function: is minimized with respect to , whereŪ ss is the modified steady-state solution after subtracting the effect of the boundary conditions similar to (20). For an arbitrary set of parameters, this modified steady-state solution for the momentum for example resembles Figure 3. The optimum way to minimize this objective function is to enrich the set of RB vectors with the modified steady-state solution. However, this solution is not known beforehand. Therefore, local basis vectors should be added to the set of RB vectors to reproduce at least the steady-state solution. To this end, the set of basis vectors introduced in Figure 4 are added. These basis vectors are designed to take into account the discontinuities in the solution, while the continuous RB vectors determine the solution between two adjacent discontinuities.
Remark 11. The selection of the local basis vectors is problem specific. For 2D and 3D cases, based on the difference between the solution in the absence and presence of the local nonlinearities, some basis surfaces or volumes should be defined accordingly.

Online stage
After obtaining the enriched basis vectorsΦ from Algorithm 3, the linear part of system (16) is reduced from Σ lin toΣ lin by a standard Galerkin projection into the space spanned byΦ while the local nonlinear dynamics Σ nl (area discontinuities and boundary equations) remain exactly as in the original model (16). A schematic view of the reduced model is shown in Figure 5. This figure illustrates that auxiliary outputs that are necessary for computing the output of Σ nl , y n W , and y n V , are provided via a feedback interconnection. Then, the outputs from Σ nl are fed into the reduced linear systemΣ lin to incorporate the effect of the boundary conditions and area variation.
Similar to (22) and by embedding the effect of discontinuous diameter, we havê ] T + C V C [ a n a n m ] T , where the effect of discontinuous diameter is incorporated by operators such asB * 1 and variablesŴ n . Recall that operators such as CC 0 ∈ R N y ×2 are of low dimension and are defined in the offline stage. Operators such as C W C 0 are defined in the intermediate stage due to their dependencies over the locations of area variations. In the next section, the notion of offline-online decomposition, including the intermediate stage taking into account the area variation, is explained.

Offline-intermediate-online decomposition
In the offline stage, no area variation is considered. It means that B * i in (16) is zero during the offline stage and the RB vectors are generated without considering the discontinuous diameter. The affine property of the model with respect to the parameters (Assumption 1) leads to the offline-online decomposition of the computational costs. 5 Since changing the locations of the diameter discontinuity changes the structure of the wave propagation inside the domain, the diameter profile of the drillstring and the annulus is not a function of the spatial coordinate x in the offline stage. If the locations of area variation change from one simulation to another in the offline stage, the snapshots and, therefore, the RB vectors experience discontinuity at different locations. Then, reconstructing the solution from these RB vectors leads to a stair-case problem. This problem also arises in many other applications. 22 To avoid this issue, no area discontinuity is considered in the offline stage. It means that the snapshots and therefore the RB vectors, which are obtained by performing a POD on the modified snapshots, are continuous. To obtain the parameters corresponding to the worst RB approximation during the Greedy algorithm, the actual error (as in line 10 of Algorithm 3) or an estimate of the error indicated in Reference 7,23 can be used.
However, as diameters are spatially dependent and should be set by the user before running the online simulation, we define an intermediate stage between the offline and the online stages. In this intermediate stage, the diameter is defined by the user. Also, the operators in (26) that are dependent on this diameter are defined in this stage. Although this stage involves computation scaled by the actual degrees of freedom of the original system, this is done only once per each change in the diameter profile. If this stage is absent, the entire expensive offline stage should be repeated for each change in the diameter profile, which would be far more expensive computationally.
In the online stage, by changing any of the varying parameters as well as the diameter profile, the solution can be obtained accurately in a computationally cheap fashion. However, if the locations of the area variation vary, the

NUMERICAL RESULTS
In this section, the numerical results are divided into two categories. First, results are obtained by using the same external control inputs (as u n c in (16) and (26)) during the offline and online stages (both stages simulate one drilling scenario). Second, without varying the offline stage, the external control inputs in the online stage are changed to simulate another scenario in drilling. The two test cases mentioned above are relevant to two real drilling scenarios, which are explained below.
Connection: The action of adding a stand of drillpipe to the drillstring to continue drilling deeper is called connection. To do so, the pump is stopped and choke is partially closed to maintain the downhole pressure in a specific range. A new drillpipe is attached to the top of the drillstring to elongate the drillstring. After this, the pump is turned on and the choke is opened up and the drilling procedure resumes. 2,24 Choke plugging: As the mud along with the cuttings travels up along the annulus, the mixture passes through the choke manifold. When some of the cuttings get stuck in the choke, due to their weight and change of flow direction, the flow area of the choke decreases if no action is taken by the operator. Any change in the choke area might affect the choke pressure and consequently the bottomhole pressure. The contingency of having some cuttings getting stuck in the choke yielding in the unplanned reduction of the flow area in the choke manifold is called choke plugging. 2,24

Training for connection scenario
In this study, the varying parameters are the mud viscosity , the well length L, the well inclination , the speed of sound c l , the reference density of the mud 0s , and the drillstring diameter d and wellbore diameter D out . The initial condition of the system (16) is its steady-state solution corresponding to a specific parameters of interest or the parameters selected during the Greedy algorithm. Then, by changing the inputs to the system (q p and z c ), the dynamics of the system are excited. In other words, the dynamic solution is a perturbation of the steady-state solution. Table 2 contains the range of the varying parameters. Due to memory constraints, only 128 equispaced samples are considered in the discrete parameter domain. There are also some fixed parameters from one simulation to another, which are listed in Table 3. In the last row of Table 2, the set of parameters selected for the online simulation o is reported, which does not lie in the discrete parameter domain  h . The diameters used in the online stage are shown in Figure 6. The number of grid cells is N = 500 and the fixed time horizon is T = 50 seconds. In each simulation, the time-step Δt is changed such that c l Δt∕Δx = 0.8, which is the Courant-Friedrichs-Lewy (CFL) number 14 associated to each simulation. The output of interest is the pressure in the last 10 percent of the domain in the annulus (ie, downhole in the well) where the pressure needs to be calculated accurately. The boundary conditions variability over time during offline stage are shown in Figure 7. This type of input simulates the connection scenario commonly performed in practice.
Remark 12. As the length of the spatial domain L varies, we scale the RB vectors on a unit length and scale back for any new given length.

F I G U R E 7 Input signals for the connection scenario
Remark 13. The interconnected modelΣ may get unstable due to the interconnection of the area discontinuity and the reduced linear subsystemΣ lin . One way to reduce the gain of this feedback interconnection is to reduce the CFL number (for CFL= 0.8, we did not encounter instability in the simulations).
The maximum of the error indicator computed in the POD-strong Greedy algorithm (line 10 of Algorithm 2) during the offline stage among all the members of the training set  h is plotted in Figure 8 where no area variation is imposed.
A comparison of the pressure profile between the FV and RB solutions is shown in Figure 9 at four different time instants during the online stage where the area discontinuity of Well 1 is considered. It should be noted due to the low compressibility of the drilling mud, that the pressure does not vary significantly over the discontinuous area and its discontinuties are not visible in Figure 9. The maximum error in the downhole pressure approximation is less than 1 bar. The dimension of the full-order model is 2000 while the dimension of the reduced-order model is 80 plus 6 basis vectors added in the intermediate stage to capture the area discontinuity (one diameter discontinuity in the drillstring for both momentum and density and two in the annulus for both momentum and density). To decrease the approximation error, the number of RB vectors should be increased. However, increasing the number of RB vectors does not resolve the so-called Gibbs phenomenon due to the hyperbolic nature of the system under investigation. 25 In addition to the pressure, Figures 10 and 11 show the comparison of the velocity profiles between the FV and the RB solution. To show the applicability of the method, the RB results are also compared with the POD results applied to the snapshots of the full-order model with discontinuity. Apparently, POD works better; however, when the location of the area discontinuity varies, the entire offline stage should be implemented again. The aforementioned Gibbs phenomenon is more visible in Figures 10 and 11. However, in the MPD context, the approximation of pressure is much more important than approximating the velocity.

F I G U R E 11
Comparison between velocity profiles obtained by the finite-volume solution, the reduced basis solution and direct proper orthogonal decomposition applied to the discontinuous solution in the annulus for the online parameter o as in Table 2, connection scenario The inaccuracies observed in Figures 10 and 11 are due to the hyperbolic nature of the system, which is beyond the scope of this paper. It can be inferred from the results that the proposed methodology for capturing the local nonlinearities works efficiently and accurately and can be a promising approach to speed up real-life drilling operations.
The time-wise comparison of the bottomhole pressure computed by the full-order and reduced-order model together with the approximation error is shown in Figure 12. Apparent from this figure, the wave reflection occurred due to the area variation has been captured accurately by the reduced-order model (see the step-wise increase in the downhole pressure over time, which is due to the wave reflection at the area discontinuities and boundaries). This phenomenon F I G U R E 12 Time-wise comparison of the bottomhole pressure between the finite-volume and reduced basis solution for the online parameter o as in Table 2 significantly changes the dominant resonance frequency of the system, which has to be captured by the reduced-order model. The time-wise change of the velocity at the choke together with the approximation error can be seen in Figure 13, which further confirms the accuracy of the reduced-order model. The CPU time allocated to solve the offline stage, the intermediate stage, the online stage and the full-order model together with the speedups are shown in Table 4 by using different numbers of RB functions. Number of RB functions for each state is the same and therefore the total number of RB functions is multiplied by four. Since there is one diameter discontinuity in the drillstring and two diameter discontinuities in the annulus, six extra local basis functions are added to enrich the previous set of RB functions.   Figure 6 tested for connection scenario up to 70 can be achieved while maintaining a high accuracy in the output of interest. The large CPU time of the offline stage (around 17 hours) is due to the high dimension of the parameter domain and the requirement to solve the full-order model (without area discontinuity) 128 times.
To further evaluate the performance of the proposed method, five random points in the parameter domain (outside the training parameter domain) are selected and new simulations for the connection scenario in the two geometries shown in Figure 6 are computed and compared with the full-order solutions. In Table 5, the following relative error indicator for the different parameter samples is reported: (27) The reported relative errors confirm the accuracy of the approach for different well geometry configuration and different parameter samples. The results also reveal that the relative error is higher if the number of area discontinuities is higher due to the many discontinuities in the state variables.

Choke plugging scenario
In the previous section, both the offline stage and the online stage are computed with the same external control inputs; only the parameters of the simulation are changed. In this section, based on the offline stage performed on connection scenario, we simulate a choke plugging scenario (with different inputs compared to the offline stage) and show the results. In a choke plugging scenario, the choke opening decreases due to cuttings getting trapped inside the choke. As there are two parallel chokes installed at the choke manifold, the plugged choke is closed and the stand-by choke is opened. Then, the personnel cleans the plugged choke. To simulate this effect, the choke opening z c is reduced for some period of time and then increased to the previous level. Throughout all these actions, the pump flow is constant. The input signals are shown in Figure 14.
Similar to the previous test case, the time-wise comparisons of the bottomhole pressure and the velocity at the choke together with their corresponding approximation errors are shown in Figures 15 and 16. The results confirm that, although the training (the offline stage) is performed for a connection scenario, it can simulate other drilling scenarios. Speedup in this case is similar to Table 4 of the previous test case.
The results presented in this section confirm that the proposed method successfully captures the effect of the local nonlinearities inside the computational domain and at the boundaries of the system during MOR.

CONCLUSION
In this paper, a new approach for addressing the issue of localized nonlinearities with the RB method has been proposed. These nonlinearities appear either at the boundaries of the system or inside the computational domain. As the boundary conditions for industrial systems are often nonlinear, it is vital to change the RB ansatz, such as the one suggested in this paper. For the nonlinearities inside the domain, a local enrichment approach bearing similarities with EIM is incorporated and its outputs are coupled with the rest of the system dynamics. Using this approach, crucial underlying physics of the phenomena such as mass conservation at a particular point and wave reflections are preserved during MOR. Any other spatially dependent local nonlinearities can be handled by this method. Locating pump stations and orifices in a pipeline can be seen as other applications which may benefit from the proposed approach. The speedup obtained by using the proposed approach can expedite, for instance, optimization for well planning, real-time simulations, and the implementation of (model-predictive and optimal) control techniques.