Analog Models of Lithospheric-Scale Rifting Monitored in an X-Ray CT Scanner

Rifting and continental break-up are fundamental tectonic processes, the understanding of which is of prime importance. However, the vast temporal and spatial scales involved pose major limitations to researchers. Analog tectonic modeling represents a great means to mitigate these limitations, but studying the complex internal deformation of lithospheric-scale models remains a challenge. We therefore present a novel method for lithospheric-scale rifting models that are uniquely monitored in an X-ray CT scanner, which combined with digital image correlation (DIC) techniques, provides unparalleled insights into model deformation. Our first models illustrate how the degree of coupling between competent lithospheric layers, which are separated by a weak lower crustal layer, strongly impacts rift system development. Low coupling isolates the upper crust from the upper lithospheric mantle layer below, preventing an efficient transfer of deformation between both layers. By contrast, fast rifting increases

However, we face various challenges when studying the evolution of rifts and rifted margins. Firstly, these tectonic systems extend over vast areas so that studying them in detail is highly challenging. Furthermore, direct access to rifts and rifted margins is often very limited since large parts of these systems are situated below sea level, and on top of that, are often covered by thick post-rift sedimentary sequences. These challenges have been somewhat overcome through the use of geophysical techniques, especially seismic methods, and the drilling of deep boreholes. However, possibly the most significant and lasting obstacle on our way to understanding rift processes is caused by the long timescales at which they occur, as it is simply not possible to directly observe the development of rift systems over millions of years.
In order to overcome this major obstacle, researchers have since long applied scaled analog models, which allow the simulation of tectonic systems and associated processes that act on vast spatial and temporal scales within a matter of hours or days in a modest laboratory, using relatively simple analog materials at limited financial costs. A large number of analog modeling studies focusing on rift tectonics have been published, and these studies have yielded a wealth of useful insights, as summarized in various reviews and overview papers (Allemand & Brun, 1991;Bahroudi et al., 2003;Beslier, 1991;Brun, 1999Brun, , 2002Brun et al., 2018;Corti et al., 2003;Corti, 2012;Koyi, 1997;McClay, 1990McClay, , 1996Michon & Merle, 2000Naylor et al., 1994;Vendeville et al., 1987;Zwaan et al., 2019;. Of these analog modeling studies, a large part focuses on crustal-scale structures, which is permissible when considering small-scale processes or very early rift stages, but these relatively simple models cannot capture the full rift system. Especially important when considering the full scope of rifting and continental break-up is the influence of the asthenosphere, that is, the isostatic compensation that asthenospheric flow exerts on a rift system. In analog modeling studies that simulate lithospheric-scale rifting in a normal gravity field, the isostatic influence of the asthenosphere is generally simulated with a dense low-viscosity fluid (syrup or honey), on which the model materials representing the lithosphere float (e.g., Allemand & Brun, 1991;Autin et al., 2010Autin et al., , 2013Benes & Davy, 1996;Benes & Scott, 1996;Beniest et al., 2018;Brun & Beslier, 1996;Capelletti et al., 2013;Molnar et al., 2017Molnar et al., , 2018Nestola et al., 2013Nestola et al., , 2015Samsu et al., 2021;Sun et al., 2009; Figure 1). The rheological layering of lithospheric-scale models in an enhanced gravity field (i.e., centrifuge models) is very similar, although the materials are generally of a higher viscosity in such experiments (e.g., Agostini et al., 2009;Corti, 2008;Corti et al., 2013).
However, all of these previous lithospheric-scale rifting models have an important disadvantage in that the internal model evolution has so far not been monitored and analyzed in real time. Where the internal deformation of crustal-scale models can in some cases be partly observed via a glass sidewall (while accounting for boundary friction effects), the internal deformation of lithosphere-scale models remains elusive since the model lithosphere generally contains different (sticky viscous) layers that hamper direct observation of the model interior through glass. These viscous layers also hamper the use of surface structures and topography (including the topography of the base of the lithosphere, e.g., Nestola et al., 2013Nestola et al., , 2015 as reliable indicators for internal deformation since these layers often decouple internal model deformation from the surface structures (e.g., Allemand et al., 1989;Molnar et al., 2017Molnar et al., , 2018Zwaan et al., 2019Zwaan et al., , 2021Zwaan, Chenin, et al., 2022; Figure 1). Sectioning is another  Allemand and Brun (1991). option to gain insight into a model's interior, but this is a destructive method that can only be used to interpret the model evolution through the reconstruction of markers, and only after completion of the model (e.g., Beslier, 1991;Brun & Beslier, 1996). Furthermore, sectioning of lithospheric-scale models containing various viscous layers is often a challenging task since the models are highly unstable and deform easily while making the sections. The same goes for removing the crustal layers to scan the topography of the base of the crust after the model run (Nestola et al., 2015).
So far the only reliable non-destructive method to truly peek into otherwise opaque analog models materials during an experimental run is X-ray CT-scanning (e.g., Colletta et al., 1991, Sassi et al., 1993Schmid et al., 2022aSchmid et al., , 2022bSchreurs et al., 2003;Holland et al., 2011;Zwaan, Schreurs, & Adam, 2018). This technique, which basically relies on the visualization of density differences within the model materials (e.g., induced by layering and dilatation along faults) has been around for a number of decades and has been applied in various studies to analyze crustal-scale rifting models (e.g., Schmid et al., 2022aSchmid et al., , 2022bZwaan et al., 2016Zwaan et al., , 2020Zwaan, Schreurs, & Adam, 2018). Yet to our knowledge CT-scanning has so far never been used to analyze the much more complex deformation in lithospheric-scale models, providing a strong incentive to push forward the state of the art in modeling.
The aim of this paper is therefore to present the unique possibilities of a newly designed analog model set-up for studying lithospheric-scale rifting, which allows for the first time the monitoring of model-internal deformation via X-ray CT scanning methods. The resulting CT imagery does not only allow a qualitative assessment of internal model evolution, but we can even quantify strain within the model through the use of digital image correlation (DIC) techniques normally applied for strain analysis of model surfaces. Furthermore, next to simple orthogonal rifting our new set-up enables us to simultaneously simulate lateral movements as well, opening many possibilities for simulating different tectonic regimes such as oblique extension and (oblique) basin inversion. The results from four rifting models presented in this paper demonstrate the general viability and potential of this new modeling procedure.

Materials
Similar to previous analog modeling studies of lithospheric-scale rifting (e.g., Allemand & Brun, 1991;Brun & Beslier, 1996;Nestola et al., 2015; Figure 1), we apply both brittle and viscous materials to simulate the brittle and ductile parts of the lithosphere, as well as the asthenosphere (Figure 2c). Assuming a stable four-layer continental lithosphere prior to deformation, the model consists from top to bottom of a strong (brittle) upper crust (UC), a (weak) ductile lower crust (LC), a strong (brittle) upper lithospheric mantle (ULM), a weak (ductile) lower lithospheric mantle (LLM) and a very weak (ductile) asthenosphere with the following relative thicknesses of the lithospheric layers from top to bottom: 2:1:1.5:2 ( Figure 2c and Table 1).
For both the brittle upper crust and brittle upper mantle layer, we use FS900S feldspar sand produced by Amberger Kaolinwerke (https://www.quarzwerke.com) ( Table 1). This angular feldspar sand as provided by the company has a dominant grain size between 100 and 250 μm (65% of grains, see grain size distribution in Willingshofer et al., 2018). In contrast to other modeling studies (e.g., Beniest et al., 2018), we do not have the means to sieve and remove the fine fraction of this FS900S sand to the 100-250 μm range prior to use. However, new ring shear tests completed at GFZ Potsdam show that not removing the finer fraction of the feldspar sand does not significantly affect the internal friction angle . The density of the feldspar sand is ca. 1,300 kg/m 3 when added into the model by means of sieving from a height of ca. 30 cm.
The ductile lower crust and ductile lower lithospheric mantle are simulated using viscous mixtures of SGM-36 silicon oil (Polydimethylsiloxane or PDMS) formerly produced by Dow Corning, now part of Dow Chemical (www.dow.com), and F120 corundum sand acquired from Carlo AG (www.carloag.ch) ( Table 1). The PDMS has a bulk density of 965 kg/m 3 (Rudolf et al., 2016), whereas the corundum grains in the corundum sand have a bulk density of 3,950 kg/m 3 . The grain size of the corundum sand lies between 88 and 125 μm. By mixing these components in different ratios, we obtain viscous materials with different (near-Newtonian) rheology and varying densities. We use a mixture with a density of ca. 1,300 kg/m 3 and a viscosity of 6 × 10 4 Pa·s (viscous mixture 1) to simulate the lower crust, while a mixture with a density of ca. 1,400 kg/m 3 and a 8 × 10 4 Pa·s (viscous mixture 2) for the lower lithospheric mantle (Zwaan, Schreurs, Ritter, et al., 2018; Table 1). These four layers representing the lithosphere float on top of a Glucosweet 44 glucose syrup layer simulating the asthenosphere (Figure 2c and Table 1). This viscous glucose syrup is acquired from ADEA (https://www. adea-srl.it/) and has a density of ca. 1,450 kg/m 3 . The higher density of the syrup with respect to the overlying model materials prevents the latter materials from sinking into the syrup. Rheometer test performed at University of Roma TRE show that the syrup has a near-Newtonian rheology with a relatively low viscosity of ca. 65 Pa·s, when compared to the viscosity of 8 × 10 4 Pa·s of the lowermost viscous layer (Schmid et al., 2022c;Zwaan, Schreurs, Ritter, et al., 2018; Table 2). This viscosity difference is representative of the viscosity difference between the lithosphere and asthenosphere in nature (Figures 1 and 2c).
In order to mark layering (both within sand layers, as well as between brittle and viscous layers) on CT imagery (Section 2.4), we add thin layers (≤1 mm) of "Liaver beads," that is, foamed glass beads produced by Liaver (www.liaver.com) (grain size between 100 and 300 μm, density = ca. 530 kg/m 3 , Warsitzka et al., 2019). The Liaver beads have very different X-ray attenuation coefficients with respect to the other analog model materials and clearly stand out on CT images. Since the Liaver beads have similar internal friction angles as the feldspar sand (Table 1), and are only used as thin layers within brittle model materials or to mark the transition between the syrup and the lower viscous layer, their use is not considered to have any major impact on model deformation. Feldspar sand properties after . b Liaver beads properties after Warsitzka et al. (2019). c Pure PDMS rheology after Rudolf et al. (2016) and Zwaan, Schreurs, Ritter, et al. (2018). d Rheology of viscous mixtures 1 and 2 after Zwaan, Schreurs, Ritter, et al. (2018). e Glucose syrup rheology after Schmid et al. (2022c). f Viscosity value holds for model strain rates <10 −4 /s. g Power-law exponent n (dimensionless) represents sensitivity to strain rate.

Set-Up
The general modeling set-up principle is similar to that of previous analog models of lithospheric-scale rifting with a lithosphere floating on a heavy fluid, with the exception that our set-up specifically allows for oblique rifting and X-ray CT-scanning (Figures 1 and 2). Our new lithospheric-scale modelling set-up is based on a conceptual design by Zwaan (2017) and has been constructed by engineers from IPEK, Ostschweizer Fachhochschule in Rapperswil, Switzerland. It includes a number of components that are added to the highly versatile NAMAZU machine, also developed by IPEK, which has been used for various previous analog modeling studies at the Tectonic Modeling Laboratory of the Institute of Geological Sciences at the University of Bern, Switzerland (e.g., Alonso-Henar et al., 2015;Fedorik et al., 2019;Klinkmüller, 2011;Schori et al., 2021;Zwaan et al., 2016Zwaan et al., , 2020Zwaan et al., , 2021Zwaan, Chenin, et al., 2022;Zwaan, Schreurs, Adam, 2018).
Our new set-up includes a large basin containing the glucose syrup and is installed on the NAMAZU machine ( Figure 2). A rectangular framework of inner sidewalls that contain the materials representing the model lithosphere is hung down in this basin (Figures 2a and 2c). These inner sidewalls are fixed to bars that on their turn are connected to computer-controlled motors. Motors Y1 and Y2 can induce orthogonal extension of the model materials contained within the inner framework, or, in combination with simple-shear motion (imposed by motor X1), can produce oblique rifting (Figures 2d-2f). The direction of (oblique) rifting is defined as angle α, which is the angle between the normal to the long model axis and the combined motion directions of the motors ( Figure 1f). Importantly, as the model is deforming, the syrup can freely flow below the inner model sidewalls, allowing for isostatic compensation at the base of the simulated lithosphere that analog models with a fixed base lack (e.g., Zwaan et al., 2019). The short inner sidewalls are equipped with hinges and consist of overlapping plates to accommodate the (oblique) rifting applied in the models, while still containing the model materials representing the lithosphere within the inner framework (Figures 2a and 2f). Importantly, the experimental set-up is constructed from X-ray transparent materials that allow for optimal X-ray CT-scanning.

Model Preparation and Parameters
Preparing the models involves a number of steps including freezing of the syrup to stabilize the model materials during model construction. All steps are described in detail in the supplementary material (i.e., in a GFZ data publication by Zwaan & Schreurs, 2023b). The details of the four models completed for this study are provided below and in Table 2. It may be noted that the models presented in this study do not represent a full systematic study, but rather aim to explore the potential provided by our new modeling approach.
Model A serves as a reference model. It involves orthogonal rifting and the total lithospheric thickness is half that of the subsequent models (i.e., the thicknesses of the upper crust, lower crust, upper lithospheric mantle and lower lithospheric mantle were 10, 5, 7.5 and 10 mm, respectively). The model does not contain any structural inheritance to localize deformation under the standard divergence velocity (10 mm/hr). The total model duration amounts to 150 min for a total of 25 mm of divergence, and the model is CT-scanned.
Model B and C also involve orthogonal rifting at the same standard divergence velocity as Model A, but have a standard model thickness of 20, 10, 15 and 20 mm for the upper crust, lower crust, upper lithospheric mantle and lower lithospheric mantle, respectively. Another contrast to Model A is the insertion of a weakness or "seed," a rod of viscous material at the base of the brittle upper lithospheric mantle layer that runs along the whole length of the model, to localize deformation. Such seeds have been used in various previous modeling studies to simulate linear weaknesses such as crustal shear zones (e.g., Le Calvez & Vendeville, 2002;Osagiede et al., 2021;Schmid et al., 2022aSchmid et al., , 2022bZwaan et al., 2016Zwaan et al., , 2021Zwaan, Chenin, et al., 2022), but we apply them in the modeled upper lithospheric mantle instead ( Figure 1c). Furthermore, Model B undergoes a total of 40 mm of divergence over a 4-hr period, where Model C, which is a rerun of Model B in the CT scanner, undergoes and extra hour of deformation for a total of 50 mm of divergence.
Finally, CT-scanned Model D involves the same standard layering as Models B and C, but initially undergoes α = 45° oblique rifting at the standard divergence velocity of 10 mm/hr for the first 195 min (32.5 mm of divergence). Subsequently, the divergence velocity is tripled to 30 mm/hr for the remaining 105 min of the model run, covering an additional 52.5 mm of divergence, so that the total divergence amounts to 85 mm.

Scaling
Applying analog modeling scaling laws (e.g., Hubbert, 1937;Ramberg, 1981;Weijermars & Schmeling, 1986) demonstrates the suitability of our models for simulating continental rifting processes (detailed scaling calculations can be found in Appendix A). Specifically, the (standard) models have a geometric scaling ratio of 6.7·10 −7 , so that 1 cm corresponds to 15 km in nature, and our 6.5 thick lithospheric layering represents a ca. 100 km thick lithosphere in nature. The standard divergence velocity in our Models B, C and D (10 mm/hr) translates to 2.6-3.1 mm/yr in nature, whereas the faster velocity in the second phase of Model D (30 mm/hr) translates to 7.8-9.3 mm/yr in nature. The 10 mm/hr divergence in Model A, which has a thinner simulated lithosphere, amounts to 10.6-12.4 mm/yr. These values are close to plate divergence velocities observed in natural rift systems (e.g., the East African Rift, Saria et al., 2014), and also dimensionless scaling ratios are very similar between model and nature, demonstrating that scaling requirements are fulfilled (Appendix A).

Model Monitoring and Analysis
We use various techniques to monitor and analyze our models and the deformation they undergo. Firstly, we use a rig containing three high-resolution 36.3 megapixel D810 Nikon cameras, one of which provides top view images, and the other two inclined (stereoscopic) images of the model surface ( Figure 2b). These cameras are computer-programmed to take pictures every 30 s, creating time-lapse series. At the same time, we apply a timer for the lighting on one side of the model to switch on and off every 30 s, so that we obtain two time series, of images with and without shade (each with 1 min intervals). Note that we apply a thin <1 mm thick, 4 × 4 cm surface grid made of corundum sand to highlight surface deformation and we sprinkle coffee powder on the surface for optimizing DIC analysis.
The time series images with shade provide visual clues about surface deformation evolution, whereas the time series images without shade are ideal for 2D (map view) displacement and strain analysis. This strain analysis is done using digital image correlation (DIC) DaVis 10.1 software from LaVision (www.lavision.de), which compares surface patterns on images (if needed after correction for their inclination) from different time steps and detects horizontal surface displacements (displacement vectors) (e.g., Adam et al., 2005;Boutelier et al., 2019). The result is a detailed and quantified analysis of surface model deformation over time, which we use in particular for generating incremental and cumulative maximum normal strain (MNS) data. These maximum normal strain data represent the deformation of the longest axis of the strain ellipse, and serves as a proxy for extensional deformation (Schmid et al., 2022a(Schmid et al., , 2022bZwaan et al., 2021;Zwaan, Chenin, et al., 2022).
Furthermore, the inclined images, preferably without shade, allow for stereoscopic reconstruction of the model topography over time through Agisoft Photoscan photogrammetry software (www.agisoft.com) (e.g., Maestrelli et al., 2020). By comparing the different distortions of the images from the stereoscopic cameras, the software can reconstruct the 3D model surface for selected time steps. These 3D surfaces are georeferenced using markers with known locations placed on the model set-up, to produce properly georeferenced digital elevation models (DEMs) over time. These DEMs are subsequently post-processed in QGIS (www.qgis.org), an open source geographic information system (GIS) software.
Next to 3D photography and the quantitative surface analyses it enables, the key innovation in terms of model monitoring and analysis in our study is the application of X-ray CT-scanning techniques (e.g., Colletta et al., 1991;Holland et al., 2011;Sassi et al., 1993;Schmid et al., 2022aSchmid et al., , 2022bSchreurs et al., 2003;Zwaan, Schreurs, & Adam, 2018, Figure 2b Finally, vertical CT section time series from selected transects in orthogonal rifting Models A and C are used for 2D DIC analysis in DaVis, providing the first-ever truly quantified insights into internal deformation of analog models of lithospheric-scale rifting. This analysis is done on time series of vertical CT sections on which, similar as done for model surface imagery, displacements are traced. We extract horizontal and vertical in-section (y-and z-) displacements, as well as vorticity data. This vorticity data indicates the degree to which displacements deviate from the nearby displacement field, revealing in-section displacement along faults and shear zones.

Orthogonal Rifting Model A
The orthogonal rifting model A involves a simulated lithosphere with half the standard thickness and no seed in the upper mantle layer (see Section 2.3). Results of this model, which is stretched at a divergence velocity of 10 mm/hr during 150 min, are presented in Figure 3. Topography analysis (Figures 3a and 3b) reveals the development of significant boundary effects in the shape of (half-)grabens along the longitudinal sidewalls, whereas no surface deformation is visible in the center of the model. Furthermore, these boundary effects are strongest toward the mid-point of the longitudinal sides.
The same general pattern of deformation is shown by the results of the DIC analysis (Figures 3c and 3d). The cumulative maximum normal strain data ( Figure 3c) highlights that normal faulting is most active along the sidewalls, as is horizontal displacement, which is illustrated by the cumulative horizontal displacement data (Dy in Figure 3d). These Dy data in Figure 3d also reveal that the regions immediately adjacent to the short ends of the model, close to where the short overlapping sidewalls moved apart, undergo some displacement. A further indication of this differential displacement is provided by the slight warping of the surface grid, visible in the background of the DIC results (Figures 3c and 3d).
CT imagery provides additional insights into the internal development of Model A (Figures 3e and 3f). The layering within the simulated lithosphere is not readily visible, but the transition between the model lithosphere and model asthenosphere can be visualized (Figure 3e). The final stage CT section of Model A reveals the boundary effects normal faulting along the longitudinal sidewalls within the model lithosphere, and the rising of the model asthenosphere ( Figure 3f).

Orthogonal Rifting Model B
Model B involves a standard lithosphere thickness, contains a seed in the upper lithospheric mantle layer, and is stretched at a rate of 10 mm/hr over a period of 4 hr ( Figure 4). We analyze the model topography and surface deformation only, since the model is not CT-scanned.
The results from the topography analysis show how the model, in stark contrast to Model A (which lacked a seed), starts developing a slight central depression at the surface at about 60 min of deformation ( Figure 4b).
After 75 min of deformation, two graben structures start to appear at the short ends of the model (Figure 4c). At this point in time, technical issues with the camera rig prevented the acquisition of time-lapse imagery for ca. 90 min, hence subsequent topography data are only available from t = 180 min on (Figure 4d). These topography data show the presence of a double graben system, flanking a central depression, as well as some deformation (boundary effects) along the longitudinal sidewalls. Both grabens are well developed, and their combined surface expression is fairly symmetrical, and remains so until the end of the model run ( Figure 4e).
The incremental normal strain data obtained from DIC analysis highlights normal faulting over the preceding 5 mm of divergence and reveals that the double graben system seen at t = 75 min in the topography analysis results (Figure 4c) is in fact already established at t = 60 min, and is initially somewhat better developed near the short ends of the model (Figure 4h). The vectors displaying horizontal displacement in the DIC maps show that the central depression between the double grabens remains in place during the model run, whereas the materials on both side of the double graben system move steadily outward (Figures 4f-4j).

Orthogonal Rifting Model C
The topography and DIC analysis results from Model C, which is a rerun of Model B monitored in the CT scanner (with 60 min extra runtime), are presented in Figure 5. In addition to CT-scanning, we also perform the same topography and DIC analysis on Model C as done for Model B, allowing for a direct comparison between both models (Figures 4 and 5). Furthermore, we use the CT imagery for a qualitative assessment of model-internal deformation ( Figures 6 and 7), as well as for a quantification of internal deformation through DIC analysis on selected CT sections (Figure 8).

Topography and Model Surface DIC Analysis
The topography analysis of Model C reveals a similar surface structure evolution as in Model B (Figures 4a-4e and 5a-5e), with an initial central depression forming, followed by the development of grabens flanking this central depression. Like in Model B, the grabens are initially best developed toward the short ends of the model, before covering the whole length of the model (Figures 4c-4e and 5c-5e). A major difference with Model B, however, is that the double graben arrangement in Model C is highly asymmetric nearer to the short ends of     Figures 4b, 4g and 5b, 5g). The normal faulting initiates at the short ends of the model before growing toward the model center, and is in fact rather symmetrical in the earlier stages of model evolution (Figure 5g). Only later on, parts of the lower graben are abandoned so that the upper graben could become the dominant structure at the model surface (Figures 5h-5j). The DIC results highlight that the bulk of deformation is accounted for by the grabens along the model's central axis (Figures 5f-5j). The additional grabens in the upper right corner, as well as the boundary effects along the longitudinal sidewall, which are well visible on topography data (Figure 5e), are shown to only localize minor amounts of normal faulting (Figures 5f-5j). An important detail is indicated by the displacement arrows: the central area between the double grabens is only (nearly) stationary (i.e., no vectors visible) where this area is flanked by two active grabens (Figures 5f-5j).

CT Analysis
The acquisition of CT imagery and additional DIC analysis on CT sections allows us to analyze the internal model evolution of Model C (Figures 6 and 7), with an emphasis on the differences between the symmetric and asymmetric rift structures previously identified by the topography and DIC analysis (Sections I and II, respectively, locations of which are shown in Figure 5).
Section I, taken near the center of Model C (Figure 5), shows the development of a symmetric double rift structure (Figures 6a-6g, and 7a, 7c). In the earliest stages of rifting (t = 30 min, Figure 6a), the seed in the simulated upper lithospheric mantle layer localizes deformation. As rifting progresses, this leads to the development of a graben structure in the upper lithospheric mantle with a central depression at the model surface (Figures 6c, 6d, and 7c). We subsequently observe the development of two upper crustal grabens flanking this central depression, which are linked to the graben in the upper lithospheric mantle layer by low-angle shear zones (LASZs) in the simulated lower crust (Figures 6e and 7c). At this point in time (t = 120 min), the upper lithospheric mantle layer is split in two parts, and also a minor graben forms away from the central double graben structure ( Figure 6e). As rifting proceeds, this minor graben does not develop significantly, in contrast to the continuously developing double grabens (Figures 6f, 6g, 7a, 7c). The rifting also leads to substantial thinning of the modeled lower crust (and of the model lithosphere in general), which is accommodated by the strong rise of the lower lithospheric mantle layer and the modeled asthenosphere from t = 120 min on, putting the lower lithospheric mantle layer and upper crustal layer in near-direct contact towards the end of the model run (Figures 6g and 7c).
Section II is taken ca. 20 cm away to the right from Section I in Model C (Figures 5a-5e), and illustrates the development of an asymmetric double rift structure (Figures 6h-6n, 7b, 7d). Similar to Section I, initial deformation localizes in the upper lithospheric mantle layer (Figures 6b and 6i). However, an early asymmetric transfer of deformation to the modeled upper crust via a low-angle shear zone in the lower crustal layer causes the development of a graben on one side of the central model axis (Figure 6j). As rifting continues and the upper lithospheric mantle is split apart, this graben continues evolving (Figures 6k-6n, 7c, 7d). In the meantime, a smaller secondary graben appears on the other side of the central model axis, whereas a curved fault develops in the center of the model (Figures 6m and 6n). Another minor graben develops away from the central model axis as well (Figure 6k-6n). Around this point in time (t = 120 min), the model asthenosphere and lower lithospheric mantle start to rise considerably, compensating the thinning of the lower mantle (and lithosphere in general) (Figure 6l-6n, 7b, 7d). As the model continues evolving, the initial graben remains dominant while the lithosphere continues thinning up to the point that the lower lithospheric mantle and upper crust are in contact (Figures 6n and 7d). It is also worth noting how much the upper lithospheric mantle to the left in Section II is strongly tilted upward when compared with the symmetric setting of Section I (Figures 7c and 7d).

DIC Analysis on CT Sections
In addition to visual inspection and interpretation of the CT imagery, we also apply DIC techniques to quantify model-internal deformation occurring in both Section I and II of Model C (Figure 8). Firstly, cumulative horizontal displacement data (Dy) reveals the development of differences in Dy values over time in both Section I and II (Figures 8a-8d). The symmetric rift structure in Section I is characterized by a stationary upper crustal block between two lithospheric domains that are diverging (Figures 8a and 8b). This arrangement is reminiscent of the stationary domain seen in DIC analysis of top view imagery in Models B and C (Figures 4f-4j and 5f-5j). By contrast, the asymmetric rift structure in Section II lacks such a stationary domain, as the diverging domains are in direct contact (Figures 8c and 8d). However, both sections have an additional low Dy domain in the lower part of the modeled lithosphere as well (Figures 8a-8d). Note that the boundaries between the different Dy domains in Sections I and II delineate the low-angle shear zones in the lower crust as identified on CT imagery (Figures 6g  and 6n), as well as shear zones deeper in the lithosphere (Figures 8a-8d). In the case of the asymmetric rift structure in Section II, we observe a large-scale shear zone running through the whole lithosphere (Figures 8c and 8d).
In addition to horizontal displacement values, we also present vertical displacement values (Dz) and cumulative vorticity data obtained through DIC analysis of CT sections (Figures 8e-8l). The Dz data reveal that the overall model lithosphere does not undergo significant vertical displacement (Figures 8e-8h), with the exception of the central rift zone, where the model surface is subsiding, and the model asthenosphere is strongly rising (as also visible in the profiles in Figure 7). The cumulative vorticity data indicate changes in displacement patterns, highlighting shear zones in the model, as well as their sense of shear (Figures 8i-8l). Although they are not highly precise, the shear zones indicated by the cumulative vorticity results are in good agreement with those previously interpreted through the CT imagery and the Dy data (Figures 6, 8a-8d, 8i-8l).

Oblique Rifting Model D
We present the topographic and surface DIC analysis of CT-scanned 45° oblique rifting Model D in Figure 9. During the first 195 min of deformation at the standard divergence rate of 10 mm/hr, the model is being stretched and a depression formed along its central axis. At t = 180 min some faint deformation zones on both sides of this depression become visible on the DIC results (Figures 9a-9c and 9g-9i). However, no clear faulting is observed, which is in clear contrast to the situation in orthogonal rifting Models B and C over the same period (Figures 4a  -4d, 4f-4i, 5a-5d, 5f-5i, 9a-9c and 9g-9i).
In order to force deformation in the model, we triple the divergence rate to 30 mm/hr in the second phase of the model run, from t = 195 min onward. As a result, faulting starts to localize along the deformation zones flanking the central depression (Figures 9d, 9e, 9j-9k). These faults are arranged in a left-stepping en echelon fashion, leading to the development of two zones of en echelon grabens along the central axis of the model, flanked by a number of additional graben structures that are mostly located in the top-right corner (Figure 9f). However, these additional grabens (as well as the boundary effects along the long ends of the model) only accommodate a minor part of the total deformation, as demonstrated by our DIC results (Figures 9j-9l) CT imagery, presented in Figure 10, provides insights into the internal deformation of Model D. The images show how during the first phase of the model run, deformation localizes along the seed, leading to the upper lithospheric mantle layer splitting apart (Figures 10a-10e). Yet, apart from some boundary effects along the long sidewalls, no clear faulting occurs in the upper crustal layer, as only the central depression seen on topography data develops (Figures 9a-9c and 10a-10e). The situation totally changes during the second model phase as faster rifting leads to the localization of faulting in the upper crustal layer, and the development of a symmetric double graben system (Figures 10f-10h). The en echelon nature of this double graben system, and its relation to the deformation deeper in the model lithosphere as well as the oblique divergence direction is revealed by 3D CT imagery (Figure 10i).
Due to the lateral displacement component in Model D, the same 2D DIC analysis on CT sections as shown for orthogonal extension Model C is not performed. However, we may assume that the 2D cumulative displacement and deformation patterns in Model D are broadly similar to those obtained for Section I from Model C, given the similarity in CT section view (Figures 6g and 10h).

Discussion
Our model results obtained through a novel combination of external and internal monitoring and analysis techniques provide insights into three key topics: (a) the general localization of faulting, (b) symmetric and asymmetric rift systems, and (c) the development of rift systems in 3D during oblique extension. We summarize these results in text and figures (Figures 11-13), and compare them with results from previous modeling studies. Here it must be kept in mind that our models do not represent a systematic parameter study due to the variations in model parameters we apply. Still we believe the results are sufficiently robust for a productive discussion, especially when linked to observations from previously published modeling work. We subsequently describe the limitations of our modeling approach and potential future implementations for our understanding of natural systems, as well as of other tectonic settings.

Localization of Faulting: Effects of Seed and Coupling
The first key observations concern the general localization of deformation in our models, as summarized in section view in Figure 11. The absence of deformation in the form of faulting in the upper crustal layer along the central axis of Model A without a seed, highlights that such a seed is required to localize a rift structure in this model set-up (compare Model A with Models B-D, Figures 3-11a, 11c). The need for a seed or weakness to localize deformation in brittle-viscous set-ups has been shown by various previous analog and numerical modelers (e.g., Oliveira et al., 2022;Zwaan et al., 2019), and is linked to the decoupling effect of the weak viscous layer representing the lower crust. This decoupling isolates the competent brittle layers in the model lithosphere, and deformation will simply focus along the sidewalls, which form the weakest part of the model, whereas the brittle layers simulating the upper crust and upper lithospheric mantle remain stationary and undeformed. However, in other analog modeling studies of lithospheric-scale rifting, researchers induce deformation by using a narrower set-up with U-shaped sidewalls (e.g., Allemand & Brun, 1991;Autin et al., 2010;Brun & Beslier, 1996;Nestola et al., 2013Nestola et al., , 2015. The edges of these U-shaped sidewalls act as velocity discontinuities (VDs), triggering deformation along the central model axis. But as pointed out by Zwaan et al. (2019), and even partiality visible at the short ends of our models as well ( Figure 5), faulting then tends to initiate at these VDs before growing toward the center of the model.
Furthermore, our models show that divergence velocity affects the strength coupling between the competent (brittle) layers and plays an important role in the localization of faulting as well. Orthogonal rifting models B and C generate double graben structures due to a transfer of deformation from the seed in the upper lithospheric mantle through the lower crust into the upper crust via low-angle shear zones ( Figure 6). However, such transfer of deformation does not occur in (the initial phase of oblique rifting of) Model D, where only a central depression forms, with some faulting along the sidewalls, whereas the upper lithospheric mantle layer splits apart due to the presence of the seed (Figures 10 and 11b). This difference is likely (partially) caused by the orthogonal divergence velocity component (Y1 + Y2, which determines the amount of new area created in the system) during the first phase of Model D being smaller than in Models B and C ( Table 2). As a result of this lower initial orthogonal divergence velocity in Model D, the strength of the lower crustal layer, which has a strain-rate dependent rheology, seems to be reduced. Therefore, the upper lithospheric mantle and upper crustal layers are more decoupled and localized deformation cannot be transferred between both these competent layers (e.g., Brun, 1999Brun, , 2002Michon & Merle, 2000Zwaan et al., 2019). Vice versa, faster rifting, as in Phase 2 of Model D, increases coupling, leading to the initiation of faulting in the upper crustal layer due to transfer of localized deformation from the seed upward, whereas faulting along the sidewalls is reduced (Figures 9, 10, and 11c). These findings are in accordance with previous modeling work showing that faster divergence allows (localized) deformation in the mantle to have more influence on upper crustal faulting, or to even overprint and suppress deformation otherwise occurring along weaknesses in the upper crust, leading to complex interactions and intricate rift structures Zwaan, Chenin, et al., 2022).

Symmetric Versus Asymmetric Rift Development
The orthogonal rifting Models B and C, both with seeds, provide important insights in the development of symmetric and asymmetric rift structures, as summarized in section view in Figure 12.
Both Models B and C show the development of symmetric rift structures as a result of symmetric transfer of deformation from the seed in the competent lithospheric mantle into the upper crustal layer (Figures 4 and 5, Section I in Figures 6-8, 12). As described in the preceding paragraph, a degree of coupling related to a moderate divergence rate is needed for this deformation transfer, which occurs along two shear zones in the lower crustal layer. The development of such double shear zones and associated double graben structures has been observed in both lithospheric-scale analog models (e.g., Brun & Beslier, 1996;Michon & Merle, 2000Nestola et al., 2015), as well as in crustal-scale brittle-viscous models, where instead of a seed the edge of a base plate (VD) represents a fault in the upper lithospheric mantle (e.g., Allemand et al., 1989;Michon & Merle, 2000Zwaan et al., 2019Zwaan et al., , 2021Zwaan, Chenin, et al., 2022). Moreover, they are also observed in numerical modeling studies (Chenin et al., 2018(Chenin et al., , 2020Dyksterhuis et al., 2007;Oliveira et al., 2022). In the case of a symmetric rift structure, two of these shear zones form simultaneously on both sides of the seed, causing the development of two equally sized grabens in the upper crustal layer, with a relatively undeformed "H-block" (i.e., "hanging wall block, " Lavier & Manatschal, 2006;Péron-Pinvidic & Manatschal, 2010) in between (Figures 12b and 12c).
By contrast, in some parts of Model C, the rift structure develops in an asymmetric fashion ( Figure 5, Section II in Figures, 6-8, 12). In this asymmetric system we find the early development of a single shear zone in the lower crust, inducing a single graben in the upper crust early on (Figure 12e). In later stages, a second shear zone may develop, but the initial shear zone remains dominant and evolves into a large shear zone crossing the whole lithosphere, whereas the lower crust is exhumed at the bottom of the dominant graben (Figure 12f). It is not clear why this asymmetry occurs in (parts of) Model C, especially since the structures in the equivalent Model B are fully symmetrical (compare Figures 4 and 5), but previous analog models show that the mode of rifting is influenced by the divergence rate (Brun & Beslier, 1996;Michon & Merle, 2003). We therefore speculate that the set-up we apply in our models is close to a tipping point, so that small differences in model preparation, which are known to affect analog modeling results (e.g., Schreurs et al., 2006Schreurs et al., , 2016, may locally shift the system from a symmetrical to an asymmetrical style, and vice versa. Furthermore, numerical modeling studies show that a variety of additional factors can affect the symmetry of a rift system (e.g., Brune et al., 2014;Huismans & Beaumont, 2002Huismans et al., 2005;Nagel & Buck, 2004), but a detailed comparison of our analog models with these numerical results is beyond the scope of this paper.
However, even though we cannot pinpoint with certainty what causes the shift in rift style in this study, our new models allow us to examine and compare both symmetric and asymmetric rifting situations. In both cases we find a sharp rise of the asthenosphere, isostatically compensating the strong thinning of the lithosphere, which in the center chiefly consists of upper crustal and lower lithospheric mantle material as both the lower crust and upper lithospheric mantle are practically split in two (Figure 6g, 6n, 7, 12c, 12f). This rise of the asthenosphere is a universal observation in other lithospheric-scale analog and numerical models undergoing strongly localized thinning (e.g., Allemand et al., 1989;Brun & Beslier, 1996;Brune et al., 2014;Chenin et al., 2018Chenin et al., , 2020Nestola et al., 2013Nestola et al., , 2015, and prevents the local "collapse" of the modeled lithosphere that may occur in advanced stages of rifting if no such isostatic compensation is included. Furthermore, there is no unrealistic regional subsidence of undeformed parts of the lithosphere due to stretching of the model and viscous thinning of lower crustal layer, previously observed in two-layer brittle-viscous models (e.g., Schmid et al., 2022aSchmid et al., , 2022bZwaan et al., 2020).

Oblique Extension and 3D Rift Structures
Apart from purely 2D considerations detailed in the previous sections, our orthogonal and oblique rifting models also provide insights into 3D rift development (summarized in Figure 13). Orthogonal divergence Models B and C generate symmetric and asymmetric rift systems with through-going faults parallel to the long axis of the models (Figures 4, 5, 13a, 13b). By contrast, oblique divergence Model D, after developing only limited faulting in the upper crust during the initial rifting phase, forms two series of obliquely oriented basins during the second, faster divergence phase (Figures 9, 10i, 13c). These series of left-stepping en echelon basins follow in fact the trace of the double grabens in the orthogonal divergence situation (Figure 13a). This indicates that, after increasing the divergence velocity in the second phase of Model D, the lower crust with its strain rate-dependent rheology is sufficiently strengthened to allow for the development of shear zones in the lower crust that transfer deformation (coupling) from the seed in the upper lithospheric mantle into the upper crustal layer (Figure 13a). Note that one needs to be careful when comparing model D with models B and C, due to the application of the second phase of faster rifting in Model D. Still, Model D illustrates what structures may be expected in oblique rift settings.
Indeed, within this general framework, the en echelon faults themselves are typical of oblique rift systems, as normal faults form perpendicular to the (local) extension direction (which itself may deviate from the plate divergence direction in oblique divergence settings, Withjack & Jamison, 1986). Examples of en echelon faulting in oblique divergence settings are found in lithospheric-scale analog and numerical models (e.g., Agostini et al., 2009;Autin et al., 2010Autin et al., , 2013Brune, 2014;Duclaux et al., 2020), as well as in crustal-scale experiments applying a brittle-viscous base plate set-up also reproduced the double graben style with en echelon faulting (e.g., McClay & White, 1995;Tron & Brun, 1991;Zwaan et al., 2021;Zwaan, Chenin et al., 2022). These findings indicate that some aspects of lithospheric-scale models (i.e., deformation in the crustal parts) can be well-represented by simpler crustal-scale models, which was previously suggested by Michon and Merle (2003) as well.
Finally, even though not observed in our study, there is a good possibility that, similar to their orthogonal divergence counterparts, oblique divergence models could develop asymmetric rift systems with a single zone of en echelon faulting ( Figure 13d). As in the case of the orthogonal divergence models, an asymmetric rift under oblique divergence conditions would lead to earlier exhumation of the lower crust than in the symmetric equivalent as all deformation would be localized in one graben. However, this lower crustal exhumation would itself be delayed compared to the lower crustal exhumation in the orthogonal divergence models, due to the smaller orthogonal divergence component (Y1 + Y2) applied in our oblique divergence set-up (Figures 2f and 13, Table 2).

Model Limitations, and Opportunities for Improvement and New Studies
Our novel application of CT-scanning as well as subsequent DIC analysis on CT imagery is a significant step forward in the field of analog modeling of lithospheric-scale rifting. It provides us with unprecedented insights into the internal deformation evolution of such lithospheric-scale models (Figures 3, 6, 7, 8, and 10). Especially the quantification of deformation over time in these models by means of DIC analysis on CT imagery (Figure 8) allows for direct comparisons with numerical modeling studies. Additional steps could involve the expansion of our 2D DIC analysis into the third dimension (so-called digital volume correlation or DVC analysis applied on CT volumes, Adam et al., 2013;Poppe et al., 2019;Schmid et al., 2022aSchmid et al., , 2022bZwaan, Schreurs, & Adam, 2018) to fully quantify the evolution of internal model deformation in 3D. However, the resolution (512 × 512 pixels) of standard medical CT scanners poses a limitation here. The imagery allows for excellent visual inspection and interpretation, but the resolution is relatively low for detailed DIC analysis. This limitation could however be mitigated to a degree by applying (industrial) CT scanners with higher resolution, by selecting a smaller scan window, or by applying thicker model layering. As a matter of fact, the 8-cm-thick lithosphere we apply is already much thicker than in most other analog modeling studies, and has been intentionally increased to improve the resolution of our results (compare the scans from Model A with those from Models B and C, Figures 3, 6 and 10). We estimate that a further doubling of this lithospheric thickness would be feasible.
The new set-up and general modeling approach itself has also proven itself to be successful. However, the procedure is somewhat challenging due to the limitation imposed by the requirements for CT-scanning: the materials need to be X-ray transparent, and the whole set-up needs to fit into a CT scanner (Figure 2b). As a result, the syrup basin is relatively small, so that adding the relatively thick lithosphere risks causing large displacement of syrup, in contrast to other modeling set-ups that either have much larger basins or much thinner lithospheres. Our solution to freeze the syrup for stability during model preparation means that completing a single model takes a week, limiting the modeling output (see details in Zwaan & Schreurs, 2023a). This lower output is however more than offset by the CT-scanning results, and it would be possible to prepare multiple models in parallel to double or triple the modeling output. Furthermore, this new set-up allows for easy simulation of different degrees of oblique divergence, and as such the study of 3D tectonic processes in rift systems (Figures 2f, 9, 10, and 13).
A further consideration is the lithospheric layering and model materials. Various authors have applied different lithospheric layering (e.g., Allemand et al., 1989), and we have only explored a very small part of the potential parameter space. Rerunning models with these different types of lithospheric layering in the CT scanner would allow for new insights. Moreover, a limitation in most analog models of lithospheric-scale rifting has always been that most materials generally do not undergo thermal effects. This may not be much of an issue in models such as those presented in this paper, which mainly aim at the earlier stages of rifting when thermal effects are considered of minor importance Zwaan, Chenin, et al., 2022). However, as the lithosphere starts necking, and the hot asthenosphere starts rising to the surface (Figures 6, 7, 8, 10, and 12), we should expect phase changes and variations in rheology to kick in. Therefore, the inclusion of model materials, of which the rheology is affected by temperature changes (e.g., Boutelier et al., 2003Boutelier et al., , 2012Boutelier & Oncken, 2011;Chemenda et al., 2002;Krýza et al., 2019), in CT-scanned models of lithospheric-scale rifting would be a promising avenue for future studies. This could be combined with parallel numerical modeling efforts, which will allow for both benchmarking and verification of both modeling methods (e.g., Brune et al., 2017;Panien et al., 2006;Zwaan et al., 2016), and for the addition of other parameters that are challenging to include in lithospheric-scale analog models (e.g., surface processes). In this context, it is also important to stress that the current method is only suitable for the modeling of magma-poor rift settings, since magmatism can strongly affect rift evolution (e.g., Buck, 2004Buck, , 2006. Finally, we see great opportunities beyond the field of rifting, as our new method could also be applied for the simulation and detailed analysis of (oblique) collisional tectonics (e.g., Calignano et al., 2015Calignano et al., , 2017Luth et al., 2010;Sokoutis & Willingshofer, 2011;Willingshofer et al., 2013;Willingshofer & Sokoutis, 2009) and lithospheric-scale basin inversion (Boutelier et al., 2003;Cerca et al., 2005;Gartrell et al., 2005).

Conclusion and Future Outlook
In this paper we present a novel method for the modeling of (magma-poor) lithospheric-scale rifting processes, which can be uniquely monitored in a CT scanner. We show how the application of CT-scanning and DIC analysis on CT imagery provides unparalleled details regarding the (internal) evolution of such lithospheric-scale models, and our first set of models provides the following first-order insights: • The degree of coupling between the competent layers in the lithosphere in the presence of a weak lower crust has an important influence on the development of the modeled rift system. Low coupling due to slow rifting isolates the upper crust from the upper lithospheric mantle layer below, preventing an efficient transfer of deformation between both layers. By contrast, fast rifting increases coupling and allows deformation in the mantle to efficiently induce deformation in the upper crust. • When sufficient coupling occurs and deformation is transferred from the mantle into the upper crust, we observe either the development of a symmetric or asymmetric (double) rift system in our models. Although the reason why the system may develop one or the other style is not fully clear, we obtain detailed insight in the evolution of either style. • Oblique divergence leads to en echelon graben arrangements, and a delayed exhumation of lower crustal material. We speculate that these structures may occur in both symmetric or asymmetric (double) rift systems.
These initial models provide an incentive to further simulate rifting processes on lithospheric scales, and to apply advanced techniques to extract as much information as possible out of these. There is indeed a broad range of opportunities for improvements to our new modeling method, for example, involving the testing of different lithospheric layering and the inclusion of new materials with temperature-sensitive rheology, for new studies within and beyond the field of rift tectonics.

Appendix A: Model Scaling
Applying analog modeling scaling laws (e.g., Hubbert, 1937;Ramberg, 1981;Weijermars & Schmeling, 1986) demonstrates the suitability of our models for simulating continental rifting processes. The parameters for our scaling calculations are provided in Table A1.
Here ρ*, h* and g* represent the density, length and gravity ratios respectively, and the equation yields a σ* of ca. 3·10 −7 . For brittle materials, scaling is relatively straightforward due to their strain-rate independent rheology. Most important is the internal friction angle of the feldspar sand we used, which is very similar that of upper crustal rocks (e.g., Byerlee, 1978; Table A1). The dynamic similarity in the brittle layers can be validated by calculating the R s ratio (Bonini et al., 2001;Mulugeta, 1988;Ramberg, 1981): R s = gravitational stress/cohesive strength = (ρ·g·h)/C. Here ρ is the density, g the gravitational acceleration, h the height and C cohesion of the brittle materials. Our calculations yield R s values for the upper crust and upper lithospheric mantle that very similar between model and nature (Table A1), confirming the adequate dynamic scaling of the brittle behavior of the materials in our models.
Scaling of viscous materials is more complex due to their strain-rate dependent rheology. In the case of viscous materials, the stress ratio (σ*) and viscosity ratio (η*) produce the strain rate ratio ( *) according to the ensuing formula: * = σ*/η* (Weijermars & Schmeling, 1986), which is in the order of 1·10 10 . Subsequently, we can obtain the velocity ratio v* and time ratios t* through the following equations: * = v*/h* = 1/t*. As a result, our standard divergence velocity of 10 mm/hr translates to 2.6-3.1 mm/yr in nature, which is similar to divergence rates observed in nature (e.g., Saria et al., 2014).
To test dynamic similarity of the viscous layers, we apply the Ramberg number R m (Weijermars & Schmeling, 1986): R m = (ρ·g·h 2 )/(η·v). The R m values obtained for the different viscous layers in our models are very similar to the R m values for their natural counterparts (Table A1). Together with the R s values the Ramberg numbers indicate that scaling requirements are reasonably fulfilled for our standard model runs.
In the case of Model A, with half the layer thickness, the scaling is still adequate, but the divergence rate of 10 mm/hr, translates to 11 mm/yr. In the case of the enhanced divergence velocity in Model D, scaling remains according to the standard model parameters in Table A1, except that the 30 mm/hr divergence velocity during the second model phase translates to 9.6 mm/yr. These divergence velocities of 11 and 9.6 mm/yr are somewhat at the high end in nature. However, in the case of Model A did not significantly affect model evolution due to the decoupling between the model layers and the lack of a seed, whereas the increased coupling is needed in Model D, in order to force localization. Note that in the case of Model D, a slower divergence (e.g., 15 mm/h, translating to 4.8 mm/yr) may also have sufficed, but we chose to somewhat force the system in order to obtain a result.

Data Availability Statement
Supplementary material is available in the form of publicly accessible data publications, stored at the servers of GFZ Data Services. New details on the material properties of the glucose syrup and feldspar sand used in these models can be found in Schmid et al., 2022c ( ). The initial inspiration for developing this new modeling approach was sparked by discussions with Yago Nestola, who was working at the University of Parma at the time. We would furthermore like to thank the engineers from IPEK Rapperswil (Dario Niggli, Peter Eichenberger, Rudolf Kamber, and Theodor Wüst) for working out our initial design and for constructing the components for the new model set-up. We are grateful to Nicole Schwendener from the Institute of Forensic Medicine (University of Bern) for helping us scan our models. Discussions with Nicolas Molnar (RWTH Aachen University) helped to improve the model design. Thanks also to Timothy Schmid from the Tectonic Modeling group at the Institute of Geological Sciences (University of Bern) for helping with the DIC analysis and the syrup testing, and to Kirsten Elger and Florian Ott from GFZ Potsdam for helping us prepare the GFZ data publications containing the supplementary data of this paper (Schmid et al., 2022c;Zwaan & Schreurs, 2023a, 2023b. We are grateful to Donna Shillington and an anonymous reviewer for providing useful feedback that helped us to improve the paper, and we thank editor Laurent Jolivet for providing additional input and guiding the overall review process. This project was funded by the Swiss National Science Foundation, Grant 200021-178731 "4D Analog modeling of oblique rifts and obliquely rifted margins" (https://data.snf. ch/grants/grant/178731). FZ was furthermore supported by a GFZ Discovery Fund fellowship. Open Access funding enabled and organized by Projekt DEAL.