On the criticality of return flows in viscous accretionary wedges and its implications for deep-crustal exhumation in subduction zones

In subduction zones, the accretionary wedges play a vital role in mediating the burial processes of incoming oceanic sediments and eventually their return pathways to the surface. A direction of the previous tectonic models invoked the standard corner flow theory, assuming a slab-parallel shear and a rigid, fixed overriding plate, to elucidate the crustal recycling processes in tectonic wedges. To deal with more complex subduction-collisional settings, where they have deformable overriding plates, and associate a horizontal slab migration (advance or rollback) component during subduction, we develop a generalized corner flow model to revisit the problem of return flow mechanics, providing a criticality analysis of the return flows as a function of the geometric, kinematic, and rheological conditions in accretionary wedges. A new set of analytical solutions is presented to evaluate the limiting conditions in which a wedge can set in significant return flows, leading to focused exhumation of the deep-crustal materials. The theoretical results suggest that, for moderate wedge-taper angles (∼30°), the viscosity ratios (μ_r) between the overriding plate and the wedge $\ge \sim {\mathrm{10}}^{\mathrm{3}}$ provide favourable tectonic settings for the return flow kinematics in accretionary wedges. Decrease in μ_r, or addition of slab roll back weakens the return flows, whereas slab advance greatly strengthens the return flows. The analytical solutions are also utilized to demonstrate reversals in the shear-sense patterns across the wedge. We expand this study by incorporating results from scaled laboratory experiments to evaluate applicability of the generalized theoretical model. It is shown from the theoretical model that the total pressure in the accretionary wedge dynamics becomes close to the lithostatic value when the rheological setting has low-viscosity (10¹⁹ Pa s) wedge materials.

1 Introduction

Accretionary wedges that generally accommodate large volumes of sediments fed by subducting oceanic plates play a critical role in modulating the mechanics of active plate interfaces in subduction zones (Agard et al., 2018; Angiboust et al., 2022; Beaumont et al., 1999; Cawood et al., 2009; Cloos and Shreve, 1988). The accreted sediments act as lubricating agents in weakening the plate interface and facilitate subduction (Behr et al., 2022; Behr and Becker, 2018; Hu et al., 2021; Pusok et al., 2022). Geophysical observations suggest that a fraction of the wedge sediments can also be dragged to deeper mantle regions along the subducting slab in the form of a narrow channel, called subduction channel (Agard et al., 2009; Gerya et al., 2002; Shreve and Cloos, 1986; Abers, 2005). Petrological investigations reveal the omnipresence of high- and ultra-high-pressure (HP-UHP) and low-temperature (LT) metamorphic rocks preferentially in the accretionary belts of both paleo- and modern subduction zones (Agard et al., 2001; Glodny et al., 2005; Platt, 1993; Warren, 2013). The pressure-temperature-time (P-T-t) paths of HP–LT rocks thus provide key fingerprints of the deep-crustal tectonic processes- beginning from burial of sediments and shallow-crustal materials, followed by their HP metamorphism at greater depths and finally their recycling back to the Earth’s surface during active subduction (Van Dinther et al., 2012; Yamato et al., 2007; Kerswell et al., 2023; Vaughan-Hammon et al., 2022). During this recycling process, the wedge materials also undergo syn-metamorphic deformations in multiple episodes, giving rise to complex tectonic structures at convergent plate boundaries (Allemand and Lardeaux, 1997; Ring et al., 2020). Understanding the kinematics of accretionary wedges is thus crucially important to interpret the mechanisms of deep-crustal exhumations and associated tectonic processes in subduction zones, which is the central theme of the present theoretical study.

Geochronological records of exhumed blueschists and eclogites suggest that high-pressure rocks can be exhumed at varied rates (<1 mm yr⁻¹ to a few cm yr⁻¹) depending upon the tectonic environments (Agard et al., 2009; Ernst, 1988; Rubatto and Hermann, 2001). Conceptual, analytical and numerical models consider a wide range of factors, such as erosion and isostatic adjustment (Brandon et al., 1998), buoyant flow of low density crustal units (Chemenda et al., 1995; Ernst et al., 1997), corner flow circulation within accretionary wedges (Burov and Yamato, 2008; Cloos, 1982; Moulas et al., 2021), extension-driven ductile thinning and normal faulting (Platt, 1986; Ring and Brandon, 2008), and several others (see Platt (1993); Warren (2013), for review) to account for such variations in the exhumation rates. Among them, the corner flow (CF) model (Batchelor, 1967), which was implemented by Cloos (1982), deals with the problem of exotic high-pressure/low-temperature (HP/LT) blocks within a mudstone matrix in the Franciscan Complex of California. According to this tectonic model, the crustal materials in the lower part of the wedge are dragged down with the subducting plate, eventually set to reverse their flow directions due to the dynamic pressure developed at the bottom of the downward-tapering wedge, as illustrated in Fig. 1. This reversal of flow, hereafter called return flow, causes syn-subduction exhumation of deep-crustal rocks back to the surface. The corner flow mechanism successfully explains the coherent occurrence of extensive HP metamorphic belts with progressively increasing P-T conditions (Agard et al., 2001). This flow kinematics accounts for sequential records of prograde high P metamorphism of the oceanic metasediments during their burial, followed by their retrograde mineral assemblages in the uplift phase (Jolivet et al., 2003). The prograde metamorphism in a lower geothermal gradient condition, as prevalent in most of the subduction zones results in formation of blueschist, or even eclogites in accretionary wedge complexes. Furthermore, entrainment of oceanic sediments or crustal fragments into the deeper subduction channel can often lead to partial melting at mantle depths (Wang et al., 2014; Hebert et al., 2009). Petrological data syntheses (e.g., Agard et al., 2001) point to, at least two different corner-flow circulations in the tectonic evolution of HP and UHP metamorphic complexes: a relatively shallower (up to 40 km) circulation within accretionary complexes that can be correlated with the formation and exhumation of blueschists and most eclogites, and a deeper circulation in subduction channels, which is responsible for the formation and exhumation of UHP eclogites (Agard et al., 2009; Burov et al., 2001; Wang et al., 2019).

The return flow kinematics has been extensively used in analytical studies of viscous flows in narrow (subduction channel) and wide (accretionary wedge) corner regions at the convergence between two rigid or deformable plates (Cloos, 1982; Grujic, 2006; Iwamori, 2003; Mancktelow, 1995; Moulas et al., 2021). Numerical CF models, including both fixed and evolving corner geometry produce return flows in accretionary wedges and subduction channels (Allemand and Lardeaux, 1997; Angiboust et al., 2022; Gerya et al., 2002; Yamato et al., 2007; Beaumont et al., 2009; Sanhueza et al., 2022). These models impose subduction-related velocity boundary conditions either at the wedge walls, or at the boundaries of a large model box. Some of them (e.g., Van Dinther et al., 2012) also consider dynamic subduction, driven by the slab pull force. Interestingly, a number of these models (see Candioti et al., 2021; Pusok et al., 2022) yield flow kinematics in the accretionary wedge different from those predicted from the classical corner flow theory (Cloos, 1982, Fig. 1), which needs further investigation. In addition, the exhumation velocities obtained from dynamic CF models (e.g., Van Dinther et al., 2012 (∼2.6 mm yr⁻¹); Yamato et al., 2007 (<6 mm yr⁻¹)) are comparable to those estimated from radiometric dating of exhumed continuous HP units (Agard et al., 2002 (1–5 mm yr⁻¹); Federico et al., 2005 (3–4 mm yr⁻¹)), but low compared to those obtained in earlier analytical and some numerical models (in the order of cm yr⁻¹ in Allemand and Lardeaux, 1997; Cloos, 1982; Marques et al., 2018). These contrasting findings stimulate a few fundamental questions on the mechanics of return flows: (1) what factors govern the onset versus the absence of return flows in downward-tapering accretionary wedges, (2) what modulates the exhumation rates of the returning wedge materials, and (3) how does the varying flow kinematics affect the deformation patterns within the wedge and in the overriding plate? These questions motivate the present theoretical study, which combines analytical and analogue modelling, with an objective to explore the conditions that facilitate the return flows kinematics and associated deformations in the deeper, viscously deforming part (between ∼20 to ∼40 km depth) of accretionary wedges (Fernández-Blanco et al., 2020; Pavlis and Bruhn, 1983). The modelling approach is, however, simplified by treating the accretionary wedge as a downward tapering closed region, filled with materials of linear viscous rheology, and by assuming no buoyancy effects on the wedge kinematics. Also, this modelling excludes sedimentation, erosion, or any entrainment of wedge materials into the mantle.

Although the classical CF model of Cloos (1982) provides an analytical solution for the velocity field in a corner region, the solution is applicable to a setting with the boundary velocity parallel to the channel's base, and that too for a fixed and rigid hanging wall of the channel. However, the assumption of a rigid hanging wall in the CF modelling may not always be geologically appropriate (Schellart, 2024), as many convergent zones accommodate significant deformations within the overriding plates. To address this limitation, Moulas et al. (2021), advanced the analytical approach by incorporating a deformable overriding plate. In addition, subduction zones often evolve through slab rollback or slab advance (Bezada et al., 2013; Cawood et al., 2009; Di Giuseppe et al., 2009; Nakakuki and Mura, 2013), which originate from a complex interplay of various forces, such as, ridge push, slab pull, and mantle convection (Ribe, 2010; Schellart, 2008; Xue et al., 2022). In such cases, subduction is associated with a non-parallel component of slab velocity. The effect of slab-oblique kinematics on the flow and deformation patterns in the viscous accretionary wedges is yet to be explored, which needs further advancement of the CF modelling.

In this article, we extend the theoretical model of Moulas et al. (2021) to develop a modified corner flow solution, incorporating the effect of slab-oblique subduction on the flow kinematics in an accretionary wedge with a deformable hanging wall. A parametric analysis is presented to show the effects of obliquity of subduction motion, wedge geometry (narrow vs. wide), the wedge and the overriding plate viscosities, and an imposed overriding plate velocity, on the flow and deformation pattern within the accretionary wedge. Additionally, a series of scaled laboratory experiments were performed in studying the kinematic evolution of wedges within deformable walls of varying viscosities under oblique subduction conditions. We evaluate the results from our theoretical CF model in the light of the analogue model findings, and find the specific settings that favour return flows (and hence exhumation) of deep-crustal materials within an accretionary wedge. Finally, we utilize the CF model to examine the possibilities of significant non-lithostatic pressure developments within the wedge, as reported in some previous studies (e.g., Gerya, 2015; Mancktelow, 1995; Raimbourg and Kimura, 2008).

Figure 1 (a) A schematic presentation of the return flow pattern in accretionary wedges. U: velocity at the wedge base. (b) Particle paths in the ideal corner flow within a wedge. The particles undergo the following sequential events: (1) burial in the time interval between T₁ and T₂; (2) high P-low T metamorphism at T₂, where the total pressure (P) exceeds the lithostatic pressure (P_lith), and finally, (3) return flow-driven exhumation in the time interval between T₂ to T₃. Note that the wedge develops opposite shear senses in the lower and upper parts. (c) Velocity profile across the wedge. Negative and positive values of the radial velocity component (u_r) indicate burial and exhumation, respectively. θ₁: wedge taper angle. r: distance from the wedge tip.

2 Theory

2.1 Basic premises

The present theoretical model is developed in the premises of downward-tapering accretionary wedges, which has been extensively used in earlier theoretical models (Cloos, 1982), analogue experiments (Luján et al., 2010) and numerical simulations (Pusok et al., 2022; Yamato et al., 2007). To simplify the theoretical treatment, the wedge is assumed to terminate at a fixed-point S, located at a depth of about 40 km (S-point boundary condition of Willett et al., 1993) (Fig. 2a). This is a limitation in our model setup, as the wedges in subduction zones can describe a narrow opening at their base, allowing a small fraction of sediments to be transported to mantle depths (Agard et al., 2009; Guillot et al., 2009). Traditionally, researchers have modelled accretionary wedges using Coulomb materials, primarily to analyze their mechanical stability as a function of the taper angle (Dahlen et al., 1984; Dahlen, 1990). This modelling approach is, however, applicable for shallow-crustal deformations. Below the brittle-ductile transition, at middle and lower crustal depths, viscous processes dominate in the wedge deformations (Mannu et al., 2016; Pavlis and Bruhn, 1983), and thereby demand a viscous rheological approach. To capture the exclusive effects of viscous wedge kinematics, our theoretical treatment excludes the additional influence of gravity that can arise due to the density contrast between the wedge and the overriding plate. In natural settings, however, buoyancy forces arising from density differences can also act as additional factor in driving the upward velocity of crustal blocks, especially at sub-crustal level (Ernst et al., 1997; Maiti et al., 2020; Schmalholz and Schenker, 2016). Also, we consider Newtonian viscous rheology to develop the present theoretical model. This linear rheological approximation is evidently simplification of the natural rheological settings that generally show a power-law creep behavior (Ranalli, 1995). In addition, their viscosity can vary with depth. However, it has been shown previously that the non-linearity of the rheology plays little role in the first-order kinematics of a corner flow (Ribe, 2015, 2018).

In the CF model of Cloos (1982), the upper bounding plate of the wedge (i.e., the overriding plate, hereafter abbreviated as OP) is considered rigid and held fixed, and the lower bounding plate (i.e., the subducting plate, hereafter abbreviated as SP) moves downward with the velocity vector oriented parallel to the wedge boundary (Fig. 1a). However, our theoretical modelling also considers a non-parallel component of slab velocity with rigid or deformable overriding plate. The subducting slab movements towards the OP and the ocean are considered to represent the slab advance and slab rollback kinematics, respectively. We simplify this theoretical approach assuming that the slab undergoes translational motion without any along slab gradient of the linear subduction vector. Our model also only considers the flow in the vertical plane perpendicular to the trench axis, and does not consider any trench-parallel flow in the wedge.

2.2 Rigid overriding plate model

The 2D accretionary wedge model is represented by a triangular domain of linearly viscous (μ₁) materials, sandwiched between a SP and an OP. The wedge walls converge to each other downward, forming a taper angle θ₁ (Fig. 2a, b). The rigid OP is held fixed, and the wedge base is subjected to a velocity U, obliquely to the wedge-subducting plate interface at an angle ϕ (positive anticlockwise; Fig. 2a, b). The slab-parallel component of U is the subduction velocity, U_s, and its horizontal component represents the slab migration velocity, U_m. In the following discussions, we use the term subduction obliquity (ϕ) to define the non-parallel component of slab velocity (measured in the trench axis-perpendicular vertical plane), which is different from the term generally used to describe the angular relationship of the convergence velocity vector with the trench line in the horizontal plane (McCaffrey, 1992). Positive and negative values of ϕ represent slab rollback and advance, respectively in this theoretical treatment. We consider no-slip boundary conditions at the two wedge walls, and also assume that the wedge undergoes no volume loss or material accretion during the deformation. A Cartesian frame xy is chosen with origin at the S-point of the wedge, with the x and y axes oriented in horizontal and vertical directions, respectively, as shown in Fig. 2b. The SP dips at angle α with respect to the horizontal surface. Considering our focus on the flow field in a triangular region, we choose a polar coordinate system for the theoretical analysis, similar to previous CF studies (Cloos, 1982; Moulas et al., 2021). The origin of the coordinate frame is set at the S-point, where r and θ represent the distance from the wedge tip (S-point) and the angular distance (positive anticlockwise) measured from the bottom wedge boundary respectively. The instantaneous velocity field within the wedge can be expressed in terms of polar coordinates (r, θ) as:

$$\begin{array}{}\text{(1a)}& {\displaystyle }{u}_{\mathrm{r}}={\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial \mathit{\theta }}}\text{(1b)}& {\displaystyle }{u}_{\mathit{\theta }}=-{\displaystyle \frac{\partial \mathit{\psi }}{\partial r}}\end{array}$$

where u_r and u_θ are the radial and tangential velocity components, respectively and ψ is the stream function. Assuming negligible inertial forces and mass conservation, the stream function ψ for a 2D viscous flow satisfies the biharmonic equation:

$$\begin{array}{}\text{(2)}& {\mathrm{\nabla }}^{\mathrm{2}}\left({\mathrm{\nabla }}^{\mathrm{2}}\mathit{\psi }\right)=\mathrm{0}\end{array}$$

Applying the boundary conditions in the present problem, we obtain the following relations:

$$\begin{array}{}\text{(3)}& {\displaystyle }\begin{array}{rl}{u}_{\mathrm{r}}& ={\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial \mathit{\theta }}}=-U\cos \mathit{\varphi },u_{\mathit{\theta }}=-{\displaystyle \frac{\partial \mathit{\psi }}{\partial r}}\\ & =-U\sin \mathit{\varphi }\text{ at }\mathit{\theta }=\mathrm{0}\end{array}\text{(4)}& {\displaystyle }{u}_{\mathrm{r}}={\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial \mathit{\theta }}}=\mathrm{0},u_{\mathit{\theta }}=-{\displaystyle \frac{\partial \mathit{\psi }}{\partial r}}=\mathrm{0}\text{ at }\mathit{\theta }={\mathit{\theta }}_{\mathrm{1}}\end{array}$$

These boundary conditions require that ψ is proportional to r everywhere in the wedge domain, including its boundaries. The stream function can be thus expressed as,

$$\begin{array}{}\text{(5)}& \mathit{\psi }(r,\mathit{\theta })=rf\left(\mathit{\theta }\right)\end{array}$$

Substituting ψ in Eq. (2), we find:

$$\begin{array}{}\text{(6)}& {\mathrm{\nabla }}^{\mathrm{2}}\left[{\displaystyle \frac{\mathrm{1}}{r}}\left(f+f^{\prime \prime }\right)\right]={\displaystyle \frac{\mathrm{1}}{r^{\mathrm{3}}}}\left(f+\mathrm{2}{f}^{\prime \prime }+f^{iv}\right)=\mathrm{0}\end{array}$$

where the primes denote the order of differentiation on f. The solution of Eq. (6) follows:

$$\begin{array}{}\text{(7)}& f\left(\mathit{\theta }\right)=A\sin \mathit{\theta }+B\cos \mathit{\theta }+C\mathit{\theta }\sin \mathit{\theta }+D\mathit{\theta }\cos \mathit{\theta }\end{array}$$

where A, B, C and D are four constants, which need be solved using four equations. Substituting f in Eq. (5) and putting the derivative expression in Eq. (1) we get:

$$\begin{array}{}\text{(8a)}& {\displaystyle }\begin{array}{rl}{u}_{\mathrm{r}}& =A\cos \mathit{\theta }-B\sin \mathit{\theta }+C\left(\sin \mathit{\theta }+\mathit{\theta }\cos \mathit{\theta }\right)\\ & +D\left(\cos \mathit{\theta }-\mathit{\theta }\sin \mathit{\theta }\right)\end{array}\text{(8b)}& {\displaystyle }{u}_{\mathit{\theta }}=-\left(A\sin \mathit{\theta }+B\cos \mathit{\theta }+C\mathit{\theta }\sin \mathit{\theta }+D\mathit{\theta }\cos \mathit{\theta }\right)\end{array}$$

Figure 2 (a) Schematic diagram of the tectonic setting considered for the present theoretical treatment. The subducting plate, dipping at angle α, moves downward obliquely at an angle ϕ to the bottom boundary of the wedge. Dashed rectangle demarcates the region of interest in this study. (b) Theoretical model setup in polar coordinate system (r, θ). The wedge and the overriding plate describe taper angles of θ₁ and (θ₂−θ₁), respectively. μ₁: wedge material viscosity. The overriding plate is considered either rigid or deformable (viscosity μ₂). See text for details.

Putting θ=0, and θ₁ in Eq. (8) and comparing the resulting expression with the boundary conditions of Eqs. (3), and (4), it follows

$$\begin{array}{}\text{(9a)}& {\displaystyle }A+D=-U\cos \mathit{\varphi }\text{(9b)}& {\displaystyle }B=U\sin \mathit{\varphi }\text{(9c)}& {\displaystyle }\begin{array}{rl}A\cos {\mathit{\theta }}_{\mathrm{1}}& -B\sin {\mathit{\theta }}_{\mathrm{1}}+C\left(\sin {\mathit{\theta }}_{\mathrm{1}}+{\mathit{\theta }}_{\mathrm{1}}\cos {\mathit{\theta }}_{\mathrm{1}}\right)\\ & +D\left(\cos {\mathit{\theta }}_{\mathrm{1}}-{\mathit{\theta }}_{\mathrm{1}}\sin {\mathit{\theta }}_{\mathrm{1}}\right)=\mathrm{0}\end{array}\text{(9d)}& {\displaystyle }A\sin {\mathit{\theta }}_{\mathrm{1}}+B\cos {\mathit{\theta }}_{\mathrm{1}}+C{\mathit{\theta }}_{\mathrm{1}}\sin {\mathit{\theta }}_{\mathrm{1}}+D{\mathit{\theta }}_{\mathrm{1}}\cos {\mathit{\theta }}_{\mathrm{1}}=\mathrm{0}\end{array}$$

The solutions of Eq. (9) are:

$$\begin{array}{}\text{(10a)}& {\displaystyle }A={\displaystyle \frac{U}{({\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}-{\sin }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}})}}\left[\sin \mathit{\varphi }\left({\mathit{\theta }}_{\mathrm{1}}+\sin {\mathit{\theta }}_{\mathrm{1}}\cos {\mathit{\theta }}_{\mathrm{1}}\right)-{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}\cos \mathit{\varphi }\right)\text{(10b)}& {\displaystyle }B=U\sin \mathit{\varphi }\text{(10c)}& {\displaystyle }C={\displaystyle \frac{U}{({\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}-{\sin }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}})}}[\cos \mathit{\varphi }\left({\mathit{\theta }}_{\mathrm{1}}-\sin {\mathit{\theta }}_{\mathrm{1}}\cos {\mathit{\theta }}_{\mathrm{1}}\right)-\sin \mathit{\varphi }{\sin }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}]\text{(10d)}& {\displaystyle }D={\displaystyle \frac{U}{({\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}-{\sin }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}})}}[-\sin \mathit{\varphi }\left({\mathit{\theta }}_{\mathrm{1}}+\sin {\mathit{\theta }}_{\mathrm{1}}\cos {\mathit{\theta }}_{\mathrm{1}}\right)+\cos \mathit{\varphi }{\sin }^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}]\end{array}$$

2.3 Deformable overriding plate model

To deal with the kinematic problem of accretionary wedges with deformable OP, we solve the velocity field in two corner regions (Anderson and Davis, 1993; Moulas et al., 2021) representing the wedge and the overriding plate (Fig. 2b). The theoretical model setup is considered same as in the rigid OP model, except that the present case incorporates an additional corner region with viscosity μ₂ in the hanging wall of the wedge, representing the viscous crust of the overriding plate (Moulas et al., 2021). The OP–wedge viscosity ratio is denoted as ${\mathit{\mu }}_{\mathrm{r}}(={\mathit{\mu }}_{\mathrm{2}}/{\mathit{\mu }}_{\mathrm{1}})$. The wedge and the overriding plate have taper angles θ₁ and θ₂−θ₁ respectively. In this analysis, we also impose a velocity (U_op) at the Moho of the overriding plate. Several authors have, however, argued that the cold upper mantle beneath the OP in proximity to the subduction zone is stationary and does not participate in the surrounding mantle corner flow (Lee and Kim, 2023; van Keken and Wilson, 2023). The presence of this cold region is attributed to mechanical decoupling at the interface between the forearc mantle and subducting slab (Abers et al., 2017; Wada and Wang, 2009). Hence, in most of our theoretical analyses, we consider U_op to be equal to zero.

Following a similar approach of Anderson and Davis (1993), and Moulas et al. (2021), we choose two stream functions: ψ₁ and ψ₂ for solving the velocity fields in the wedge and the OP crustal domains, respectively. Similar to the theoretical treatment in Sect. 2.2.1, the two stream functions are expressed in the form:

$$\begin{array}{}\text{(11a)}& {\displaystyle }& {\displaystyle }{\mathit{\psi }}_{\mathrm{1}}(r,\mathit{\theta })=r{f}_{\mathrm{1}}\left(\mathit{\theta }\right)\text{(11b)}& {\displaystyle }& {\displaystyle }{\mathit{\psi }}_{\mathrm{2}}(r,\mathit{\theta })=r{f}_{\mathrm{2}}\left(\mathit{\theta }\right)\end{array}$$

where,

$$\begin{array}{}\text{(12a)}& {\displaystyle }& {\displaystyle }{f}_{\mathrm{1}}\left(\mathit{\theta }\right)=A_{\mathrm{1}}\sin \mathit{\theta }+B_{\mathrm{1}}\cos \mathit{\theta }+C_{\mathrm{1}}\mathit{\theta }\sin \mathit{\theta }+D_{\mathrm{1}}\mathit{\theta }\cos \mathit{\theta }\text{(12b)}& {\displaystyle }& {\displaystyle }{f}_{\mathrm{2}}\left(\mathit{\theta }\right)=A_{\mathrm{2}}\sin \mathit{\theta }+B_{\mathrm{2}}\cos \mathit{\theta }+C_{\mathrm{2}}\mathit{\theta }\sin \mathit{\theta }+D_{\mathrm{2}}\mathit{\theta }\cos \mathit{\theta }\end{array}$$

A_i, B_i, C_i, D_i (i=1, 2) represent eight constants. To find the velocity fields in the wedge and OP regions, we need to evaluate the eight constants, which evidently demands eight boundary conditions. At the subducting plate-wedge interface, the two velocity boundary conditions in terms of polar coordinates (r, θ) is:

$$\begin{array}{}\text{(13a)}& {\displaystyle }& {\displaystyle }(u_{\mathrm{r}}{)}_{\mathit{\theta }=\mathrm{0}}=-U\cos \mathit{\varphi }\text{(13b)}& {\displaystyle }& {\displaystyle }(u_{\mathit{\theta }}{)}_{\mathit{\theta }=\mathrm{0}}=-U\sin \mathit{\varphi }\end{array}$$

Velocity boundary conditions at the Moho in the overriding plate give,

$$\begin{array}{}\text{(14a)}& {\displaystyle }& {\displaystyle }(u_{\mathrm{r}}{)}_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{2}}}=-U_{\mathrm{op}}\text{(14b)}& {\displaystyle }& {\displaystyle }(u_{\mathit{\theta }}{)}_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{2}}}=\mathrm{0}\end{array}$$

Assuming no slip at the wedge–OP interface, the velocity continuity across this interface leads to:

$$\begin{array}{}\text{(15a)}& {\displaystyle }& {\displaystyle }\left(\right(u_{\mathrm{r}}{)}_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{1}}}{)}_{\mathrm{wd}}=\left(\right(u_{\mathrm{r}}{)}_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{1}}}{)}_{\mathrm{op}}\text{(15b)}& {\displaystyle }& {\displaystyle }\left(\right(u_{\mathit{\theta }}{)}_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{1}}}{)}_{\mathrm{wd}}=\left(\right(u_{\mathit{\theta }}{)}_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{1}}}{)}_{\mathrm{op}}\end{array}$$

where the subscripts “wd” and “op” refer to the wedge and the overriding plate. The total normal stress as well as the shear stress must be continuous across the wedge–OP interface. The deviatoric stress tensor τ_ij is expressed in polar coordinates (r, θ) as,

$$\begin{array}{}\text{(16)}& \begin{array}{rl}& \left[\begin{array}{cc}{\mathit{\tau }}_{\mathrm{r}r}& {\mathit{\tau }}_{\mathrm{r}\mathit{\theta }}\\ {\mathit{\tau }}_{\mathit{\theta }\mathrm{r}}& {\mathit{\tau }}_{\mathit{\theta }\mathit{\theta }}\end{array}\right]\\ & =\left[\begin{array}{cc}\mathrm{2}{\mathit{\mu }}_i\frac{\partial u_{\mathrm{r}}}{\partial r}& {\mathrm{\mu }}_i\left(\frac{\mathrm{1}}{r}\frac{\partial u_{\mathrm{r}}}{\mathit{\delta }\mathit{\theta }}+\frac{\partial u_{\mathit{\theta }}}{\mathit{\delta }r}-\frac{u_{\mathit{\theta }}}{r}\right)\\ {\mathrm{\mu }}_i\left(\frac{\mathrm{1}}{r}\frac{\partial u_{\mathrm{r}}}{\partial \mathit{\theta }}+\frac{\partial u_{\mathit{\theta }}}{\partial r}-\frac{u_{\mathit{\theta }}}{r}\right)& \mathrm{2}{\mathrm{\mu }}_i\left(\frac{\mathrm{1}}{r}\frac{\partial u_{\mathit{\theta }}}{\partial \mathit{\theta }}+\frac{u_{\mathrm{r}}}{r}\right)\end{array}\right]\end{array}\end{array}$$

where μ_i (i=1, 2) denotes the viscosities of the two corner regions. The total stress tensor σ_ij is given by:

$$\begin{array}{}\text{(17)}& {\mathit{\sigma }}_{ij}=-P{\mathit{\delta }}_{ij}+{\mathit{\tau }}_{ij}\end{array}$$

where P is the total pressure, given by the sum of static pressure (P_s) and dynamic pressure (P_d). As P_s has no implications in the present corner flow solution, P=P_d, which can be determined from the stress-equilibrium equation in the radial direction,

$$\begin{array}{}\text{(18)}& -{\displaystyle \frac{\partial P_{\mathrm{d}}}{\partial r}}+{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial {\mathit{\tau }}_{\mathrm{r}\mathit{\theta }}}{\partial \mathit{\theta }}}=\mathrm{0}\end{array}$$

Using the boundary condition P_d=P₀ at r=∞, we find the integration constant, which appears after solving the differential equation in Eq. (18). P₀ is set as 0 in the given problem, assuming that the dynamic pressure becomes zero at infinite distance away from the S-point of the wedge. The choice of P₀, however, does not qualitatively change the velocity or strain fields either in the wedge or the overriding plate.

The normal components (τ_rr and τ_θθ) of the deviatoric stress tensor (Eq. 16) in the present setting are zero as there is no extension along as well as across the radial direction, i.e., ${\dot{\mathit{\epsilon }}}_{\mathrm{r}r}={\dot{\mathit{\epsilon }}}_{\mathit{\theta }\mathit{\theta }}=\mathrm{0}$. This wedge kinematics results because the velocity components at any point (r, θ) within the wedge are independent to r and hence, the radial directions are the directions of the flow apophyse across the wedge. Consequently, the total normal stress (Eq. 17) perpendicular to the wedge-overriding plate interface equals the dynamic pressure P_d. The stress continuity across the interface requires (Moulas et al., 2021),

$$\begin{array}{}\text{(19)}& \left(\right(P_{\mathrm{d}}{)}_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{1}}}{)}_{\mathrm{wd}}=\left(\right(P_{\mathrm{d}}{)}_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{1}}}{)}_{\mathrm{op}},\left(\right({\mathit{\tau }}_{\mathrm{r}\mathit{\theta }}{)}_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{1}}}{)}_{\mathrm{wd}}=\left(\right({\mathit{\tau }}_{\mathrm{r}\mathit{\theta }}{)}_{\mathit{\theta }={\mathit{\theta }}_{\mathrm{1}}}{)}_{\mathrm{op}}\end{array}$$

Equations (13), (14), (15) and (19) provide the eight boundary conditions, which are utilized to solve the eight unknown constants in Eq. (12). Their expressions are provided in the Appendix A. Using these constants we determine the two stream functions (ψ₁ and ψ₂) and find the velocity fields in the polar coordinates (r, θ) at the two corner regions from Eq. (1). It is noteworthy that for given U and U_op, the velocity field is dependent on θ₁, ϕ, θ₂ (=180° − α), and μ_r (=μ₂μ₁) in an accretionary wedge with deformable OP (see Eqs. A6–A13 in Appendix A). This is in contrast to the case with rigid and fixed OP, where the velocity field is entirely dependent on the first two of the aforementioned factors, as discussed in Sect. 2.1. The velocity components determined in the polar coordinate system is converted to the Cartesian frame xy using the following transformation:

$$\begin{array}{}\text{(20)}& \left[\begin{array}{c}u\\ v\end{array}\right]=\left[\begin{array}{cc}\cos \mathit{\alpha }& -\sin \mathit{\alpha }\\ \sin \mathit{\alpha }& \cos \mathit{\alpha }\end{array}\right]\left[\begin{array}{cc}\cos \mathit{\theta }& -\sin \mathit{\theta }\\ \sin \mathit{\theta }& \cos \mathit{\theta }\end{array}\right]\left[\begin{array}{c}{u}_{\mathrm{r}}\\ u_{\mathit{\theta }}\end{array}\right]\end{array}$$

3 Theoretical results

3.1 Return flux

Return flow in a downward-tapering wedge is considered to be one of the most effective exhumation mechanisms to explain the high P- and low T rock assemblages in the accretionary wedges or subduction channels (Agard et al., 2009; Cloos and Shreve, 1988; Yamato et al., 2007). For a quantitative analysis of the return flows in a wedge, we introduce a non-dimensional parameter, termed as return flux (F_R), which is defined as the ratio of the volumetric rate of materials being returned to that of materials undergoing burial at an instant within the wedge. This parameter accounts for burial and exhumation within a wedge driven by slab-parallel subduction, as well as additional contributions arising from wedge squeezing (or unsqueezing) in cases of non-parallel subduction. In the polar coordinate system under consideration, positive and negative values of the radial component (u_r) of the velocity vector indicate the rates of return and burial of wedge materials, respectively (Fig. 1c). For the present simplified 2D problem, we assume the velocity field remaining uniform along the trench, and express F_R as:

$$\begin{array}{}\text{(21)}& F_{\mathrm{R}}={\displaystyle \frac{\text{volumetric rate of return}}{\text{volumetric rate of burial}}}={\displaystyle \frac{{\sum }_i^{N_{\mathrm{1}}}(\text{positive }{u}_{\mathrm{r}}{)}_i}{\left|{\sum }_i^{N_{\mathrm{2}}}(\text{negative }{u}_{\mathrm{r}}{)}_i\right|}}\end{array}$$

where N₁ and N₂ represent number of areal elements of regions with positive and negative radial velocity (u_r) within the wedge. High F_R implies a condition of strong return flows (hence exhumation), whereas zero or low F_R indicates either no return flows or a kinematic condition in which the return flows are deflected towards the overriding plate (for deformable overriding plate). A F_R=1 condition indicates that the wedge maintains volume conservation by allowing return flows to occur in volumes equal to those by burial-driven flows at any instant. For slab-parallel subduction, F_R can vary between 0 and 1. We, thus, choose a value of 0.5 (middle between the two extremities) as a reference to qualitatively express the return flow strength (strong vs. weak) in wedges. An addition of bulk horizontal shortening in the wedge due to an oblique U condition would result in an extra volume of wedge material extrusion (return), leading to a F_R>1 condition. For the F_R calculations, we discretize the wedge into 1000 equal small arc elements at r=10 km. Each arc covers an angle of ${\mathit{\theta }}_{\mathrm{1}}/\mathrm{1000}$. Since the velocity field in the corner flow solution is independent of r, F_R does not depend on the choice of r. The infinitesimal strain field at a point inside the wedge or the OP can be expressed by a stretching-rate tensor (Means et al., 1980) in polar coordinates (r, θ) as:

$$\begin{array}{}\text{(22)}& \left[\begin{array}{cc}{\dot{\mathit{\epsilon }}}_{\mathrm{r}r}& {\dot{\mathit{\epsilon }}}_{\mathrm{r}\mathit{\theta }}\\ {\dot{\mathit{\epsilon }}}_{\mathit{\theta }\mathrm{r}}& {\dot{\mathit{\epsilon }}}_{\mathit{\theta }\mathit{\theta }}\end{array}\right]=\left[\begin{array}{cc}\frac{\mathit{\delta }{u}_{\mathrm{r}}}{\mathit{\delta }r}& \frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\mathrm{1}}{r}\frac{\mathit{\delta }{u}_{\mathrm{r}}}{\mathit{\delta }\mathit{\theta }}+\frac{\mathit{\delta }{u}_{\mathit{\theta }}}{\mathit{\delta }r}-\frac{u_{\mathit{\theta }}}{r}\right)\\ \frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\mathrm{1}}{r}\frac{\mathit{\delta }{u}_{\mathrm{r}}}{\mathit{\delta }\mathit{\theta }}+\frac{\mathit{\delta }{u}_{\mathit{\theta }}}{\mathit{\delta }r}-\frac{u_{\mathit{\theta }}}{r}\right)& \frac{\mathrm{1}}{r}\frac{\mathit{\delta }{u}_{\mathit{\theta }}}{\mathit{\delta }\mathit{\theta }}+\frac{u_{\mathrm{r}}}{r}\end{array}\right]\end{array}$$

The eigenvalues and eigenvectors of the stretching-rate matrix (Eq. 22) give the values and the orientations of maximum and minimum instantaneous stretching axes, ISA_max and ISA_min respectively. The vorticity w in polar coordinates (r, θ) can be written as:

$$\begin{array}{}\text{(23)}& w={\displaystyle \frac{\mathrm{1}}{r}}\left[{\displaystyle \frac{\mathit{\delta }\left(r{u}_{\mathit{\theta }}\right)}{\mathit{\delta }r}}-{\displaystyle \frac{\mathit{\delta }{u}_{\mathrm{r}}}{\mathit{\delta }\mathit{\theta }}}\right]={\displaystyle \frac{\mathrm{2}(D\sin \mathit{\theta }-\mathrm{C}\cos \mathit{\theta })}{r}}\end{array}$$

For wedges with deformable OPs, the vorticity in the wedge and the OP is determined by replacing C, D of Eq. (18) by C₁, D₁ and C₂, D₂ respectively. Equation (18) shows that the vorticity is inversely proportional to the distance from the tip of the wedge. For our analysis, the vorticity at any point (r, θ) is normalized by the factor $U/(r{\mathit{\theta }}_{\mathrm{1}}$). Positive and negative values of the normalized vorticity parameter ($\dot{w}$) indicate anticlockwise and clockwise rotation, respectively. As mentioned in the preceding Sect. 2.2, the wedge kinematics and that in the overriding plate depend on several factors (e.g., ϕ, μ_r). We, however, present a parametric analysis to show the impact of the most influential factors in controlling the velocity and strain fields in the wedge. The present analysis covers a depth range between 20 and 40 km (Fig. 2a). Unless mentioned otherwise, we discuss in the foregoing sections the results for the case of accretionary wedges with a taper angle of 30° and viscosity of 10¹⁹ Pa s (Behr and Becker, 2018), keeping the OP fixed (i.e., U_op=0). The OP-wedge viscosity ratio (μ_r) is varied between 1 to 10⁴, implying that our modelling covers a weak (∼10¹⁹ Pa s) to an extremely strong (∼10²³ Pa s) overriding plate condition in the present analysis. Such a wide range of the OP viscosity is chosen considering that the strength of crustal materials can vary widely depending on the composition and thermal structure, among several other factors. According to available estimates, a cold, anhydrous, mafic overriding mass is found to be strong with viscosity $>\sim {\mathrm{10}}^{\mathrm{22}}$ Pa s, whereas a warm, hydrous, quartz-rich one is significantly weak with viscosity $<\sim {\mathrm{10}}^{\mathrm{20}}$ (Shinevar et al., 2015). Our modelling considers a subduction dip angle of 30°, which is within the ranges of mean shallow dips determined from several subduction zone transects (Hu and Gurnis, 2020; Lallemand et al., 2005). The wedge base is subjected to a velocity (U) of 3 cm yr⁻¹, as observed in many modern-day subduction zones (Doglioni et al., 2007). In natural settings, however, faults and shear zones may localize at the subducting slab-wedge contact, and U can thus effectively vary from a few millimetres to centimetres per year.

3.2 Parametric analysis

3.2.1 Role of subduction obliquity

Accretionary wedges with SP movements parallel to the basal wedge boundary (ϕ=0°) against a fixed rigid OPs localize return flows close to the upper boundary (Fig. 3a-i). F_R calculated from Eq. (21) is 1, implying that the crustal materials undergoing burial are potentially recycled back to the surface in the same quantities at any instant. The analytical solutions indicate that the degree of subduction obliquity strongly controls the flow patterns in wedges. Slab advance kinematics (i.e., negative ϕ) greatly facilitates the return flux, e.g., F_R=4.56 when $\mathit{\varphi }=-\mathrm{10}$°. The wedge materials extrude at high velocities (∼U) along the inclined wedge roof (Fig. 3b-i), and may even locally exceed the subduction velocity for high negative ϕ ($<\sim -\mathrm{15}$°). In contrast, the slab rollback kinematics (i.e., ϕ is positive) significantly reduces F_R, dropping it to zero for considerable rollback rates. At extremely slow rollback ($\mathit{\varphi }<\sim {\mathit{\theta }}_{\mathrm{1}}/\mathrm{3}$), the wedge still produces weak return flow. However, at relatively higher rollback velocities, all the wedge materials are buried at high rates, which may locally exceed the subduction velocity (Fig. 3c-i, d-i).

Figure 3 Velocity and vorticity fields within wedges with rigid walls in tectonic settings with varying subduction obliquities (ϕ). θ₁=30°, and α=30°. U_mag: flow velocity magnitude and $\dot{w}$: a normalized vorticity parameter. Short red lines mark the orientations of the minimum instantaneous stretching axes (ISA_min). Note that addition of a slab advance component substantially multiplies the return flux (F_R), whereas slab rollback lowers F_R.

The subduction obliquity also controls the viscous strain distributions in the wedge beneath a rigid OP. When ϕ has negative, zero or extremely low positive values ($\mathit{\varphi }<\sim {\mathit{\theta }}_{\mathrm{1}}/\mathrm{3}$), the wedge flows show a vorticity sense reversal, from negative (clockwise, thrust-sense shear) in the basal part to positive (anticlockwise, normal-sense shear) in the roof (Mandal et al., 2002). The upper part of the wedge is characterized by vertical to sub-vertical ISA_min axes (Fig. 3a-ii, b-ii). The field of normal-sense shear in the roof region progressively shrinks with increasing ϕ, and vanishes for high slab rollback velocities, corresponding to $\mathit{\varphi }>\sim {\mathit{\theta }}_{\mathrm{1}}/\mathrm{3}$. For a relatively large ϕ value ($>\sim \mathrm{2}{\mathit{\theta }}_{\mathrm{1}}/\mathrm{3}$), the basal part of the wedge produces normal-sense shear with sub-vertical ISA_min axes (Fig. 3d-ii). Conditions in the range $\sim {\mathit{\theta }}_{\mathrm{1}}/\mathrm{3}<\mathit{\varphi }<\sim \mathrm{2}{\mathit{\theta }}_{\mathrm{1}}/\mathrm{3}$, in contrast, produces thrust-sense of shear in the entire wedge domain (Fig. 3c-ii).

3.2.2 Control of overriding plate $/$ wedge viscosity ratio

Our analytical solutions indicate that, in case of a deformable OP, the viscosity ratio, μ_r can strongly influence the flow pattern in the wedge. μ_r≥10³ yields the first order flow pattern similar to that observed in the rigid OP setting. For μ_r=10³, the ϕ=0° subduction condition gives rise to a high return flux (F_R=0.7), which is, however, less than that (F_R=1) in the equivalent rigid OP setting (Fig. 4a-i). Negative ϕ leads to a condition of exceptionally strong return flux (e.g., F_R=3.43 when $\mathit{\varphi }=-\mathrm{10}$°), where the return flows attain high velocities (0.7U), comparable to the subduction rates (Fig. 5a-i). However, lower μ_r ($<\sim {\mathrm{10}}^{\mathrm{3}}$) weakens the confinement of viscous flows within the wedge, allowing the flows to become distributed in the OP region. In such conditions, the materials lying beneath the wedge-roof flow towards the overriding plate, instead of returning back towards the surface (Figs. 4b-i, c-i, 5d-i). This kinematic transformation results in significant weakening of the return flux within the wedge. For ϕ=0 and −10°, a decrease in μ_r from 10³ to 10² drastically reduces F_R from 0.7 to 0.06, and 3.43 to 0.55, respectively. At μ_r=10, the wedge hardly develops any return flux (Fig. 4c-i). For considerable slab rollback velocity (positive ϕ), similar to the case with rigid overriding plate, F_R drops to zero for all the viscosity ratios (Fig. 5b-i, c-i, e-i, f-i). The slab rollback process facilitates the burial rates of wedge materials, which may locally exceed the subduction velocity. The condition of relatively low μ_r (≤10²) sets in trench-ward movement of the overriding plate during the slab rollback (Fig. 5e-i, f-i).

Figure 4 Velocity and vorticity fields in wedges with deformable wall settings undergoing slab-parallel subduction. θ₁=30°, and α=30°. Decrease in wall to wedge viscosity ratio (μ_r) lowers F_R.

In a mechanical condition with μ_r=10³, negative, zero or small positive ϕ values result in reversal of shear across the wedge, showing a transition of shear sense from thrust in the basal part to normal beneath the wedge roof (Figs. 4a-ii, b-ii, 5a-ii, d-ii). For a given ϕ, decreasing μ_r causes the normal-sense shear field to shrink, and finally disappear when μ_r is extremely low (∼10), and the whole wedge is dominated by thrust-type shear (Fig. 4c-ii, see also Sect. S2 in the Supplement). For ${\mathit{\mu }}_{\mathrm{r}}\le \sim {\mathrm{10}}^{\mathrm{2}}$, the OP undergoes deformations in response to the shear stress generated at the wedge-slab interface (Figs. 4b-ii, c-ii, 5d-ii, f-ii). The addition of slab advance kinematics induces strong sub-horizontal shortening, combined with thrust shear in the OP. On the other hand, the slab rollback, as shown earlier for the wedge with rigid OP, the entire wedge is affected by thrust shear when ϕ is neither extremely low nor high (Fig. 5b-ii, e-ii). In contrast, higher ϕ conditions develop normal shear with sub-vertical ISA_min axes in the trench-ward side of the wedge (Fig. 5c-ii, f-ii). For moderate to low viscosity ratios (${\mathit{\mu }}_{\mathrm{r}}\le \sim {\mathrm{10}}^{\mathrm{2}}$), the slab rollback is associated with considerable sub-horizontal extension (ISA_max) in the overriding plate, away from the wedge (Fig. 5f-ii).

Figure 5 Effects of subduction obliquity (ϕ) and viscosity ratio (μ_r) on the velocity and vorticity fields in wedges. θ₁=30°, and α=30°.

3.2.3 Influence of taper angle

Accretionary wedges can be narrow or wide, depending on the sediment influx and outflux of sediments into or from the wedge (Clift and Vannucchi, 2004; Cloos and Shreve, 1988) and the effects of other factors, such as subduction geometry (Pusok et al., 2022). Our analytical solutions show that increasing θ₁ strengthens the return flux in a wedge (Fig. 6a-i, b-i, c-i). For ϕ=0° and μ_r=10², an increase in θ₁ from 15 to 120° increases F_R from 0 to ∼1. Narrow wedges with ${\mathit{\theta }}_{\mathrm{1}}<\sim \mathrm{15}$° do not produce any return flows (Fig. 6a-i) when μ_r is low to moderate ($\le \sim {\mathrm{10}}^{\mathrm{2}}$). In such conditions, the flows fails to redirect wedge materials back to the surface as they are strongly deflected towards the OP, which becomes more pronounced when the plate convergence involves slab advance (see Sect. S1). This kinematics implies that the wedge–OP interface will progressively migrate rearward and result in distortion of the initial wedge geometry. Such an “unstable” state of wedges with low taper angles, and weak walls has been discussed in detail by Moulas et al. (2021).

Large-taper wedges (θ₁=120°) develop thrust shearing in a bi-vergent configuration, as illustrated in Fig. 6c-ii. The oppositely verging thrust-shear zones are separated by a zone of normal-sense shearing. In contrast, narrow wedges with moderate to low μ_r ($<\sim {\mathrm{10}}^{\mathrm{2}}$) never produce a condition of shear reversal, but allow the OP to accommodate considerable deformations (Fig. 6a-ii).

3.2.4 Role of the OP velocity

The overriding plates often experience a trench-ward velocity (U_op) due to far-field tectonic stresses or drag of the underlying mantle flows (Cerpa and Arcay, 2020; Jarrard, 1986; van Dinther et al., 2010). We use the present analytical solutions to investigate the effects of U_op on wedge flows. In case of no slab rollback, addition of U_op strengthens the convergence within the wedge, resulting in extrusion of the wedge materials from its roof regions to the surface (Fig. 6d) (Nettesheim et al., 2018). For relatively lower rollback velocities (e.g., ϕ=15°), the wedge shows curvilinear flow trajectory with upward convexity (Fig. 6e), as observed in some of the numerical models (e.g., Candioti et al., 2021), implying that crustal materials under a steady state kinematic state of the wedge would be extruded for a small extent, and then recycled back to deeper regions. High rollback velocities, in contrast, facilitate the burial process of wedge materials with the subducting slab, allowing their little or no subsequent uplifts (Fig. 6f).

Figure 6 (a–c) Influence of taper angle (θ₁) on the kinematics of wedges in settings undergoing slab-parallel subduction. (d–f) Effect of trench-ward moving overriding plate movement ($U_{\mathrm{op}}=U/\mathrm{3}$) on the flow patterns of wedges: μ_r=10², and α=30°.

3.2.5 Effects of subduction dip

For a constant dip of the wedge–OP interface, decreasing slab dips give rise to increasing taper angles, which in turn facilitate the return flux (F_R) within the wedge. A steep dip of the wedge–OP interface yields higher extrusion rates (higher vertical component of velocity in Cartesian xy-frame) of the extruding wedge materials (see also Sanhueza et al., 2022). However, steeply dipping slabs generally exhibit high rollback velocities (hence increase in ϕ) due to slab pull force (Stegman et al., 2010), which in turn affect the wedge kinematics. From a dynamic point of view, lower subduction dips allow the convergent system to transmit higher stresses to the wedge and the overriding plate, and give rise to intense deformations in the forearc region (Chemenda et al., 2000; Guillaume et al., 2009).

4 Observations in analogue experiments

The analytical solutions described above provide the instantaneous kinematic behaviour of a wedge for a given geometry and under a specific set of boundary conditions. To validate the analytical results, and evaluate the wedge kinematics in finite timescale, a series of scaled laboratory experiments were performed. The experiments allow us to decipher the evolution of velocity fields in wedges with deformable walls and varying boundary kinematics, as described in the preceding sections.

4.1 Model setup

The model accretionary wedge was made up of polydimethylsiloxane (PDMS), a Newtonian viscous fluid (viscosity ∼10⁴ Pa s) at strain rates below 10⁻² s⁻¹, and a density of ∼965 kg m⁻³ (Weijermars, 1986). The wedge was ∼ 5.6 cm thick, and had a taper angle of 45° Fig. 7a). One set of models had overriding plates of high-viscosity (∼ 10⁶ Pa s) modelling clay-PDMS mixture (4:1 ratio), while the other set had overriding plates of a relatively low-viscosity (∼ 10⁵ Pa s) mixture of the same materials (1:4 ratio). The wedge had a basal dip of 45°, and a vertical boundary with the overriding plate. The model had a length of 20 cm in y direction, and its upper surface was set free. This configuration develops a deformation regime close to plane strain condition with negligible deformation along the y direction, especially in the earlier stages of the experiments. The deformation box walls were lubricated with soap to reduce the friction between the viscous materials and glass walls.

In a reference frame fixed to the base of the deformable wall, a velocity was applied either parallel or oblique to the wedge base to simulate slab-parallel subduction and slab advance/rollback kinematics, respectively (Fig. 7b, c) (details of the kinematic boundary conditions provided in Sect. S3). During an experimental run, photos of the model vertical section were taken at regular time intervals, which were post-processed to obtain incremental displacement fields using the PIVLab (Thielicke and Stamhuis, 2014), an open-source image correlation software run in MATLAB.

Figure 7 (a) Laboratory model setup used for wedge experiments (details provided in the text). Dashed rectangle shows the region of interest. (b) Velocity conditions imposed at the model boundaries to simulate slab-parallel subduction, slab advance, and slab rollback in the models. (c) Boundary velocity conditions with respect to a fixed basal buttress.

4.2 Scaling

Analog models were scaled to natural settings to obtain length, dynamic, rheological, and kinematic similarities (Schellart and Strak, 2016; Ramberg, 1981). A length ratio l_r ($l_{\mathrm{model}}/l_{\mathrm{nature}}$) of $\mathrm{1.43}\times {\mathrm{10}}^{-\mathrm{6}}$ was used in the model. Under this scaling, the 5.6 cm thick wedge model represented a ∼40 km thick viscous accretionary wedge. We, however, limit our analysis to the lower 20 km (equivalent to ∼3 cm in model) part of the wedge as natural accretionary wedges mostly exhibit brittle deformation in the upper crustal regime. The model viscous wedge was made thicker to minimize the mechanical influence of the model's upper “free surface” in the region of interest (i.e., between a depth of 20 to 40 km), especially during the early stages of experiments. The experimental results of the whole model wedge is provided in Sect. S4.

The subduction in all the models was driven by moving a plastic sheet at a velocity (V_model) of 0.7 mm min⁻¹ beneath wedge base, which ultimately passed through a slot beneath the model (Fig. 7). Considering a basal velocity (V_nature) of 1 cm yr⁻¹ in geological settings, we obtain a velocity scaling ratio (V_r) of ∼36 792, and a time ratio ($t_{\mathrm{r}}=l_{\mathrm{r}}/V_{\mathrm{r}}$) of $\sim \mathrm{3.89}\times {\mathrm{10}}^{-\mathrm{11}}$. This relatively lower value of V was chosen as a part of the subduction velocity may be partitioned in faults or shear zones at the wedge-subducting plate contact. However, the analytical results suggest that the choice of V hardly influences qualitatively the flow patterns in CF settings. Assuming a depth-averaged density of 2800 kg m⁻³ for accretionary wedge rocks, the density scaling ratio ρ_r was set at ∼0.36. As the experiments were run under normal gravity conditions, the scaling ratio for gravitational acceleration (g_r) is 1. These scaling factors yields a viscosity scaling ratio, ${\mathit{\mu }}_{\mathrm{sr}}={\mathit{\rho }}_{\mathrm{r}}\times l_{\mathrm{r}}\times t_{\mathrm{r}}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{17}}$. Our model wedge is thus scaled to a natural accretionary wedge with a relatively higher depth-averaged viscosity of $\sim \mathrm{5}\times {\mathrm{10}}^{\mathrm{20}}$ Pa s. However, for a given OP–wedge viscosity ratio, this absolute wedge viscosity (μ₁) hardly influences the corner flow pattern, as extrapolated from the analytical solutions. μ₁ can only play a role in determining the intensity of gravity-driven flows that occur during later stages of wedge deformation.

4.3 Experimental results

4.3.1 Slab-parallel subduction (ϕ=0°)

This model setting (ϕ=0°) develops strong return flows in the wedge when the wall to wedge viscosity ratio is large (μ_r≈10²) (Fig. 8a-ii). Materials close to the moving wall are dragged down, while materials in the rear side are uplifted. During early stages of the experiments, the maximum vertical velocity (v_max) across the wedge is about 0.5U, where U is the basal velocity. A lower viscosity ratio (μ_r≈10) results in the return flow to strongly deflect towards the overriding plate (Fig. 9a-ii), and substantially lowers v_max within the wedge to 0.25U. The deflected flow induces considerable deformations and uplift (at rates ∼0.23U) in the wall close to the wedge. These flow patterns are similar to the analytically predicted velocity fields (Figs. 8a-i, 9a-i), especially close to the wedge corner region. Away from the corner, the velocity vectors in the return flow become slightly steeper than those predicted by the theoretical solution. This small difference arises from the influence of the “free upper surface”, which is located at a finite distance above the corner in laboratory models, whereas the analytical model assumes a wedge that extends infinitely upward from the corner. In μ_r≈10 experimental models, the deflected return flows progressively distort the wedge geometry by forming a bulge into the wall near the corner region (Fig. 9a-ii). This bulging results from dynamic pressure localization at the wedge corner, which is strong enough to deform the relatively weak overriding plate (see Sect. S1). However, the wedge grossly retains its initial wedge geometry in models with relatively strong overriding plates (μ_r≈100). These experimental findings are in agreement with the numerical models of Moulas et al. (2021), showing that the corner geometry becomes unstable with lowering of μ_r.

.

Figure 8 Evolution of wedge kinematics in μ_r≈10² experimental models: (a) slab-parallel subduction (ϕ=0°), (b) slab-advance ($\mathit{\varphi }=-\mathrm{8}$°), and (c) slab rollback (ϕ=30°) settings. Left panel: Instantaneous velocity field predicted from the analytical solutions with similar setups shown as a reference. Central (ii), and right (iii) panels: Model flow patterns (in a reference frame fixed to the base of the deformable wall) at two different stages of the experiments. The solid and dashed red lines demarcate the present boundaries of the wedge and their initial positions, respectively

Figure 9 Evolution of wedge kinematics in μ_r≈10 experimental models: (a) slab-parallel subduction (ϕ=0°), (b) slab-advance ($\mathit{\varphi }=-\mathrm{8}$°), and (c) slab rollback (ϕ=30°) settings. Left panel: Instantaneous velocity field predicted from the analytical solutions with similar setups shown as a reference. Central (ii), and right (iii) panels: Model flow patterns (in a reference frame fixed to the base of the deformable wall) at two different stages of the experiments. The solid/dashed red lines same as in Fig. 8.

With progressive subduction, the model wedge develops an upward-convex topography at the wedge surface with steep slope at the wedge front (Figs. S3a-ii, S4a-ii in the Supplement). The topography stops growing in the later stages of the experiments, and eventually gives rise to gravity-driven flow from the rear region toward the wedge front. In this stage, the concurrence of extrusion flows at depth, and their subsequent collapse in the shallow region lead to a closed flow path in the upper parts of the wedge (Figs. S3a-ii, S4a-ii), as observed in several previous analytical (e.g., Emerman and Turcotte, 1983; Platt, 1993), and numerical models (e.g., Van Dinther et al., 2012; Allemand and Lardeaux, 1997; Yamato et al., 2007). The effect of this closed flow pattern extends even to the bottom 20 km of the wedge in the later stages of experiments (Figs. 8a-iii, 9a-iii). For low μ_r, the trenchward gravity-driven flow within the wedge also drags the overriding-plate material, resulting in change in the return-flow kinematics from earlier stages of the experiment. This effect is especially evident away from the wedge corner, where the previously sideways-deflected return-flow vectors reorient toward a more vertical orientation (compare Figs. 9a-ii and 9a-iii). At this stage, the wedge–OP interface migrates rearward at depth close to the corner, but migrates towards the wedge front away from the corner, resulting in progressive reduction in taper angle of the wedge (Fig. 9a-iii).

In natural settings, surface erosion and infilling of sediments that can smoothen out the topographic slope in the wedge, and are thus expected to reduce the late-stage gravity-driven flow. We ran an additional set of experiments with μ_r≈10 by adding new materials to the evolved wedge surface after 45 km shear movement, restoring grossly its horizontal initial orientation (see Sect. S5, Fig. S5). Further shearing in this modified model setting gives rise to return flows strongly deflecting towards the overriding wall, similar to that observed in the initial stage of the experimental run.

4.3.2 Slab advance ($\mathit{\varphi }=-\mathrm{8}$°)

Slab advance kinematics facilitated the uplift process at rapid rates in the rear part of the wedge, as predicted by the analytical solutions (Figs. 8b, 9b). During the initial stages, the maximum extrusion velocity v_max within the wedge attains values of ∼0.8U, and ∼0.4U in the strong and weak overriding plate models, respectively (Figs. 8b-ii, 9b-ii). For lower μ_r, the strongly deflected return flows induce significant deformations and uplift (at rates ∼0.16U) in the overriding plate.

Progressive slab advancement leads to bulk horizontal squeezing of the wedge, driving the process of vertical extrusion. After 50 km of shear and 10.5 km of shortening, wedge materials at depth thus experience uplift at relatively high rates of ∼0.6U and ∼0.4U for high and low μ_r experiments, respectively (Figs. 8b-iii, 9b-iii). In the later stages of experiments, the advancing slab comes in direct contact with the overriding plate (Fig. 8b-iii), forcing uplift (at rates ∼0.13U) and deformation even in the stronger wall (μ_r≈10²) of the wedge. For low μ_r, the wedge–OP interface is strongly deflected sidewise towards the overriding plate, and distorted near the corner (Fig. 9b-iii).

4.3.3 Slab rollback (ϕ=30°)

Slab Rollback facilitates burial of materials in the entire wedge, leaving return flows completely absent in the system. This experimental finding agrees with the analytical solutions (Figs. 8c, 9c). The downward flows during rollback cause progressive vertical subsidence, accompanied by bulk horizontal extension of the wedge. Low μ_r (≈10) model induces this extension to occur even in the weak overriding wall (Fig. 9c-ii).

With progressive rollback, the initial narrow geometry of the wedge undergoes significant modifications (Figs. 8c-iii, 9c-iii), but materials continue to flow down with maximum velocity in proximity to the moving wall. This flow pattern gives rise to a steep surface slope close to the moving wall (see Figs. S3, S4), which eventually triggers additional gravity-driven sub-horizontal flows of wedge materials from the rear to the frontal side in later stages of the experimental run. Effect of this sub-horizontal flow is also observed at the bottom 20 km of the wedge. This results in shallowing of the previously steep velocity vectors during the later stages (Figs. 8c-iii, 9c-iii). That effect is more prominent in models with lower μ_r, where the magnitude of negative dynamic pressure, which drives the corner-flow kinematics during rollback, is relatively less (see Sect. S1). In natural subduction zones, however, opening at the base of the wedge during slab rollback can lead to entrainment of the wedge materials into the mantle. In such cases, the flow pattern at the bottom parts of the wedge will differ from the results observed in our models with a rigid base.

4.4 Analogue and analytical model results: a synthesis

Both analogue experiments and analytical solutions suggest that the following parameters: slab obliquity (ϕ), overriding plate-wedge viscosity ratio (μ_r), wedge taper angle (θ₁), and overriding plate velocity (U_op) largely determine the intensity of return flows in viscous accretionary wedges.

For slab-parallel subduction (ϕ=0°), a high μ_r ($>\sim {\mathrm{10}}^{\mathrm{2}}$) condition favours deep-crustal material extrusion, forming strong return flows in the rear side of the wedge (Figs. 4a, 8a). Lowering of μ_r ($<\sim {\mathrm{10}}^{\mathrm{2}}$) weaken the return flows, while facilitates lateral deflection of the extruding materials, resulting in intense deformations in the overriding plate (Figs. 4c, 9a). In contrast, slab advance (negative ϕ) leads to bulk horizontal shortening, resulting in enhanced vertical extrusion (Figs. 5a, d, 8b). Analogue model results show that after finite advance, the subducting slab comes directly in contact with the overriding plate, inducing significant uplift and deformation in the overriding plate, irrespective of μ_r. On the other end, slab rollback ($\mathit{\varphi }>\sim {\mathit{\theta }}_{\mathrm{1}}/\mathrm{3}$) promotes the burial process across the entire wedge without any return flow (Figs. 5b, c, e, f, 8c, 9c). Under a low μ_r condition, this rollback process can also induce extensional deformations in the overriding plate. Analogue model results further indicate that that gravity-driven flow may become important with progressive deformation of the low-viscosity wedge materials. Analytical solutions reveal that for low to moderate μ_r, wedges with small taper angles suppress return flows, while wider wedges facilitate it. Trenchward movement of the overriding plate also enhances the upward extrusion process.

The subduction settings with slab-parallel subduction, slab advance, or highly oblique ($\mathit{\varphi }>\sim \mathrm{2}{\mathit{\theta }}_{\mathrm{1}}/\mathrm{3}$) rollback results in reversal of shear sense across the wedge for moderate to high μ_r (Figs. 4a, 5a, c, d, f). However, the sense of shear remains unchanged in cases of low μ_r, or less-oblique rollback (Figs. 4c, 5b, e).

5 Discussion

5.1 Conditions for return flows in viscous accretionary wedges

The analytical, and analogue model results discussed above suggest that the return flows in accretionary wedges, which play a crucial role in deep-crustal exhumation processes can occur under specific conditions. The role of each factor in setting the return flows is summarized in Fig. 10. It can be readily seen that appreciable return flows (F_R>0.5) occur in the following conditions: little or no slab rollback (ϕ<5°), high overriding plate-wedge viscosity ratio (${\mathit{\mu }}_{\mathrm{r}}\ge \sim {\mathrm{10}}^{\mathrm{3}}$), and moderate to high taper angle (${\mathit{\theta }}_{\mathrm{1}}\ge \sim \mathrm{30}$°). Wedges with large taper angles produce strong return flows, independent of the OP viscosity. Narrow wedges, i.e., low θ₁ (e.g., 15°) become almost devoid of any focused return flows, unless the OP is rigid or has high viscosity (μ_r>10³). A strong overriding plate accommodates little strains, forcing the viscous flows to localize within the wedge. Increase in μ_r thus promotes the return flux (Fig. 10b). The return flow is also facilitated by the slab advance (ϕ<0°), but countered by the slab rollback process in subduction zones (Fig. 10a). In presence of a stronger OP (${\mathit{\mu }}_{\mathrm{r}}\ge \sim {\mathrm{10}}^{\mathrm{3}}$), additional squeezing of the wedge during slab advance can lead to a F_R>1 condition. These contrasting flow kinematics originate from the competing dynamics of two processes: dominance of a downward flow, induced by the subducting slab drag, and an upward flow, driven by a negative upward dynamic pressure gradient along the downward-tapering wedge (Fig. 1c). The first process enables the burial mechanism of the wedge materials, whereas the second process forces the materials to extrude up, setting in return flows. The analytical results suggest that deformable OP settings with high μ_r, low ϕ, and large θ₁ develop strong negative upward dynamic pressure gradients along the wedge (Fig. S1), which leads to the generation of return flows. For small θ₁, and low μ_r, dynamic pressure gradient across the wedge forces the return flow to deflect towards the wall (see Sect. S1). Progressive increase in ϕ lowers this pressure gradient and then reset its trend in the reverse direction (Fig. S2, see also Fig. 13). This finding predicts the absence of return flows in wedges where the slab rollback is active. Trench-ward OP motion also facilitates the extrusion process in a wedge (Fig. 6d). Other geologically relevant factors, which are excluded in the present study, such as presence of buoyant rocks, normal faulting, and erosion can promote return flows of deep-crustal materials in accretionary wedge materials (see Warren, 2013, for a review).

Figure 10 Analytically determined plots of return flux (F_R) as a function of (a) subduction obliquity (ϕ), and (b) viscosity ratio (μ_r), for varying taper angles (θ₁) of the wedge. Note that for low to moderate taper angles (θ₁≤30°), significant return flux (F_R>0.5) occurs when μ_r is large (≥∼10³), and there is no slab rollback (ϕ0° or negative).

5.2 Exhumation of high P metasediments by return flows: model vs. nature

Subduction zones generally host distinct types of high-pressure low-temperature (HP–LT) metamorphic rocks: metamorphosed oceanic sediments, oceanic crust (+ mantle), and continental crust, each of them records its own exhumation process (Agard et al., 2009, 2001; Burov et al., 2001; Gerya et al., 2002; Angiboust et al., 2012). The exhumation of oceanic sediments in accretionary wedges usually occurs as a syn-convergent steady process operating on tens of million years, whereas oceanic crustal fragments record episodic exhumation events corresponding to different stages of subduction (Husson et al., 2009; Yamato et al., 2007; Wang et al., 2019). They also show varying modes of occurrences in the subduction settings. High pressure metasediments occur as extensive, continuous metasedimentary outcrops with smooth, progressively increasing P-T conditions (Agard et al., 2001; Ernst, 1993; Vaughan-Hammon et al., 2022). On the other hand, exhumed oceanic crustal bodies form distinct slices of units, often associated with serpentinites, and exhibit contrasting P-T peaks, as reported from the Western Alps (Agard et al., 2009; Schwartz et al., 2001; Angiboust et al., 2012). Previous studies have used the corner flow mechanism to explain the exhumation of high-pressure metasedimentary rocks in terms of return flows in accretionary wedges (Allemand and Lardeaux, 1997; Gerya et al., 2002; Shreve and Cloos, 1986). The present theoretical analyses suggest that accretionary wedges can produce such return flows only in specific conditions determined by a number of parameters, such as subduction obliquity, OP-wedge viscosity ratio, and taper angle of the wedge, as shown in Fig. 10.

The relative volume of exhumed high-pressure metasedimentary rocks varies significantly in the accretionary wedges of convergent settings (Kerswell et al., 2023). Agard et al. (2009) have provided estimates recycled volume budgets of subducted sediments, ∼80 %–90 % in the Cascades, whereas 30 %–50 % in the Western Alps. However, the estimates yield remarkably low values, <1 % in the Chile and Franciscan complex. Although low preservation potential of blueschist facies rocks or differential erosional effects are possible factors for such variations, the flow kinematics in the accretionary wedge seems to be one of the most influential factor in determining the amounts of uplifted masses. Our theoretical analyses reveal that extensive return flows can occur in specific kinematic, rheological and geometrical conditions. For example, in the Cascades and the Western Alps, the ratio of returned to subducted material volumes at any instant within a steady state accretionary wedge (F_R) has been estimated high (>0.5) (Agard et al., 2009). The present theoretical results suggest the following: (1) subduction with little or no slab rollback, (2) high relative strength of the overriding plate (${\mathit{\mu }}_{\mathrm{r}}\ge \sim {\mathrm{10}}^{\mathrm{3}}$), and (3) large wedge taper angles (θ₁) as potential conditions favoring such high return flow kinematics. For small taper angles (${\mathit{\theta }}_{\mathrm{1}}\le \sim \mathrm{30}$°) and weak overriding plates (${\mathit{\mu }}_{\mathrm{r}}\le \sim {\mathrm{10}}^{\mathrm{2}}$), the settings hardly return materials within the wedge, the quantity of which can be further reduced by slab rollback (Figs. 4b-i, 5b-i, 9a, c). For relatively weak overriding plates (OP), the return flows originating from deeper regions eventually deflect towards the overriding plate, resulting in exhumation of both the wedge and OP materials in the retroside regions (Butler et al., 2011), at extremely slow rates (<0.1U for a setting with μ_r=10², and θ₁=30°) (Fig. 11a). However, such sluggish uplift processes do not support them to travel much in the upward direction. Extremely weak (μ_r<10) OP settings, however, produce flows closely similar to those formed in a single (i.e., μ_r=1) corner model with a large taper angle, enabling the crustal materials to uplift at considerable velocities and accumulate in the retroside of the wedge. Slab rollback at a relatively higher velocity (keeping the OP fixed) suppresses the return flows, instead it facilitates leaking of crustal materials through a narrow passage into the mantle (Figs. 3c-i, d-i, 4b-i, c-i, e-i, f-i). However, the same tectonic setting favours exhumation in other rocks, such as oceanic crustal units, low-density continental blocks and mantle rocks (Brun and Faccenna, 2008; Husson et al., 2009; Wang et al., 2019), as the slab rollback and associated upper plate extension create free spaces within the lower- or upper plate-side of the fore-arc, setting a dynamic condition for the exhumation of (U) HP rocks from mantle depths. On the other hand, our model results suggest that a period of slab advance during subduction can facilitate the extrusion process in the wedge, resulting in exhumation of deep-crustal materials (Figs. 3b-i, 5a-i, d-i, 8b).

Figure 11 Vertical-velocity profile (at y=10 km) across the wedge and the overriding plate in analytical models with different θ₁ (a–d), μ_r, and ϕ. Positive and negative values of the vertical velocity indicate uplift/exhumation and subsidence/burial of wedge materials, respectively.

The calculated extrusion velocities of HP-LT rocks in accretionary belts vary on a wide range, from <1 mm yr⁻¹ to >1 cm yr⁻¹ (see Table 1b of Agard et al., 2009, and references therein), but they mostly lie in the order of mm yr⁻¹. In contrast, previous CF model studies yield extrusion velocities much higher than 1 cm yr⁻¹ (Allemand and Lardeaux, 1997; Cloos, 1982; Gerya et al., 2002), which perhaps occur due to the geometrical and kinematic boundary conditions chosen in the models. The present theoretical analysis, however, suggests that the return flows can be significantly slow due to specific mechanical and geometric conditions: low relative viscosities of the overriding plates (${\mathit{\mu }}_{\mathrm{r}}<\sim {\mathrm{10}}^{\mathrm{3}}$) and small taper angles (${\mathit{\theta }}_{\mathrm{1}}\le \sim \mathrm{30}$°) (Fig. 11a, b). For θ₁=30°, accretionary wedge models with μ_r=10³ and 10² yield maximum vertical uplift velocities (v_max) about one-fourth and one-eighth of the wedge basal velocity (U), respectively (Fig. 11a). As the material particles track curvilinear paths, they would be extruded mostly at velocities less than v_max, and their time-averaged velocity would be less than v_max. Our theoretical estimates are consistent with those (1–5 mm yr⁻¹) calculated from the Schistes Lustres complex in the Western Alps ($U=\sim \mathrm{1}$–2 cm yr⁻¹, Yamato et al., 2007) and other accretionary belts (Brandon et al., 1998; Glodny et al., 2005; Ring et al., 1999). The theoretical results further indicate that syn-subduction slab rollback (ϕ>0°), even small in magnitudes can again weaken the return flow kinematics, reducing the exhumation velocities to a large extent (<1 mm yr⁻¹; Fig. 11a–d). Additionally, discrete slip along the wedge-subducting plate interface can also decrease the exhumation velocity (v), as v is directly proportional to the velocity (U) at the base of the wedge. Other mechanisms, such as gravitational spreading at shallow crustal level (Van Dinther et al., 2012), deformation localization, syn-subduction sediment accumulation can also contribute to a reduction in the exhumation velocity. Therefore, in contrast to the prevailing notion that the extremely slow exhumation rates result solely from erosion (Brandon et al., 1998; Platt, 1993) and other reasons discussed above, our model findings suggest that corner flow processes can also modulate the kinematics of extrusion of deep HP rocks to the surface at slow rates.

High-pressure rocks in many convergent belts record unusually high uplift rates. As an example, the Ampelos/Dilek nappe of the Cycladic blueschist unit in the eastern Aegean underwent extrusion at rates 3–3.5 cm yr⁻¹ during incipient collision of the Anatolian microcontinent with Eurasia (Ring et al., 2007). Our model results suggest that mechanically strong overriding plates (${\mathit{\mu }}_{\mathrm{2}}>\sim {\mathrm{10}}^{\mathrm{22}}$ Pa s) and/or significant slab advance (ϕ<0°) kinematics form potential convergent settings for high-rate extrusion in wedges with moderate tapering angles (Figs. 8b, 11a, b). In addition, large crustal scale normal-sense shear zones (or faults) within the wedge or at the wedge-overriding plate contact can exhume high P rocks at fast rates (Bento dos Santos et al., 2021; Searle, 2015). Presence of buoyant crustal blocks and erosion can further facilitate the extrusion kinematics (England and Holland, 1979; Ernst et al., 1997; Godin et al., 2006). The present theoretical estimates also account for relatively higher extrusion rates (∼1 cm yr⁻¹) in wider wedges with large taper angles (Fig. 11c, d).

5.3 Deformation patterns

Accretionary wedge deformations have generally been explained using the critical taper theory in the framework of Coulomb rheology (Dahlen, 1990; Dahlen et al., 1984; Buiter, 2012). However, this rheological consideration may not provide good approximations for wedges thicker than ∼15 km as thermally activated ductile deformation mechanisms become significant in controlling the rheology of the wedge materials at depths (Pavlis and Bruhn, 1983; Platt, 1986). The depth of brittle-ductile transition is strongly controlled by the prevalent thermal regime (Stewart et al., 2000). Below this depth, viscous flow regimes govern the deformation patterns across the wedge (Fig. 2a).

The model results presented in the preceding section indicate that the deformation styles in convergent settings strongly depend on a number of geometric and kinematic parameters: slab motion, overriding plate's relative viscosity (μ_r), and wedge taper angle. Three distinct modes of deformation can be recognized in the field diagram of Fig. 12, each with a characteristic vorticity pattern. The central part of the diagram defines a field of no shear reversal, i.e., the whole wedge is characterized by the same shear sense, as imposed by the SP. The central field is bounded by two lines, the locus of which varies with μ_r. To the left and right of this field, the shear sense reverses in the upper (i.e., near the OP) and lower (i.e., near the SP) regions of the wedge, respectively. The variation in the wedge kinematics is essentially a consequence of the interplay between the slab-driven Couette flow and the dynamic pressure-driven Poiseuille flow (see Sect. S2 for a detailed discussion). For moderate to high μ_r ($\ge \sim {\mathrm{10}}^{\mathrm{2}}$), subduction, with or without slab advance (i.e., zero or negative ϕ), the wedge settings produce two synchronous domains of shear deformations with opposite shear senses: thrust shear in the lower part and normal shear in the upper part of the wedge (Figs. 4a-ii, b-ii, 5a-ii, b-ii) . The shear reversal is accompanied by a systematic variation of the strain fields, with sub-horizontal principal shortening axes close to the subducting plate, which reorient to become sub-vertical towards the overriding plate. Wider wedges, characterized by oppositely-dipping walls transforms the shear patterns into doubly-vergent thrust sense shears, separated by a zone of normal-sense shear in their axial regions (Deville, 2023; Luján et al., 2010) (Fig. 6c-ii). However, the shear reversibility disappears if μ_r is low ($\sim <{\mathrm{10}}^{\mathrm{2}}$) and the slab retreat (i.e., positive ϕ) becomes active. For a moderately tapering wedge, the later can give rise to thrust-sense shearing in the entire wedge (Figs. 3c-ii, 4c-ii, 5b-ii, e-ii), with the exceptions that high slab rollback velocities develop normal-sense (with subhorizontal extension) and thrust-sense shearing close to the subducting and the overriding plate, respectively (Figs. 3d-ii, 5c-ii, f-ii). Similar variations in the shear kinematics across the wedge are also evident from the deflection of marker lines after finite deformation in the experimental models (Fig. 12i–iv). The analogue models also suggest that after considerable deformation, the topographic slope generated at the upper surface of relatively low viscosity wedge sets in gravity-driven flow that leads to horizontal extension in the rear part of the wedge, as also reported from previous analogue experiments (Chattopadhyay and Mandal, 2002; Luján et al., 2010).

Figure 12 Fields of “shear reversal” and “no shear reversal” across the wedges in a parametric space defined by subduction obliquity (ϕ), taper angle (θ₁), and viscosity ratio (μ_r). To the left and right sides of the central “no shear reversal” field, shear sense reverses differently, as shown in the analytically determined vorticity regimes at various positions (I–V small red squares) of the field diagram for μ_r=10³. Short red dashes in I–V mark the orientations of ISA_min. Red solid lines within the wedges (i–iv) indicate deflection of initial straight markers (red dashed lines) after 10 km of shearing in experimental models (μ_r=10²) for varying subduction obliquities (see small blue squares).

The shear reversal phenomenon has been reported from extruded rocks in several subduction zones (e.g., Hellenides, Ducharme et al., 2022; Western Alps, Agard et al., 2001; Sambagawa Belt, Osozawa and Wakabayashi, 2015) and collisional belts (e.g., the Cycladic blueschist unit in the Eastern Aegean, Ring et al., 2007; Himalayan mountain belts, Beaumont et al., 2001; Grujic et al., 2002; Dutta and Mukherjee, 2021, 2019; the Porto–Viseu Metamorphic Belt in the Central-Iberian Zone, Bento dos Santos et al., 2021). Considering a steady state flow kinematics, crustal materials following their burial to greater depths are recycled back to the surface by return flows and experience shear deformations in both lower and upper domains of the wedge. Consequently, they can record thrust-sense shear structures in them, overprinted by normal-sense shear structures during exhumation, as observed in the Cycladic Blueschist Unit at Fabrika (Ring et al., 2020). Such corner-flow driven complex tectonics can produce deformational structures under sub-horizontal shortening in the burial phase, which are superposed by structures related to sub-vertical shortening. Multiple structural fabrics, and fluid-assisted deformations develop in the course of a progressive tectonic event, without any change in the external kinematic frame (Richter et al., 2007; Muñoz-Montecinos and Behr, 2023). The strain fields at the time of exhumation indicate sub-vertical crustal shortening and tectonic thinning (Feehan and Brandon, 1999; Ring and Brandon, 2008; Kryl et al., 2021), manifested in sub-horizontal foliations, which are superposed on earlier tectonic structures, such as folds and fabrics in several paleo-accretionary wedges e.g., San Juan-Cascade nappes, coastal accretionary belt of Chile (Richter et al., 2007; Muñoz-Montecinos et al., 2020). Crustal scale extensional (normal) shear zones that localize in the upper part of the wedge largely control the exhumation of high grade metamorphic rocks beneath low grade rocks in subduction and collisional settings (Searle, 2015). The extensional shearing in the deeper region gives way to normal faulting at the shallow crustal level in the hinterland as the return flows result in vertical crustal uplift in the rear of the wedge (Figs. 11, S3, S4). This uplift develops normal faulting to retain the stability (Platt, 1986). The normal faults further assist in the exhumation of high P materials, especially at shallow crustal depths (Platt, 1993; Ring et al., 1999).

In subduction zones, the subducting oceanic slabs can induce strong deformations in their overriding plates, often manifested in characteristic topographic developments (Buiter et al., 2001; Dasgupta and Mandal, 2018; Hampel and Pfiffner, 2006; Shijie Zhong and Gurnis, 1994). Previous studies suggest that the forces responsible for such OP deformations originate in several ways, viscous drag at the base of the overriding plate, suction effects between the subducting and the overriding plate, and shear generated at the slab interface (Chen et al., 2017; Schellart, 2024, and references therein) A detailed discussion of their independent role is beyond the scope of this study. However, our theoretical analysis suggests that the shear force generated by the subducting plate can induce large crustal deformations in the overriding plate only under specific conditions constrained by the overriding plate viscosity, and the subduction kinematics (slab advance vs. rollback). High μ_r ($>\sim {\mathrm{10}}^{\mathrm{3}}$) allows the overriding plates to accommodate little shear exerted by the wedge flows. Consequently, they undergo little deformations (${\dot{\mathit{\epsilon }}}_{xy}<{\mathrm{10}}^{-\mathrm{16}}$ s⁻¹) and hardly produce topographic uplifts. This condition rather facilitates the return flows in localizing high uplifts within the wedge rear. Moderate to low μ_r favours the overriding plate to transmit the wedge-induced shear and undergo deformations with considerable uplift rates (Figs. 9a, 11a, b). This condition thus results in a widely uplifted region in the forearc. Such fore-arc uplifts are more pronounced in wedges with large taper angles (which is equivalent to μ_r=1) (Fig. 11c, d). This interpretation explains the occurrence of broad forearc highs in many long-lived, large accretionary complexes, such as the Lesser Antilles, the Alaskan, the Makran, and the Cyprus-Anatolian margin (Pavlis and Bruhn, 1983). Several authors have presented similar interpretations for the tectonics of fore-arc highs, emphasizing the role of upward viscous flows in the lower crust (Fernández-Blanco et al., 2020; Pavlis and Bruhn, 1983). The uplift velocity of the fore-arc region increases when the slab advance accompanies subduction (Figs. 8b, 9b, 11). The slab rollback, in contrast, forces weaker overriding plate (${\mathit{\mu }}_{\mathrm{r}}<\sim {\mathrm{10}}^{\mathrm{2}}$) to flow towards the wedge, resulting in fore-arc subsidence, as also observed in the physical experiments (Fig. 9c-i, ii). These results are in agreement with the previous experimental results where the subduction setup produced fore-arc subsidence and sub-horizontal extension during the trench retreat (Chen et al., 2017; Xue et al., 2022).

5.4 Can accretionary wedges generate significant tectonic overpressures?

Pressure estimates from geobarometry of exhumed metamorphic rocks are generally equated with the depth (h) of metamorphism and reconstruction of the tectonic history, assuming that the lithostatic pressure, p=ρgh, where ρ is the average rock density, and g is the gravitational acceleration. However, this condition applies essentially to a static fluid condition. A direction of studies accounts for the effects of dynamic pressure, calculated as the difference between the total pressure (or mean stress) and lithostatic pressure, i.e., excess pressure (Mancktelow, 1995; Petrini and Podladchikov, 2000). This can be hundreds of MPa or even several GPa higher (tectonic overpressure) or lower (tectonic underpressure) than the corresponding lithostatic pressure. Several numerical model studies claim the possibility of tectonic overpressure (TOP) in subduction-collisional zones (Li et al., 2010; Marques et al., 2018; Raimbourg and Kimura, 2008; Reuber et al., 2016). It has been shown theoretically that in the brittle regime, faulting in compression occurs under pressures of tens to hundreds of MPa, where the total pressure reaches up to twice the lithostatic pressure in dry rocks (Mancktelow, 2008; Petrini and Podladchikov, 2000). Mancktelow (1995), Moulas et al. (2021) analytically determined tectonic overpressure in downward tapering viscous subduction channels bounded by rigid and deformable wall(s), respectively. We use our analytical solution to evaluate the maximum dynamic pressure (P_dmax) in a wedge, which is found to be a function of several parameters: subduction velocity vector (U: both magnitude and direction), overriding plate velocity (U_op), the taper angle (θ₁) and viscosity (μ₁) of the wedge, and overriding plate-wedge viscosity ratio (μ_r). The following conditions: high μ₁, μ_r, U, U_op and low ϕ (negative), θ₁ values give rise to higher values of P_dmax (Fig. 13).

Figure 13 Plot of the total pressure (P), normalized to the lithostatic pressure (P_lith) along the acute bisectors of the wedges for varying μ₁, θ₁, μ_r, and ϕ. P_lith is calculated assuming the overburden density as 2800 kg m⁻³.

The magnitude of dynamic pressure in the downward-tapering wedge increases towards the tip and theoretically becomes infinity at the wedge tip (S-point), a singularity point. Natural accretionary wedges, however, are likely to possess a finite width at its base instead of terminating at a single point, as some volumes of the incoming sediments leak through the bottom of the accretionary wedges and form the deeper subduction channel (Agard et al., 2009). Hence, we find it reasonable to not consider the anomalously high dynamic pressure at the bottom 2 km of our model wedge. This gradient in dynamic pressure drives the exhumation of deep-crustal HP rocks (Fig. 1b). The dynamic pressure scales linearly with the viscosity (μ₁) of the wedge materials. The choice of viscosity values is a critically important issue in the dynamic pressure calculations. Field data, compiled from a range of exhumed subduction complexes and experimental rock deformations suggest μ₁ in the ranges of ∼10¹⁸ to 10²⁰ Pa s at the subduction interfaces (Abila et al., 2025). For our calculations, we thus consider μ₁=10¹⁹ Pa s for metasedimentary rocks in the wedge (Behr and Becker, 2018). Our analytical results show that for a moderate subduction velocity of 3 cm yr⁻¹ and a moderate taper angle of 30°, maximum TOP is less than 8 % of the lithostatic pressure in a wedge with μ₁=10¹⁹ Pa s (Fig. S1). The magnitude of P_dmax reduces further to <4 % with a decrease in the overriding plate viscosity, μ₂≤10²¹ Pa s (Fig. 13a, b, c). The magnitude of dynamic pressure can be of one order higher in wedges with higher viscosity (μ₁≥10²⁰ Pa s) (Fig. 13d). Hence, some high-pressure rocks record extremely fast exhumation rates, as calculated from lithostatic pressure approximations might reflect the influence of TOP at relatively shallow depths ($<\sim \mathrm{40}$ km) (Marques et al., 2018). Presence of high viscosity units can also lead to higher TOP in less competent materials with lower deviatoric stress, due to force-balance requirements (Schmalholz et al., 2014; Schmalholz and Podladchikov, 2013). TOP also increases for narrower wedges with lower taper angles (Fig. 13b). Extremely small taper angles, however, does not produce significant dynamic pressure, unless the overriding plate is relatively very strong (viscosity $>\sim {\mathrm{10}}^{\mathrm{22}}$ Pa s). A weak overriding plate cannot sustain a narrow wedge, and with progressive time, the wedge develops a rearward bulge close to the corner (Fig. 9a-iii). This causes further drop in dynamic pressure at the wedge base. Tectonic overpressure becomes extremely low (<2 %) in wider wedges (Fig. 13c). In addition, our theoretical results indicate that the process of slab advance can elevate the TOP in the subduction zone (Figs. 13, S2). On the contrary, slab rollback creates a condition of tectonic underpressure (TUP) within the wedge, the magnitude of which increases towards depth (Fig. 13).

To summarize, several previous studies have emphasized the importance of tectonic overpressure in exhumed HP units. The present study suggests that for accretionary wedges comprising of low viscosity sediments, the nonlithostatic pressure component is small, generally in the order of tens of MPa (Reuber et al., 2016; Yamato et al., 2007). Significant tectonic overpressure can develop if any of these following conditions are satisfied: (1) wedges consist of materials of high viscosity, (2) the wedge is narrow with a strong hanging wall, and (3) the subduction involves slab advance.

6 Model limitations

The theoretical modelling approach used in this study is simplified in the fact that it assumes an already existing homogeneous viscous wedge above the subducting plate, and also excludes the complex processes (e.g., accretion, sedimentation, erosion) in the accretionary wedges (Angiboust et al., 2022; Clift and Vannucchi, 2004; Simpson, 2010). Moreover, the models consider the whole accretionary wedge and the overriding plate crust as Newtonian viscous materials of constant viscosities and equal densities, and thereby do not account for the effects of non-linear, or depth-dependent rheology of the crust (Ranalli, 1995). In natural wedges, presence of crustal heterogeneities, or discrete shear zones can alter the velocity field predicted from our analytical models with a continuum approximation (Borderie et al., 2018; Grujic, 2006; Jaquet et al., 2018). However, considering that small-scale (with respect to the bulk wedge dimension) perturbations, such as shear zones, are averaged out, we analyze the first-order velocity field using the continuum approximation. The effect of buoyancy forces in orogenic wedges, arising from density contrast between the wedge materials and the surrounding (Burov et al., 2014; Candioti et al., 2021) is also neglected in the models. Another limitation is the assumption of a rigid upper mantle, represented by the rigid model base in experiments, which do not allow any wedge material to flow across the Moho. In natural subduction zones, however, part of the incoming sediments may leak through the accretionary wedge bottom in the form of a deep subduction channel (Clift and Vannucchi, 2004; Guillot et al., 2009). In spite of these limitations, the present theoretical and analogue models provide significant insights into the effects of the various parameters (geometric, kinematic and rheological) on the first order kinematics of a mature, steady-state wedge (Yamato et al., 2007). Our analytical model gives us an instantaneous velocity field under a given wedge geometry. For a deformable overriding plate, or slab-oblique kinematics, the wedge, however, undergoes geometrical modifications with progressive subducting slab movement. Particularly, progressive slab rollback replaces the S-point in the tapered geometry by a separation of finite width on the base. The analytical results are applicable if the separation remains small in the course of slab migration. Analogue models, run for significant finite movements, however, allow us to investigate the temporal evolution of wedge kinematics. Experimental results suggest that the velocity pattern predicted from analytical solutions remain qualitatively unchanged even after considerable finite deformation, until gravity-driven flows begin to capture the wedge kinematics at shallow depths.

7 Conclusions

A generalized analytical solution of the corner flow theory has been derived to study the kinematics of crustal flows in viscous accretionary wedges, with a special emphasis on the criticality of the conditions required for return flows and thereby exhumation of deep-crustal rocks in convergent tectonic boundaries. The analytical solution, supported by scaled physical laboratory model experiments, suggest that the subduction obliquity (ϕ), overriding plate-wedge viscosity ratio (μ_r) and the taper angle (θ₁) of the wedge are principal factors in determining the velocity and strain fields within the accretionary wedge. The principal findings of this study are concluded in the following: 1) Accretionary wedges can produce focused return flows specifically in settings with ${\mathit{\mu }}_{\mathrm{r}}\ge \sim {\mathrm{10}}^{\mathrm{3}}$, ${\mathit{\theta }}_{\mathrm{1}}\ge \sim \mathrm{30}$°, and subduction with little or no slab rollback. The slab rollback (ϕ>0°) process drastically weakens the return fluxes in the wedges. (2) Return flow velocities in the deeper regions of wedges can range from <1 mm yr⁻¹ to a few cm yr⁻¹. Slab advance (ϕ<0°) during the subduction promotes the vertical flow velocities, attaining maximum values comparable to the subduction velocity, provided the overriding plate is strong ($\sim {\mathit{\mu }}_{\mathrm{2}}\ge \sim {\mathrm{10}}^{\mathrm{22}}$ Pa s) and wedge taper angle is not low (${\mathit{\theta }}_{\mathrm{1}}\ge \sim \mathrm{30}$°). (3) The shear sense may reverse across the wedge under a critical combination of the aforesaid parameters. (4) Weak overriding plates (${\mathit{\mu }}_{\mathrm{r}}<\sim {\mathrm{10}}^{\mathrm{2}}$) accommodate significant amounts of global deformations, resulting in vertical uplifts and subsidence for slab-parallel subduction or advance, and rollback, respectively. (5) Accretionary wedges comprising low-viscosity (∼10¹⁹ Pa s) sediments are unlikely to develop large tectonic overpressures (or underpressures).

Appendix A

We show the detailed expressions of the eight equations used for deriving the eight constants A_i, B_i, C_i, D_i (i=1, 2). These eight equations are derived from the eight boundary conditions. The velocity conditions at the wedge base (Eq. 13) follow:

$$\begin{array}{}\text{(A1a)}& {\displaystyle }& {\displaystyle }{A}_{\mathrm{1}}+D_{\mathrm{1}}=-U\cos \left(\mathit{\varphi }\right)\text{(A1b)}& {\displaystyle }& {\displaystyle }-B_{\mathrm{1}}=-U\sin \left(\mathit{\varphi }\right)\end{array}$$

At the Moho in the overriding plate (OP), the velocity boundary conditions (Eq. 14) are expressed as:

$$\begin{array}{}\text{(A2a)}& {\displaystyle }\begin{array}{rl}{A}_{\mathrm{2}}\cos \left({\mathit{\theta }}_{\mathrm{2}}\right)& -B_{\mathrm{2}}\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)+C_{\mathrm{2}}\left(\sin \right({\mathit{\theta }}_{\mathrm{2}})+{\mathit{\theta }}_{\mathrm{2}}\cos ({\mathit{\theta }}_{\mathrm{2}}\left)\right)\\ & +D_{\mathrm{2}}\left(\cos \right({\mathit{\theta }}_{\mathrm{2}})-{\mathit{\theta }}_{\mathrm{2}}\sin ({\mathit{\theta }}_{\mathrm{2}}\left)\right)=-U_{\mathrm{op}}\end{array}\text{(A2b)}& {\displaystyle }\begin{array}{rl}-\left(A_{\mathrm{2}}\sin \right({\mathit{\theta }}_{\mathrm{2}})& +B_{\mathrm{2}}\cos \left({\mathit{\theta }}_{\mathrm{2}}\right)+C_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & +D_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}\cos \left({\mathit{\theta }}_{\mathrm{2}}\right))=\mathrm{0}\end{array}\end{array}$$

Using the continuity condition for the two velocity components across the wedge–OP interface (Eq. 15), we obtain the following relations:

$$\begin{array}{}\text{(A3)}& {\displaystyle }\begin{array}{rl}{A}_{\mathrm{1}}\cos \left({\mathit{\theta }}_{\mathrm{1}}\right)& -B_{\mathrm{1}}\sin \left({\mathit{\theta }}_{\mathrm{1}}\right)+C_{\mathrm{1}}\left(\sin \right({\mathit{\theta }}_{\mathrm{1}})+{\mathit{\theta }}_{\mathrm{1}}\cos ({\mathit{\theta }}_{\mathrm{1}}\left)\right)\\ & +D_{\mathrm{1}}\left(\cos \right({\mathit{\theta }}_{\mathrm{1}})-{\mathit{\theta }}_{\mathrm{1}}\sin ({\mathit{\theta }}_{\mathrm{1}}\left)\right)\\ & =A_{\mathrm{2}}\cos \left({\mathit{\theta }}_{\mathrm{1}}\right)-B_{\mathrm{2}}\sin \left({\mathit{\theta }}_{\mathrm{1}}\right)+C_{\mathrm{2}}\left(\sin \right({\mathit{\theta }}_{\mathrm{1}})\\ & +{\mathit{\theta }}_{\mathrm{1}}\cos \left({\mathit{\theta }}_{\mathrm{1}}\right))+D_{\mathrm{2}}(\cos \left({\mathit{\theta }}_{\mathrm{1}}\right)-{\mathit{\theta }}_{\mathrm{1}}\sin \left({\mathit{\theta }}_{\mathrm{1}}\right))\end{array}\text{(A4)}& {\displaystyle }\begin{array}{rl}-\left(A_{\mathrm{1}}\sin \right({\mathit{\theta }}_{\mathrm{1}})& +B_{\mathrm{1}}\cos \left({\mathit{\theta }}_{\mathrm{1}}\right)+C_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{1}}\sin \left({\mathit{\theta }}_{\mathrm{1}}\right)\\ & +D_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{1}}\cos \left({\mathit{\theta }}_{\mathrm{1}}\right))\\ & =-\left(A_{\mathrm{2}}\sin \right({\mathit{\theta }}_{\mathrm{1}})+B_{\mathrm{2}}\cos ({\mathit{\theta }}_{\mathrm{1}})\\ & +C_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\sin \left({\mathit{\theta }}_{\mathrm{1}}\right)+D_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\cos \left({\mathit{\theta }}_{\mathrm{1}}\right))\end{array}\end{array}$$

Lastly, the continuity of normal and shear stress components across the wedge–OP interface (Eq. 19) leads to:

$$\begin{array}{}\text{(A5a)}& {\displaystyle }\begin{array}{rl}& \left(\mathrm{2}{\mathit{\mu }}_{\mathrm{1}}\right(D_{\mathrm{1}}\cos \left({\mathit{\theta }}_{\mathrm{1}}\right)+C_{\mathrm{1}}\sin \left({\mathit{\theta }}_{\mathrm{1}}\right)\left)\right)/r\\ & =\left(\mathrm{2}{\mathit{\mu }}_{\mathrm{2}}\right(D_{\mathrm{2}}\cos \left({\mathit{\theta }}_{\mathrm{1}}\right)+C_{\mathrm{2}}\sin \left({\mathit{\theta }}_{\mathrm{1}}\right)\left)\right)/r\end{array}\text{(A5b)}& {\displaystyle }\begin{array}{rl}& -\left(\mathrm{2}{\mathit{\mu }}_{\mathrm{1}}\right(C_{\mathrm{1}}\cos \left({\mathit{\theta }}_{\mathrm{1}}\right)-D_{\mathrm{1}}\sin \left({\mathit{\theta }}_{\mathrm{1}}\right)\left)\right)/r\\ & =\left(\mathrm{2}{\mathit{\mu }}_{\mathrm{2}}\right(C_{\mathrm{2}}\cos \left({\mathit{\theta }}_{\mathrm{1}}\right)-D_{\mathrm{2}}\sin \left({\mathit{\theta }}_{\mathrm{1}}\right)\left)\right)/r\end{array}\end{array}$$

By solving the eight equations above, we find the expressions of the eight constants:

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}{A}_{\mathrm{1}}& =-\left(\mathrm{2}U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}\cos \right(\mathit{\varphi })-U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}\cos ({\mathit{\theta }}_{\mathrm{2}})-U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}\cos (\mathit{\varphi })\\ & +\mathrm{2}U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}^{\mathrm{2}}\cos \left(\mathit{\varphi }\right)+\mathrm{2}U{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}\cos \left(\mathit{\varphi }\right)+U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \left(\mathit{\varphi }\right)\\ & +U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \left({\mathit{\theta }}_{\mathrm{2}}\right)+U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}\cos \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\right)\cos \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & -U{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}\sin \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\right)\sin \left(\mathit{\varphi }\right)+U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}\sin \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\right)\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & -\mathrm{2}U{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\sin \left(\mathit{\varphi }\right)+\mathrm{2}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)-\mathrm{4}U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}\cos \left(\mathit{\varphi }\right)\\ & -\mathrm{4}U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\cos \left(\mathit{\varphi }\right)-U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\right)\cos \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & +U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\sin \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\right)\sin \left(\mathit{\varphi }\right)\\ & -U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\sin \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{2}}\right)\sin \left(\mathit{\varphi }\right)-U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\sin \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\right)\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & +U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}\cos \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\right)\cos \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{2}}\right)\cos \left(\mathit{\varphi }\right)+U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}\sin \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\right)\\ & \cdot \sin \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{2}}\right)\cos \left(\mathit{\varphi }\right)+\mathrm{2}U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\sin \left(\mathit{\varphi }\right)\\ & -\mathrm{2}U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}\sin \left(\mathit{\varphi }\right)\\ & -\mathrm{2}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)+\mathrm{2}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & +\mathrm{4}U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\cos \left(\mathit{\varphi }\right)\\ & -U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\right)\cos \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{2}}\right)\cos \left(\mathit{\varphi }\right)\\ & -U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\sin \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\right)\sin \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{2}}\right)\cos \left(\mathit{\varphi }\right))\\ & /(\mathrm{2}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}+\mathrm{2}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}^{\mathrm{2}}+\mathrm{2}{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}+{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}})\\ & +{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-\mathrm{2}{\mathit{\theta }}_{\mathrm{2}})-{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}-{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}\\ & -\mathrm{4}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}-\mathrm{4}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\\ & -{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\right)+{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{2}}\right)\\ & -{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-\mathrm{2}{\mathit{\theta }}_{\mathrm{2}})+\mathrm{4}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}})\end{array}\end{array}$$

$$\begin{array}{}\text{(A7)}& B_{\mathrm{1}}=U\sin \left(\mathit{\varphi }\right)\end{array}$$

$$\begin{array}{}\text{(A8)}& \begin{array}{rl}{C}_{\mathrm{1}}& =\left({\mathit{\mu }}_{\mathrm{2}}\right(U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)-U{\mathit{\mu }}_{\mathrm{2}}\sin \left(\mathit{\varphi }\right)+U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & -U{\mathit{\mu }}_{\mathrm{1}}\sin (\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{1}})+U{\mathit{\mu }}_{\mathrm{1}}\sin (\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{2}})\\ & +U{\mathit{\mu }}_{\mathrm{2}}\sin (\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{1}})\\ & -U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}\sin (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-{\mathit{\theta }}_{\mathrm{2}})+U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}\sin (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-{\mathit{\theta }}_{\mathrm{2}})\\ & -\mathrm{2}U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{1}}\cos \left(\mathit{\varphi }\right)+\mathrm{2}U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\cos \left(\mathit{\varphi }\right)+\mathrm{2}U{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\cos \left(\mathit{\varphi }\right)\\ & +\mathrm{2}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{1}}\cos \left({\mathit{\theta }}_{\mathrm{2}}\right)-\mathrm{2}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\cos \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & -\mathrm{2}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\cos \left({\mathit{\theta }}_{\mathrm{2}}\right)\left)\right)\\ & /(\mathrm{2}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}+\mathrm{2}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}^{\mathrm{2}}+\mathrm{2}{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\\ & +{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}\cos \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\right)+{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-\mathrm{2}{\mathit{\theta }}_{\mathrm{2}})-{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}-{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}\\ & -\mathrm{4}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}-\mathrm{4}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\\ & -{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\right)+{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{2}}\right)-{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\\ & -\mathrm{2}{\mathit{\theta }}_{\mathrm{2}})+\mathrm{4}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}})\end{array}\end{array}$$

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}{D}_{\mathrm{1}}& =\left({\mathit{\mu }}_{\mathrm{2}}\right(U{\mathit{\mu }}_{\mathrm{1}}\cos (\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{1}})-U{\mathit{\mu }}_{\mathrm{1}}\cos (\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{2}})\\ & -U{\mathit{\mu }}_{\mathrm{2}}\cos (\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{1}})-U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-{\mathit{\theta }}_{\mathrm{2}})\\ & +U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-{\mathit{\theta }}_{\mathrm{2}})\\ & +U{\mathit{\mu }}_{\mathrm{2}}\cos \left(\mathit{\varphi }\right)+U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}\cos \left({\mathit{\theta }}_{\mathrm{2}}\right)-U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}\cos \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & +\mathrm{2}U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{1}}\sin \left(\mathit{\varphi }\right)-\mathrm{2}U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\sin \left(\mathit{\varphi }\right)-\mathrm{2}U{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\sin \left(\mathit{\varphi }\right)\\ & -\mathrm{2}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{1}}\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & +\mathrm{2}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)+\mathrm{2}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)\left)\right)\\ & /(\mathrm{2}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}+\mathrm{2}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}^{\mathrm{2}}+\mathrm{2}{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}+{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}})\\ & +{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-\mathrm{2}{\mathit{\theta }}_{\mathrm{2}})-{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}-{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}-\mathrm{4}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}-\mathrm{4}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\\ & -{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\right)+{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{2}}\right)\\ & -{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-\mathrm{2}{\mathit{\theta }}_{\mathrm{2}})+\mathrm{4}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}})\end{array}\end{array}$$

$$\begin{array}{}\text{(A10)}& \begin{array}{rl}{A}_{\mathrm{2}}& =-\left(\right(U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}\cos (\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}+\mathrm{2}{\mathit{\theta }}_{\mathrm{2}}))/\mathrm{2}-U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}\cos ({\mathit{\theta }}_{\mathrm{2}})\\ & -\left(U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}\cos \right(\mathit{\varphi }\left)\right)/\mathrm{2}+\left(U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}\cos \right(\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\left)\right)/\mathrm{2}\\ & -\left(U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}\cos \right(\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{2}}\left)\right)/\mathrm{2}\\ & +\left(U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}\cos \right(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}+{\mathit{\theta }}_{\mathrm{2}}\left)\right)/\mathrm{2}+\left(U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}\cos \right(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-{\mathit{\theta }}_{\mathrm{2}}\left)\right)/\mathrm{2}\\ & +U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\sin (\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{2}})-U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}\sin (\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{1}})\\ & +\mathrm{2}U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}^{\mathrm{2}}\cos \left(\mathit{\varphi }\right)\\ & +\mathrm{2}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}\cos \left({\mathit{\theta }}_{\mathrm{2}}\right)+\mathrm{2}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}\cos \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & +\left(U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \right(\mathit{\varphi }\left)\right)/\mathrm{2}+U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & -\left(U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \right(\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}+\mathrm{2}{\mathit{\theta }}_{\mathrm{2}}\left)\right)/\mathrm{2}\\ & -U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}\sin (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-{\mathit{\theta }}_{\mathrm{2}})-\left(U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \right(\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\left)\right)/\mathrm{2}\\ & +\left(U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \right(\mathit{\varphi }+\mathrm{2}{\mathit{\theta }}_{\mathrm{2}}\left)\right)/\mathrm{2}-\left(U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \right(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}+{\mathit{\theta }}_{\mathrm{2}}\left)\right)/\mathrm{2}\\ & +U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\sin \left(\mathit{\varphi }\right)-U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}\sin \left(\mathit{\varphi }\right)+U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & -\left(U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \right(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-{\mathit{\theta }}_{\mathrm{2}}\left)\right)/\mathrm{2}-U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\sin (\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{2}})\\ & +U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}\sin (\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{1}})-\mathrm{4}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}\cos \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & -\mathrm{2}U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\cos \left(\mathit{\varphi }\right)-\mathrm{2}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\cos \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & +U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}\sin (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-{\mathit{\theta }}_{\mathrm{2}})\\ & -U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\sin \left(\mathit{\varphi }\right)-U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}\sin \left(\mathit{\varphi }\right)\\ & +U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)+\mathrm{2}U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\cos \left(\mathit{\varphi }\right)\\ & +\mathrm{2}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\cos \left({\mathit{\theta }}_{\mathrm{2}}\right))\\ & /(\mathrm{2}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}+\mathrm{2}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}^{\mathrm{2}}+\mathrm{2}{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}+{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}})\\ & +{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-\mathrm{2}{\mathit{\theta }}_{\mathrm{2}})-{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}-{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}-\mathrm{4}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}-\mathrm{4}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\\ & -{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\right)+{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{2}}\right)\\ & -{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-\mathrm{2}{\mathit{\theta }}_{\mathrm{2}})+\mathrm{4}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}})\end{array}\end{array}$$

$$\begin{array}{}\text{(A11)}& \begin{array}{rl}{B}_{\mathrm{2}}& =\left(\right(U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}\sin (\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}+\mathrm{2}{\mathit{\theta }}_{\mathrm{2}}))/\mathrm{2}-U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}\sin ({\mathit{\theta }}_{\mathrm{2}})\\ & -\left(U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}\sin \right(\mathit{\varphi }\left)\right)/\mathrm{2}-\left(U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}\sin \right(\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\left)\right)/\mathrm{2}\\ & +\left(U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}\sin \right(\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{2}}\left)\right)/\mathrm{2}\\ & +\left(U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}\sin \right(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}+{\mathit{\theta }}_{\mathrm{2}}\left)\right)/\mathrm{2}-\left(U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}\sin \right(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-{\mathit{\theta }}_{\mathrm{2}}\left)\right)/\mathrm{2}\\ & +\left(U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\sin \right(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-{\mathit{\theta }}_{\mathrm{2}}\left)\right)/\mathrm{2}\\ & +\mathrm{2}U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}^{\mathrm{2}}\sin \left(\mathit{\varphi }\right)+\mathrm{2}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & +\mathrm{2}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)+U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-{\mathit{\theta }}_{\mathrm{2}})\\ & -\left(U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\sin \right(\mathit{\varphi }\left)\right)/\mathrm{2}+U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & -\left(U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\sin \right(\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}+\mathrm{2}{\mathit{\theta }}_{\mathrm{2}}\left)\right)/\mathrm{2}\\ & +\left(U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\sin \right(\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\left)\right)/\mathrm{2}+\left(U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\sin \right(\mathit{\varphi }+\mathrm{2}{\mathit{\theta }}_{\mathrm{2}}\left)\right)/\mathrm{2}\\ & -\left(U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\sin \right(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}+{\mathit{\theta }}_{\mathrm{2}}\left)\right)/\mathrm{2}-U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\cos \left(\mathit{\varphi }\right)\\ & +U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}\cos \left(\mathit{\varphi }\right)-U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}\cos \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & +U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\cos (\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{2}})-U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}\cos (\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{1}})\\ & -U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\cos (\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{2}})+U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}\cos (\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{1}})\\ & -\mathrm{4}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & -\mathrm{2}U{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\sin \left(\mathit{\varphi }\right)-\mathrm{2}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & -U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-{\mathit{\theta }}_{\mathrm{2}})+U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\cos \left(\mathit{\varphi }\right)\\ & -U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}\cos \left(\mathit{\varphi }\right)+U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}\cos \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & +\mathrm{2}U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\sin \left(\mathit{\varphi }\right)+\mathrm{2}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\sin \left({\mathit{\theta }}_{\mathrm{2}}\right))\\ & /(\mathrm{2}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}+\mathrm{2}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}^{\mathrm{2}}+\mathrm{2}{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}\\ & +{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}+{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}\cos \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\right)+{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-\mathrm{2}{\mathit{\theta }}_{\mathrm{2}})-{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}\\ & -{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}-\mathrm{4}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}-\mathrm{4}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\\ & -{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\right)+{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{2}}\right)\\ & -{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-\mathrm{2}{\mathit{\theta }}_{\mathrm{2}})+\mathrm{4}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}})\end{array}\end{array}$$

$$\begin{array}{}\text{(A12)}& \begin{array}{rl}{C}_{\mathrm{2}}& =\left({\mathit{\mu }}_{\mathrm{1}}\right(U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)-U{\mathit{\mu }}_{\mathrm{2}}\sin \left(\mathit{\varphi }\right)+U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & -U{\mathit{\mu }}_{\mathrm{1}}\sin (\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{1}})+U{\mathit{\mu }}_{\mathrm{1}}\sin (\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{2}})\\ & +U{\mathit{\mu }}_{\mathrm{2}}\sin (\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{1}})-U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}\sin (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-{\mathit{\theta }}_{\mathrm{2}})\\ & +U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}\sin (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-{\mathit{\theta }}_{\mathrm{2}})-\mathrm{2}U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{1}}\cos \left(\mathit{\varphi }\right)\\ & +\mathrm{2}U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\cos \left(\mathit{\varphi }\right)+\mathrm{2}U{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\cos \left(\mathit{\varphi }\right)\\ & +\mathrm{2}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{1}}\cos \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & -\mathrm{2}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\cos \left({\mathit{\theta }}_{\mathrm{2}}\right)-\mathrm{2}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\cos \left({\mathit{\theta }}_{\mathrm{2}}\right)\left)\right)\\ & /(\mathrm{2}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}+\mathrm{2}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}^{\mathrm{2}}+\mathrm{2}{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}+{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}})\\ & +{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-\mathrm{2}{\mathit{\theta }}_{\mathrm{2}})-{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}-{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}-\mathrm{4}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}\\ & -\mathrm{4}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}-{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\right)\\ & +{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{2}}\right)-{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-\mathrm{2}{\mathit{\theta }}_{\mathrm{2}})\\ & +\mathrm{4}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}})\end{array}\end{array}$$

$$\begin{array}{}\text{(A13)}& \begin{array}{rl}{D}_{\mathrm{2}}& =\left({\mathit{\mu }}_{\mathrm{1}}\right(U{\mathit{\mu }}_{\mathrm{1}}\cos (\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{1}})-U{\mathit{\mu }}_{\mathrm{1}}\cos (\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{2}})\\ & -U{\mathit{\mu }}_{\mathrm{2}}\cos (\mathit{\varphi }-\mathrm{2}{\mathit{\theta }}_{\mathrm{1}})-U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-{\mathit{\theta }}_{\mathrm{2}})\\ & +U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-{\mathit{\theta }}_{\mathrm{2}})+U{\mathit{\mu }}_{\mathrm{2}}\cos \left(\mathit{\varphi }\right)\\ & +U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}\cos \left({\mathit{\theta }}_{\mathrm{2}}\right)-U_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}\cos \left({\mathit{\theta }}_{\mathrm{2}}\right)+\mathrm{2}U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{1}}\sin \left(\mathit{\varphi }\right)\\ & -\mathrm{2}U{\mathit{\mu }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\sin \left(\mathit{\varphi }\right)-\mathrm{2}U{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\sin \left(\mathit{\varphi }\right)\\ & -\mathrm{2}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{1}}\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)+\mathrm{2}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)\\ & +\mathrm{2}{U}_{\mathrm{op}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}\sin \left({\mathit{\theta }}_{\mathrm{2}}\right)\left)\right)\\ & /(\mathrm{2}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}+\mathrm{2}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{2}}^{\mathrm{2}}+\mathrm{2}{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}+{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}+{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}})\\ & +{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-\mathrm{2}{\mathit{\theta }}_{\mathrm{2}})-{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}-{\mathit{\mu }}_{\mathrm{2}}^{\mathrm{2}}\\ & -\mathrm{4}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}-\mathrm{4}{\mathit{\mu }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}-{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}\right)\\ & +{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos \left(\mathrm{2}{\mathit{\theta }}_{\mathrm{2}}\right)-{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}\cos (\mathrm{2}{\mathit{\theta }}_{\mathrm{1}}-\mathrm{2}{\mathit{\theta }}_{\mathrm{2}})\\ & +\mathrm{4}{\mathit{\mu }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{2}}{\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}})\end{array}\end{array}$$

Substituting the expressions of these eight constants in Eqs. (12) and (11), the two stream functions (ψ₁ and ψ₂) are calculated. These two stream functions (Eq. 11) are then utilized to find the velocity field in the wedge as well as the overriding plate.