The effect of slab gaps on subduction dynamics and mantle upwelling

Gaps within a subducting plate can alter the surrounding mantle flow field and the overall subduction zone dynamics by allowing hot sub-slab mantle to flow through the gaps and into the mantle wedge. This through-slab flow can produce melting of the slab gap edges as well as significant upwelling that can lead to anomalous alkaline volcanism and/or dynamic uplift in the overriding plate, while the altered mantle flow patterns affect the trench evolution. Numerous geodynamic models have investigated the processes that form slab gaps, but few studies have examined the dynamics of slab gap-altered mantle flow, its effects on trench morphology and kinematics, or the controlling parameters on these processes. Here, laboratory subduction models with a pre-cut gap in a subducting silicone plate are used to explore how slab gap size, and slab gap depth influence the surrounding mantle flow field and trench dynamics. Results suggest that both the vertical extent and the depth of the top (trailing edge) of the slab gap are crucial parameters for modulating overall subduction dynamics. They show that a slab gap, which occurs near the surface and initially comprises 30% of the subducting plate width, can extend enough vertically in the slab to produce significant vertical flow through the gap. Changes to the trench geometry and kinematics are also evident in the models, such that doubleand triple-arc geometries are formed during subduction of a shallow slab gap. All of these results are consistent with observations of slab gaps and their induced surface expressions, or the lack thereof, in Eastern Anatolia, East Java, Italy, and Argentina.


Introduction
A variety of driving and resisting forces act on the subducting lithosphere as it descends into the mantle, leading to a vast array of slab geometries and behaviors (e.g. Billen, 2008;Forsyth and Uyeda, 1975;Jarrard, 1986). Variations in slab age, upper plate structure, and rheology affect all these forces, and hence, the overall subduction dynamics. For instance, variation in the effective viscosity of the slab can modify the bending resistance, whereas changes in density alter the driving slab pull force . The change in slab pull force can result in modification of slab dip, with angles ranging from nearly vertical to nearly horizontal, and slab folding in the mantle transition zone (Billen, 2010;Stegman et al., 2010). In addition to variations in the plate age-controlled thermal buoyancy of the slab, compositional heterogeneities (e.g., oceanic plateaus, aseismic ridges, or micro-continents) within the down-going plate can affect the buoyancy of the subducting plate and contribute to a resistance to further subduction (Cloos, 1993;Martinod et al., 2005). This can lead to deformation of the subducting plate and the formation of a slab gap manifested as a tear, break, or hole within the slab (Fig. 1A) (Arrial and https://doi.org/10.1016/j.tecto.2020.228458 Received 13 November 2019; Received in revised form 14 April 2020; Accepted 19 April 2020 Billen, 2013;Cloos, 1993;Mason et al., 2010).
Numerical and analog subduction models have demonstrated several mechanisms responsible for the formation of slab gaps including the subduction of oceanic plateaus (Arrial and Billen, 2013;Mason et al., 2010), subduction of aseismic ridges (Hu and Liu, 2016), subduction of spreading ridges (Burkett and Billen, 2009;Guillaume et al., 2010), the collision of continental blocks (Fernández-García et al., 2019; van Hunen and Allen, 2011) and by interactions between mantle plumes and the slab (Betts et al., 2012). In most of these cases, the presence of a buoyant feature on the subducting plate (e.g. an oceanic plateau, aseismic ridge, young oceanic lithosphere, or continental crust) contributes to a resistance to subduction. Continued pull from the subducted, negatively buoyant slab can yield tearing and potential detachment of the deeper slab in the upper mantle (Burkett and Billen, 2009;Mason et al., 2010; van Hunen and Allen, 2011;Hu and Liu, 2016). These studies illustrate that the formation of a slab gap depends upon a number of variables such as rheology of the slab, rheology of the surrounding mantle, eclogitization of the buoyant feature, as well as the size and initial density of the buoyant feature. Similarly, in the case of an impinging plume, the formation of a slab gap is entirely dependent on the slab being extremely weak (Betts et al., 2012).
While most of these studies have focused on the cause and evolution of the slab gap rather than the consequences of the gap on the surrounding mantle flow field, few analog and numerical models have examined the effect of a slab gap on subduction dynamics and mantle flow. Regardless of the mechanism for gap formation, these studies reveal complex flow patterns in the mantle (Fig. 1B;Jadamec, 2016;Guillaume et al., 2010;Király et al., 2018;Liu and Stegman, 2012). From these studies it is apparent that scenarios without a plume and where the slab rolls back, mantle flow through the slab gap locally reverses the poloidal flow field and adds an additional toroidal component to the flow around the gap. However, the relationship between slab gap parameters and subduction dynamics remains poorly explored. The presence of a gap has been invoked to explain anomalous observations in volcanism, dynamic uplift, subsidence, and mantle flow in different settings (Ferrari, 2004;Hall and Spakman, 2015;Liu and Stegman, 2012;von Blackenburg and Davies, 1995;Wortel and Spakman, 2000). The common thread in such studies is that during lateral motion of the slab through the mantle, either due to slab steepening or trench rollback, a gap might allow for the upwelling of hot sub-slab asthenosphere through the slab gap and into the mantle wedge (Jadamec, 2016), which can induce melting and subsequent volcanism as well as broad dynamic uplift in the overriding plate and changes in trench geometry (Faccenna et al., 2014b;Király et al., 2018).
Here we seek to explain the range of surface expressions of slab gaps in a generic setting by quantifying the effects of two key features of slab gaps: (i) depth, and (ii) size. We use analog subduction models to investigate the role of slab gaps on mantle flow and trench kinematics, and the implied consequences for dynamic uplift and volcanism. These results are then applied to end-member case studies that span the range of observable surface responses to slab gaps in Eastern Anatolia, East Java, Italy, and Argentina.

Methodology
Physical laboratory experiments are widely used in geodynamics and tectonics because their simplicity captures the most important spatial and temporal features of large-scale processes. In the past, analog models have helped to parameterize different aspects of subduction processes, such as trench motion (e.g., Funiciello et al., 2008), slab and trench curvature (e.g., Schellart, 2010), slab-dip variations (Guillaume et al., 2009), and mantle flow patterns disturbed by a slab gap (Guillaume et al., 2010;Király et al., 2018). Analog models have the advantages of being inherently three-dimensional (3D) and running under a natural gravitational field, which together yield a realistic physical response to the initial experimental setup. Here, we build upon previous laboratory subduction models of Király et al. (2018) by modeling single subduction of a fixed oceanic plate with a pre-cut hole, representing a slab gap, in its center. The models vary the width of the slab gap as well as the width of the subducting plate to explore how variations in slab gap size and subducting plate width control perturbations in the local mantle flow field. The effects of variations in slab gap depth are assessed through the time evolution of the models from the moment the slab gap enters to the subduction zone until it entirely reaches the bottom of the tank.

Model setup
The analog models presented here were run in the Laboratory of Experimental Tectonics, at the University of Roma Tre. These models were based on a preliminary set of experiments during the 2017 CIDER program at the University of California, Berkeley. The experimental setup consisted of a Plexiglas tank filled with glucose syrup, sized 80x80xll cm 3 , overlain by a silicone putty plate (sized 45x30x1.5 cm 3 or 45x20x1.5 cm 3 ). The physical properties of these two materials were tuned to represent the Earth's upper mantle and oceanic lithosphere, respectively (Funiciello et al., 2006;Király et al., 2018;Strak and Schellart, 2014). Glucose syrup, a high density (1430 kg/m 3 ), low Á. Király, et al. Tectonophysics 785 (2020) 228458 viscosity (~200 Pa·s at room temperature), transparent Newtonian fluid, served as the analog of Earth's upper mantle. Oceanic lithosphere was modeled with silicone putty, which consisted of Polydimethylsiloxane (PDMS) silicone with iron powder mixed in to calibrate its density (1505 kg/m 3 ) and viscosity (~40,000 Pa·s). This mixture is visco-elastic, displaying Newtonian behavior under experimental strain rates Weijermars and Schmeling, 1986). These materials simplify the Earth's rheological profile to a system with two Newtonian layersa lithosphere and the underlying upper mantle. Following the standard scaling procedure for stresses created by the negative buoyancy of the subducting lithosphere and the resistance of the mantle (e.g. Cruden, 2008, Funiciello et al., 2003, or see a review in Schellart and Strak, 2016), scaled for length, density, and viscosity in a natural gravity field (g model = g nature ) (Davy and Cobbold, 1991;Hubbert, 1937;Weijermars and Schmeling, 1986), the models have a time scale ratio which translates 1 min of experimental time to~1 Myr in nature. Model parameters and the scaling relationships of a reference experiment are listed in Table 1. The experiments were monitored by two cameras placed above and on the side of the experimental tank, taking top and lateral view photos for capturing the evolution of the experiment. Mantle flow was determined by analyzing the sequence of photographs using the PIVLab MATLAB toolbox (Thielicke and Stamhuis, 2014), which uses Particle Image Velocimetry (PIV) to resolve the velocity field between two successive photos. In order to understand the effect of a slab gap on the vertical mantle flow, two light sources were used to highlight a~1 cm wide trench-normal vertical section in the center of the model domain from the top and the bottom of the box (a white cobra light from the bottom and a red laser from the top; Fig. 2). This light illuminated micro-bubbles in the glucose syrup (each highlighting 2-5 pixels i.e. 0.2-0.5 mm on the side-view photos) that acted as the particles used in the PIV analysis with PIVLab. PIV analysis was also done on the topview photos; however, our choice of lighting did not allow for adequate signal-to-noise ratio to properly resolve the horizontal mantle flow.
The results of four selected models of time-dependent, buoyancy driven subduction with a slab gap are summarized in Table 2. The first three models held the plate width constant and varied the width of the imposed slab hole. The fourth model varied the slab width using the intermediate gap size of 6 cm. The two subducting plate widths tested were 20 cm (which scales to 1200 km, Table 2) and 30 cm (1800 km) ( Table 2), both leaving enough space between the edges of the slab and the sides of the box (> 20 cm) to minimize the effects of the walls on subduction. The length of the subducting plate was set to 45 cm (2700 km), which was long enough that the wall to which the plate was fixed did not affect the flow during the slab gap subduction (~15-20 cm between the wall and the trench at the end of the experiments). The slab gap was cut 12 cm (720 km) from the western (leading) edge of the plate and~30 cm from the eastern side of the Plexiglass box, where the plate was fixed (Fig. 2), minimizing any effects from the edges of the Plexiglass box. Slab gap widths of 3 cm (180 km), 6 cm (360 km), and 9 cm (540 km) were pre-cut into the 30 cm plate, making up 10%, 20%, and 30% of the plate width, respectively (Table 2). A slab gap of 6 cm was pre-cut into the 20 cm plate width, comprising 30% of the subducting plate width. Hence, we tested a range of slab-width and slab gap-width ratios as well as the absolute size of the slab gap. The length of each slab gap was initially cut to 2 cm (120 km) (Fig. 2).

Table 1
Physical and scaling parameters of the laboratory experiments. With these parameters, the geometric, kinematic, and dynamic similarity criteria (Hubbert, 1937) are fulfilled. A Reynolds number which is ≤≤ 1 represents a laminar flow. The experiments are scaled for stresses using the buoyancy number Fb = Δϱgh/η (U/L) (e.g., Faccenna et al., 1999), which compares the gravity induced stresses to the viscous resistance (Δϱ being the density difference between the lithosphere and the asthenosphere, g the gravitational acceleration, h the thickness if the lithosphere, η the viscosity of the mantle U the characteristic velocity, and L a characteristic length in the system).  Lights (a white cobra light and a laser light) were aligned into 1 viewing plane one sourced from below the box and one from the top. The side view camera is placed close to (and parallel to) the South border of the box. Note that plate width varies between 3, 6, and 9 cm depending on the experiment, while plate width varies between 20 and 30 cm (see Table 2). Green line marks the position of the trench. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Á. Király, et al. Tectonophysics 785 (2020)

Assumptions
In this section we describe several key approximations: 1. The viscosity and density of each layer are assumed to be constant, and therefore, serve as the average effective viscosity and average effective density. The assumption of a constant Newtonian viscosity, rather than non-Newtonian viscosity, for the upper mantle would prevent localized reductions in asthenosphere viscosity in regions of high strain rate and the associated increased flow rates (Jadamec and Billen, 2010;Jadamec, 2016). However, the overall pattern in the dynamics is captured. 2. The bottom of the tank represents an impermeable boundary. This approximates the Earth's upper-lower mantle boundary at the 660km discontinuity in the case of a large (> 1 order of magnitude) viscosity ratio between the upper and lower mantles. This is supported by studies which infer a large viscosity (~30-100 times) increase at the upper and lower mantle boundary (Mitrovica and Forte, 1997;Ringwood and Irifune, 1988). 3. In order to image the mantle flow around the slab, the experiment does not include any lateral or overriding plates. Hence, it is assumed that both the overriding and lateral plates as well as the subduction interface have the same properties as the upper mantle. This assumption affects the rate of subduction, but allows the model to capture the overall behavior of the system ( van Dinther et al., 2010;Yamato et al., 2009). 4. Subduction is internally-driven, as it is solely forced by slab pull and does not include any external forcing, such as ridge push or mantle wind. This has the advantage of naturally conserving energy in a closed system (Schellart and Strak, 2016). 5. The subducting plate is fixed to the wall at its trailing edge, which can be thought of as modeling an oceanic plate attached to a large, immobile continent. This maximizes the trench retreat rate (Funiciello et al., 2008) and also gives a maximum estimate for the mantle flow. 6. The slab gap is modeled as a pre-cut hole on the subducting plate that interacts with the trench after a~12 cm long section of the plate is subducted (Fig. 2). The formation of the slab gap is not modeled, so the effect of the slab gap is introduced at the surface. While the hole is initially cut in a rectangular shape, the shape of the hole (being weaker than the plate) is altered by stretching (bỹ 50-60% in length from 2 cm to~3.25 cm) as the plate subducts prior to the gap entering the trench.

Analog modeling results
Initially, the general behavior of each model is similar to that of previous single-slab subduction models (e.g., Funiciello et al., 2006;Strak and Schellart, 2014). Prior to being affected by the slab gap, subduction undergoes three phases: 1) an initial free sinking phase with accelerating subduction velocity, 2) slab interaction with the impermeable base of the tank and a short-term deceleration of subduction velocity, and 3) a steady state phase with slab rollback and a relatively constant trench retreat rate. In the models presented here, the steady state phase is interrupted by the subduction of the pre-cut slab gap.
Before reaching the trench, the hole in the lithosphere expands in size as a result of its lower viscosity (with respect to the surrounding plate) and the pull force from the subducted slab down-dip of the gap. When the trench interacts with the gap, the geometry and kinematics of the trench as well as the mantle flow pattern change. In the following sections, we describe these features (trench geometry, trench kinematics, and mantle flow) as observed in the best imaged experiment (w20W6; wplate width, Winitial gap width in cm). This model is then compared to other models with varying slab gap sizes and a larger plate width (Table 2). Readers are referred to Fig. 2 for the geographic directions (north, east, south, west) that are used when discussing the modeling results.

Trench geometry and kinematics
The evolution of trench morphology is illustrated in 5-minute increments in Fig. 3A as well as in the top view supplementary videos. During slab rollback, when the trench reaches the pre-cut hole in the oceanic plate (representing slab detachment of the central segment of the plate), the trench geometry naturally changes in the center. Because the process of slab gap formation is not modeled here, we can only outline the border between the mantle and the oceanic lithosphere materials, which results in a trench evolution that jumps to the easternmost edge of the pre-cut hole in the central part of the plate, leaving two separate subduction arcs on the northern and southern sides (Fig. 3A, blue lines). After these two segments reach the eastern edge of the gap, the trench forms a "triple-arc" geometry with two actively subducting (northern and southern) segments and a central non-subducting segment (Fig. 3A, green lines). The two, active segments slowly consume the central passive section as the slab gap sinks, dragging down lithosphere at their edges. When the edges of the two active sections meet, there is a brief period with a double-arc trench geometry ( Fig. 3A, yellow lines). This phase is short-lived because the center, where the two subducting trench edges meet, quickly accelerates eastward as it is pulled by the other two segments. The original singlearc trench geometry is then quickly recovered (Fig. 3A, red lines).
The same phases are observed in the trench kinematics. Fig. 3B shows the trench retreat velocity through time in model w20W6 at the center, north, and south segments of the trench (shown by arrows on Fig. 3A). The background of the plot is color-coded according to the changes in the trench morphology. Before the trench reaches the edge of the pre-cut hole, the retreat velocity is even along the trench. As previously observed (e.g., Funiciello et al., 2006;Strak and Schellart, 2014), the trench retreat velocity increases during the initial free sinking phase and stabilizes after interaction with the 660 km discontinuity (bottom boundary of the model domain). The steady state velocity in this model is~2 mm/min. As the trench approaches the edge of the gap, the retreat rate in the center decreases. This happens because the gap represents a weak zone which promotes plate motion towards the trench. When the pre-cut hole enters into the subduction zone, and creates a slab gap, the trench retreat slightly increases on the north and south sides of the gap until these two separated subduction zones reach the eastern edge of the gap at the surface. During the triplearc trench geometry (green shaded area on Fig. 3B), the northern and southern arcs show a decreasing retreat rate, due to the lack of slab pull Table 2 List of laboratory models run. Model names are given by wplate width and Winitial gap width. For all models, the initial length of the slab gap was set to 2 cm (120 km) and the length of the plate was 45 cm (2700 km). Values in cm are the experimental values and those given in km are the nature equivalent values. The recovery time represents the time between the start of slab gap subduction and the peak retreat velocity in the center.  Király, et al. Tectonophysics 785 (2020) 228458 at the center, which means that they have to drag the central, nonsubducting part. The eastern edge of the gap transitions from the initial advancing style to an increasingly fast retreat. The retreat rate in the center peaks when the entire slab gap is consumed at the surface and the subduction of a single continuous arc resumes. At the same time, the retreat rate of the two sides also increases. When the trench reassumes a single arc geometry, the trench retreat velocity re-stabilizes at the steady state velocity of~2 mm/min along the entire trench.

Slab gap and mantle flow
Prior to subduction of the pre-cut hole, the flow field in the highlighted vertical plane, as well as at the top of the model on a horizontal plane, exhibits a pattern that is consistent with 3D models of subduction with a single slab (e.g., Funiciello et al., 2006;Strak and Schellart, 2014). PIV analysis shows clearly defined toroidal cells around the edges of the slab, from the top view photos, and poloidal flow in the vertical plane with flow towards the slab and around the slab tip in the mantle wedge region (Fig. 4A). In the sub-slab region, the poloidal flow has a predominantly horizontal component away from the slab, due to the slab rollback, and small horizontal velocities towards the slab just below the lithospheric plate (Fig. 4A).
PIV analysis shows a slight delay between the opening of the slab gap (gap interaction with the trench; Fig. 4B) and the time when the mantle finally flows through the gap (Fig. 4C). This happens when the leading edge of the gap reaches a depth of~4 cm, which is~240 km when scaled to the Earth. The top view PIV reveals that the slab gap induces a small-scale toroidal flow at the edges of the hole, which promotes the formation of the two separate arcs to the north and to the south (as shown also in Király et al., 2018). In the vertical plane, the disturbance of the flow increases as the slab gap opens and becomes larger. At t = 52 min in the experiment, the leading edge of the slab gap reaches a depth of 7 cm (Fig. 4C), while the trailing edge of the gap remains at the surface. This opens up enough space to create an upward flow that is strong enough to overprint the general pattern of corner flow. In the sub-slab area, a general reduction in velocities is also observed (Fig. 4C). The same upward flow above the slab persists for > 20 min until the trailing edge of the gap enters the trench (Fig. 4D). As subduction resumes at the surface, the trailing edge of the slab gap sinks deeper (Fig. 4E). The upward flow pattern remains in the mantle wedge within a small-scale poloidal cell, but the pattern now occurs deeper in the mantle. Meanwhile, the near-surface mantle flow returns to its undisturbed pattern.
Normalized vertical/horizontal velocity ratio is calculated to quantify the upward flow in front of the slab gap (Fig. 5). To determine the normalized velocity ratio, we extract the PIVlab velocity vectors (e.g. Fig. 4) in a 2.36 cm depth layer (140 km) between 13.55 cm west of the trench and 1.77 cm east of the trench (total length equivalent of 920 km surrounding the trench) for each time step. For each velocity vector, a minimum value of ± 0.25 mm/min is applied to the horizontal component to avoid biasing by dividing with small numbers. The ratio of vertical to horizontal velocity is then calculated at each location. Finally, the velocity ratio is scaled by the magnitude of adjusted velocity vector to avoid corruption of the signal by low-magnitude velocities. This can be described as follows: where v is the vertical component of velocity, u is the horizontal component of velocity, and w is the minimum value of 0.25. Fig. 5A shows this ratio across the slab gap in 10 min increments, while Fig. 5B shows the peak of this ratio through time. The normalized velocity ratio develops a peak in front of the trench (to the west) from 50 to 60 min whereas the profile otherwise remains relatively flat (Fig. 5A). This peak is transient, appearing when the leading edge of the slab gap reaches at least a half-mantle depth (5.5 cm; 330 km) and prior to the trailing edge sinking below~2 cm (120 km) (where the ratio is measured) (Fig. 5B).

The effect of slab gap size
In this section, the results of the models with a plate width of 30 cm (1800 km) and a slab gap varying from 3 to 9 cm (180 to 540 km; Figs. 6 and 7) in width are compared to the previously described model with 20 cm plate width and 6 cm initial slab gap width (w20W6).
The models with 6 cm slab gap width and 20 vs 30 cm plate width show that increasing the size of the plate results in a quicker recovery time of the single arc geometry following slab gap subduction (Fig. 6). Otherwise, trench velocity measurements show similar trends for the w20W6 and the w30W6 experiments, showing a small increase in the Á. Király, et al. Tectonophysics 785 (2020) 228458 eastward trench velocity on the northern and southern parts of the trench when the slab gap enters the subduction zone, followed by a slow-down period (Fig. 6C). In w30W6, the northern and southern subduction zones (which are each 5 cm wider than in the w20W6 model) consume the remnants of the slab gap (i.e. the central nonsubducting part) more quickly, which also results in a lower amplitude peak in the retreat velocity at the center of the trench (Fig. 6C). With respect to model w20W6, the PIV analysis of model w30W6 illustrates a short-lived, small magnitude upwelling (supplementary videos w20W6_lat.mp4 and w30W6_lat.mp4; Fig. 7B). However, due to poor lighting during experiment w30W6, additional quantitative comparisons between these two models could not be performed and is something that should be conducted in future studies. Experiments w30W3 and w30W9 represent two end-members with respect to tested slab gap size (initially 3 cm vs 9 cm wide gaps) with the same slab width (30 cm). The evolution of trench morphology is very Á. Király, et al. Tectonophysics 785 (2020) 228458 distinct in these two models ( Fig. 6A and B): w30W3 (small 3 cm wide gap size) goes through the previously described "double-" and "triplearc" geometries in~20 min of experimental time, while w30W9 (9 cm wide slab gap) displays a long-lived "triple-arc" geometry with two highly curved trenches to the north and south of the central non-subducting section. When the remnants of the slab gap are consumed in the center, the trench retreat rate increases the most in model w30W9 and the least in w30W3 (Fig. 6C).
With regards to the mantle velocity field, w30W3 does not produce any significant mantle upwelling (Fig. 7A), though a small disturbance in the flow is observed when the gap reaches the lower part of the upper mantle (supplementary video w30W3_lat.mp4). In w30W9 (Fig. 7C), the gap results in similar mantle upwelling to the w20W6 model except that upwelling lasts for a longer period of time and occurs at a greater distance from the trench (supplementary video w30W9_lat.mp4).

Discussion of analog model results
Regardless of the assumptions involved in the modeling methods, the four time-dependent models show that a slab gap can have a significant impact on the subduction dynamics, manifested as changes in the 3D flow field in the mantle and the geometry of both the slab and trench. As evidenced by numerical models (e.g., Conder and Wiens, 2007;Jadamec and Billen, 2010;Jadamec, 2016) and simple Poiseuille flow calculations (Appendix A), it is expected that the use of a non-Newtonian rheology will increase flow through the slab gap relative to the Newtonian case and allow for a smaller gap size to produce an appreciable amount of vertical flow. The replacement of glucose syrup with a silicone overriding plate that can resist deformation is also expected to reduce the amount of trench curvature (Meyer and Schellart, 2013), and decrease the velocities of trench retreat and slab rollback (Holt et al., 2015;Meyer and Schellart, 2013;Sharples et al., 2014). The model setup promotes trench retreat, which increases sub-slab pressure as well as mantle flow through the slab gap and around the slab edges. It is important to note that the majority of subduction zones are classified as retreating (Faccenna et al., 2007;Heuret and Lallemand, 2005;Schellart et al., 2008). Overall, these critical subduction parameters serve to control the magnitude of mantle flow and trench motion, while remaining independent of the general relationships described above. Thus, while testing the effects of a non-Newtonian rheology, overriding plate, and trench motion polarity on slab gap subduction is beyond the scope of this study, our approach to focus on relative motion and flow patterns rather than precise size or depth values is supported, providing useful mechanical insights into how a slab gap can affect regional mantle flow and trench kinematics. A further limitation in our model setup is the pre-cut hole that models the slab gap. This enforces the slab gap to open at the surface and instantly, without consuming energy on breaking the lithospheric plate. It can be considered an endmember, representing an idealistic plate configuration where the boundary between the oceanic and continental plate is extremely weak, and hence the continental block detaches at the surface. The model results for the trench geometry and kinematics should therefore be applied with care when considering natural settings.
Analysis of the flow field in the vertical plane of our experiments reveals two key findings on how slab gaps affect mantle flow. First, the ambient mantle flow field surrounding a sinking slab is only perturbed in experiments with a larger pre-cut hole (> 20% of slab width). This is seen clearly in Fig. 7, as well as the supplementary videos, which illustrate that strong, near-surface vertical flow adjacent to the slab is only achieved in models w20W6 and w30W9 (pre-cut hole 30% of slab width). In contrast, model w30W3 (pre-cut hole 10% slab width) primarily exhibits horizontal flow near the surface throughout the duration of the experimental run. Unfortunately, the poor image quality of experiment w30W6 prevents a precise determination of vertical mantle flow through the slab gap. However, our PIV analysis suggests a less significant, short-lived upwelling in this model. These results imply that slab gaps do not inherently perturb mantle flow, but instead they must be above a certain size to perceptibly alter mantle flow (Fig. 8). Further, Figs. 4 and 5 show that mantle upwelling only occurs when the vertical extent of the gap is larger than 5-6 cm (half-mantle depth). When the initial gap size is small (3 cm~180 km wide and 2 cm~120 km long), the slab gap does not become large enough to allow for significant flow through the gap.
Calculations that support the observation that a larger gap size is required to alter mantle flow are presented in Appendix A. We approximate the slab gap as a channel and use the Poiseuille flow solution, where an increase in the pressure difference facilitates flow through the channel (Appendix A). Comparison of the Newtonian and non-Newtonian solutions indicate that non-Newtonian flows achieve a greater flow velocity for a smaller gradient in dynamic pressure and/or smaller hole size (Appendix A). While this simple calculation assumes an average value for shear stress at the base of a lithospheric plate Á. Király, et al. Tectonophysics 785 (2020) 228458 (Melosh, 1977), the result is consistent with 2D and 3D numerical models of subduction induced mantle flow that compare a Newtonian versus non-Newtonian viscosity (Jadamec and Billen, 2010;Jadamec, 2016;Király et al., 2017), corroborating a characteristic of non-Newtonian viscosity such that the flow is focused, which would promote stronger through-gap mantle flow.
The second key finding is that the alterations to the mantle flow pattern induced by the slab gap persist for the duration of gap sinking after the gap initially exceeds the necessary size. However, the nearsurface flow pattern is only briefly perturbed (Figs. 4 and 8) while the gap is connected to the surface. In agreement with Király et al. (2018), this temporal pattern also applies to the near-surface horizontal flow. The horizontal surface flow is only perturbed when the slab gap is close to the surface. Once the central section resumes subduction, the mantle flow towards the newly formed and rapidly retreating trench overprints the previous small-scale toroidal pattern. A mantle upwelling in the vertical plane is also observed in experiments where the slab gap extends from the surface to a greater depth of at least 5.5-6 cm, or 330-360 km in Earth's mantle. This indicates that gap-induced vertical mantle flow is closely tied to the depth of the gap itself. Nearsurface flow perturbations are transient and only occur when a gap forms at or near the surface (Fig. 8). The overall pattern of paired sets of toroidal mantle flow through a submerged slab gap is consistent with previous 3D analog (Guillaume et al., 2010) and numerical (Jadamec, 2016;Magni, 2019) models of slab gap subduction. Magni (2019) also notes that the toroidal flow through a slab tear does not necessarily have a vertical component. This is due to the larger depth of gap formation (200 km) and the orientation of the tear that is perpendicular to the retreating slab (Magni, 2019). These results compare well with our experimental results which show that the gap position must be close to the surface and the gap size has to be large enough for upwellings to occur, as well as with our Poiseuille flow calculations which show the relationship between flow and pressure gradient.
Previous geodynamic models of slab gap subduction (Guillaume et al., 2010), emergent slab holes (Liu and Stegman, 2012), and mantle flow around a single slab's edges Jadamec and Billen, 2010;Király et al., 2017;Strak and Schellart, 2014) predict localized mantle upwelling adjacent to the slab edge, which would facilitate decompression melting and the formation of anomalous volcanism on the surface. The models here extend the implications from these previous studies by systematically demonstrating the change in mantle flow as a function of gap size and testing the effect of slab width on the time-dependent evolution of slab-gap related mantle upwelling. While the aforementioned studies as well as our experimental results provide insight on the patterns of mantle flow within the vicinity of a slab gap or around a slab edge, future work quantifying the relative pressure gradients associated with the release of sub-slab pressure following the opening of a slab gap versus the sub-slab pressure release around a lateral slab edge could provide a stronger link between decompression melting and associated anomalous volcanism.
Slab gaps that form close to the surface also affect trench geometry and kinematics. As Fig. 6 shows, a small slab gap of only 10% of the slab width has a short-lived effect on trench geometry, and hence less effect on the trench kinematics. A larger slab gap (30% of the slab width, models w20W6 or w30W9) leads to a longer period of a "triple-arc" geometry with two separated trenches. This is caused by two factors: i) a greater decrease of the slab pull force, and ii) a wider central (nonsubducting) segment in between the two arcs (in the central part), which has to be consumed before a single-arc geometry is restored. Consequently, when subduction finally starts in the central, non-subducting segment of the slab, the trench retreats quickly over a large distance to form a "single-arc" system. This quick retreat rate could induce anomalous extension in the overriding plate. The recovery time is not directly related to the ratio of slab gap width to plate width. For example, the model with a 9 cm (540 km) wide gap on a 30 cm (1800 km) wide plate takes much longer to recover than the model with 6 cm (360 km) wide gap on a 20 cm (1200 km) wide plate (Table 2).

Natural occurrence of slab gaps
The experiments presented here demonstrate that two key parameters (depth of the trailing edge of the gap and its vertical extent) regulate whether slab gaps may produce observable consequences (Fig. 8). We review four case studies (Fig. 9) that highlight the geological observations attributed to slab gaps and discuss how these observations are consistent with our results.

Eastern Anatolia
Anatolian tectonics are dominated by the convergence of the African and Arabian plates with Eurasia. While subduction is ongoing beneath Central Anatolia along the Cyprus trench, Eastern Anatolia has been undergoing continent-continent collision since mid-Miocene slab breakoff (Keskin, 2003;Şengör et al., 2003). In Eastern Anatolia, evidence of the termination of subduction through closure of the southern branch of the Neo-Tethys ocean is preserved in the Bitlis-Zagros Suture Zone (BZSZ) (Şengör and Yilmaz, 1981). The collision of Arabia with Eurasia led to the formation of the high elevation East Anatolian Plateau, which is characterized by an average elevation of~2 km. Today, continued convergence is largely accommodated by northeast-and southeast-striking conjugate strike-slip fault systems (Şengör and Yilmaz, 1981). Despite continued shortening, a crustal thickness of 35-45 km indicates that the high elevation of the plateau is isostatically undercompensated . Furthermore, relatively slow seismic velocities in the uppermost mantle indicate that there is either thin or missing mantle lithosphere beneath the plateau (Delph et al., 2015), suggesting that high elevations may instead be supported by hot, shallow asthenosphere (Şengör et al., 2003). Slow seismic velocities extend throughout the upper mantle and provide little evidence for subducted Neo-Tethyan lithosphere above the mantle transition zone east of the imaged Cyprus slab (e.g. Biryol et al., 2011;Portner et al., 2018).
Voluminous Miocene-to-recent volcanism is present in the East Anatolian Plateau (Keskin, 2003;Pearce et al., 1990) and post-dates the onset of plateau uplift (~12 Ma), with the first small eruptions occur-ring~11 Ma (Pearce et al., 1990). The apparent north-to-south progression of volcanism across this region from Miocene to present, together with the southernmost and the youngest volcanic provinces that exhibit the weakest arc-related geochemical signature, suggests that plateau uplift is related to slab steepening and detachment processes (Keskin, 2003;Reid et al., 2019;Schleiffarth et al., 2018). The relative contribution of slab steepening and detachment to plateau uplift and volcanism is a subject of debate. Geodynamic models suggest that mantle lithosphere delamination may have a more direct control on plateau development (Göğüş and Pysklywec, 2008;Bartol and Govers, 2014;Memiş et al., 2020). Because our models do not include delamination, when comparing our results to Eastern-Anatolia we disregard the effect of this process, while we acknowledge that delamination may overprint or even control the observed uplift and volcanism of the plateau.
The combination of seismic, geological, and geochemical evidence show that the East Anatolian Plateau is a high elevation, isostatically undercompensated plateau that is covered with young volcanics with variable geochemical signatures. Detachment and continued sinking of the subducted Neo-Tethyan lithosphere adjacent to active subduction beneath Central Anatolia has opened a growing slab gap in the upper mantle (Portner et al., 2018;Reid et al., 2019). A relative absence of fast seismic velocities in the upper mantle east of the Cyprus slab ( Fig. 9A) indicates that the plateau is underlain by a slab gap of > 400 km down-dip (Portner et al., 2018). Slow seismic velocities within the gap are often interpreted as an indication of hot upwelling asthenosphere through the slab gap, which might have provided a dynamic support for high elevation of the Eastern Anatolian Plateau and might have caused recent volcanism (Biryol et al., 2011;Keskin, 2003;Portner et al., 2018;Şengör et al., 2003).
In this example, slab detachment created a significant slab gap, which is both large (> 400 km down-dip) and shallow (< 100 km depth), allowing for significant perturbations to the ambient mantle flow field, including near-surface vertical flow. This flow may have contributed directly to the observed Miocene-recent uplift and volcanism on the overriding Eurasian plate (Şengör et al., 2003). Further, though they do not rule out the potential for delamination to contribute to or cause these observations, our results are consistent with gap formation as a contributor. On the other hand, the complex history over the last~10 Myrs (after the slab-gap opening) in the area does not allow us to compare other surface features, such as trench geometries, with the models.

East Java
The Sunda Arc is the result of continuous subduction of the Indo-Australian plate beneath the Sundaland block of the Eurasian plate since at least~45 Ma (e.g. Hall, 2002). Volcanism along the arc has not Á. Király, et al. Tectonophysics 785 (2020) 228458 (Bird, 2003). The red triangles are volcanoes active in the Holocene (Siebert and Simkin, 2011). In each cross-section, black arrows summarize the interpreted mantle flow surrounding each slab gap. Gray circles are earthquakes from the Comprehensive Catalog (ComCat) of the USGS National Earthquake Information Center. The top 100 km of each cross-section is shaded due to relatively poor resolution of the shallow mantle in teleseismic tomography. Gray lines are %dVp contours of +0.5 in Eastern Anatolia, East Java, and the Apennines, and +0.8 in Argentina.
The tomography models plotted are (A) ANA2_P_2018 (Portner et al., 2018), (B) SAM4_P_2017 , (C and D) UU-P07 (Amaru, 2007). Abbreviated labels are as follows: BZSZ, Bitlis-Zagros Suture Zone; EAP, East Anatolian Plateau. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Á. Király, et al. Tectonophysics 785 (2020) 228458 been continuous throughout this time, with at least two notable periods of quiescence (Smyth et al., 2008). In the East Java arc segment, the recent volcanic history is punctuated by transient Pleistocene back-arc volcanism (Edwards et al., 1994). The back-arc volcanism is characterized by a range of geochemical compositions including small volumes of potassic alkaline lavas that indicate heterogeneous mantle source with an enriched mantle component (Edwards et al., 1994). In light of patterns of Wadati-Benioff zone seismicity and seismic tomography, this anomalous volcanism has been attributed to a descending slab gap (Hall and Spakman, 2015). An absence of seismicity in the Wadati-Benioff zone at 250-500 km depth is coincident with a large gap in the fast seismic velocity anomaly attributed to the subducting Indo-Australian lithosphere ( Fig. 9B; Amaru, 2007;Widiyantoro et al., 2011). In the model of Hall and Spakman (2015) this gap formed when buoyant material in the subducting slab arrived at the trench and was impeded from passing further into the subduction zone. The arrival of buoyant material impeded subduction and created a tear in front of the buoyant block. As subduction continued around the opening gap, convergence has been continuous which caused thrusting in Java. The model suggests that subduction resumed behind the plateau causing a hole in the sinking slab. The resultant gap then passed underneath the region of the backarc where the potassic alkaline volcanic deposits were erupted. The presence of the gap contributed to an enriched mantle source required for potassic alkaline volcanism. As subduction continued and the gap sank along with the slab, the high-K alkalic back-arc volcanism ceased (Edwards et al., 1994;Hall and Spakman, 2015). Given that this gap has a moderately large size (~300 km down-dip) this history is consistent with a sinking slab gap. The thrusting phase in Java might correspond to the advance of the leading edge of the gap in the models. However, as the slab gap is now at greater depth, the trench has a normal, "singlearc" geometry. The gap may have perturbed mantle flow through the duration of its history, but it was only able to affect near surface flow and induce volcanism during its shallow, initial stage. A detailed study of the local mantle flow field is necessary to confirm the predicted perturbations to modern flow.

Apennines
The Italian peninsula is a narrow (~200 km wide) mountain range including the Apennines and the Calabrian arc (Fig. 9C). The Apennines-Calabrian chain is characterized by significant along-strike variations in topography, uplift, crustal thickness, and seismicity (Faccenna et al., 2014b). The central Apennines has the highest relief (2.9 km) as well as the highest mean elevation (~0.6 km), which is 0.3-0.5 km higher than the mean elevation of the northern Apennines and Calabria. This region also has the thinnest crust (35-40 km) and the shallowest seismicity within the whole area. While the Wadati-Benioff zone extends continuously to a depth of~400 km below the Calabrian arc, the deepest earthquakes in the northern Apennines are at 100-150 km depth. There are no deep or intermediate depth earthquakes in the central Apennines and seismicity is restricted to shallow crustal earthquakes with extensional focal mechanisms. P-wave seismic tomography (Amaru, 2007;Piromallo and Morelli, 2003) reveals a slow seismic velocity anomaly beneath the central-south Apennines that spans from the surface to a depth of~400 km, which is interpreted as a slab window (Faccenna et al., 2014b;Magni et al., 2014;Piromallo and Morelli, 2003). Dynamic topography calculations using P-wave seismic tomography models suggest that up to 0.8 km of topography above the slab window is dynamically supported, that is consistent with isostatically non-supported (residual) topography calculations (Faccenna et al., 2014b). While the slab detached in the Central Apennines, subduction continued to the south along the Calabria arc and probably transitioned to subduction-style delamination in the Northern Apennines (Chiarabba et al., 2014). The central area has a complex pattern of seismic anisotropy indicating that mantle flow is oriented roughly perpendicular to the subduction zone. On the contrary, measurements of seismic anisotropy from the Calabrian arc suggest orientation of the mantle flow parallel to the subduction zone (refer to the compilation in Faccenna et al., 2014a). The seismic anisotropy pattern in the area is in agreement with the horizontal mantle flow around a slab window and in between two oppositely dipping slabs, as shown in Király et al. (2018).
These observations support our prediction of the perturbed mantle flow by an extensive slab gap and if the gap is shallow, it is a possible contribution to dynamic uplift of the overriding plate. Additionally, the observed Apennines-Calabria "triple-arc" geometry with more arcuate mountain chains in the Northern-Apennines and in Calabria connected with the high topography area in the Central-Apennines ( Fig. 9C; Király et al., 2018), agrees with the "triple-arc" geometry produced in the analog models with a shallow slab gap. However, the geodynamic setting in the Central-Mediterranean is further complicated by the opposite subduction of the eastern side of the Adria plate. The presence of two subduction zones might enhance the mantle flow through the slab gap as the sub-slab pressure is increased under the Adria plate (Király et al., 2018).

Argentina
Teleseismic tomography reveals a~200 km wide gap in the Nazca slab beneath central Argentina at~300 km depth ( Fig. 9D; Portner et al., 2017). Shear-wave splitting analysis suggests a notable degree of mantle flow through the slab gap, in contrast to the adjacent mantle flow field (Lynner et al., 2017). The slab gap is located down-dip of the Pampean flat slab, where the subducted segment of the Juan Fernández Ridge may increase the buoyancy of the slab and contribute to shallow subduction (van Hunen et al., 2002). Geodynamic modeling of aseismic ridge subduction agrees with the observations (Hu and Liu, 2016), indicating that the slab gap formation may be a direct consequence of differential sinking rates between the buoyant ridge and the surrounding slab.
There is no clear surface expression of upwelling through the slab gap. The closest surface features to the gap are the eastern Sierras Pampeanas in central Argentina, a series of thick-skinned basementcored uplifts of Paleozoic age (Goddard et al., 2018). The latest, Pliocene uplift of the mountains was predated by small volumes of volcanism that expanded west-to-east into the back arc with trend in composition that reflects an increase in melting depths with thickening of the crust until magmatism terminated~5 Ma (Ramos et al., 2002). However, the eastern Sierras Pampeanas are only the easternmost extent of the Sierras Pampeanas, which is a larger structural domain that extends westward to the Andean Cordillera. Deformation and volcanism within the Sierras Pampeanas have long been attributed to migration of the Juan Fernández Ridge and flattening of the Pampean slab segment (Ramos et al., 2002 and references therein), suggesting the easternmost portion is unlikely to be related to the presence of the imaged slab gap.
The observed~300 km deep slab gap beneath Argentina is likely to have been formed deeper than the gaps in the other described case study regions. While this~200 km wide slab gap is enough for initiation through-slab flow (Lynner et al., 2017), the depth of the gap prevents near-surface vertical flow and, consequently, surface expressions such as volcanism and/or uplift within the overriding plate or the deformation of the subduction trench. This further supports our conclusion of the analog modeling that despite the importance of the slab gap size, the most significant upwelling through a slab gap occurs adjacent to the gap and generates a surface expression only when the gap is near the surface.

Conclusions
Laboratory subduction models show how the presence of a gap in a slab influences the surrounding mantle flow field and overall Á. Király, et al. Tectonophysics 785 (2020) 228458 subduction zone dynamics. Our experiments show that slab gaps can induce a small-scale toroidal flow pattern around the edges of the slab gap, upwelling of mantle through the slab gap, and a double-or triplearc geometry with varying trench kinematics. In our results these effects are the most pronounced if the slab gap size is sufficiently large and if the gap is shallow (near the surface). The model results support observations of slab gap-induced surface expressions found in Eastern Anatolia, Pleistocene East Java, and Italy, and explain the lack of slab gap-induced surface observables when gaps are deep, such as the ones found in modern East Java and Argentina. Thus, the analog modeling results are consistent with observations of natural subduction zones that contain slab gaps. Results can also be used to evaluate anomalous volcanism, dynamic uplift/subsidence, and/or mantle flow patterns that have previously been interpreted as slab gap-related. Supplementary data to this article can be found online at https:// doi.org/10.1016/j.tecto.2020.228458.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.