Rotational Extension Promotes Coeval Upper Crustal Brittle Faulting and Deep‐Seated Rift‐Axis Parallel Flow: Dynamic Coupling Processes Inferred From Analog Model Experiments

The lower parts of warm, thick continental crust can flow in a ductile fashion to accommodate thinning of the upper brittle crust during extension. Naturally occurring continental rifts with a rift‐axis parallel deformation gradient imply an underlying rotational component. In such settings, rift‐parallel crustal flow transports material perpendicular to the direction of rifting. We use analog experiments to investigate rotational rifting and coeval crustal flow. To test the effect of rift‐axis parallel flow on rift evolution, we use different gravitational loads resulting in a range of horizontal pressure gradient magnitudes which drive horizontal lower‐crustal flow. The use of (three dimensional) 3D Digital Volume Correlation techniques on X‐ray CT data combined with 3D Digital Image Correlation techniques applied to topographic stereo images provides detailed insights on the contemporaneous evolution of ductile flow patterns and brittle rift structures, respectively. Our results depict a complex flow field in the ductile lower crust during rotational rifting with: (a) extension‐parallel horizontal inward flow and vertical upward flow that compensates thinning of the brittle upper crustal layer; (b) rift‐axis parallel lateral flow, that compensates greater amounts of thinning further away from the rotation axis; and (c) different degrees of mechanical coupling between the brittle and viscous layers that change during rift propagation. Our analog experiments provide insights into ductile lower crustal flow patterns during rift evolution. The results emphasize the three dimensionality of rifting, which is an important effect that should be considered when estimating the amount of crustal extension from two dimensional (2D) cross sections.

ductile deformation coupled? To answer these questions, this study presents brittle-ductile 3D analog models of crustal-scale, rotational continental rifting with a particular focus on the role of rift-axis parallel ductile flow. While the previous work of T. C.  analyzed rift propagation and near-surface deformation in detail, this study provides a quantitative analysis of internal deformation. By combining 3D stereo Digital Image Correlation (3D stereo DIC) and Digital Volume Correlation (DVC), we present a novel approach to quantify surface and internal model deformation and gain a comprehensive understanding of crustal flow and its role during rotational continental rifting with an along-strike divergence velocity gradient.

Methods
We use an experimental apparatus which simulates rotational extension around a pivoting point ( Figure 2). An identical setup has been used and described in previous studies of T. C. , , and . Two short, curved and two mobile long sidewalls confine the model domain ( Figure 2a). The mobile sidewalls rotate around a vertical axis which separates the model box into 25 and 65 cm long contractional and extensional model domains, respectively ( Figure 2b). The applied maximum divergence velocity is defined along a circular segment at distance 65 cm with respect to the rotation axis. Simultaneously with the opening motion in the extensional domain, shortening occurs in the contractional domain. The initial model surface is defined by a quasi-rectangular area of 90 cm length and 31 cm width. For our analyses, we leave out the curved model parts close to the short sidewalls, resulting in a rectangular area of interest of about 80 by 31 cm. The model overlies a 5 cm thick foam base that is initially compressed and homogenously expands during the experiment runs. This setup allows for a gradually decreasing basal velocity profile with maximum divergence velocity at the mobile sidewalls and zero velocity near the rotational axis ( Figure 2c). In the contractional domain, divergence velocities switch polarity and the foam undergoes further contraction. Rigid plexiglass bars sandwiched between the basal foam provide additional support to prevent the overlying model from sinking in. This basal setup has been successfully tested and applied for rotational rift models in T. C. .
The mechanical two-layer model consists of a crustal scale, brittle-viscous setup in which a brittle layer simulates the upper crust and the underlying viscous layer the lower crust ( Figure 2c). In the extensional domain, rifting in the brittle layer is localized by applying a viscous seed on top of the viscous layer (Figure 2a). This seed (a semi-cylindrical polydimethylsiloxane [PDMS]/corundum sand mixture rod with radius of ca. 0.5 cm) acts as a pre-existing weak zone along which deformation initiates. We deliberately place the seed along the entire model length to control surface deformation. Note that the absence of a seed results in deformation localization solely along the mobile sidewalls. The applied divergence velocity gradient causes rift propagation toward the rotation axis over time with most mature rift stages developing furthest away from the rotation axis. We place an additional sand layer on top of the contractional domain in some of our models before starting the experiments. This additional gravitational load (see Table 1 for layer (a) Cut-out view of the experimental apparatus confined by two curved and fixed, short, and two long movable sidewalls, respectively. The vertical rotation axis separates an extensional from a contractional domain which contains-in certain experiments-an additional sand layer to induce a pressure-gradient driven flow in the viscous layer. (b) Top view of the experimental apparatus and the region of interest (green rectangle). The pivoting motion around the rotation axis causes a divergent velocity gradient with decreasing velocities toward the rotation axis in the extensional domain. In the shortening domain the pivoting motion induces a convergent velocity gradient with decreasing velocities toward the rotation axis. (c) Model set up cross section and schematic strength profile. The two-layer crustal-scale model consists of a brittle sand layer on top of a viscous mixture with a total thickness of 6 cm. To ensure symmetric and homogenously distributed extension, the model sits on top of a foam base, supported by interlayered plexiglass bars, which expands homogeneously after compression prior to the model run. thickness) induces a horizontal pressure gradient and thereby, horizontal flow in the viscous basal layer in the opposite direction of rift propagation (Figures 2 and 3). Prior to the model runs, we sieve a fine 4 by 4 cm grid of corundum sand on the model surface for visual aid and sprinkle coffee powder to create a "salt and pepper" pattern for DIC.

Analog Model Setup
We conducted a series of crustal scale brittle-viscous analog experiments with variable additional gravitational load-induced pressure gradients causing lateral flow in the viscous layer. For models with an identical brittle-ductile thickness ratio T BD (or identical brittle-ductile strength ratio, S BD ), we compare the deformation evolution in models with an initially flat model surface (i.e., no additional gravitational load) with corresponding models with an additional gravitational load. This thickness ratio T BD is defined as the brittle layer thickness divided by the viscous layer thickness. The corresponding brittle-ductile strength ratio S BD is listed in Table 1. The overall thickness of the two-layer model was 6 cm for all experiments, representing a 30 km thick, continental crust (McKenzie et al., 2000). However, the thickness ratio T BD , between the brittle and ductile layers varies and defines three different model types with T BD = ½, 1, and 2 simulating different crustal configurations. All experiments were run with a maximum divergence velocity of 10 mm/hr (symmetrically 5 mm/hr each side), which linearly decreases toward the rotation axis. This results in a maximum extension of 40 mm or about 13% maximum extension after 4 hr of each model run. Similarly, maximum shortening of 15.4 mm or 5% is reached at the final stage at the distant circular segment in the contractional domain. Lastly, we induced a pressure-gradient driven horizontal flow in the viscous layer by adding an additional quartz sand layer on top of the contractional domain in some models prior to the experiment runs. Thus, the viscous layer flows in the opposite direction to that of rift propagation. For our model naming convention, Mx-y*, the "x" and "y" denotes the brittle-ductile thickness ratio T BD and the height of the additional quartz sand layer on top of the contractional domain, respec tively. Models which have been analyzed by XRCT are further denoted with an asterisk. As a reference setup, we use the models M1-x and M1-x* with a brittle-ductile thickness ratio of 1. The complete model series consists of 10 experiments listed in Table 1.
The viscous basal layer consists of a mixture of corundum sand and PDMS with a mixing ratio 1:1 to achieve a density of 1,600 kg/m 3 . The mixture has a quasi-linear viscosity of 1 × 10 5 Pa and a stress exponent of 1.05 (Zwaan, Schreurs, Ritter, et al., 2018). For the upper brittle layer of the analog models, we use dry quartz sand with a bulk density of 1,560 kg/m 3 to simulate the upper brittle crust. The desired density is achieved by sieving sand into the model box from a height of 30 cm. The brittle-viscous models thus have a density gradient increasing with depth avoiding density instabilities and spontaneous upwelling of the viscous lower layer. Dilation and localized deformation in the quartz sand locally lowers the peak friction coefficient over time, resulting in strain softening of about 16% (i.e., the difference between the coefficient of peak friction and dynamic friction coefficient normalized by the coefficient of peak friction (Panien et al., 2006) from mutual two-point regression analysispresented in Santimano et al., 2015).
For XRCT scanned experiments, we mix small quantities (weight ratio 1:50) of Zirshot ceramic microbeads with the quartz and corundum sands which enhances volumetric patterns in CT scans to facilitate DVC analysis (Adam et al., 2013). Additionally, we sieve a thin layer of corundum sand from 30 cm height on top of the viscous mixture before adding quartz sand to enhance the X-ray absorption contrast between the viscous and brittle domains. The influence of the ceramic microbeads and the thin corundum sand layer on the overall mechanical properties is negligible (Klinkmüller, 2011;Panien et al., 2006;T. Schmid et al., 2020aT. Schmid et al., , 2020b 2020a. All material properties are listed in Table 2. Note. Densities denoted with * are achieved by sieving from 30 cm height.

Scaling
We apply standard scaling equations from Hubbert (1937) and Ramberg (1981) to ensure proper scaling of our models with respect to nature. For brittle Mohr-Coulomb type materials and viscous materials, dynamic similarity is given by the stress ratio * = * * ℎ * , and strain-rate ratio * = * ∕ * , respectively.
We ensure dynamic and kinematic similarity between nature and the experiments by calculating the Smoluchowski number S m , and the Ramberg number R m for the brittle and viscous layers, respectively. S m is defined as the ratio between gravitational stress and cohesive strength (Ramberg, 1981) S m = ρgh/(C + μρgh), where ρ, h, C, and μ are the density, thickness, cohesion and friction coefficient, respectively. Since upper crustal rocks have a cohesion values of 50 MPa and internal friction coefficients of ∼0.6 (Byerlee, 1978), this yields S m ∼1 for both our experiments and nature. For our reference model setup with a brittle-ductile thickness ratio T BD = 1, the Ramberg number R m = ρgh 2 /ηv yields values of 51 and 53 for our experiments and nature, respectively. For models with T BD = 2, R m values are 23 and 24 for the experiments and nature, respectively. For models with T BD = 0.5, R m ∼90 and ∼95 for the experiments and nature, respectively. The Reynolds number R e = ρvh/η is defined as the ratio between inertial forces and viscous forces and is <<1 for all of our experiments and nature.
Based on the above scaling laws, the material properties and similar non-dimensional numbers for the models and nature, we consider our analog experiments to be properly scaled. Model parameters and non-dimensional numbers are given in Table 3.

Contractional Boundary Condition and Induced Pressure Gradient
During opening of the extensional domain of the apparatus, shortening in the contractional domain induces a lateral pressure gradient, resulting in rift-axis parallel flow of the lower viscous layer from the contractional to the extensional model domain (Figure 3a). With subsequent deformation localization in the contractional domain, differential topography increases the lateral pressure gradient, enhancing rift-axis parallel flow ( Figure 3b). Additionally, we simulate different additional gravitational loadings by using sand layers of 1 and 3 cm on top of the Scaling ratios 1 × 10 −6 0.55 1 2 × 10 −6 1 × 10 −6 6 × 10 9 2 × 10 −16 1 × 10 4 2 × 10 −10 a Note that for calculating R m the maximum velocity is used. contractional model domain. To analyze the net effect of such additional gravitational load on the rift-axis parallel flow in the lower viscous layer, the flow component induced by contracting the sidewalls must be subtracted from the total flow. However, in our models the contribution of the contracting sidewalls is negligible.
Based on the above scaling relationships, a 3 cm thick additional sand layer translates to 15 km in nature, which exceeds realistic elevation values. However, we use such an exaggerated gravitational load to enhance and visualize the displacement dynamics, which are present in all our models, even without an additional gravitational load ( Figure 3). The additional sand load causes a pressure gradient from the contractional to the extensional model domain, which drives rift-parallel flow in the viscous layer. If the relative displacement between the bounding layers (i.e., foam base and overlying sand layer) is zero, such a pressure-gradient driven flow is known as Poiseuille flow and its analytical 1D solution is given by: (e.g., Turcotte & Schubert, 1982). Here, u(z) denotes the lateral flow velocity at the vertical position z in a channel with height h. Further, μ is the viscosity of the fluid in the channel and dp/dx the lateral pressure gradient along the horizontal length of the channel. In some of our experiments, we induce a lateral pressure gradient by an initial topographic elevation with pressure p topo given by (Clark & Royden, 2000), where ρ s is the density of the additional brittle sand layer (1,560 kg/m 3 ), and g, and T(x) are the gravitational acceleration and topographic elevation along the model length x, respectively. To measure the net horizontal flow due to the pressure gradient (Figure 3c), we XRCT scan a model with initially flat topography (M1-0*) and quantify the flow velocity due to contracting sidewalls ( Figure 3a). Note that the flow velocity is diverging and decreases toward the extensional far end. Along the moving sidewalls the flow velocity is zero due to friction (i.e., no-slip boundary condition).

Deformation Monitoring
We use an automated light and camera setup to monitor surface deformation by means of top view images for qualitative description as well as sets of stereo images for further 3D displacement analysis using 3D stereoscopic DIC (Adam et al., 2005). For the qualitative surface evolution description, we use a Nikon D200 (10.2 Mpx) DSLR camera mounted above the experiment apparatus. Two Nikon D810 (36 Mpx) DSLR cameras are aligned on a horizontal bar above the model surface with an angle of ca. 30° with respect to each other providing a setup for stereoscopic images and subsequent 3D DIC analysis.
For imaging the internal model deformation, we run models with a reference thickness ratio T BD = 1, that is, a 3 cm thick sand layer overlying a 3 cm thick viscous layer, with initial flat topography (M1-0*), and models with additional quartz sand layers of 1 cm (M1-1*), and 3 cm (M1-3*) on top of the contractional domain in a medical XRCT scanner (64 slice Siemens Somatom Definition AS X-ray CT-scanner). The models are scanned with a 20-min interval (corresponding to 3.3 mm divergence or 1% maximum extension). Additionally, we scanned an experiment with T BD = 1 with no material (i.e., upper, and lower crust and additional sand layer) placed on the compressional part of the apparatus (M1-0E*). The resolution of the scans is 512 by 512 pixels for each slice with a size of 37 by 37 cm (∼0.72 by 0.72 mm/px). For every time step, models must be halted and scanned three times for quantitative DVC. Each individual scan has a duration of ca. 1 min with a time gap of ca. 30 s between the three scans. The CT scanner required such cool-down periods which were kept at a minimum by reducing the horizontal scan area (68 by 31 cm; Figure 2b) to avoid ongoing deformation in the time-dependent viscous layer.

Data Analysis and Post Processing
For the quantitative analysis of surface deformation, we use the StrainMaster module from the commercial DaVis image correlation software (Ver. 8.4, LaVision). Since image correlation yields time-series incremental displacement fields on a fixed Eulerian grid, we use Lagrangian summation, to obtain finite deformation. Based on displacement fields, the infinitesimal strain tensor E is obtained from: with the following components in 2D: E xx and E yy are the normal strain components in x and y directions, respectively and E xy and E yx are the shear strain components. The maximum principal stretching axis (largest eigenvalue of the strain tensor) defines the magnitude of the maximum normal strain, E max , where: Since E max is invariant to the coordinate system, we use it to describe deformation in our experiments where rigid body rotation occurs around a fixed axis.
Unlike cutting wetted analog models at their final stage, CT scanning provides a non-destructive manner to gain insights into transient internal deformation during experiment runs. As in our surface strain analysis, we use commercial DaVis software (Version 10.2, LaVision) for quantifying internal 3D deformation captured in CT scan images. The DVC module of this software computes 3D displacement fields by cross correlating intensity patterns from subsequent volumetric (CT scanned) data sets collected over time (Adam et al., 2013;Poppe et al., 2019;Zwaan, Schreurs, & Adam, 2018). The data volume is divided into smaller voxel sub-volumes to determine local displacement vectors by identifying similar intensity patterns in subsequent volume data sets. In contrast to 2D correlation, peaks are represented by a spherical blob in 3D and for each sub-volume the best matching position is associated with displacement vector components dx, dy, and dz from the initial to the deformed volume. As each sub-volume yields only one displacement vector, spatial vector resolution is increased by placing sub-volumes next to each other at a distance d, smaller than the voxel edge length L to create a sub-volume overlap. Making use of a multi-pass cross correlation, DaVis first applies a computationally efficient Fast Fourier Transform (FFT) algorithm for coarse sub-volumes to predict displacement fields for subsequent iterations which are done using a direct correlation (DC) algorithm.
Using the smallest XRCT slice thickness of 0.6 mm/px the volumetric data consists of voxels with dimensions of 0.72 by 0.72 by 0.6 mm/px. Since DVC requires cubic voxels with the same edge length along all three axes, the images must be re-processed before conducting quantitative analysis. In our anisotropic voxels the y axis (slice thickness of 0.6 mm/px) provides the best resolution in the volumetric data set, and we adjust pixel sizes (0.72 mm by 0.72 mm) in the xz-plane to obtain isotropic 0.6 mm sided voxels. We apply a nearest neighbor resampling algorithm in ImageJ to increase the pixel number in the xz-plane ( Figure 4) without creating artificial pixel intensity values due to interpolation. As the models are scanned multiple times per time step, we further stack the voxel data to increase intensity values. In addition, this step reduces noise created by the high X-ray reflectivity Zirshot ceramic beads and improves the image quality. For DVC analyses, we use sub-volume sizes of 128 pixels with an overlap of 50% for the first step (FFT). Further iterations (DC) subsequently decrease sub-volume sizes to 64, 32 and 16 pixels with overlaps of 50% and 75% for the final step. Where sub-pixel displacements are not captured by the discrete correlation peak, fitting of a Gaussian curve restores accuracy. The final sub-volume displacement vectors are assembled to construct incremental (here dt = 20 min) 3D displacement fields. Incremental displacements are in the order of 1-2 mm for all three displacement components with an uncertainty of 0.03 and 0.06 mm in planes parallel and normal to the CT scan slices, respectively.
Processing of 3D DIC data follows the protocol of T. C. . Like DIC data, we process data from DVC analyses in MATLAB (Version 18b) and Python (Version 3.7) to quantify deformation.

Results
The complete model series consists of 10 experiments (see Table 1). Before quantifying surface deformation, we describe the first-order evolution based on top view images and CT scans, to provide qualitative insights into the deformation evolution. First-order deformation features such as rift propagation show strong similarities in all models. Hence, we use a model with the reference brittle-ductile thickness ratio T BD = 1 and an additional sand layer of 3 cm on top of the contractional domain, to describe the overall deformation evolution. Note, that an identical setup with an initially flat topography has been previously described in  and T. C. .

General Rift Evolution for Models With Gravitational Load
Differences in strain evolution among models with or without additional gravitational load, and with varying brittle-ductile ratios are inconspicuous and hence, difficult to assess. This highlights the necessity to further analyze deformation quantitatively. However, we briefly describe characteristic stages of rift evolution under the influence of an additional gravitational load. Figure   . Positions of all extracted 2D slices from Digital Volume Correlation analyses. X slices refer to vertical, longitudinal transects. Y slices refer to cross section transects. As an example, in both slices the colors refer to the vertical displacement component Dz and white arrows indicate the total displacement, projected into the plane. All transects, including horizontal z slices can be found in an additional data publication on the GFZ data repository (T. C. Schmid, Rudolf, et al., 2022).
as the rotation axis, where bulk extension is zero (Figures 5c and 5f). Formerly isolated intra-rift fault segments have partially connected and consist of longer fault segments propagating toward the rotation axis.
The intra-rift faults show a less advanced stage of propagation when compared to models with an identical crustal configuration but flat topography (e.g., M1-0*; T. C. . Figure 6 compares rift morphologies of models M1-0* and M1-3* at the final stage. The rift morphology in model M1-0* shows three distinct pairs of conjugate normal faults with successively younger faults toward the rift center (Figures 6a, 6b and 6d). Each fault generation bounds a graben structure resulting in increasing subsidence toward the rift center.
In contrast, rifting in model M1-3* shows a less advanced stage with only two generations of conjugate normal faults present (Figures 6e, 6f and 6h). Note that thrusting in the contractional domain is entirely suppressed in model M1-3*. XRCT scan imagery from the corresponding cross section reveals that the gravitational load entirely inhibits thrust development (Figures 6c and 6g).

Propagation of Rift Boundary Faults
We investigate lateral rift propagation by tracking the position of the rift boundary fault tips with respect to the rotation axis over time. For this, we set a threshold of 10% of cumulative surface strain identical to T. C. . Figure 7 shows the mean rift tip position (i.e., the mean position of the two conjugate rift boundary faults) for all 10 experiments. For all models, rift propagation follows a distinct trend with rapid rift lengthening in an early stage and slower propagation in a second stage. All models show onset of rift propagation at around 1.3% of maximum extension. However, rift onset in models with an initially flat topography tends to occur earlier compared to models with gravitational loading. For models with a brittle-ductile thickness ratio T BD = 2, the mean rift propagation path is nearly identical but with a decreasing brittle-ductile ratio, the onset of rift propagation is successively retarded for models with a gravitational load.

Inward Migration of Fault Activity
The delayed rift propagation in experiments with gravitational loading is clearly visible in surface strain-rate maps obtained from DIC analysis ( Figure 8) which we use as a proxy for fault activity. T. C. Schmid, Schreurs, and Adam (2022)   Presumably, the rift structure in model M1-3* will become as mature as that in model M1-0* at a later time.

Vertical Surface Motions
The cumulative vertical displacement, Dz (Figure 9), is used to compare the surface topography in the final stage of all experiments with an additional gravitational load (Mx-3 and Mx-1) with the corresponding experiments with an initially flat topography (Mx-0). Note that these motions only indicate net uplift and subsidence and do not necessarily represent the absolute topography. The cumulative horizontal displacement vector field is also indicated in Figure  The zero-elevation line (black line) separates zones of uplift and subsidence which, in the case of initially flat topography models, coincides with domains separating contraction from extension. Note that the extensional domain shows patterns of diffuse subsidence outside the rift. This regional subsidence is due to the homogeneous thinning of the lower viscous layer, which is placed on top of rigid basal setup that does not allow isostatic compensation. The maximum rift and pop-up structure width depends on the absolute thickness of the brittle layer. Maximum subsidence values decrease with lower T BD ratio.
In models with an additional gravitational load (Figures 9f-9j), parts of the extensional domain, that are located farther away from the rotation axis show identical subsidence patterns when compared to models with identical  Solid lines correspond to models with different brittle-ductile thickness ratios, T BD , but initial flat topography. Dashed lines correspond to models with different T BD and an additional gravitational load of either 1 or 3 cm. Note that a decreasing T BD leads to subsequent rift propagation delay in models with an additional gravitational load but has no effect on models with initial flat topography. and has a sharp boundary to the adjacent uplift bulge, where maximum uplift values are around 5 mm. The uplift bulge continues far into the extensional domain and eventually uplift is large enough to suppress net subsidence in the rift near the rotation axis ( Figure 9j). Note that the horizontal vector field in models with an additional gravitational load deviates increasingly from the reference state (i.e., models with initial flat topography) as T BD decreases.

Internal Deformation Analysis
We now focus on quantification of the internal deformation in XRCT scanned models with an intermediate brittle-ductile thickness ratio T BD = 1 (M1-0E*, M1-0*, M1-1*, and M1-3*). Specifically, we analyze the flow field in the viscous domain by means of cumulative displacement vector fields and focus on the horizontal component of rift-parallel flow (i.e., along the model y-axis). Figure 4 gives an overview of the data slices extracted from DVC analyses of all four XRCT scanned models. For each model, we extract individual 3D displacement components on five x-slices (i.e., yz-planes; x 1 − x 5 ) and five y-slices (i.e., xz-planes; y 1 − y 5 ). Additionally, we extract displacement vectors from horizontal z-slices in both the brittle and viscous domains (i.e., xy-planes; z 1 and z 2 ). For all slices, the displacement vectors represent the total 3D displacement, projected into the plane of the corresponding CT slice. The entire slice compilation is available in an additional data publication (T. C. Schmid, Rudolf, et al., 2022). Here, we present a representative selection of displacement data to document quantitative 3D deformation features in our XRCT-scanned models.

Displacement Components in Vertical Rift Perpendicular Transects y 1 and y 5
Apart from differences in the absolute values of displacements, similar observations are made in the y-slices of all XRCT scanned models. Therefore, we only show the y 1 and y 5 slices of model M1-0* in Figure 10  The magnitudes of the maximum horizontal displacement component Dx are smaller than expected based on the maximum divergence velocity (20 mm/hr on each side). This results partly from the fact that the applied maximum divergence velocity is defined along a circular segment. However, this deviation is on the order of 1% and therefore negligible. More importantly, the DVC analysis does not capture the entire model length (see Figure 2b) and is also limited by a fixed data mask width along the y axis and can therefore not account for the increasing differences in model width along the y axis over time. Therefore, the DVC analysis cannot capture maximum displacements near the mobile sidewalls (i.e., model widths, −150 and 150 mm) which results in smaller maximum values of the horizontal displacement component Dx.
In the contractional domain (Figures 10a-10c), Dx values range between −5 and 5 mm (Figure 10a). The left part (i.e., from x = −150 to 0 mm) moves to the right (indicated by light red colors) while the right part (i.e., from Generally, the horizontal displacement components Dx and Dy are more pronounced in the extensional domain (Figures 10d-10f). The Dx component shows diffuse stretching in the viscous layer, whereas deformation in the brittle layer localizes above the viscous seed, resulting in a graben system (Figure 10d). Note that the displacement polarity changes at the brittle-ductile interface which results in pockets of reversed viscous flow below the rift axis (from ca −50 to 50 mm). The change in the horizontal displacement direction at the brittle-ductile interface is coupled with the vertical displacement Dz (Figure 10f). In the viscous layer below the rift zone, the viscous seed material rises, while the surrounding viscous and brittle material subsides due to gravitational unloading in the rift center. The combination of horizontal Dx, and vertical Dz displacement causes a convection-like cell indicated by the displacement vectors. The horizontal displacement component Dy (Figure 10e) has maximum values in the viscous layer below the rift. In contrast to the contractional domain, the zone of pronounced displacement is restricted to the model center.

Rift-Axis Parallel Horizontal Displacement Dy in the x 3 Longitudinal Transect
Longitudinal transects of all XRCT scanned models show that the rift-axis parallel horizontal displacement (Dy) is focused in the viscous domain ( Figure 11). D y displacement values range between −5 and 5 mm for models M1-0E* ( Figure 11a In experiments with an additional gravitational load, a continuous channel of rift-axis parallel displacement in the viscous domain is subsequently localized below the topographic step (Figures 11c and 11d). In Model M1-1*, maximum displacement values of 15 mm are located below the topographic step at the rotation axis and gradually decrease to 0 at ∼200 mm distance from the rotation axis (Figure 11c). This is also valid for model M1-3* (Figure 11d). However, in this experiment maximum values of rift-axis parallel displacement double and reach up to 30 mm.

Vertical Displacement Dz in the x 3 Longitudinal Transect
The increased rift-axis parallel displacement in models with an additional gravitational load also affects the vertical displacement component Dz (Figures 12c and 12d). In the brittle layer, negative values as large as −8 mm reflect subsidence within the rift. For models without an additional gravitational load (M1-0E*, M1-0*; Figures 12a and 12b), subsidence of the brittle layer occurs across the entire extensional domain and reaches maximum values near the far end at ∼500 mm from the rotation axis. In the viscous layer, maximum values of Dz = 8 mm occur near the far end at ∼500 mm and indicate upward flow simultaneously with the thinning of the brittle layer (Figures 12a and 12b; see also Figure 10f). In model M1-0* (Figure 12b), subsidence in the brittle layer decreases toward the rotation axis and changes to uplift in the contractional domain. With an additional

Discussion
In this study, we used DVC in combination with 3D DIC to depict and understand the coupling of near-surface rift development and deep-seated crustal flow in a rotational rift setting. To simulate different initial crustal configurations, we used a series of brittle-ductile thickness ratios T BD , where increasing values indicate increasing crustal strengths. Our results show that DVC applied to XRCT images provides a unique opportunity to gain insights on the time-dependent internal deformation of analog models. This method allows for a fully quantitative description of both the surface deformation and the deformation of internal layers representing lower parts of the crust.
A previous study by Zwaan, Schreurs, and Adam (2018) reports DVC analyses performed on crustal scale analog rift models undergoing orthogonal extension. They showed that the interaction between two segmented rift branches caused substantial horizontal displacements in the viscous layer (lower crust analog), parallel to the direction of rift propagation. This shows that slight deviations (i.e., an offset between two parallel rift branches) from a cylindrical setup can result in considerable out of plane displacements. In the case of a rotational model setup, rift-axis parallel displacement components are even more pronounced.  used the same experimental apparatus as this study and gave preliminary insights on rift-axis parallel displacement in a rotational system. Although they applied 2D DIC on a single longitudinal XRCT transect (Figure 10b in , rift-axis parallel displacement patterns in our model M1-0* agree well with their findings. We also observe a similar enhanced horizontal Dy displacement in our model M1-0E* (see Figures 11a , 11b and 14b) from which we conclude that the contracting boundary sidewalls in our experimental apparatus do not influence displacement patterns in the extensional domain. Rift-axis parallel flow has also been documented in lithospheric scale numerical models (e.g., Le Pourhiet et al., 2018;Van Wijk & Blackman, 2005), highlighting the importance of out of plane displacement components even in models with reasonably cylindrical boundary conditions.

The Influence of Rift-Axis Parallel Horizontal Flow on Rift Propagation
The presence of a horizontal pressure gradient (resulting from a topographic load at the contractional end of the model domain) influences the evolution of the rift structure. Figures 7 and 8 illustrate the delay in rift evolution in terms of rift tip propagation, the timing of the abandonment of rift boundary faults and successive initiation of faulting activity along intra-rift faults (inward migration). To understand the cause of such delayed evolution, we consider deformation in the viscous layer. Figures 11 and 12 show the rift-axis parallel horizontal displacement Dy, and the vertical displacement component Dz, respectively. The subsiding brittle layer in the extensional domain is coupled with an upward flow of the viscous material underneath the rift axis ( Figure 12). For increasing additional gravitational loads, the Dz upward flow component appears to be perturbed and overprinted by a strong rift-axis parallel horizontal flow component Dy (Figure 12d). T. C.  propose that intra-rift faults initiate early as antithetic faults but remain largely passive while slip is accommodated along the rift boundary faults. It is only after they reach the brittle-ductile interface, where shear stresses are high (Corti et al., 2010), that intra-rift faults accommodate considerable displacement and concomitant fault activity along the rift boundary faults decreases. Hence, the upward flow of the viscous material increases and facilitates further activation of intra-rift faults. Zwaan, Schreurs, and Adam (2018) showed the inhibiting effect of syn-rift sedimentation on the ascent of viscous material below the rift structure. Although we do not consider the effect of sedimentation, the presence of a large additional gravitational load induces a dominant rift-axis parallel horizontal displacement component Dy, which is strong enough to overprint or inhibit vertical upward flow of the viscous material. This effect delays rift evolution, as expressed by retardation of rift tip propagation, inward fault migration, and activation of intra-rift faults.

Effect of Deep-Seated Deformation on Surface Topography and Rift Evolution
Surface analyses show that increasing the brittle-ductile thickness ratio T BD has no recordable influence on rift evolution in models with an initial flat topography (Figures 7 and 9). However, in models with an additional gravitational load, we observe enhanced vertical motions close to the rotation axis ( Figure 9) as well as horizontal motions. Such horizontal motions at the surface are especially prominent in models with an intermediate or low brittle-ductile thickness ratio T BD (Figures 9g-9j). For XRCT scanned models with T BD = 1, this observation is consistent with the internal displacement fields (Figures 10-12) and is likely governed by increased flow in the viscous layer of the model and, the additional gravitational load. This becomes evident when rewriting Equation 1 in terms of the volumetric flow rate where μ is the viscosity that is required to maintain horizontal flow in a channel with thickness h, driven by a lateral pressure gradient dp/dx, along the horizontal length of the channel to minimize lateral variations in thickness (e.g., Kruse et al., 1991;Turcotte & Schubert, 1982). The proportionality μ∼h 3 in Equation 6 requires lower viscosities for smaller channel thicknesses to maintain a constant volumetric flow rate Q, when other parameters are kept constant. Since we use models with constant viscosity, a larger channel thickness (i.e., lower T BD ) results in enhanced horizontal flow and hence, enables more surface deformation. We must emphasize therefore that it is the absolute thickness of the viscous layer (in addition to viscosity and pressure gradient), rather than the brittle-ductile thickness ratio, that is the principal determinant for the magnitude of lower crustal flow.
This enhanced lower crustal channel flow is ultimately expressed by increased vertical motions in models with an intermediate to low T BD (Figures 9 and 13). Regardless of the brittle-ductile thickness ratio, local contraction results in a zone of uplift in models with flat initial topography (Figures 9a-9e and 13a). However, in models with an additional gravitational load (Figures 9f-9j and 13-13e), the intensity of vertical motion at the surface increases in experiments with thick viscous layers (i.e., T BD = 0.5). While models with thin viscous layers (i.e., T BD = 2) still yield net uplift in the contractional domain (Figures 9f and 13b), models with viscous layers of intermediate thickness (i.e., T BD = 1) show a net subsidence (Figures 9g-9j and 13c-13e). In these models the channel thickness is sufficient to maintain a rift-parallel, horizontal channel flow component which transports material away from the contractional domain, counteracting the tectonic uplift caused by contraction (Figures 9g-9j  and 13c-13e). This results in a net subsidence in the contractional domain and, where the horizontal flow component decreases, the excess material in the lower crust enables synchronous uplift during extension while rifting occurs. In the extreme case (i.e., thick viscous layer; T BD = 0.5), uplift raises the propagating rift tip above its original elevation, resulting in net uplift of the rift system when crossing the uplift barrier (Figures 9j and 13e).

Conceptual Model for Rift-Axis Parallel Flow
The  ∼30 mm (Figure 14a) 1 predicts 0 mm. Interestingly, the differences between measured and theoretical values are systematically larger with increasing additional gravitational loads. This discrepancy can be attributed to two causes: First, the rotational opening of the experimental apparatus provokes rift-axis parallel flow away from the rotation axis even without an additional gravitational load (Figures 10, 11 and 14). Second, the strong rift-axis parallel displacement component measured at the model surface (gray arrows Figures 9g-9j) indicates that the contact between the brittle and viscous layers is not a fixed boundary (with respect to the viscous layer) and therefore allows higher flow velocities in the viscous domain than expected in the case of pure Poiseuille flow. Figure 14b shows Dy profiles through the brittle and viscous layers at distinct positions along the rift axis. The effect of the pressure gradient (i.e., parabolic shape of the displacement profile) is evident in the viscous layer near the topographic step and decreases farther away from the rotation axis in the extensional domain. The brittle layer shows uniform Dy values through the entire layer thickness indicating that the brittle layer moves as a rigid block. It is only in a thin layer at the brittle-ductile interface, where a displacement gradient in the granular material occurs. This transition is likely due to energy dissipation and is expressed by enhanced shear strain, which reduces displacement values in the brittle layer compared to the viscous layer ( Figure 14c). The displacement gradient dDy/dz reveals anticlockwise rotation at the brittle-ductile transition indicating a top-to-left shear sense (red arrows; Figure 14c).
The general equation for the velocity, u, in a planar Couette flow of a linear viscous fluid in a channel with thickness h, viscosity μ, and an applied horizontal pressure gradient dp/dx is: where u 0 is the relative velocity between channel and bounding plates (Turcotte & Schubert, 1982). In the case of Poiseuille flow, the bounding plates are fixed and u 0 = 0, such that the lateral pressure gradient produces a parabolic velocity profile with maximum values in the center of the channel (Figure 14d, left panel). We have shown above that rift-axis parallel viscous flow in our experiments is not confined by fixed bounding plates (i.e., the brittle layer moves as well), which ostensibly agrees with a Couette flow component in addition to pressure-gradient driven Poiseuille flow. However, there is a fundamental difference between the flow described in Equation 7 and the rift-axis parallel flow observed in our experiments. For Couette flow, the upper boundary moves and drives flow in the viscous layer below. In the case of our experiments, however, it is the viscous layer that exerts a drag force on the brittle layer, resulting in energy dissipation due to shearing at the brittle-ductile interface (Figure 14d, middle panel). It is this shearing, which ultimately causes a difference in the Dy values in the viscous and brittle layer (Figure 14d, right panel).

Coupling Between Brittle and Viscous Layers During Rift-Axis Parallel Displacement
The horizontal displacement field at the surface of models with an additional gravitational load is systematically different to those of all models with an initial flat topography (Figure 9). Such modified horizontal displacements indicate a rift-axis parallel overprint, similar to those observed in longitudinal DVC transects (Figures 11 and 14). However, as discussed above, the rift-axis-parallel displacement component at the model surface is reduced compared to the displacement component in the viscous part of the model suggesting a degree of brittle-viscous decoupling between the two model layers (sensu Zwaan et al., 2019).
To further investigate the mechanical coupling of brittle and ductile deformation during rifting, we use the quantitative results for models with T BD = 1 and compare the evolution of the rift-axis parallel mean values in the brittle ( Figure 15a) and viscous domain (Figure 15b). We restrict the investigation to a 10 cm wide central area (i.e., close to the X 3 longitudinal transect; Figures 11 and 14) in the extensional domain to exclude displacement effects due to rotational opening.
The rift-axis parallel displacement in the brittle domain (Figure 15a) shows positive and negative trends for models with (i.e., M1-1* and M1-3*) and without (i.e., M1-0* and M1-0E*) an additional gravitational load, respectively. However, for model M1-0*, this negative trend is negligible, within the error. Note that positive and negative values indicate mean displacements toward and away from the location of maximum opening, respectively. In the viscous domain, mean values all show positive trends indicating rift-axis parallel flow toward the location of maximum opening. Plotting the evolution of the ratio between mean values from the brittle and viscous domain (Figure 15c) reveals the degree of brittle-viscous coupling between the two model layers with respect to rift-axis parallel displacements. If the mean Dy values of the brittle and viscous layer in one experiment are totally coupled, this ratio should be 1. Any ratio below 1 indicates an increasing contrast between the rift-axis parallel displacement in the brittle and viscous layer. Positive ratios <1 indicate displacement directions in the viscous domain are larger than those in the brittle domain. Any ratio <0 indicates opposing displacement directions in the brittle and viscous layers. Note that absolute values of mean displacements in the viscous domain are always larger than in the brittle domain and hence, all ratios plot between −1 and 1.
In the case of an additional gravitational load (i.e., M1-1* and M1-3*), both displacement directions are away from the rotation axis toward the extensional domain. The proportional increase in both brittle and viscous Dy values is expressed by a constant ratio of ∼0.2. This suggests that for both experiments, rift-axis parallel displacements in the brittle domain are reduced by a factor of 5 in contrast to the viscous flow below the rift. Interestingly, this ratio only increases by a small amount for a larger gravitational load. The increased flow in the viscous domain seemingly affects the rift-axis parallel displacement in the brittle domain equally, regardless of the gravitational load.
For models without an additional gravitational load, negative ratios indicate opposing displacement directions in the brittle and viscous layers. However, the evolution documents an early phase (up to ca 100 min) with continuously increasing displacement contrasts suggesting a first phase of increased decoupling (Figure 15c). After this early phase, ratios increase for both models as Dy mean values in the viscous domain increase faster than Dy mean values in the brittle layer decrease (i.e., become more negative). We speculate that in this later phase the intensified viscous flow (away from the rotation axis) is due to the increased differential topography along the rift axis which counteracts the opposing Dy displacement in the brittle domain. Eventually, the viscous flow likely diminishes the negative Dy mean values in the brittle layer such that it flips the direction and hence, increases the brittle-viscous coupling.   . In our experiments M1-0E* and M1-0* the corresponding rift tip propagation does not deviate substantially from those of other experiments (Figure 7). It therefore appears that the increase in coupling between Dy mean values in the brittle and viscous layers (after ca 120 min) occurs too late during the model runs to significantly influence lateral rift propagation. Nevertheless, the evolution of the coupling between Dy mean values in the brittle and viscous layers in models M1-0E* and M1-0* highlights that it is not a static but rather dynamic process which may change during progressive rotational extension.

The Role of Rift-Axis Parallel Crustal Flow in Natural Settings
Our XRCT scanned models show a substantial amount of rift-axis parallel displacement in the viscous domain ( Figures 10, 11, 14 and 15). Also in the absence of a gravitational load, we observe rift-parallel, central flow channels with maximum Dy values up to 5 mm ( Figure 11) and mean values of ∼1 mm ( Figure 15) during 40 mm of rift extension. Based on our scaling, this translates to maximum and mean flow velocities of 2 and 0.4 mm/a in nature, respectively. Such values may seem minor when compared to the orthogonally oriented divergence velocity (i.e., 8.8 mm/a) but are they negligible? Various studies have documented, an out-of-plane flow component in the lower crust and its implications on rift evolution (e.g., Amato et al., 2004;Gautier et al., 2008;Little et al., 2011;MacCready et al., 1997;P. D. Clift, 2015).
Based on gravitational potential and buoyancy considerations, Clift (2015) proposes rift-parallel lower crustal flow from the continental crust toward the oceanic rift propagator in the Woodlark Basin as well as in the South China Sea, similar to observations in our analog models. In addition, Little et al. (2007Little et al. ( , 2011Little et al. ( , 2013 conclude for the Woodlark basin that, based on structural data from the D'Entrecasteaux Islands, regional, lower crustal flow parallel to the rift-axis must have been largely decoupled from the general NNW-SSE extension direction between the Woodlark microplate and the Australian plate (Little et al., 2011). The exhumed gneiss domes on the D'Entrecasteaux Islands show depth-dependent trends in stretching lineation orientation, indicating rift-axis parallel lower crustal flow toward the Papuan Peninsula in contrast to a dominant rift-axis normal stretching lineation in upper crustal levels (Little et al., 2011). The proposed flow direction (based on shear fabrics) for rift-parallel flow is opposite of what we observe in our rotational models and may be explained by asthenospheric inflow into the thinned crust (Abers et al., 2002;Little et al., 2011;Mondy et al., 2018), which we do not consider in our models. The rotational experiment by Mondy et al. (2018) documents a prominent rift-axis parallel flow component in the asthenosphere following the tip of the propagating rift toward the rotation axis. Simultaneously, vertical upward flow compensates for the thinning crust and impedes rift-axis parallel flow in the lower crust away from the rotation axis (Mondy et al., 2018). Interestingly, in early stages of the experiment presented in Mondy et al. (2018), the rift-axis parallel flow component of the asthenosphere is directed toward the extensional far end. However, the authors do not state its effect on lower crustal flow. In our models, we use a rigid basal setup and therefore cannot account for the effect of the asthenosphere on the general model evolution.
In this regard, one could argue that our models represent an early stage of continental rifting, where the role of the asthenospheric flow is negligible.
Despite these two contrasting suggestions from Clift (2015) and Little et al. (2011) for the rift-axis parallel lower crustal flow direction, field data from gneiss domes in core complexes on the D'Entrecasteaux Islands is key for uncovering coeval but orthogonal extension directions at varying crustal levels (i.e., rift-axis parallel lower or mid crustal flow and rift-axis perpendicular stretching in the upper crust). The documentation of two orthogonally oriented stretching lineations in such gneiss domes also led other authors to propose the existence of rift-axis parallel lower crustal flow, orthogonal to the regional extension direction (e.g., Amato et al., 2004;Gautier et al., 2008;MacCready et al., 1997). MacCready et al. (1997) document rift-axis parallel northward flow of the middle crust in the Ruby Mountains, Nevada as a response to a rift-axis parallel extension gradient with increasing upper crustal stretching values toward the north. Similarly, Gautier et al. (2008) document in the Nigde Massif, Turkey, rift-axis parallel, lower crustal flow (perpendicular to the extension direction) and determine a northerly flow because of an extension gradient with increasing values toward the northern end of the basin. In both cases, the driving mechanism of regional flow seems to be an extension gradient in which rift-axis parallel lower crustal flow compensates greater thinning of the upper crust. Both settings have kinematic boundary conditions that cause rift-axis parallel lower crustal flow toward areas with higher extension, identical to what we document in our models. Gautier et al. (2008) describe such rift-axis parallel flow as regional-scale lateral channel flow, in contrast to a more local-scale lateral inward flow toward the rift axis, which compensates thinning of the upper crust parallel to the extension direction (Block & Royden, 1990). Similarly, we observe rift-perpendicular ( Figure 10d) and rift-axis parallel (Figures 10e and 11) horizontal flow components in addition to vertical upward flow (Figures 10f and 12) in our models. Eventually, the combination of all displacement components (as documented in nature as well as in our models) give rise to complex 3D flow patterns which can transport a considerable amount of material in or out of a 2D plane. This effect clearly must be considered when estimating crustal extension from 2D rift-perpendicular cross sections (P. D. Clift, 2015).

Conclusion
We performed a series of brittle-viscous analog models of rotational rifts with an extension gradient in which gravitational loading causes a pressure-gradient driven rift-axis parallel horizontal flow in the lower viscous layer, which is synchronously with rift development in the upper brittle layer. From these experiments we gain insights on the interplay of brittle and ductile deformation in natural rifts that have formed under the influence of along strike extension gradients. The analysis of ductile (deep-seated) and brittle (near-surface) deformation by means of DVC and 3D stereo DIC, respectively provides a quantitative and comprehensive picture of progressive deformation in the experiments. This approach allowed us to detect mechanisms in the viscous model layer that influence the near-surface brittle deformation and eventually cause delayed rift propagation. This results in important implications for the role of crustal flow in natural rift settings, summarized below.
• When compared to models with an initial flat topography, gravitational loading influences the timing of inward migration of faults and retards activation of intra-rift faults. Hence, models with an additional gravitational load show less advanced rift maturity after an identical run time when compared to flat topography models. However, this does not affect lateral rift propagation, regardless of the brittle-ductile thickness ratio T BD . • Models with initial flat topography develop rift structures with subsiding basins in the extensional domain and uplifted areas confined by thrusts in the contractional domain. These patterns are identical for all brittle-ductile thickness ratios. Additional gravitational loading, however, enhances vertical displacements and changes patterns of subsidence and uplift. Subsidence due to gravitational loading drives horizontal flow in the viscous model layer, which results in uplifted areas in the extensional model domain. • For thicker lower viscous layers, flow channels are larger and hence provide more space for enhanced horizontal ductile channel flow. If an additional gravitational load is applied, this creates more pronounced uplift in the extensional domain. • Enhanced rift-axis parallel flow in the viscous model layer overprints vertical upward flow near the rift axis.
Upward flow originally helps to activate intra-rift faults and drives inward fault migration. Hence, viscous rift-axis parallel delays inward migration of activity along faults and enables boundary faults to be active for a longer time. This results in a comparatively less mature rift stage at a given time. • Rift-axis parallel flow in the viscous model domain is driven by a pressure gradient due to differential topography (Poiseuille flow). Viscous flow drags the brittle upper layer as a rigid block. Shearing occurs along the brittle-ductile interface where energy dissipation results in reduced Dy displacement values compared to Dy displacement values in the viscous layer. • DVC depicts enhanced horizontal rift-axis parallel flow patterns in the viscous domain for models with additional gravitational loads. Correlation between growing mean values of rift-axis-parallel displacement components in the viscous and brittle layer, shows different degrees of coupling between these two domains. For models with an additional gravitational load, stable mechanical coupling (of some degree) between the brittle and ductile layer is expressed by identical displacement directions toward the extensional domain. In models without additional gravitational loads, mechanical coupling is transient and changes during the model run, showing that this is a dynamic process, rather than a static one. • Rift-axis parallel ductile flow, documented in rift settings with an extension gradient, results in material transport out of the rift transect. When scaled to nature, displacements from our models confirm substantial outflow of lower crustal material parallel to the rift axis. This mechanism may explain discrepancies between subsidence estimates based on upper crustal extension and subsidence in profiles and emphasizes the three dimensionality of rifting.

Data Availability Statement
Rheological measurements of the brittle and viscous materials used in this study are available in the form of open access data publications via GFZ Data Services (T. Schmid et al., 2020aSchmid et al., , 2020bZwaan, Schreurs, Ritter, et al., 2018, respectively). An additional open access data publication on the GFZ Data Service (T. C. Schmid, Rudolf, et al., 2022) provides additional images and movies of our models. Links to these data sets are provided in the reference list.