Asperity Size and Neighboring Segments Can Change the Frictional Response and Fault Slip Behavior: Insights From Laboratory Experiments and Numerical Simulations

Accurate assessment of the rate and state friction parameters of rocks is essential for producing realistic earthquake rupture scenarios and, in turn, for seismic hazard analysis. Those parameters can be directly measured on samples, or indirectly based on inversion of coseismic or postseismic slip evolution. However, both direct and indirect approaches require assumptions that might bias the results. Aiming to reduce the potential sources of bias, we take advantage of a downscaled analog model reproducing megathrust earthquakes. We couple the simulated annealing algorithm with quasi‐dynamic numerical models to retrieve rate and state parameters reproducing the recurrence time, rupture duration and slip of the analog model, in the ensemble. Then, we focus on how the asperity size and the neighboring segments' properties control the seismic cycle characteristics and the corresponding variability of rate and state parameters. We identify a tradeoff between (a–b) of the asperity and (a–b) of neighboring creeping segments, with multiple parameter combinations that allow mimicking the analog model behavior. Tuning of rate and state parameters is required to fit laboratory experiments with different asperity lengths. Poorly constrained frictional properties of neighboring segments are responsible for uncertainties of (a–b) of the asperity in the order of per mille. Roughly one order of magnitude larger uncertainties derive from asperity size. Those results provide a glimpse of the variability that rate and state friction estimates might have when used as a constraint to model fault slip behavior in nature.

Supporting Information: Supporting Information may be found in the online version of this article.
In nature, a, b, and D c can be estimated indirectly from geodetic and seismological observations, in combination with numerical modeling of earthquake sequences (e.g., Frank et al., 2017;Johnson et al., 2006;Thomas et al., 2017).Rate and state friction parameters of velocity weakening areas (i.e., asperities) are inferred mainly using numerical models based on the theoretical framework of rate-state friction.These models are powerful tools, as they allow testing how different values of friction parameters influence the seismic cycle behavior (e.g., Barbot, 2019).(a-b) and D c are successfully identified when there is a good match between numerical model behavior and a variety of observables such as the recurrence time of a characteristic earthquake of a given magnitude (e.g., Thomas et al., 2017), or interseismic, coseismic, and postseismic deformation patterns (e.g., Barbot et al., 2012;Premus et al., 2022).
On the contrary, for velocity strengthening patches the product (a−b)σ can be determined knowing the static Coulomb stress change ΔCFF from a coseismic slip model and the postseismic slip evolution (e.g., Perfettini & Avouac, 2004): where V pl is the plate convergence velocity, and V + is the sliding velocity on the fault immediately after the earthquake.V + is constrained by the inversion of the geodetic data (e.g., Lin et al., 2013;Perfettini et al., 2010;Thomas et al., 2017) or with a dense aftershock catalog (Frank et al., 2017).Successively, (a-b) can be derived assuming that σ is known (e.g., Zhao et al., 2022).
den Hartog et al. (2021) found that geodetically derived-and laboratory measurements of (a-b) match particularly well for different areas of the Longitudinal Valley Fault (Taiwan), indicating good agreement across scales.This evidence corroborates the concept that laboratory-derived frictional properties of fault samples are essential to predict fault behavior on real Earth.Such achievement, if supported by additional instances (e.g., different faults and tectonic environments), has fundamental implications for producing realistic rupture scenarios that eventually lead to solid hazard assessments.However, adding more instances is challenging as it requires both (a) having access to fault samples from a large fault that can display a variety of slip behaviors (i.e., creeping zones and locked patches) and (b) the availability of long geodetic time series that ideally include multiple seismic cycles and lluminate postseismic slip.
Here we build on this topic, retrieving rate and state friction parameters that reproduce the observed "seismic" behavior of a downscaled laboratory megathrust.Our controlled experimental approach circumvents some uncertainties inherent natural systems, for example due to fault geometry (Moreno et al., 2009), temperature (Den Hartog et al., 2021), effective normal load (Frank et al., 2017;Thomas et al., 2017), slip and rake uncertainties associated to rupture models used for estimating ΔCFF (e.g., Wang et al., 2018).We profit from a novel seismotectonic analog model mimicking the first-order characteristics of a megathrust seismic cycle under well-constrained boundary conditions (Mastella et al., 2022).We focus on success and pitfalls when trying to reproduce fault behavior with a simple single asperity configuration.We show how the frictional properties and 10.1029/2023JB026594 3 of 16 dimensions of areas surrounding seismic asperities contribute to tuning the seismic cycle characteristics.We also show how the size of the asperity controls model behavior and how the corresponding rate and state parameters should be tuned to take into account the observed variations across different asperity sizes.

Foamquake: A Downscaled Seismogenic Subduction Megathrust
We use a novel seismotectonic analog model to reproduce multiple seismic cycles of a generic subduction megathrust in a scaled (i.e., model to nature length dimensions scaling factor ∼3 × 10 −6 ) and realistic way (Mastella et al., 2022).The reliability of this model is ensured by a variety of analog earthquake source parameters proportionalities (e.g., seismic moment-duration, moment-area) as compared with those derived from large subduction earthquakes (Murotani et al., 2013), together with the emergence of rupture cascades and superimposed cycles, similar to space-time rupture patterns observed at subduction megathrusts in nature (Philibosian & Meltzner, 2020).The analog model features a foam rubber (here the name Foamquake; Figure 1a), wedge-shaped plate (145 × 90 × 20 cm 3 ) that is underthrusted by a 10° dipping conveyor belt (analog of the subducting plate) driven at constant velocity V pl = 0.01 cm/s.A 1 cm thick layer (the subduction megathrust analog) embeds a rectangular asperity made of rice grains surrounded by a matrix made of quartz sand along the contact of the wedge and the belt (Figure S1 in Supporting Information S1).Rice displays velocity weakening behavior (i.e., (a-b) = −0.026)while sand is nearly velocity neutral (i.e., (a-b) ∼ 10 −5 ), as verified by ring shear test measurements (Mastella, Corbi, Funiciello, Rosenau, et al., 2021).We use data from five models with different asperity lengths (Lasp) in the 10-80 cm range.All models have a downdip asperity width (Wasp) of 30 cm and an updip limit at 10 cm distance from the trench.The downdip extent and depth range of the asperity are set to scale with the global average (Heuret et al., 2011).
Data (available open access at Mastella, Corbi, Funiciello, & Rosenau, 2021) consist of time series of model surface velocity fields produced using the Particle Image Velocimetry (PIV) technique (Sveen, 2004).The stick-slip behavior of Foamquake is captured by PIV as alternating phases of "slow" surface velocities toward the back end of the wedge (corresponding to stick phases, the analog of interseismic periods showing landward motion) and velocity peaks directed toward the wedge tip (corresponding to slip phases, the analog of coseismic periods showing trenchward motion; Figure 1b).Models are monitored at 50 frames per second for 2 min.Those settings allow the detection of tens of seismic cycles and determining rupture durations as slip episodes are captured by multiple frames.From the velocity field, we (a) sample a line-time time series to visualize the along-strike variability of model behavior (Figure 1d) and (b) calculate the following relevant parameters at the center of the section (i.e., a reference measurement point located at the center of the surface projection of the asperity; Figure S1 in Supporting Information S1): the recurrence time, the rupture duration and the fault slip.Recurrence time R t is defined as the time between consecutive velocity peaks, rupture duration R dur is inferred as the time span during which velocity exceeds a threshold of 0.1 cm/s (Figure 1e), and the fault slip D as the total displacement (which is computed from the velocity multiplied by time between consecutive frames) over each event (Figure 1c).We assume that parameters measured at the model surface mirror their corresponding equivalents along the analog megathrust.This assumption implies a minor (i.e., ∼7%; Mastella et al., 2022) underestimation of D and negligible impact on R t and R dur .

Numerical Modeling: General Information and Model Setup
We use multi-cycle numerical simulations to explore how rate and state parameters affect the seismic cycle characteristics (i.e., R t , R dur , and D; measured at the center of the asperity as for laboratory experiments) to identify the combination of rate and state parameters (i.e., a, b, and D c ) that allows reproducing the observed Foamquake behavior.Since Foamquake displays relatively low (∼0.1-1.0 cm/s) coseismic velocities, we use the quasi-dynamic approximation.This approach utilizes the radiation damping term to approximate the effect of inertia and represent a simplified and computationally efficient alternative to fully dynamic simulations (Thomas et al., 2014).In particular, we use QDYN (Luo et al., 2017), an open-source software that has been used to model a variety of features of the seismic cycle (e.g., Luo & Ampuero, 2017).
Our numerical models are designed to mimic the laboratory scale Foamquake experiment except for a dimension reduction which in the numerical models is set as 2D, that is a linear fault embedded in a 2D elastic medium loaded perpendicular to the cross-section, while Foamquake is 3D.We select thrust fault type with state variable evolving according to the aging law (Equation 2).Numerical models are designed to represent the same section sampled in Foamquake (Figure 1a).We opt for 2D models (e.g., Kaneko et al., 2010) for two reasons: (a) being computationally less demanding than 3D modes, they are suitable for extensive parameter exploration and (b) we focus exclusively on model sensitivity to changes applied along the strike of the fault.Recently, Li et al. (2022) studied the advantages and limitations associated with dimensionality reduction in numerical models of the seismic cycle, reporting a good match between 2D and 3D models for high aspect ratio asperities (Lasp/Wasp ∼ 5) when the investigated parameter is the recurrence time, while lower aspect ratio asperities (Lasp/Wasp ∼ 1) provide a better fit with coseismic characteristics.To disentangle potential bias introduced by 2D models, we run complementary 3D benchmark models to check for consistency of the seismic cycle characteristics retrieved with 2D modeling (Figure S2 in Supporting Information S1).Numerical models include a central velocity weakening patch (i.e., (a-b) vw < 0; subscript vw indicates velocity weakening) that is confined laterally by nearly velocity neutral segments (i.e., (a-b) vs = 10 −5 unless otherwise specified; subscript vs indicates velocity strengthening).The fault is loaded from the sides by steady motion V pl = 0.01 cm/s, mimicking the subducting plate velocity of Foamquake.Finally, numerical models share with Foamquake the same characteristics: total fault length (145 cm), asperity length Lasp (varying between 10 and 80 cm), shear modulus G (13.6 kPa), shear wave speed V S (30 m/s) and normal load profile characterized by 50 Pa above the asperity and 13 Pa at its sides (Figure 2a).

Estimating the Frictional Parameters of the Velocity Weakening Patch: The Reference Model
We selected a Foamquake experiment with Lasp = 20 cm as the reference model.This configuration represents the optimum in terms of scaling similarity with respect to nature.In particular, the reference model produces analog earthquakes with scaled to nature magnitudes Mw in the 8.5-9.4 range and a series of proportionalities between rupture parameters (e.g., moment-area, moment-slip) matching closely with the corresponding proportionalities observed in large megathrust earthquakes (Mastella et al., 2022).
This reference model is characterized by velocity peaks with maximum amplitude ∼0.45 cm/s (mean ∼0.25 cm/s; Figure 1b) and quasi-periodic recurrence behavior, as highlighted by the coefficient of variation C v = 0.36 (C v is defined as the standard deviation divided by the mean of the recurrence time); C v < 0.5 is generally associated to quasi-periodic behavior (e.g., Moernaut, 2020).To further characterize the model, we calculate the median Me and the standard deviation SD of the inferred parameters.R t shows a multi-modal distribution with peaks at ∼0.8 s, ∼1.5 s and ∼2.2 s, while R dur = 0.03 ∓ 0.01 s (Me ∓ SD) and D = 6.8 × 10 −5 ∓ 2.8 × 10 −5 m (Me ∓ SD) are normally distributed (Figure S3 in Supporting Information S1).Notice that mean values would be equivalently representative of R dur and D but the median seems more appropriate for R t , where it coincides with the group with R t ∼ 1.5 s which is the most representative of the asperity seismic cycle (see Section 5.2).
Compared to Foamquake, numerical simulations appear more regular (Figures 2b-2d).In particular, they have periodic behavior (C v < 0.01) and sequences of events equal to each other for duration, slip and lateral extent (Figures 2c and 2d).
To constrain the frictional parameters of the velocity weakening patch of Foamquake we first followed the same strategy adopted in previous studies (e.g., Barbot et al., 2012;Thomas et al., 2017).The procedure consists of two steps.In step one, we select an initial (guess) model where D c is set taking into account that the critical nucleation size h* = (π/2) GbD c /(b-a) 2 σ (Rubin & Ampuero, 2005) must be smaller than the asperity, a or b are set to a given value (e.g., a = 0.01) and the parameter (a-b) is constrained knowing the recurrence time T of a "characteristic" event, that is, one that saturates the velocity weakening region according to (Barbot et al., 2012): where W is the asperity size, V refers to velocity and subscripts co, int and r refer to coseismic, interseismic and long-term, respectively.After checking model behavior, in step two, (a-b) is fine tuned to reach the desired magnitude and recurrence time.We proceeded by setting a = 0.01 and computing T for various couples of V r and b (Figure 3a).b is allowed to vary in the 0.01 to 0.03 range.Since V int is difficult to constrain in our models (due to the small displacement between consecutive frames), we let it vary from 10 −4 m/s (i.e., the conveyor belt's velocity) to 10 −15 m/s, spanning the range from fully-to virtually zero creeping.We set the other parameters V co = 0.2 × 10 −2 m/s, V r = 10 −4 m/s, σ = 50 Pa, W = 0.2 m as in the laboratory.We identify the region of the b -V int space where the computed T is within the observed range of Foamquake.The corresponding values of b are in the 0.015-0.025range.Then, we compute h* for the same range of values of b, and D c comprised between 10 −4 m and 10 −7 m (Figure 3b).We observe that, depending on the chosen value of b, D c should be smaller than 5 × 10 −5 m for satisfying the h* < W condition and obtaining instabilities that nucleate and propagate in velocity weakening regions (e.g., Barbot, 2019).
Using the constraints above for b and D c , we run a first set of numerical simulations by selecting sparse values in the phase plot of Figure 3b and for each of them we compute R t , R dur , and D (Figures 3c-3e).This set of simulations is characterized by R t which increases from ∼0.4 to ∼1.5 s with increasing D c .Therefore, our initial choice of rate and state parameters allow us to properly mimic the recurrence time of Foamquake.However, simulations matching the recurrence time display longer slip duration and larger slip as compared to Foamquake.The numerical simulation with D c = 10 −5 m and b = 0.025 exemplifies this issue, being characterized by a recurrence time perfectly matching Foamquake's median and by ∼73% longer duration and ∼64% larger slip (filled diamond in Figures 3c-3e).Similarly, simulations with b = 0.02 and D c < ×10 −6 m display a slip that matches that of Foamquake, but shorter durations and recurrence times.None of the simulations from this initial set satisfy the condition of matching simultaneously R t , R dur , and D as observed in the laboratory: they reproduce only individual aspects of the seismic cycle of Foamquake.We implement the Metropolis-Hastings rule for accepting or rejecting a simulation (e.g., Passarelli et al., 2010), that is the probability of accepting a solution with larger error than the previous one decreases progressively during the search.As a consequence, initial iterations scan a broad region of the parameters space which likely includes-or can be close to the global minimum that is, the best combination of the parameters a, b and D c .We run 2,000 simulations and observe that our sampler converges to a minimum error after ∼1,900 iterations (Figure 4a).
The best numerical simulation identified with this procedure provides R t , R dur , and D that are closely matching the corresponding median of Foamquake.From the whole set of simulations we selected only those that simultaneously satisfy the condition of being in the Me ∓ SD/2 interval of Foamquake model parameters, hereafter named successful simulations.For ease of visualization, successful simulations are compressed in the D-time plane and compared against Foamquake's median (Figure 4b).Distributions of rate and state parameters of successful simulations and corresponding R t , R dur , and D are reported in Figure 4, panels c-h.Based on these distributions we infer (a-b) vw = −0.012∓ 4.6 × 10 −4 (Me ∓ SD/2).We note that this value of (a-b) is ∼2 times smaller than the results inferred from a ring shear tester for the same material but with ∼4 times higher normal loads.Such a normal-load dependency is consistent with theoretical considerations and expectations from rock testing (e.g., Scholz, 1998) and has been verified for our analog material (rice) experimentally (Rosenau et al., 2009).

Estimating (a-b) of the Velocity Strengthening Lateral Segments
The rate and state friction parameters search performed for the velocity weakening patch relies on the assumption that D c and b are constant along the entire fault length and that neighboring segments are nearly velocity neutral (i.e., (a-b) vs = 10 −5 ).The latter assumption is consistent with ring shear tests, which show minimal changes of sand friction for shear rate that varies three orders of magnitude (Mastella, Corbi, Funiciello, Rosenau, et al., 2021;Rosenau et al., 2017).Again, the normal load used in the ring shear tester is generally larger than in Foamquake and a normalload dependency of (a-b) can be reasonably expected (e.g., Scholz, 1998).Additionally, fault geometries in the ring shear tester and in Foamquake are different.Therefore, we here aim at an independent assessment of (a-b) vs in the experiment.Consequently, we run an additional set of numerical simulations where we systematically varied the value of (a-b) of both the velocity weakening patch (a-b) vw and neighboring segments (a-b) vs with the dualistic objective of (a) testing how their mutual combination control the seismic cycle characteristics; and (b) sensing how our previous rate and state parameters estimation for the velocity weakening patch is affected by the assumption of velocity neutrality of lateral segments.We set b = 0.0242 and D c = 1.05 × 10 −5 m constant along the entire fault length as for the best numerical simulation identified by means of SA and explore values of (a-b) vs in a wide range, from ideally neutral (i.e., (a-b) vs = 10 −5 ) to very strengthening (i.e., (a-b) vs = 10 −1 ).For the asperity we vary (a-b) vw between −0.02 and −0.01.Such a wide range of (a-b) vw creates slower than threshold-to seismic slip velocities.
First, we analyze the amplitude of velocity peaks, being the parameter crucial for event detection (Figure 5a).The amplitude of velocity peaks increases both with increasing values of (a-b) vw and decreasing values of (a-b) vs .Velocity peaks become progressively smaller, up to overtaking the velocity threshold (i.e., no event is detected) for (a-b) vw < −0.0105 and (a-b) vs < 10 −4 .When (a-b) vs > 10 −4 the transition to smaller than threshold peak velocities occurs for (a-b) vw > −0.011, with the transition that shifts toward (a-b) vw > −0.017 for (a-b) vs = 0.1.We also observe that, for a given value of (a-b) vs , larger (a-b) vw values promote longer R t , short R dur , and large D (Figures 5b-5d).On the contrary, for a given value of (a-b) vw , larger (a-b) vs values reduce R t , R dur , and decrease D.
The vast majority of numerical simulations with velocity peaks above the threshold are characterized by R dur of a few tens of ms as for Foamquake.Next, we analyze whether a given (a-b) vw − (a-b) vs combination provides a model behavior that falls within the observed Foamquake's reference experiment Me ∓ SD/2 (Figure S4 in Supporting Information S1).This procedure is applied for each of the model parameters (i.e., R t , R dur , and D).
Simulations with R t comprised between Foamquake's Me ∓ SD/2 are primarily confined in the (a-b) vs < 10 −2 range and (a-b) vw spanning from −0.02 to −0.01.Simulations with D between the Me ∓ SD/2 of Foamquake have  Our results highlight that there is not a single (a-b) vw -(a-b) vs combination that allows reproducing all parameters of Foamquake's reference experiment in ensemble, but rather a group of combinations that require mutual tuning of both parameters.However, the low (a-b) vs values obtained from ring shear tester measurements suggest values of (a-b) vw ∼ −0.012 as preferred and physically realistic.

Estimating the Frictional Parameters of the Velocity Weakening Patch: The Effect of Varying Asperity Size
After investigating the friction parameters of the reference experiment we now analyze the impact of asperity size on the inversion procedure.We base our analysis on five experiments which systematically explore the along-strike length of the asperity Lasp (i.e.,10,20,40,60,80 cm).
In Foamquake experiments, the amplitude of velocity peaks increases with increasing Lasp (Figures 6a-6e).
In particular, we observe analog earthquakes with velocity peaks generally smaller than 0.5 cm/s for Lasp up to 20 cm and ∼ two times faster velocity peaks for Lasp ≥ 40 cm.Similarly, R t , R dur , and D increase with Lasp.Also for those parameters, we notice a sharp transition between different model behaviors for experiments with Lasp ≤ 20 cm and the others.Compared to experiments with Lasp ≤ 20 cm, experiments with Lasp ≥ 40 cm display R t , R dur , and D larger by factors ∼1.7, ∼1.7, and ∼2.5, respectively (Figures 6f-6h).
Using the same SA algorithm implemented for the reference Foamquake experiment, we retrieve rate and state friction parameters that allow mimicking the seismic cycle characteristics of each Foamquake experiment (Figures 7a-7e).As for the reference model, we assume (a-b) vs = 10 −5 .
This procedure allows understanding how rate and state parameters should be tuned to properly reproduce the effect of asperity size on seismic cycles.We obtain a good fit for all experiments, indicating that the identified rate and state parameters allow numerical simulations to experience the same seismic behavior as in the laboratory.The SA generally achieves a good match for the observed versus simulated R t , D, and R dur .The parameter a remains about constant and comprised between 0.012 and 0.015 for models with different Lasp (Figure 7f).With respect to a, b decreases from ∼0.04 for the smallest Lasp to ∼0.02 for models with the largest Lasp (Figure 7g).This dependence of b to Lasp causes (a-b) to increase with Lasp, with values of (a-b) increasing from −0.026 to −0.007 (Figure 7h).Estimates of D c are between ∼7 × 10 −6 and ∼1.2 × 10 −5 m, without a clear trend with respect to Lasp (Figure 7i).Our results suggest that modeling the seismic cycle of a fault of a given length with asperities of different sizes implies tuning (a-b), in a way that larger asperities are associated with lower (a-b) values (i.e., less velocity weakening), while keeping D c constant.Our results suggest also that both a and b should be tuned to obtain the desired seismic cycle characteristics, instead of keeping one of the two parameters fixed at a given value.

Discussion
Analog models and numerical simulations allowed us to explore in a complementary way how the seismic behavior and empirical parameters of the rate-state law are tuned by different physical properties of a given fault.Analog models allowed investigating the effect of asperity length while numerical simulations have been used to invert rate and state parameters of both the asperity and neighboring segments from the experimental time series.
Here we discuss how our results frame into a wider context focusing on model dependence on asperity length and similarities and differences within the two complementary approaches.

Rate and State Parameters Variability With Asperity Length
Our results show that (a-b) vw apparently increases with asperity length while D c remains relatively constant.This scale dependency comes as a surprise since a and b are believed originally to be (scale-independent) material properties (Scholz, 1998).Our initial expectation was to retrieve values of (a-b) vw similar to that obtained with ring shear tester (i.e., (a-b) = −0.026)and possibly constant for experiments with different asperity lengths but with the same normal stress and loading rate.On the contrary, we observed about constant (i.e., variations <2‰) values of a and decreasing values of b for increasing asperity lengths, resulting in (a-b) that increases of a factor 3.8 over the investigated asperity length range.Interestingly, the retrieved (a-b) value for the smallest Lasp matches particularly well the ring shear tester measurement consistent with their comparable size (i.e., the cross-section of the measurement cell of the ring shear tester is 5 cm wide), but models with larger Lasp require larger (a-b) values to fit the observed laboratory behavior.To explain this behavior we run complementary numerical simulations where we kept a, b and D c constant and equal to those retrieved for the reference Foamquake model and varied exclusively the asperity length.We observe that recurrence time, slip duration and slip increase with asperity length for both 2D and 3D numerical simulations, with slip showing the greatest dependence on asperity size in the case of 2D simulations (Figure S2 in Supporting Information S1).The observed increase of (a-b) vw with asperity length therefore is required to contrast this behavior: larger (a-b) vw values promote relatively shorter recurrence times, smaller slip and shorter slip durations.Those observations suggest that a and b might depend on the size and dimensionality of the seismogenic patch.D c is a key friction parameter that controls various earthquake characteristics and processes (Rabinowicz, 1951;Scholz, 1988).D c assumes different meanings for the different constitutive laws in which it is introduced as the critical slip distance or slip weakening distance (Marone & Kilgore, 1993;Ohnaka, 2000).A scaling relationship between D c and Mw has been proposed based on theory and empirical observations (Cocco et al., 2023;Gallovič & Valentová, 2023;Ohnaka, 2000).Typical laboratory measurements of D c from cm-sized samples are in the 10 −6 to 10 −5 m range, while field-based estimates and numerical simulations suggest larger values of D c in the 10 −3 to 1 m range (e.g., Gallovič & Valentová, 2023;Johnson et al., 2006;Marone & Kilgore, 1993;Tinti et al., 2021).The physical significance of D c is debated (Scholz, 1998).In fact, D c has been suggested to represent the asperity contact properties (Rabinowicz, 1951) or to reflect different geometric fault properties such as roughness (e.g., Ohnaka, 2003;Scholz, 1988), gauge particle size (Dieterich, 1981) and thickness (Marone & Cox, 1994), thickness of the active shear zone (Marone & Kilgore, 1993) and size of the breakdown zone (Gallovič & Valentová, 2023;Ohnaka & Yamashita, 1989;Rice, 1980).Our experiments and numerical simulations do not seek discriminating between different natures of D c .However, the retrieved value of D c ∼ 10 −5 m matches closely that of geometric surface irregularities (wavelength of roughness) of rice implemented for our asperity (Rosenau et al., 2009).Upscaling D c from the analog model to natural prototypes results in D c in the order of decimeters (Rosenau et al., 2009) similar to what has been reported for nature (e.g., Johnson et al., 2006;Marone & Kilgore, 1993;Ohnaka, 2000).The fact that all Foamquake experiments are run with the same materials, together with the minor laboratory earthquake magnitude variations (average Mw variation ∼0.4 over the investigated Lasp range) observed for experiments with different asperity sizes, explain why the retrieved D c is about the same for different experiments.These observations support a modeling strategy considering that a few times larger asperities (e.g., <4x) should be modeled with smaller (a-b) rather than larger D c values.

Non-Dimensional Parameters, Rupture Behavior and Influence of the Neighboring Segments
In the frame of earthquake physics and numerical modeling of the seismic cycle, some non-dimensional parameters have been proposed to control different aspects of the seismic cycle (e.g., Barbot, 2019).In those studies, simple configurations with a central asperity (e.g., Barbot, 2019) or alternating velocity weakening and velocity strengthening patches are implemented (e.g., Luo & Ampuero, 2017).
Two of these parameters describing asperity properties are the Dieterich-Ruina-Rice number (Ru) and the velocity neutral proximity number (Rb).Ru is defined as (Barbot, 2019): where the first term relates dynamic weakening to rigidity (i.e., tendency for instability) and the second term represents the size of the asperity to the critical slip distance.Ru can be interpreted as the "criticalness" of the asperity.Small Ru values (i.e., <2) are associated with slow slip behavior.Increasing Ru, the system moves toward faster, "seismic" slip velocities and complexity in the space-time rupture history.
Rb is defined as (Barbot, 2019): Ru and Rb values computed using rate and state parameters that characterize our experiments range between ∼1-∼4 and ∼0.3-∼0.6,respectively, depending on Lasp (Figure 7 panels l and m).According to previous parametric studies based on numerical simulations, those values of Ru and Rb generate a slow slip and characteristic earthquake behavior (Barbot, 2019;Cebry et al., 2022).The slow (with respect to real earthquakes) slip velocities observed in our experiments (about mm/s) and the small Ru retrieved with our analysis fit well with previous numerical computations, where peak slip velocities slower than 10 −2 m/s, are expected for small Ru, of few units (Barbot, 2019).Additionally, for the retrieved ranges of Ru and Rb, numerical simulations predict periodic behavior.Our experiments display quasi-periodic recurrence times (C v in the 0.25-0.36range depending on Lasp) and therefore fit nicely in this framework.
Any deviations from the perfectly periodic behavior in our experiments can be explained by the inherent variability of material properties or by interactions of the asperity with laterally unconfined model edges, that reduce stiffness and enhance instability (Cebry et al., 2022;Rosakis et al., 2007).The latter, preferred hypothesis explains (a) why the observed clustering around two principal recurrence times in our reference experiment disappears with increasing Lasp (i.e., the asperity becomes the principal source for recurrence time rather than small events that propagate from model edges; Figure S5 in Supporting Information S1) and (b) the complexity that characterizes the space-time rupture history.In fact, numerical simulations indicate that small Ru produces sequences of symmetrical ruptures nucleating at the center of the asperity (Figures 8a and 8b; Figure S6 in Supporting Information S1).On the contrary, in our experiments we distinguish a variety of rupture types including ruptures that start from model edges or from the asperity and can both rupture only a fraction-or the whole lateral extent of the model (Figures 8c-8f).The failure of numerical simulations to reproduce the space-time rupture variability as observed in the laboratory is not surprising and is most likely due to the single reference point analysis framing.Future efforts to include this feature in numerical simulations could be achieved by increasing the number of model parameters (e.g., a, b, and D c variable along strike) as well as the number of constraints describing model behavior (e.g., recurrence time measured at different locations and frequency magnitude distribution).
The influence of lateral segmentation of frictional properties on seismic cycle behavior suggests the introduction of  = ( − )vw∕( − )vs; (7) another non-dimensional parameter that describes the relative strength between different segments of our models (Luo & Ampuero, 2017).Our rate and state parameters search was based on the assumption of velocity neutrality of lateral segments.This assumption creates high strength contrast between the asperity and lateral segments with α in the order of thousands (Figure 5f).Complementary simulations with variable lateral segment strengths suggest that simulations with lower (i.e., up to a hundred) α create the same behavior as the reference Foamquake model.The fact that a wide range of α works equally well for given asperity properties (i.e., 0.8 < Ru < 1.1; Figure 5e) highlights the non-uniqueness of successful simulation parameters and a threshold for α above which no impact on model behavior is expected.
We conclude that, even in the case of a simple fault, modelers face some degrees of non-uniqueness when setting rate and state parameters.Here, we highlighted the tuning between asperity-and neighboring segments' strength ratio but other parameters such as normal load and D c variations along strike deserve interest for future studies.
Additional efforts for better constraining model parameters include SA approaches to retrieve a, b and D c on both the asperity and neighboring segments with multiple target points along the fault length.

Ranking Source of Uncertainties When Assessing (a-b) of the Asperity
Our simulations shed light on uncertainties involved when trying to infer rate and state variables from observations.Assuming neighboring segment properties and asperity size are both known, (a-b) vw variations in the order per mille are required to explain Foamquake's variability.This is the range of (a-b) that we retrieved for the reference Foamquake model (Figure 7h).Roughly the same range of variability in (a-b) vw is expected if neighboring segment properties are poorly constrained but asperity size is known (Figure 5e).Interestingly, (a-b) vw varies most significantly (up to 10 −2 ) and systematically depending on the size of the asperity (Figure 7h).This observation suggests that having access to asperity size (e.g., with geodetic locking data or seismological constraints) is pivotal for inferring (a-b) vw and, in turn, producing realistic rupture scenarios for seismic hazard assessment.Importantly, defining the size of asperities in nature represents a challenge as they might be much smaller than previous estimates based on rupture extents (Burgmann et al., 2005;Zhao et al., 2022) or locked area (shadowing effect, e.g., Herman et al., 2018).The general offshore location of asperities or locked patches poses an additional uncertainty when based on geodetic observations, although the along-strike length might be one of the parameters rather well constrained by land-based geodesy (e.g., Kosari et al., 2020).
An additional source of uncertainty derives from the effective basal loading rate, which, in turn, affects the recurrence time.In our numerical models, we imposed the same loading rate (i.e., convergence velocity) at the base of the wedge as in Foamquake.However, the effective loading rate which is transferred into the foam wedge, could be smaller in the experiment due to partial coupling between the foam wedge and the analog slab.This effect is difficult to quantify in the laboratory but explains the missing successful simulations with R t > median for the

•
Rate and state friction parameters of a downscaled, laboratory subduction megathrust are constrained • Numerical models show a tradeoff between (a-b) of the asperity and (a-b) of neighboring segments • To fit the behavior of laboratory experiments with different asperity lengths it is necessary to vary (a-b) and D c

Figure 1 .
Figure 1.Foamquake's seismic behavior.(a) Oblique view of the experimental setup.(b) Timeseries of the trench orthogonal component of the velocity field sampled at the reference point.Analog earthquakes are represented by velocity peaks with progressive numbering.(c) Time evolution of cumulative slip measured at the reference point.(d) Line-time evolution of the trench orthogonal component of the velocity field sampled along the red line reported in panel (a).Analog earthquakes are represented by vertical black lines.Red rectangles shown in panels (b-d) highlight the fraction of the timeseries that is zoomed in panels (e-g), respectively.Green and red circles shown in panel e mark the intersections between the velocity threshold (blue line) and the time series, defining the onset and end of analog earthquakes, respectively.

Figure 2 .
Figure 2. Seismic behavior of a representative numerical model with the same panel organization as in Figure 1.(a) Schematic representation of the numerical 1D fault with parameters distribution.Timeseries of velocity measured at the reference point, cumulative slip and line-time ruptures distribution (panels b-d, respectively).The width of vertical lines in panels (d and g) is not meant to scale with rupture duration.Parameters as follow: a = 0.01, b = 0.02, D c = 5 × 10 −6 m.

Figure 3 .
Figure 3. Setting rate and state parameters for numerical simulations: recurrence time and nucleation size computed for various parameters combinations (panels a and b, respectively).Red curves in panel (a) constrain the range of recurrence time observed in Foamquake.The red curve in panel b marks the transition between nucleation lengths larger and smaller than the size of the asperity used in the reference Foamquake experiment.Open symbols in panel (b) highlight couples of b and D c used in a set of numerical simulations whose model behavior (i.e., recurrence time, slip duration and slip) is shown in panels (c-e).Red horizontal lines and background shading represent Foamquakes' Me ∓ SD/2. a = 0.01 is assumed constant along the fault length.Filled diamonds highlight a simulation discussed in Section 4.1.

Figure 4 .
Figure 4. Simulated annealing method and route toward the preferred rate and state parameters combination.(a) Evolution of the error E during the parameter search.The black solid line in panel a highlights the path of accepted simulations.Successful simulations that fall within Foamquake's Me ∓ SD/2 are highlighted by the black dots and reported in the following panels.Panels (b) shows displacement-time curve for successful numerical simulations (black) compared against Foamquake median (red) distribution of model parameters (panel c-e) and distribution of rate and state parameters (panels f-h) for successful simulations.Red circles and horizontal lines in panels (c-e) represent Foamquake's Me ∓ SD/2.White stars represent the simulation with lowest error.A line-time representation of this model is reported in Figures 8a and 8b.

Figure 5 .
Figure 5. Simulations sensitivity to (a-b) vs and (a-b) vw .Circles represent simulations with b = 0.0239, D c = 9.46 × 10 −6 m and given (a-b) vs and (a-b) vw .Circles are color-coded by peak slip rate, recurrence time, rupture duration and slip (panels a-d respectively).Open red circles highlight simulations with recurrence time, slip duration and slip comprised between the corresponding Foamquake's Me ∓ SD/2.The field of the (a-b) vs -(a-b) vw plane where no circle is represented corresponds to simulations with velocity peaks slower than the threshold, that is, no event detected.Panel (e) compresses information from panels (b-d) in a single plot highlighting with open circles simulations that successfully reproduce Foamquake's recurrence time, slip duration and slip, in ensemble.Panel (f) shows the same successful simulations of panel (e), but represented in the Ru-α plane (see Equations 5 and 7).

Figure 7 .
Figure 7. Rate and state parameters for Foamquake experiments with different Lasp.Panels (a-e) show Foamquake's median and successful numerical simulations reported in the slip-time plane (red and black lines, respectively).Panels (f-n) show a summary of the retrieved rate and state parameters as a function of Lasp visualized as violin plots (the width of each violin is proportional to the distribution of each parameter for successful simulations).White circles represent the median.

Figure 6 .
Figure 6.Foamquake's seismic behavior for models with different asperity lengths Lasp.Timeseries of the trench orthogonal component of velocity measured at the reference point for five experiments with Lasp specified in each panel (panels a-e).The red box of panel b highlights the reference experiment with Lasp = 20 cm.Recurrence time, slip duration and slip observed in each experiment are summarized with boxplot (red circles represent the median, thin and thick lines constrain the 90% and 50% of the distribution) in panels f-h, respectively.

Figure 8 .
Figure 8. Line-time representation of numerical (panel a and zoom in panel b) and analog (panels c-f) earthquakes.The background shading represents slip rate.Red contours represent the velocity threshold 10 −3 m/s used to detect slip phases.Horizontal blue lines constrain the central asperity.
) and reflects the relative importance of transient weakening (b) over the direct effect (a).In case of velocity neutrality, a = b and Rb becomes 0. If b dominates (strong velocity weakening) Rb approximates 1.For velocity strengthening (a > b) Rb becomes negative.