MCMC inversion of the transient and steady-state creep flow law parameters of dunite under dry and wet conditions

The rheology of the upper mantle impacts a variety of geodynamic processes, including postseismic deformation following great earthquakes and post-glacial rebound. The deformation of upper mantle rocks is controlled by the rheology of olivine, the most abundant upper mantle mineral. The mechanical properties of olivine at steady state are well constrained. However, the physical mechanism underlying transient creep, an evolutionary, hardening phase converging to steady state asymptotically, is still poorly understood. Here, we constrain a constitutive framework that captures transient creep and steady state creep consistently using the mechanical data from laboratory experiments on natural dunites containing at least 94% olivine under both hydrous and anhydrous conditions. The constitutive framework represents a Burgers assembly with a thermally activated nonlinear stress-versus-strain-rate relationship for the dashpots. Work hardening is obtained by the evolution of a state variable that represents internal stress. We determine the flow law parameters for dunites using a Markov chain Monte Carlo method. We find the activation energy 430±20\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$430\pm 20$$\end{document} and 250±10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$250\pm 10$$\end{document} kJ/mol for dry and wet conditions, respectively, and the stress exponent 2.0±0.1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.0\pm 0.1$$\end{document} for both the dry and wet cases for transient creep, consistently lower than those of steady-state creep, suggesting a separate physical mechanism. For wet dunites in the grain-boundary sliding regime, the grain-size dependence is similar for transient creep and steady-state creep. The lower activation energy of transient creep could be due to a higher jog density of the corresponding soft-slip system. More experimental data are required to estimate the activation volume and water content exponent of transient creep. The constitutive relation used and its associated flow law parameters provide useful constraints for geodynamics applications.


Introduction
Our current understanding of the rheology of the upper mantle comes mainly from rock deformation studies in the steady-state creep regime, where olivine and olivinebearing rocks are deformed in the relevant temperature and pressure conditions (e.g., Carter and Avé Lallemant 1970;Paterson 1981, 1984;Karato et al. 1986;Karato and Wu 1993;Mei and Kohlstedt 2000a). However, the rapid deformation following large earthquakes often cannot be explained solely by the steady state response of the viscoelastic substrate (e.g., Pollitz 2003Pollitz , 2005Freed et al. 2010;Qiu et al. 2018;Tang et al. 2019;Weiss et al. 2019;Tang et al. 2020) and requires the activation of transient creep, a hardening phase of the rheological behavior associated with higher strain rates than the steady-state response at similar loading conditions. Geodynamic studies rely on the constitutive formulations and parameters derived from laboratory studies to estimate and refine a thermo-mechanical model of the lithosphere-asthenosphere system and explain various observations (e.g., Freed and Bürgmann 2004;Rousset et al. 2012;Rollins et al. 2015;Masuti et al. 2016;Muto et al. 2019;Agata et al. 2019;Barbot 2020). More insights into the constitutive behavior of transient creep from laboratory data is needed to better constrain the physics of rock deformation.
Olivine at high temperature deforms mainly by dislocation creep, diffusion creep, and dislocation-accommodated grain-boundary sliding (DisGBS). Dislocation creep is characterized by a nonlinear relationship between stress and strain rate, and does not depend on grain size. Diffusion creep is characterized by a linear relation between stress and strain rate and is highly sensitive to grain size. DisGBS has the characteristics of both dislocation creep and diffusion creep with a nonlinear relation between stress and strain rate and a dependence on grain size (e.g., Hirth and Kohlstedt 1995;Hansen et al. 2011;Wang et al. 2010;Ohuchi et al. 2015). The dominant mechanism at steady state can be determined from the deformation map based on stress and grain size (Ashby 1972). Dislocation creep is usually considered the dominant mechanism under asthenospheric conditions based on deformation experiments (e.g., Karato and Wu 1993), seismic anisotropy (e.g., Montagner and Tanimoto 1990;Silver and Holt 2002;Long and Silver 2009), and mantle xenolith microstructure (e.g., Mercier 1985;Park and Jung 2015). However, the deformation of olivine under both laboratory and asthenospheric conditions is close to the boundaries between dislocation creep, diffusion creep, and DisGBS, and the physical conditions may vary either due to grain-size evolution at geological time scales (Karato 1989;Rozel et al. 2011) or due to varying loading conditions during the seismic cycle (Lambert and Barbot 2016;Barbot 2018).
Grain boundary sliding accommodated by the motion of dislocation was first observed in fine-grained olivine mineral by Hirth and Kohlstedt (1995), who proposed grain boundary sliding (GBS) as an additional degree of freedom to relax the von Mises (1928) criterion that requires the activation of five independent slip systems for the homogeneous deformation by strain partitioning of a polycrystalline mineral. Additional experiments have been conducted to determine the flow law parameters of DisGBS under different conditions. Hansen et al. (2011) conducted experiments on San Carlos olivine aggregates with different grain sizes under dry conditions and argued that the seismic anisotropy observed in the asthenosphere is not purely due to dislocation creep and that DisGBS could be one of the dominant mechanisms of deformation (see also Wang et al. 2010). In addition, high-pressure ( ∼ 6.7 GPa) experiments on olivine aggregates using a deformation-DIA apparatus, a type of multi-anvil deformation apparatus, under dry and wet conditions have been conducted, and the flow law parameters for the DisGBS mechanism at steady state were determined (Ohuchi et al. 2015).
The flow law parameters of upper-mantle rocks for dislocation creep, diffusion creep, and DisGBS under dry and wet conditions in the steady state regime are well known (e.g., Karato et al. 1986;Mei and Kohlstedt 2000a, b;Karato and Jung 2003;Hirth and Kohlstedt 2003;Hansen et al. 2011;Wang et al. 2010;Ohuchi et al. 2015). In contrast, the transient creep flow law parameters for all these mechanisms under dry or wet conditions are unknown. A few studies analyzed data in the transient creep regime (Post 1977;Chopra 1997;Hanson and Spetzler 1994;Thieme et al. 2018), but the transient creep flow law parameters were not constrained uniquely. Post (1977) and Chopra (1997) performed experiments on dunites under constant stress conditions using a solid-medium and a gas-medium apparatus, respectively, and concluded that transient creep and steady-state creep are controlled by similar physical mechanisms. Hanson and Spetzler (1994) conducted further experiments on olivine single crystals and argued that transient creep could be anisotropic, as they observed strain softening instead of strain hardening in the case of the [110] c orientation. Hanson and Spetzler (1994) and Chopra (1997) found that the dislocation density evolves during the transition from transient creep to steady-state creep, in agreement with theoretical models that assume an evolving internal microstructure (e.g., Karato 2008). Recently, Thieme et al. (2018) conducted a detailed study of the transient creep of olivine aggregates under anhydrous conditions at high temperature (1273 and 1473 K) and 300 MPa of confining pressure using a gas-medium apparatus. Dislocation microstructures were analyzed at different cumulative strain in the transient regime, but evidence for a systematic evolution of the microstructure was inconclusive.
Several mechanisms have been suggested for the transient creep of olivine aggregates (e.g., Karato 1998;Sherburn et al. 2011;Masuti et al. 2016;Holtzman et al. 2018;Hansen et al. 2020;Karato 2021) and for ice (e.g., Ashby and Duval 1985). Based on the similarity of hardening behavior between single crystal and olivine aggregate at low temperature, an intragranular mechanism is thought to be dominant in the transient regime, where the microstructural changes during transient creep are related to the evolution of backstress due to long-range elastic interaction among dislocations (Hansen et al. 2019). In contrast, based on the plastic anisotropy observed in olivine and ice, an intergranular mechanism has been proposed to be dominant in the transient regime, where the strain incompatibility and hence the internal stress arises due to the difference in strengths of different slip systems (e.g., Karato 1998;Masuti et al. 2019;Ashby and Duval 1985;Karato 2021).
The goal of this study is to establish the validity and determine the parameters of a constitutive framework for the transient creep of dunites. The microstructure evolution during transient creep could be modeled with an empirical flow law that involves the evolution of a state parameter. Sherburn et al. (2011) employ a formulation derived from viscoplasticity theory (Nemat-Nasser and Hori 1999;Nemat-Nasser 2004), using isotropic and anisotropic kinematic hardening. Since transient creep corresponds to an initial rapid displacement followed by a slower response, work hardening is appropriate. However, the formulation of Sherburn et al. (2011) contains a large number of free parameters that cannot be fully constrained with available laboratory data. The laboratory data on transient creep can be explained by a simple formulation that captures the evolution of effective viscosity throughout the experiment, showing systematic hardening towards an asymptotic steady-state value, with a ratio of transient to steady-state viscosity of 0.17 to 0.67 (Chopra 1997). These constraints can readily be assimilated in linear spring-dashpot assemblies for geophysical modeling (e.g., Sabadini et al. 1985;Yuen et al. 1986) and in studies of postseismic deformation (e.g., Pollitz 2005;Hoechner et al. 2011;Wang et al. 2012). However, a linear constitutive law is not compatible with our corporate understanding of the rheology of the upper mantle at steady state, which involves a non-linear stressversus-strain-rate relationship with explicit dependencies on temperature, water content, and grain size. Overall, there is no consensus in the mineral physics community on the exact formulation of the transient creep of olivine aggregates.
Here, we consider a formulation with few degrees of freedom that extends the rheological laws that are well established for steady state (e.g., Hirth and Kohlstedt 2003;Karato and Jung 2003). Taking geodynamic applications into account, we consider the constitutive relationship proposed by Masuti et al. (2016), whereby the mechanical behavior of the representative volume element is captured by a spring-dashpot assembly (Fig. 1). The model is compatible with the inter-granular transient dislocation creep theory proposed by Karato (1998), whereby the small strain during the initial transient stage  of deformation is accommodated by the soft-slip system and, as the deformation continues, the hard-slip system takes over to homogenize the stress (see also Ashby and Duval 1985). The nonlinear Burgers formulation has proven useful to model postseismic deformation following giant earthquakes, including the 2012 M w 8.6 Indian Ocean (Masuti et al. 2016), the 2011 M w 9.1 Tohoku-Oki (Muto et al. 2019;Agata et al. 2019), and the 2010 M w 7.2 El Mayor-Cucapah (Tang et al. 2020) earthquakes, especially when multiple cycles of earthquakes, slow-slip events, and other stress perturbations are involved (Shi et al. 2020;Barbot 2020).
In this work, we present an analysis of laboratory data for transient and steady state creep to establish the effect of grain size, stress, and temperature for dry and wet conditions. The manuscript is organized as follows. In the next section, we describe a curated data set that allows us to parameterize the effects of grain size, stress, and temperature. Then, we describe the nonlinear constitutive formulation and a numerical method to simulate forward models of deformation at constant strain rate. We present a Monte Carlo method to estimate the constitutive parameters from nonlinear inversion of laboratory data. Next, we determine the flow law parameters for the steady state and transient creep of dunites in dry and wet conditions. We verify the applicability of the constitutive framework on San Carlos olivine and some other dunite data. This is followed by a discussion of the physics underlying transient and steady state creep and a conclusion.

Experimental data and observations
The evolution of transient creep is most easily identified in experiments conducted at constant strain rate. However, only a few experiments conducted thus far in the community provide sufficient temporal resolution of stress evolution to constrain the underlying constitutive law. Available data sets for dunites using a gas-medium apparatus include the constant strain-rate experiments of Chopra and Paterson (1981) and Chopra and Paterson (1984) and the constant stress experiments of Chopra (1997); other data sets using a solid-medium apparatus include the work of Post (1977). The solid-medium apparatus, such as Griggs and D-DIA, has been valuable to determine the flow law parameters of steady-state creep for various rocks (e.g., Karato and Jung 2003;Hirth and Kohlstedt 2003;Bürgmann and Dresen 2008;Ohuchi et al. 2015). However, it cannot be used to constrain the transient creep flow law parameters without correction of the frictional stress between the deformation piston and the confining media (e.g., Kido et al. 2016;Holyoke III and Kronenberg 2010). Whether this correction is precise enough to study transient creep remains to be tested.
To mitigate potential bias from the Griggs apparatus, we do not include the data from Post (1977) in our analysis. Unfortunately, most of the data from Chopra (1997) do not include a clearly identifiable transient creep in the strain-versus-time curves (Additional file 1: Figure S1 ) and are also excluded from our inversion.
Previous work on transient creep often involved fitting individual strain-versus-time curves against a constitutive model prediction (e.g., Post 1977;Chopra 1997). This approach is limited because of the possible trade-offs between model parameters that can fit the data equally well. In addition, each strain-versus-time data set may require inconsistent parameter values. We, therefore, consider the data from Chopra and Paterson (1981) and Chopra and Paterson (1984) that cover a wide range of thermodynamic conditions and showcase a clearly identifiable phase of transient creep. They conducted constant strain rate experiments on Åheim dunites with a grain size of 900 µ m and on Anita Bay dunites with a grain size of 100 µ m under hydrous and anhydrous conditions using a gas-medium apparatus. The experimental conditions span a wide range of temperatures from 1273 K to 1673 K and strain rates from 10 −6 to 10 −4 s −1 (Tables 1 and 2), but were conducted at the same confining pressure of 300 MPa. We use their full stress-versus-strain curves that are corrected for apparatus distortion, except for runs number 2706 and 4392, which were conducted at 1273 K and did not reach steady state. Due to powerlaw breakdown at lower temperature (e.g., Katayama and Karato 2008), the deformation in these experiments could be dominated by the Peierls mechanism instead of dislocation creep. Some of the experimental data from Chopra and Paterson (1981) have been published as part of Chopra (1986), including runs number 4313, 4381, 4485, 4306, 4370, and 4371. Runs number 4709, 4699, and 4705 from Chopra (1986) are not included in our analysis as the deformation in these experiments indicated diffusion creep and not dislocation creep. Chopra and Paterson (1984) determined the flow law parameters for the steady-state creep of dunites. For the anhydrous case, they found the stress exponent to be 3.6 ± 0.1 and the activation energy to be 535 ± 15 kJ/ mol (Table 4). Under anhydrous condition, dunites do not show any grain size dependency. However, under hydrous condition the deformation of dunites is grainsize sensitive. Chopra and Paterson (1981) analyzed only the steady-state data under hydrous conditions and provided the flow law parameters for Åheim and Anita Bay dunites separately. They reported the stress exponent 3.35 ± 0.09 and 4.48 ± 0.16 , and the activation energy 444 ± 12 and 498 ± 19 kJ/mol for Anita Bay and Åheim, respectively. They also discussed the observation of grain size dependency using two different theoretical models, and concluded that either intra-or inter-granular mechanisms, or both, could be contributing to deformation. Here, we exploit these data further by introducing grain size dependency explicitly in the constitutive formulation.
In addition to the data discussed above, there are some data sets available for the San Carlos olivine aggregates (e.g., Karato et al. 1986;Thieme et al. 2018;Demouchy et al. 2013;Mei and Kohlstedt 2000a). The composition of San Carlos olivine is different from dunite and hence these data sets cannot be combined together to determine a unique set of flow law parameters. In addition, the experimental setup across most laboratories varies, requiring additional random stochastic variables that complicate the inversion (e.g., Korenaga and Karato 2008). As the experimental setup for San Carlos olivine aggregates from Karato et al. (1986) is similar to the one of Chopra and Paterson (1981) and Chopra and Paterson (1984), we use these data to check the applicability of the nonlinear Burgers model and to verify the similarities with the dunite flow law parameters. Although a majority of the data from Chopra (1997) cannot be used, we do test our model performance against those few data sets, where the transient and steady-state creep are clearly identified.

Flow law for transient creep
Several attempts have been made to formulate a rheological law for transient creep, in particular for metals (e.g., Andrade 1910;Voce 1948;Hart 1970;Poirier 1980;Chinh et al. 2004), but these formulations do not satisfy a number of fundamental principles, such as time frame invariance and reference frame indifference, that allow extrapolation for geodynamics applications. The same deformation should be predicted for any clock and any coordinate system. Recent developments in numerical modeling of geodynamic processes (e.g., Van Dinther et al. 2013;Herrendörfer et al. 2015Herrendörfer et al. , 2018Sobolev and Muldashev 2017) or seismic cycles (e.g., Barbot 2018;Shi et al. 2020;Barbot 2020) require a constitutive framework amenable to continuously evolving stress with multiple perturbations. This precludes using time explicitly in the constitutive formulation, unless as a derivative.
We consider a Burgers spring-dashpot assembly, which consists of a Maxwell element in series with a Kelvin element (Fig. 1), where the dashpots obey a nonlinear stress-versus-strain-rate relationship. The Kelvin and the Maxwell elements represent the soft-and hard-slip systems, respectively, capturing the transient and steady-state creep mechanisms. The spring in the Kelvin element represents the internal stress within the representative volume element. The spring in the Maxwell element represents the macroscopic stress on the whole aggregate. The effective stress on the soft-slip system may be different from the macroscopic stress due to internal stress. During the initial stages of deformation, the internal stress evolves continuously and the soft slip-system dominates in the transient regime.
Once the internal stress is approximately equal to the macroscopic stress, the deformation becomes dominated by the hard-slip system. In this idealized model, the Maxwell and Kelvin elements are in series, and the corresponding strain rates are additive. Hence, the thermodynamically irreversible part of the strain rate is given by where ε K is the strain rate in the Kelvin element, ε M is the strain rate in the Maxwell element, and ε i represents the total inelastic strain rate.
The strain rate in the Maxwell element is given by a nonlinear stress-versus-strain-rate relationship compatible with previous findings (e.g., Karato and Jung 2003;Hirth and Kohlstedt 2003): where A M is the pre-exponential factor, σ is the deviatoric stress, n M is the stress exponent, E M is the activation energy, R is the universal gas constant, T is the temperature, and d and p M are the grain size and its exponent. The formulation captures the behavior of various microphysical mechanisms of deformation at steady-state depending on the choice of parameters (e.g., Bürgmann and Dresen 2008). For example, diffusion creep corresponds to n M = 1 , dislocation creep corresponds to p M = 0 and n M > 1 , and DisGBS corresponds to p M > 0 and n M > 1.
In general, the strain rate in the Kelvin element is a function of the effective stress: where ε K is the cumulative strain in the Kelvin element and G K is a work-hardening coefficient. The effective stress drives plastic deformation in the Kelvin element and is modulated by G K ε K , which represents the internal stress within the representative volume element. The strain rate in the Kelvin element can then be written ε K = f K (q K ) , where f K is a zero-crossing at the origin and monotonically increasing function. The exact form of f K is unknown, with possibly linear, nonlinear, or exponential forms. For consistency with the functional form of the stress-versus-strain rate at steady state, we adopt the following power law, thermally activated formulation: where the form x ||x|| n−1 is to allow for possible reverse deformation in the Kelvin element in some conditions, for example, to simulate cyclical deformation (e.g., Hansen et al. 2019). When the effective stress is negative, expressions of the form x n would give rise to complex numbers when n is not an integer and are, therefore, inadequate. The parameters A K , n K , p K , and E K play the same roles as in equation (2), but for the Kelvin element. Masuti et al. (2016) describe a tensor equivalent of equation (4), but a scalar formulation is sufficient to model laboratory data. The formulation is appropriate to represent transient and steady-state creep consistently. For example, steady state implies ε K = 0 , such that the strain rate in the Maxwell element represents steady-state plastic flow. During transient creep, both the Kelvin and Maxwell components contribute to the total inelastic strain. Therefore, rapid deformation of the Kelvin element is associated with the phase of transient creep, even though the Maxwell component also contributes strain during that stage (Fig. 1).
The olivine flow law parameters and their uncertainties at steady state are well known for dislocation and diffusion creep (e.g., Bürgmann and Dresen 2008;Korenaga and Karato 2008). However, the flow law parameters and their uncertainties for transient creep are unknown. In this study, we focus on the experimental data conducted at the confining pressure of 300 MPa. The effective activation energy, therefore, represents an activation enthalpy with a small pressure effect. We restrict ourselves to the determination of the transient creep parameters A K , G K , n K , and E K for the anhydrous case, assuming that the steady-state flow law parameters are known from previous studies, and to A K , G K , n K , p K , E K , A M , n M , p M , and E M for the hydrous case.

Modeling approach
Unlike for the steady-state creep parameters, which can be determined using any of the least-squares methods that solve a linear inverse problem, determining the transient creep parameters requires the modeling of the time series of strain or stress under constant stress or constant strain rate conditions, respectively, and the resulting inverse problem is nonlinear. We model the stress-versus-strain evolution with the nonlinear Burgers model described above as follows.
We consider a spring-dashpot model (Fig. 1), where the total strain rate is the sum of the elastic and inelastic components (e.g., Andrews 1978; Barbot and Fialko 2010): where ε e and ε i are the elastic and inelastic strain rates, respectively. Substituting Hooke's law ε e =σ /G M , where G M is Young's modulus, we obtain Equations (2), (3), (4), and (6) form a set of coupled nonlinear ordinary differential equations associated with an initial value problem.
To simulate deformation data given constitutive parameters and physical conditions, we solve the governing equations numerically using a predictor-corrector method with a fixed step size of t = 0.5 s. We recast the governing equations described above in canonical form as with the state vector and where the function f (y, t) combines equations (2), (3), (4), and (6). In the predictor step, the initial guess is obtained using where the subscript refers to integer time steps. In the corrector step, the initial guess is improved with The constant strain-rate experiments may consist of two phases. The first stage is the hardening part where a strain rate is applied. In the relaxation parts, a sudden change of strain rate, lower than the hardening part, is imposed. We solve the hardening and the relaxation part of the stress-versus-strain curves separately with appropriate initial conditions. The initial condition for the hardening part is straightforward with y = 0 . The transition for the strain rate stepping experiments follows the continuity condition, such that the initial condition for the simulation at a different strain rate is the value of the state vector y at the end of the previous phase. The nontrivial initial conditions for each experiment are given in Table 1. The total time of integration is the maximum strain divided by the imposed strain rate. For dry and wet dunites, the maximum time of integration needed is of the order of 20,000 s.
We perform a preliminary exercise, where we model each stress-versus-strain curve for dry dunites with an independent set of parameters within the nonlinear Burgers rheology (Additional file 1: Table S2; Figure S3). Each curve can be predicted with a great degree of precision, suggesting the suitability of the proposed model. However, the inversion is poorly constrained and the resulting constitutive parameters are not unique (Additional file 1: Table S2). The aim of our study is to determine the transient creep flow law parameters uniquely within error bars, which requires inverting all the stress-versus-strain data sets with a unique set of constitutive parameters.
To do so, we first solve the governing equations using a numerical solver for a given set of parameters for all thermodynamic conditions. Then, we determine the flow law parameters that best describe all the laboratory observations simultaneously and their uncertainties using a Markov chain Monte Carlo (MCMC) method. We use Gibbs sampling, a type of MCMC technique, to explore the model space, as it provides a simple and efficient technique to search for the best-fit parameters and their uncertainties, even for nonlinear problems (Geman and Geman 1984;Casella and George 1992;Hang et al. 2020). For example, Gibbs sampling has been used to constrain the steady-state flow law parameters of olivine considering multiple mechanisms operating at the same time (Korenaga and Karato 2008;Jain et al. 2019). For a given set of physical parameters, we simulate the mechanical response, i.e., the stress-versus-strain series, for all the experimental conditions. We then compare the synthetic stress-versus-strain series to all the experimental data. We compute the global misfit between all the experimental and simulated data, assuming a uniform error in stress of ±2 MPa (Hirth and Kohlstedt 1995). The global misfit is reduced in each iteration of Gibbs sampling, and the procedure, therefore, represents a global optimization of all the model parameters. For the anhydrous case the free parameters are A K , G K , n K , E K , and G M . For the hydrous case, they are A K , G K , n K , p K , E K , and G M .
Gibbs sampling requires only a few tuning parameters: the number of samples to be generated in computing the conditional distribution for each parameter and the number of burn-in and step samples used to reduce the effect of artificial correlation and initial guess. We benchmark the Gibbs sampling inversion using synthetic data generated under similar thermodynamic conditions to those of the experimental data (Table 3), and we verify that we can recover the true values accurately (Fig. 2). With Gibbs sampling, all the proposed samples are accepted, leading Pre-exponential factor (log 10 A K ) Hardening coefficient (log 10 G K ) Stress exponent (log 10 n K ) Activation energy (log 10 E K )  Fig. 2 Gibbs sampling benchmark. We generate synthetic data using similar thermodynamic conditions to those of the experimental data using A K = 10 5 s −1 MPa −n K , G K = 10 4.602 MPa, n K = 10 0.477 , and E K = 10 5.6 kJ/mol. We generate 2 × 10 5 samples and discard the initial 10 4 samples. The histograms are formed using every 100 th samples of the remaining set. The mean of the marginal distribution matches the true value (dashed black line) from which the synthetic data were generated to correlated successive samples. The autocorrelation of Gibbs samples becomes negligible after a lag of 100 samples for all the optimizations considered in the study. We follow the same procedure for all cases described below: we generate 2 ×10 5 samples with conditional evaluations of 1000 samples in each dimension. We then discard the initial 10 4 samples of the burn-in process and we retain one in every hundred of the remaining samples.

Flow law parameters for dunites
In this section, we determine the transient creep flow law parameters and their uncertainties for dunites in the hydrous and anhydrous cases using the Gibbs sampling inversion procedure.

Dry dunites
We first apply the MCMC approach to analyze the experimental data of dunites under anhydrous conditions, assuming the steady state flow law parameters from Chopra and Paterson (1984). The thermodynamic conditions used for the dry case are listed in Table 1.
Using the retained samples from our MCMC simulation, we generate the marginal and joint distributions of the corresponding model parameters (Fig. 3). For all cases considered, we find similar values for the mean and the mode of the marginal distributions, so we do not discuss the differences. We find a strain-hardening coefficient G K = 33 ± 0.3 GPa, significantly smaller than Young's modulus G M = 92 ± 2 GPa. We find the effective activation energy E K = 430±20 kJ/mol and stress exponent n K = 2.0 ± 0.1 for the transient creep, lower than those of the steady-state creep previously established by Chopra and Paterson (1984), i.e., E M = 535 ± 15 kJ/mol and n M = 3.6 ± 0.1.
Using the mean of the marginal probability distributions as input parameters, we simulate the deformation, and compare it to the stress-versus-strain experimental curves (Fig. 4). The data at fixed strain rate are explained well. The misfit is higher than when each curve is optimized individually (Additional file 1: Figure S3) due to the tradeoffs that naturally occur when all data are optimized simultaneously. The model also explains the strain rate stepping experiments, but the fit to the relaxation phase, at a lower strain rate, is Table 1 Experimental setting for dry dunite deformation from Chopra and Paterson (1984). All the stress-versus-strain series are available from their supplementary material. The initial conditions used to solve the ordinary differential equation for the MCMC inversion and the best-fit model are provided here poorer than the initial hardening phase. This could be due to change in the internal stress of the mineral, or permanent changes in the microstructure due to the exceedingly large strain reached, affecting transient creep more than the steady state behavior. In other words, the internal stress might have a separate evolution equation.

Run number T (K)ε in s
The overall fit to run 4668 is less satisfying and is the only run, where the steady state was achieved much earlier than other experiments. The total strain of experiment 4668 is also low compared to all other experiments. It is possible that model parameters are more strongly influenced by other data, where, in general, the total strain is much larger. This is not a fundamental limitation of the nonlinear Burger, as the 4668 data can be fit perfectly in the individual curve fitting exercise (Additional file 1: Figure S3). In fact, small variations in the model parameter can improve the fit to certain data as can be seen from the confidence interval of model predictions (Fig. 4).
We also explore the potential of the internal state variable (ISV) constitutive framework of Sherburn et al. (2011) that involves 20 parameters, including the isotropic and anisotropic hardening coefficients. With this many degrees of freedom and no data currently available to constrain all the parameters, the ISV model may provide a superior predicting capability. Sherburn et al. (2011) used a subset of the data of Chopra and Paterson (1984) to estimate the best-fit parameters, but the resulting ISV model performs unsatisfactorily against the remaining data (Fig. 5). In contrast, the nonlinear Burgers model optimized using the same data set as for the ISV model is fully constrained, provides a better fit to the data, and involves only 8 parameters for transient and steady-state creep. While the model of Sherburn et al. (2011) may be a valid formulation for the deformation of olivine, the principle of parsimony and Occam's razor suggest that the simpler model that we suggest, providing even a superior fit to the data, is more meritorious.

Steady-state creep flow law parameters of wet dunites
Our main focus in this study is to determine the transient flow law parameters for dunites under hydrous and anhydrous conditions, assuming the steady-state flow law parameters from previous studies. In the case of anhydrous conditions, steady-state flow laws for the dunites are known from Chopra and Paterson (1984). However, for the hydrous condition, the steady-state flow laws of Chopra and Paterson (1984) did not account for the effect of grain size. Here, we explicitly incorporate the grain-size dependency into our model assuming a steadystate flow law of the following form: where all the symbols have been defined previously.
We determine the flow law parameters at steady state using a linear least-squares procedure and the results are summarized in Table 5 and in Fig. 6. We find that the activation energy and stress exponent of steady-state creep under hydrous condition are E M = 450 ± 60 kJ/mol and n M = 3.6 ± 0.5 , respectively, in agreement with the previous results for DisGBS (e.g., Hansen et al. 2011;Ohuchi et al. 2015;Hirth and Kohlstedt 2003). The resulting grain-size exponent for wet dunites is p M = 0.6 ± 0.2 , indicating that the grain size effect is small compared to that of stress on the deformation of wet dunites. Similar values of the grain size exponent of 0.7 ± 0.1 and activation energy of 445 ± 20 kJ/mol have been reported for San Carlos olivine aggregates under anhydrous conditions in the DisGBS regime (Hansen et al. 2011) and a grain-size exponent of ∼ 1 and activation energy of 423 ± 56 kJ/mol have been reported for DisGBS under dry and wet conditions (Ohuchi et al. 2015). Our estimate of the grain size exponent also agrees with the theoretical model of Langdon (1994). However, a higher value of the grain-size exponent around 2.0 under anhydrous conditions in the DisGBS regime has been reported (e.g., Hirth and Kohlstedt 2003;Wang et al. 2010).
The effect of water content on the strength of dunites was difficult to quantify in early experiments due to partial melting of specimens under experimental conditions ( Van der Wal et al. 1993). The experiments on Åheim dunite likely included more partial melt than with the Anita Bay samples. Due to preferential partitioning of water into the melt and different subgrain spacing between wet Anita and wet Åheim dunites, the water content in Åheim could be lower than in the Anita Bay dunite (Hirth and Kohlstedt 1996). Hence, the grain size effect discussed here may be biased by a different water content between the Anita Bay and Åheim dunites. Further experiments in the DisGBS regime with a wider range of grain sizes are needed to quantify the grain size dependency of the constitutive flow law more definitely.

Transient creep flow law parameters of wet dunites
Using the steady-state flow law parameters determined above and the data shown in Table 2, we conduct an MCMC simulation to determine the transient creep flow law parameters for wet dunites. Similar to the anhydrous case, marginal and joint distributions of the samples are generated with Gibbs sampling (Fig. 7). The distribution of the model parameters approximately follows a multivariate Gaussian distribution. The parameter values and uncertainties are summarized in Table 5. Using the mean from the marginal distributions, we plot the data against the model fit (Fig. 8). Unlike the anhydrous case, no strain-rate stepping experiment data are present for the hydrous case. The constitutive framework explains the 14 stress-versus-strain curves satisfactorily (Fig. 8).
The hardening coefficient G K = 13 ± 0.3GPa for the wet case is much smaller than for the dry case and also an order of magnitude smaller than the corresponding Young's modulus ( G M = 115 GPa). We find an activation energy for the transient creep of wet dunites of 250 ± 12 kJ/mol, which is lower than the activation energy of the transient creep of dry dunites of 430 ± 20 kJ/mol, and also lower than the activation energy of the steady-state creep of wet dunites of 450 ± 60 kJ/mol. We find the same stress exponent of transient creep in both hydrous and anhydrous cases, i.e., n K ∼ 2.0 ± 0.1 . The grain-size exponent p K ∼ 0.7 for the transient creep is similar to the steady-state value, suggesting the same underlying mechanism.

Validation with other experimental data
Since dunites contain more than 94% of olivine, we test how our model performs against the data of dry and wet San Carlos olivine aggregates from Karato et al. (1986). For the wet data, we use the transient and steady-state flow law parameters determined from our study that includes the effect of grain size. For the dry data, we use the transient creep parameters from our study and the steady-state parameters from Karato et al. (1986). Considering that the compositions of dunites and San Carlos olivine aggregates are slightly different (Additional file 1: Tables S3, S4), we only optimize the material-dependent parameters A K and A M . Our model performs well against both the dry and wet data, including for the strain-rate stepping experiment, using A K = 10 7.60 s −1 MPa −n K and A M = 10 5.47 s −1 MPa −n M (Fig. 9a& 9b). We also compare the model Table 2 Experimental conditions for the deformation of wet dunites from Chopra and Paterson (1981). All the experiments were reported in Chopra and Paterson (1981). However, some stress-versus-strain series were published as part of Chopra and Paterson (1984) and Chopra (1986) Chopra and Paterson (1984) Stress exponent n M 3.6 ± 0.1 - Chopra and Paterson (1984) Activation energy E M 535 ± 15 -kJ/mol Chopra and Paterson (1984)   Gibbs sampling. We generate 2 × 10 5 samples, discard the initial 10 4 samples and from the remaining samples, we select every 100 th samples to produce the univariate and bi-variate marginal distributions. The probability density of the model parameters approaches that of a multi-variate normal distribution with the constant stress experiment of Chopra (1997) on wet dunite, where the transient phase is most visible. The data and model predictions agree well using A K = 10 2.26 s −1 µ m −p K MPa −n K and A M = 10 4.8488 s −1 µ m −p M MPa −n M despite some long-wavelength oscillations (Fig. 9). These results suggest that the proposed constitutive framework may explain transient creep and steady state creep of upper-mantle rocks for a broad range of conditions.

Discussion
Recent advancements in the modeling techniques of postseismic deformation have provided new insights into the mechanics of deformation of the upper mantle. To further assist these endeavors, we provide a deformation map at conditions relevant to the asthenosphere using the flow parameters for DisGBS derived from this study (Fig. 10). The deformation map indicates that Dis-GBS may be an important deformation mechanism in the asthenosphere. Other investigators have indeed Stress-versus-strain series data (black dots) for dunites with best-fit models. The red line represents the best-fit model prediction with the nonlinear Burgers constitutive law. The grey band represents the confidence interval of the forward model, which is plotted by randomly sampling the posterior distribution taking the covariance into consideration argued that the DisGBS mechanism may be the dominant deformation mechanism under upper mantle conditions (Wang et al. 2010;Hansen et al. 2011;Ohuchi et al. 2015;Bollinger et al. 2019). If so, our study determines the flow law parameters for both the transient and steady-state creep of the wet DisGBS regime and may allow more consistent models of mantle dynamics. Nevertheless, our constitutive formulation and associated parameters have known limitations. With the available experimental data, we are not able to estimate the water content exponent nor the activation volume for the transient creep of dunites. The flow law parameters for olivine single crystals could be used with some caution, as the slip systems controlling the transient and steadystate creep may change as a function of water content, temperature, and depth or pressure (Masuti et al. 2019). Furthermore, several points for the application of our proposed flow law to postseismic data analysis must be Stress-versus-strain data and best-fitting predictions of the ISV and nonlinear Burgers models for dunites in anhydrous conditions. a-c) The model fits to the data used by Sherburn et al. (2011). All the constants from C 1 to C 20 are from Table 1 of Sherburn et al. (2011), except for C 13 = 0.83 , which corrects a typographic error. The reduced chi-square value for the ISV and nonlinear Burgers models are 41.1 and 13.9, respectively, indicating that the nonlinear Burgers model performs better than ISV. d-f Comparison with data not used in the optimization. Reduced chi-square values for the ISV and nonlinear Burgers models on these data are 411.2 and 105.1, respectively. The nonlinear Burgers model, which uses fewer parameters compared to the ISV model, performs better considered. First, the pre-exponential factors A K and A M incorporate a water contribution and need to be adjusted depending on the actual water content and its power exponent. Second, the activation energies E K and E M represent apparent activation energy and an appropriate activation volume needs to be chosen to account for the pressure effect. It has been argued that the pressure effect on the DisGBS mechanism under wet conditions is small (Ohuchi et al. 2015), but more experimental studies are needed to quantify the pressure effect on transient creep and steady-state creep. The transient creep parameters offer new insights into the possible mechanisms of deformation of olivine aggregates in the upper-mantle conditions of pressure, temperature, and water content. We explain the low activation energy of transient creep by extending a theoretical model for plastic anisotropy, whereby the transient and steady-state creep represent the activation of the soft-slip system and hard-slip system, respectively (e.g., Karato 1998;Ashby and Duval 1985;Masuti et al. 2019). We consider a model, where the dislocation climb rate depends not only on the diffusion coefficient but also on the jog density, i.e., the number of jogs per unit length (Karato 2010): where c j is the jog density and D is the relevant diffusion coefficient (e.g., Karato 2008;Poirier 1985). The activation energy of the climb-controlled dislocation creep E creep is the sum of the activation energy of diffusion and the activation energy of jog formation, following where E D and E j are the energies of diffusion and jog formation, respectively, and, in the case of olivine, E D is for the diffusion of Si (For details, see Weertman 1968;Kohlstedt 2006;Karato 2010, and references therein). Following other insights on the water effect on plastic anisotropy of olivine single crystals (Masuti et al. 2019), we associate a high jog density to the soft-slip system and a low jog density of the hard-slip system. Starting from a high jog density, few new jogs are required for deformation to occur and the jog formation energy is low. As  Table 5. We use the steady state data at various temperatures and strain rates from Table-II and  Table-III of Chopra and Paterson (1981). We use all the data except the experimental run, where the steady state was not achieved a result, the activation energy of transient creep may be approximated by the activation energy of diffusion, i.e., E K ≈ E D . In contrast, steady-state creep is controlled by the hard-slip system, where the jog density is low, and more jogs need to be created for deformation to occur. For steady state creep, jog formation energy is higher, so we have E M > E K . Fei et al. (2013) found an activation energy for the diffusion of Si in forsterite and olivine in both dry and wet conditions of E D ≈ 430 ± 20 kJ/mol, in striking agreement with our results of E K = 430 ± 20 kJ/ mol for the dry case, indicating a negligible activation energy for new jog formation during transient creep.
The activation energy of transient creep for the wet case ( 250 ± 12 kJ/mol) cannot be directly compared with results of diffusion Si in olivine, as the rate controlling mechanism for wet DisGBS could be different. We now discuss differing views and interpretations. Based on his results from constant stress experiments, Chopra (1997) argued that the microstructure evolves during the transient phase and that the mechanical steady-state is achieved when the final microstructural steady-state is reached. In this conception, the transient and steady state creep represent the evolving deformation of a unique, continuously unfolding process, possibly glide and climb of dislocations. His interpretation is based on the similarity of his inferred activation energy for transient creep and steady state creep, i.e., 380 ± 100 and 498 ± 38 kJ/mol, respectively, and the similarity of the stress power exponents. Using constant strain rate experiments, we find a lower activation energy and a lower stress exponent for transient creep in both hydrous and anhydrous conditions, with smaller uncertainties. Based on these results and considering the theoretical model of jog density described above, we suggest that the physical mechanisms behind the transient and steady state creep are different, i.e., one requires a higher jog density than the other.
Our results and interpretation for the mechanics of transient creep hold for olivine aggregates at the hightemperature conditions relevant to the upper mantle,     (Hansen et al. 2020). As the rate constant value of [011] c orientation at high temperature is within the uncertainty of the rate constant for olivine aggregates at low temperature, they propose the strong slip system activated by the [011] c orientation as the dominant slip system for the deformation of olivine aggregates in the transient regime. One of the possible reasons for the discrepancy with our interpretation is a change of dominant slip system with pressure (e.g., Raterron et al. 2009), as the dominant slip system is a complex function of water content, temperature, and pressure (see Fig. 4

Conclusion
We explain the experimental deformation data of dunites under hydrous and anhydrous conditions using a nonlinear Burgers assembly. The proposed rheological framework consistently captures the transient creep and steady-state creep of dunites. The activation energy and the stress exponent for transient creep are lower than those of steady state creep. This suggests that the physical mechanism behind the transient creep of olivine aggregates may be different from that of the steady-state creep, possibly due to a higher jog density of the soft-slip system. Consequently, the transient creep of olivine aggregates may be controlled by the soft-slip system, which is saturated with jogs, and the steady-state creep may be controlled by the hard-slip system, which is undersaturated with jogs. Currently available diffusion data for Si, Validation with other data. a, b Stress-versus-strain series (black dots) for San Carlos olivine aggregates from Karato et al. (1986) with best-fit model. The red line represents the best-fit model with the nonlinear Burgers rheology. The grey band represents the confidence interval. All the model parameters used are same as Tables 4 and 5, except for A K and A M (Additional file 1: Tables S3 and S4). c Constant stress experimental data of strain-versus-time series (black dots) for Åheim dunite from Chopra (1997) and best-fit model (red profile) the slowest diffusing species in olivine, support our jog density model in anhydrous conditions. Contrarily to our findings, some recent experimental studies indicate that the hard-slip system controls the deformation in the transient regime. This discrepancy may be caused by a change of dominant slip system as a function of pressure, but more experimental results are needed to resolve whether the soft-or hard-slip system dominate during the transient creep of olivine aggregates in the thermodynamic conditions relevant to the upper mantle. For hydrous dunites, we use a grain-size dependent nonlinear Burgers model to explain the observations, indicating a Dis-GBS mechanism. The effect of grain size under hydrous conditions is the same for both the transient creep and steady-state creep of dunites. Under upper mantle conditions, DisGBS may be an important mechanism of deformation and we provide the flow law parameters for both transient creep and steady-state creep. The constitutive framework and the associated physical parameters provide a reference to model deformation in the lithosphere and upper mantle at the time scales of the seismic cycle and post-glacial rebound.
Additional file 1: Figure S1. Strain time series for wet dunite from Chopra (1997). We use simple basis functions to separate the transient and steady-state components. The overall fit is shown in red.The separated transient and steady-state signals are shown in blue and black lines, respectively. Figure S2. Strain time series for wet dunite from Chopra (1997) and the best fit model (red). Each strain time series was considered separately. Figure S3. Stress versus strain data and best-fitting predictions (red) for dry dunites (Chopra and Paterson 1984), where each curve is fitted separately. Table S1. Model parameters of wet dunites determined from the constant stress experimental data of Chopra (1997). Table S2.
Flow law parameters of dry dunites determined from the constant strain rate experimental data. Each stress versus strain curve is fit individually. Table S3. Flow law parameters of dry San Carlos (SC) olivine aggregates. Steady-state parameters (n M and E M ) for the San Carlos olivine are from Karato et. al. (1986). The best-fitting values for dry dunites are given for comparison. Table S4. Flow law parameters of wet San Carlos (SC) olivine aggregates. Only A K and A M are given and the other parameters are the same as in Table 5. The values for dunite are given for comparison. The San Carlos olivine best-fitting values are within the uncertainties of our wet dunites estimation. Deformation mechanism maps drawn using the flow law parameters from this study for wet DisGBS. Under asthenospheric conditions, DisGBS can be dominating and the flow law parameters determined in this study may be applicable for postseismic studies. We use the water content exponent and activation volume of 0.35 and 8 cm 3 /mol corresponding to [110] c orientation. Under wet condition and at high temperatures, the [100](010) slip-system controls the steady-state creep (Masuti et al. 2019). For the dislocation and diffusion creep, we use the flow law parameters from Hirth and Kohlstedt (2003)