Investigating Compaction by Intergranular Pressure Solution Using the Discrete Element Method.

Intergranular pressure solution creep is an important deformation mechanism in the Earth's crust. The phenomenon has been frequently studied and several analytical models have been proposed that describe its constitutive behavior. These models require assumptions regarding the geometry of the aggregate and the grain size distribution in order to solve for the contact stresses and often neglect shear tractions. Furthermore, analytical models tend to overestimate experimental compaction rates at low porosities, an observation for which the underlying mechanisms remain to be elucidated. Here we present a conceptually simple, 3-D discrete element method (DEM) approach for simulating intergranular pressure solution creep that explicitly models individual grains, relaxing many of the assumptions that are required by analytical models. The DEM model is validated against experiments by direct comparison of macroscopic sample compaction rates. Furthermore, the sensitivity of the overall DEM compaction rate to the grain size and applied stress is tested. The effects of the interparticle friction and of a distributed grain size on macroscopic strain rates are subsequently investigated. Overall, we find that the DEM model is capable of reproducing realistic compaction behavior, and that the strain rates produced by the model are in good agreement with uniaxial compaction experiments. Characteristic features, such as the dependence of the strain rate on grain size and applied stress, as predicted by analytical models, are also observed in the simulations. DEM results show that interparticle friction and a distributed grain size affect the compaction rates by less than half an order of magnitude.

Introduction

Intergranular pressure solution (IPS, sometimes referred to as (dis)solution‐precipitation creep or even “chemical compaction”) is known to be an important time‐dependent deformation mechanism in the wet regions of the Earth's crust, under conditions of relatively low stress and temperature, down to depths of 20 km (Elliott, 1973; Lehner, 1990; Rutter, 1976, 1983; Spiers et al., 2004). Pressure solution is highly relevant for compaction of sedimentary rocks during diagenesis (Croizé et al., 2010; Tada & Siever, 1989; Yang, 2000) and also for shear deformation of dense rocks under greenschist metamorphic conditions (Elliott, 1973; Stöckhert et al., 1999), and in controlling the frictional behavior of faults (Bos et al., 2000; Hickman & Evans, 1992; Sleep & Blanpied, 1992).

Numerous analytical models exist that describe the microphysical processes operating during pressure solution and the kinetics of creep, especially for the case of compaction by IPS (e.g., Lehner, 1990; Shimizu, 1995; Spiers & Schutjens, 1990). These have been tested against a large number of experiments (e.g., Dewers & Hajash, 1995; Raj, 1982; Spiers & Schutjens, 1990). However, these models are limited in their ability to capture realistic packing geometries by the necessity to assume single‐sized grains arranged in a regular packing (cubic or hexagonal). This simplification is required to analytically express grain contact area, grain contact stress distribution, and pore wall area as functions of strain or porosity. In addition, many models assume that no shear tractions exist on a grain contact, that is, that the friction between grains is negligible. However, friction may be significant in compaction experiments on granular aggregates. Finally, these models tend to overestimate the compaction rates of low‐porosity aggregates by up to several orders of magnitude (e.g., Niemeijer et al., 2002; Schutjens, 1991). This exposes a gap in our understanding of the microscale processes operating at these conditions, which limits the confidence with which experimental results can be extrapolated to natural conditions, particularly for those cases that involve relatively low porosities (<15%). Within a geological context, making predictions of the long‐term compaction behavior of natural faults systems, sedimentary basins, and hydrocarbon reservoirs based on laboratory data may therefore not be warranted.

Classically, much of the (numerical) modeling work done on IPS is related to the compaction of aggregates and the upscaling of laboratory results to nature. Several numerical procedures have been proposed that address the stability of dissolving interfaces (Ghoussoub & Leroy, 2001; Koehn et al., 2006), the initiation of stylolitization (Dewers & Hajash, 1995), self‐organization of compacting systems (Ortoleva, 1994), and compaction of sedimentary basins (Dewers & Hajash, 1995; Renard et al., 1999; Yang, 2000). In particular, basin compaction provides scientific challenges with great relevance for industrial applications, for it is well known that the (terminal) porosities observed in sedimentary formations cannot be readily explained by analytical compaction models for pressure solution. Instead, it is often found that geological reservoirs exhibit porosities that are directly related to their maximum burial depth, rather than geological age and burial history (Giles et al., 1998). This (apparent) equilibrium state is in strong contrast with the nature of pressure solution, which is inherently a nonequilibrium process (e.g., Lehner & Bataille, 1985). Attempts have been made to explain this equilibrium state by consideration of plastic (Palciauskas & Domenico, 1989) and viscoplastic (Revil, 2001) creep that would render the instantaneous porosity a function of effective stress or by consideration of suprahydrostatic fluid pressures resulting from densification (Yang, 2000). Alternatively, the existence of a “yield” stress for pressure solution has been proposed on the basis of a vanishing driving force to explain the observed depth dependency of porosity (Revil, 2001; Van Noort et al., 2008). Furthermore, the presence of cations such as Al3+, e.g., resulting from mica dissolution, can dramatically reduce the rates of quartz dissolution (and likely also precipitation, Iler, 1973; Spiers et al., 2004).

In addition to driving the compaction of sedimentary basis, IPS is well known to operate in gouge‐filled faults under upper‐crustal and midcrustal conditions (Chester & Chester, 1998; Evans & Chester, 1995). Densification and creep of faults by IPS have been put forward as a mechanism for long‐term evolution of fault strength during interseismic periods (e.g., Angevine et al., 1982; Sleep & Blanpied, 1992) and have been employed to investigate the frictional stability and rheology of faults (Chen & Spiers, 2016; Den Harto & Spiers, 2014; Niemeijer & Spiers, 2007). In spite of the similarities between sedimentary basin compaction and fault densification, the operation of IPS in fault gouges differs in that pervasive cataclastic processes need to be considered, particularly at shear deformation rates that are relevant for earthquake nucleation (Niemeijer & Spiers, 2006). For instance, the existence of an equilibrium state, as seen in geological reservoirs, can be questioned when grain‐grain contacts are continuously rejuvenated by granular flow.

Simulating creep of a granular aggregate by IPS using the discrete element method (DEM) has several advantages over conventional analytical models: first of all, DEM does not require the assumptions of a single‐valued grain size and regular packing, as the stresses are calculated per particle. The positions, velocities, and forces acting on an individual particle can be evaluated at any given time during the simulation. This grants insights into the distribution of deformation occurring within a sample, which is very difficult to access in real‐world experiments. Second, the governing physics is implemented at the particle scale and various deformation mechanisms can operate simultaneously (for example, pressure solution and grain boundary sliding). By combining or excluding different mechanisms, the relative importance of a single mechanism can be tested. Lastly, and perhaps most importantly, individual parameters such as the interparticle friction coefficient, kinetic constants, and particle stiffness can be readily modified to investigate the role that these parameters may play in the overall deformation behavior. These features make DEM a suitable method for studying the dynamics of granular aggregates, particularly in the context of deformation of fault gouges where granular flow, cataclasis, and fluid‐rock interactions are prevalent.

So far, only a few studies have attempted to simulate compaction by IPS with DEM (Bernabé & Evans, 2014; Zheng & Elsworth, 2012, 2013). Zheng and Elsworth (2012) related the convergence rate of particles in a stressed aggregate to the undersaturation of the solute in the pore fluid through a sophisticated two‐dimensional DEM scheme, but comparison with compaction data on quartz sand reported by Niemeijer et al. (2002) showed that the compaction rates were underestimated by the model. It was argued by the authors that this reflected a difference in grain size distribution between the model and experiment, but this hypothesis was not tested further. Comparison with the open‐system compaction experiments of Yasuhara (2003) showed better agreement. Bernabé and Evans (2014) used a three‐dimensional DEM model to investigate long‐term compaction behavior of quartz with explicit coupling of the interface reaction and diffusion kinetics through the pore fluid saturation state, but the resulting compaction rates were not compared to experimental results. In order to simulate compaction over long time‐scales (up to 10,000 years), Bernabé and Evans (2014) assumed that the DEM aggregate was always in mechanical equilibrium, so that an iterative relaxation technique could be employed for efficient time integration of the aggregate. While the assumption of quasi‐static compaction is believed to be valid in these simulations, it is overly restrictive for the simulation of aggregate deformation in the dynamic regime, which is relevant for e.g. fault and earthquake mechanics. Furthermore, to simplify the implementation process, Bernabé and Evans (2014) chose to omit interparticle friction from the model, which leaves open questions to the role that it might play in the compaction behavior. For direct comparison with experimental results (where grain‐grain and grain‐wall friction may play a significant role), interparticle friction in the DEM model has been included in the current study.

The aim of this work is to present a conceptually and numerically simple way of modeling IPS in DEM, so that it can be readily incorporated in future DEM studies of the mechanics of fault gouges. The relevant grain contact input parameters are the ones used in analytical IPS models (e.g., by Spiers et al., 1990), and elasticity at the contacts is accounted for using a standard DEM approach. To validate the resulting (three‐dimensional) DEM model, we conduct short‐term compaction experiments on granular halite and directly compare the experimental compaction rates with those of the simulations. In this way, we test if the DEM model is capable of reproducing experimental compaction data without prescribing the evolution of the particle contact area with bulk strain, as required by analytical models. We explore if the proposed DEM approach preserves the grain size and stress exponents, as theoretically predicted for frictionless particles arranged in a regular close packing, for a nonzero interparticle friction. We then proceed to investigate the effects of interparticle friction and of a polydisperse grain size distribution on the overall compaction behavior to assess if these effects could explain the discrepancy between compaction rates predicted by analytical models and those observed in experiments. Since granular flow is not yet fully understood, we only consider compaction in the absence of pervasive granular flow in this work, which we leave for future study.

In summary, the aims of this work are to (1) present a method for modeling pressure solution creep using DEM; (2) assess the validity of the model by comparison with dedicated compaction experiments, with theoretical predictions; (3) investigate the effect of grain contact friction on the macroscopic compaction rates; and (4) investigate the effects of a nonuniform grain size distribution on the long‐term compaction behavior. Since analytical models fail to explain the reduction in strain rates at low porosities (Niemeijer et al., 2002; Schutjens, 1991), there is potential for the discrete element method to shed new light onto this problem by systematically investigating parameters such as grain contact friction and grain size distribution. These again will also be highly relevant for the mechanics of fault gouges.

Previous Work on Compaction via Pressure Solution

Constitutive Relations for Pressure Solution

It is widely accepted that the thermodynamic driving force for pressure solution in chemically closed systems is provided by gradients in chemical potential, induced by gradients in stress around grain surfaces. Much work has been done on deriving theoretical models (e.g., Lehner, 1995; Pluymakers & Spiers, 2014; Rutter, 1976; Shimizu, 1995) that relate the rate of deformation to the applied stress, temperature, and nominal grain diameter. These models assume that an interconnected network of fluid is present within the grain boundaries, either as a dynamically stable, island‐channel network (Raj, 1982; Spiers & Schutjens, 1990) or as an adsorbed thin film (Robin, 1978; Rutter, 1983). By dissolving mass at the highly stressed grain contacts and reprecipitating it at lower stress interfaces or at stress‐free pore walls, mass is transferred through the grain boundary fluid and creep is achieved perpendicular to the direction of the locally acting normal stress, producing deformation parallel to the local normal stress direction (i.e., parallel to the contact normal). This assumes that the resistance to grain boundary sliding is negligible, or alternatively that the dissolution/precipitation is only possible normal to grain surfaces (Lehner, 1990; Shimizu, 1995). In the situation where it is assumed that there is no mass transport out of the system, the process of pressure solution is a series of three consecutive serial processes: (1) dissolution of material at the highly stressed grain contacts, (2) diffusion of the dissolved solid out of these contacts into the pore space, and (3) precipitation on interfaces supporting lower stress, including the pore walls. Assuming (near) steady state transfer, the kinetic process that proceeds at the lowest rate will then control the overall rate of deformation (Lehner, 1990; Raj, 1982). Numerous authors have derived constitutive equations for pressure solution at single circular grain contacts (e.g., Pluymakers & Spiers, 2014; Rutter, 1976; Shimizu, 1995; Spiers & Schutjens, 1990). For compaction of a porous aggregate idealized as an array of spherical grains arranged in a regular packing, essentially equivalent equations were obtained, which take the following form for both isostatic and uniaxial compaction (Spiers et al., 2004): Here the subscripts s, d, and p denote dissolution‐, diffusion‐, and precipitation‐controlled kinetics, respectively, is the rate of deformation, A is a geometric constant, d is the mean grain size, σ is the effective stress (applied stress minus pore fluid pressure), Ω is the molar volume of the solid, R is the universal gas constant, T is the absolute temperature, D is the effective diffusion coefficient of the ionic species in the grain boundary, C the equilibrium solubility of the solid, S is the (mean) thickness of the grain boundary zone, and f(ϕ) is a dimensionless function of porosity (ϕ) that accounts for the evolution of contact area and pore wall area during compaction. The rate constants for grain boundary dissolution (I ), precipitation (I ), and for grain boundary diffusion (D C S) vary with temperature and are expected to exhibit an Arrhenius dependence on temperature, that is, to take the form: for I and I (Van Noort & Spiers, 2009), and for (D C S) (e.g., Spiers et al., 1990). Here ΔH is the apparent activation energy associated with each process. For highly soluble materials at room temperature, such as alkali metal salts, dissolution and precipitation are relatively fast, and diffusion along the grain boundary controls the creep rate, even for small grain sizes. By contrast, low solubility solids, such as quartz, often show interface‐reaction‐controlled kinetics (Spiers et al., 2004). The linear form of equation (1a) holds well for small departures from equilibrium, that is, for small contact stresses. However, when the driving force for IPS is large, a nonlinear relation between strain rate and effective stress is more appropriate (see e.g., Niemeijer et al., 2002; Pluymakers & Spiers, 2014).

Spiers et al. (1990) systematically investigated the effect of temperature, mean grain size, and applied normal stress on the deformation (creep) rate of granular halite (NaCl), in the range of 20–90°C, 100–275 μm, and 0.5–2.2 MPa, respectively. Grain sizes were controlled to lie in narrow ranges by sieving (near monodispersed). From the grain size dependence of the creep rate, Spiers et al. (1990) inferred that (at room temperature) diffusion of mass along the grain boundaries was the rate controlling process and derived the values for the constitutive parameters and as appearing in equation (3). These parameters were derived by fitting equation (1b) to the laboratory data set and assuming that the aggregate can be represented by a regular packing of single‐sized spheres, and that the values of (D C S)0 and ΔH are independent of temperature, stress, and grain size. Since the values of these kinetic or constitutive parameters are relatively well constrained for halite (De Meer et al., 2002; Hickman & Evans, 1995; Schutjens & Spiers, 1999; Spiers & Schutjens, 1990; Urai et al., 1986), creep experiments conducted on granular halite will be used here for comparison with the DEM simulation outcomes and for validation of the DEM model.

Effect of a Distributed Grain Size

It is well known that various sedimentary and transport processes can produce a wide range of grain size distributions in sedimentary rocks (McLaren & Bowles, 1985; Visher, 1969). Similarly, cataclastic fault rocks are commonly characterized by a power law grain size distribution (Storti et al., 2003). However, for mathematical convenience, analytical models typically assume a uniform grain size. In order to compare the analytical models with experiments, the mean grain size of the experiment is commonly chosen as the reference grain size. To address the potential effects of a distributed grain size on IPS, Niemeijer et al. (2009) conducted compaction experiments on granular halite aggregates with various mean grain sizes and distributions. They observed a strong correlation between the coefficient of variation (standard deviation normalized with respect to the mean) of the grain size distribution and the measured volumetric strain rate. Specifically, the strain rate was found to increase markedly with increasing coefficient of variation (i.e., with a widening distribution). Niemeijer et al. (2009) attributed this behavior to the presence of the smaller grains within the wider distribution, which allowed for enhanced creep rates. Since the theory of diffusion‐controlled IPS predicts an inverse cubic dependence of strain rates on grain size (equation (1b)), it is expected that the mass transfer rate will be higher for smaller grains if the contact stresses are the same. However, in a porous aggregate where the majority of the load is carried by a network of larger grains, effectively shielding smaller grains, the stresses acting on the smaller grains might be much lower than those acting on the larger grains, thus leading to lower macroscopic sample strain rates. Which of these competing effects dominates the overall compaction behavior of a porous aggregate is not clear.

To investigate the effect of a distributed grain size on creep of polycrystalline materials, Ter Heege et al. (2004) proposed a composite flow law for grain size sensitive and grain size insensitive creep, taking into account the possibility of a lognormally distributed grain size. This flow law was derived for two end‐member assumptions: that of sample‐homogeneous stress or of sample‐homogeneous strain rate. When either end‐member assumption is applied to IPS, it is found that the predicted creep rate decreases with increasing coefficient of variation. It should be noted that the model was derived for a bulk crystalline material, rather than a porous granular aggregate: the assumption that either the stress or the strain rate is homogeneously distributed over the sample is not valid for high‐porosity aggregates. However, as a granular sample densifies, that is, porosity decreases, it is expected that the assumptions of the model by Ter Heege et al. (2004) are more closely met, and that the dependence on the width of the grain size distribution changes. This hypothesis is difficult to test experimentally but can be tested numerically using the discrete element method.

Numerical Methods and Procedure

Basic DEM Approach

The discrete element method is a widely used numerical method for modeling the mechanics of discontinuous (granular) media and aggregates. While the method as originally proposed by Cundall and Strack (1979) has remained essentially unchanged, ongoing advances in computational performance are allowing for larger, more complex simulations that resemble the geometry and scale of the object of study increasingly well. The DEM method involves identifying which particles interact with one another, calculating the forces associated with each particle‐particle contact and hence acting on each particle, and integrating the accelerations (as calculated from forces through Newton's second law of motion) to obtain individual particle velocities and positions. By prescribing the way each particle interacts with its neighbors at the microscale (e.g., Hookean or Hertzian, elastic frictional), DEM aims to realistically capture the macroscopic material response to dynamic or kinematic boundary conditions imposed on the system. In the present study, we employ the open‐source software package granular LAMMPS (Landry et al., 2003; Plimpton, 1995), as modified by Marketos (2013) and further modified for pressure solution.

For simplicity, the particles considered here are perfectly spherical, indestructible, and unbonded, in the sense that no tensile forces are allowed and no moments of force are transmitted between them. The particles can be of different size, so that a realistic grain size distribution can be simulated. A size‐independent (constant) coefficient of friction is used to describe the sliding resistance at grain contacts. All the stresses (and forces) present in the system are assumed to be effective stresses (total stress minus the pore fluid pressure), and local excess pore pressures are assumed to dissipate instantaneously, as porosities are high. When two particles are in contact, the elastic restoring force acting in a direction perpendicular to the contact plane is typically calculated through a force‐displacement law, for example, through equation (4) below for a linear contact spring: where F is the normal contact force, k is the stiffness, and Δx el is the elastic center‐to‐center displacement or conceptual overlap between the particles (see Figure 1). A similar equation applies for the shear forces. The elastic shear displacement Δy is incremented each time step by the relative particle displacements, and the maximum shear force is constrained by a frictional slider. Since most natural materials feature irregular, indented or even sutured grain surfaces, it is unlikely that stressed contacts in such a material exhibit perfect Hertzian elasticity. Instead, many materials seem to display stress‐strain relationships that are approximately linear (e.g., Nakata et al., 1999) or fall between idealized Hookean and Hertzian behavior (e.g., Cole & Hopkins, 2016). Furthermore, the operation of pressure solution rapidly grows and flattens the grain contacts, so that the elastic response will quickly become sub‐Hertzian. In the absence of a unifying elastic theory that is practical for implementation in DEM, we simply adopt a linear force‐displacement law.

Figure 1

The elastic‐frictional interaction between two particles sharing a single contact is represented by a linear spring in the normal direction, and a spring and frictional slider in the shear direction. The force perpendicular to the plane of contact is calculated as , where k is the normal contact stiffness and Δx el is the elastic overlap in the normal direction, Δx is the total overlap, and δ is the nonelastic overlap (i.e., dissolved contact layer thickness). The elastic restoring force tangential to the plane of contact is calculated incrementally as , as it is typically done in DEM, where k is the shear stiffness and Δy is the lateral relative displacement caused by relative motion tangential to the plane of contact or by relative rotation. The shear force is limited by the frictional slider to F =μ F , with μ being a (constant) coefficient of friction.

In a system where IPS is possible, besides the elastic and frictional interactions usually considered in DEM, as soon as a stress acts on a particle contact, pressure solution will be activated at a rate V ps that is proportional to the effective stress on the contact. In order to mimic the approach of grain centers and overall densification of the system in DEM, we propose here that the center‐to‐center displacement or total overlap Δx should be reduced by the total nonelastic overlap (i.e., the dissolved particle contact layer thickness) δ, so that the resulting, elastically maintained contact force after dissolution becomes

In each time step, δ is then incremented by V ps×Δt. When the forces applied at the boundaries of the system are held constant (i.e., simulating a creep test), the average force on the contacts, and hence Δx el, remains constant. Continuous dissolution at the contact (i.e., increasing δ) will then result in an increasing total overlap Δx , producing compaction of the system over time. Alternatively, if the positions of the boundaries are kept fixed after some initial force is applied (i.e., simulating a stress relaxation test), the average total overlap Δx remains constant, but δ increases with ongoing dissolution. This causes F to decrease and the sample stresses to relax.

The rate of pressure solution V (i.e., the rate of change of δ) follows directly from the theory of pressure solution for a single contact (Lehner, 1995; Pluymakers & Spiers, 2014; Rutter, 1976) and can be written in the form

Here the subscripts i and d represent interface‐reaction controlled (dissolution or precipitation) and diffusion controlled, respectively; V ps is the rate at which the overlapping zone “dissolves” (i.e., at which δ increases); Z is a kinetic constant representing the temperature‐dependent kinetics of dissolution, precipitation (units: m Pa−1 s−1), or diffusion (units: m3 Pa−1 s−1); and A is the area of the circular intersection between the two spherical bodies, which is assumed to be fluid filled. The pore fluid itself is not modeled directly and is assumed to be homogeneous in composition (solute concentration) and present throughout the system. Because pressure solution on the contacts effectively acts as a linear viscous dashpot, no additional viscous damping is used in the simulations, except for when a new sample is created (as will be discussed later).

It is well known that the rate of interfacial ion detachment/attachment in highly soluble salts, such as NaCl, is sufficiently rapid that mass transfer rates are controlled by diffusion, even over short length scales. By comparison of end‐member constitutive relations for diffusion‐ and interface‐reaction‐controlled IPS, and by assuming values for (D C S) and the reaction rate constant of 10−19m3 s−1 (Spiers et al., 1990) and 10−4 m s−1 (Alkattan et al., 1997; Scrutton & Grootscholten, 1981), respectively, it can be shown that IPS is rate limited by diffusion for length scales larger than about 60 nm. In practice, this means that the kinetics of IPS of halite will always be controlled by diffusion (Spiers et al., 2004), even for small grain contact areas that are characteristic for the early stages of a compaction experiment. Hence, in this study we only consider diffusion‐controlled kinetics, governed by equation (6b) in the simulations. Note that equation (6a) could have easily been implemented as well.

While pressure solution operates (i.e., V ps>0), the total overlap Δx of the two spherical particles continuously increases. The total area of contact A is calculated each time step as the area of intersection of two spheres, which increases nonlinearly with Δx . Hence, it follows from equation (6a) that pressure solution is a self‐limiting process: with ongoing pressure solution, A increases and V ps decreases. This makes our DEM formulation self‐consistent, without further need for a porosity function (see equation (1a)) that accounts for the evolution of the contact area with increasing compaction. Furthermore, at constant F , the contribution of elastic deformation to the total overlap remains constant at Δx el=F /k, whereas the contribution from dissolution by pressure solution (δ) continuously increases. For bulk strains larger than the average elastic strain, pressure solution quickly dominates in determining the size of the particle contact (i.e., δ ≫ Δx ). Following this line of reasoning, it is expected that even for modest compaction strains, the choice for the elastic parameters plays a second‐order role on the rates of deformation (as will be discussed in section 3.3).

For the entire duration of the simulations, the diameter of the particles remains constant, and the particle shape does not deviate from spherical. In our current DEM approach, reprecipitation of the dissolved mass onto the pore walls is not modeled explicitly, in the sense that mass redistribution by pressure solution does not result in local reduction of the pore volume. Hence, the area of the dissolving contacts is not affected by reprecipitation. For relatively small strains (i.e., high porosities), the total contact area is small compared to the total pore wall area and this assumption holds. For large strains, however, reprecipitation should be taken into account, and a method accounting for the mass redistribution during the DEM simulation (like that employed by Zheng & Elsworth, 2012) should be considered. Preliminary calculations show that for a porosity of 10% the relative error in contact area is less than 10%. Note that for the porosity calculations we only consider the total solid mass present in the system, which, for a chemically closed system (with respect to the solid mass), is constant in time, so that our estimates of bulk porosity are not affected. Additionally, since the fluid phase is not modeled explicitly (i.e., no advective transport is considered), mass conservation is implicit in our DEM approach.

Numerical Stability and Scaling

To ensure numerical stability, the time step employed during the simulation is constrained to Δt = t /10. The critical time step t is defined as , where m is the mass of a particle and chosen such that a mechanical perturbation of a given particle cannot travel farther than its direct neighbors and is proportional to the highest eigenperiod of the system (Cundall & Strack, 1979). For realistic particle sizes, of the order of 10–100 μm, the limiting time step reduces to values that are computationally not feasible for simulating the timescales of real‐life experiments (hours to days). Therefore, in order to simulate experimental procedures, either the kinetics of pressure solution have to be nondimensionalized (i.e., scaling time) or the limiting time step has to be artificially increased by increasing the mass density ρ of the particles (i.e., density scaling), a technique commonly used in DEM (e.g., Cui & O'Sullivan, 2006; Ng, 2006; O'Sullivan et al., 2003). Here we increased the time step by using particle sizes near unity while keeping the mass density at 2,100 kg m−3. Then, we chose to nondimensionalize time and strain rate as where Z is the diffusion‐controlled kinetic constant, and σ is the externally applied axial stress. Following numerous studies (e.g., Pluymakers & Spiers, 2014; Rutter, 1976), the diffusion‐controlled kinetic constant can be expressed as where the product (D C S) is related to temperature according to equation (3). Note that the above equation only consists of thermodynamic quantities and material properties and is therefore independent of the assumed packing geometry, grain size, and stress distribution within the sample. Hence, Z is a constant and can be used in scaling. The (minimum) particle radius r scales all spatial dimensions (i.e., ), which makes all dimensions of the simulation domain proportional to the particle size. Since mass is not rescaled, this approach is essentially equivalent to density scaling, while being more convenient for obtaining the optimal degree of scaling. When scaling in this way, particle accelerations and inertial effects need to be closely considered, so as to remain within the quasi‐static regime of deformation. A suitable degree of scaling was found in small‐scale sensitivity tests to achieve reasonable computation times, while minimizing the effects of inertia on the macroscopic compaction behavior as a result of the chosen degree of scaling.

Validation Tests

The validation of the model was done by comparing experimental compaction rates (see section 4.1) with those from DEM simulations. These simulations will from here on be referred to as simulation set A. The following procedure was adopted in the validation DEM model run: the DEM sample was rendered by generating 20,000 particles at random locations within a cylindrical simulation domain of similar (scaled) diameter as the real‐life counterpart, after which they were allowed to settle by gravity. This mimics the way lab samples were prepared, where the powdered sample material was poured into a cylindrical, 1‐D compaction vessel. During the settling of the particles, particle body forces were artificially reduced (“damped”) to lower the time for the aggregate to come to rest. After settling, gravity and artificial damping were turned off, and the DEM sample was precompacted “dry” by applying an axial compressive load cycled between 1 and 10 times the load used during “wet” compaction. Load cycling was repeated until the porosity did not change by more than 0.1% (absolute value) between consecutive cycles. All simulations were conducted in a zero lateral strain (i.e., uniaxial strain) geometry, and the load was imposed by downward migration of the planar wall bounding the top of the simulation domain, the area of which was constant throughout the simulation. A proportional controller algorithm was used to maintain the stress (force divided by the wall area) exerted by the top wall at a prescribed value (within 0.01% of the target value), as described by Marketos (2013). This essentially sets the boundary velocity to be proportional to the difference between a target stress and the instantaneous stress supported by the wall. By adopting this procedure of load cycling, a realistic initial porosity of 40.43% was obtained for this particular sample. After precompaction, a constant axial load was applied to simulate creep testing conditions. Pressure solution was then activated numerically, after the sample had reached equilibrium, simulating the addition of a saturated solution phase at atmospheric pressure, and the sample was allowed to compact. The parameters used in this simulation, as well as their real‐world counterparts (if applicable), are reported in Table 1.

Table 1

The Relevant Parameters as Used in the Simulations Aimed to be Compared to and Validated by Experimental Data (as Reported in Figure 2)

		Value (range)
Parameter	Symbol	Experiment	Simulation	Units
Grain radius	r	33.5 (31.5–35.5) ×10−6	1.06 (1.00–1.12)	m
Mass density	ρ	2.1 × 103	2.1 × 103	kg m−3
Compaction vessel radius	‐	1.075 × 10−3	34.127	m
Axial stress	σ n	1.08 × 106	1.08 × 106	Pa
Particle stiffness	k	‐	2.65 × 109	N m−1
Particle friction coefficient	μ p	‐	0.6	‐
Wall friction coefficient	‐	‐	0.1	‐
Kinetic constant	(D C S)0	2.79 (1.39–4.19) ×10−15	‐	m3 s−1
Activation energy	ΔH d	2.45 (2.36–2.55) ×104	‐	J K−1
Molar volume	Ω	2.69 ×10−5	‐	m3 mol−1
Pressure solution constant	Z d	9.58 (4.77–14.4) ×10−27	10−13	m3 Pa−1 s−1
Number of particles	‐	‐	20,000	‐
Time step	Δt	‐	1.82 × 10−4	s

Note. Only those symbols are given that were reported in the text. The parameter values that were not known or not applicable to the experiments are not provided. The values of (D C S)0 and ΔH, which are used to calculate Z (at 40°C), are taken from Spiers et al. (1990).

Figure 2

Experimental compaction results for wet halite (applied stress σ =1.08 MPa; starting grain size 63 < d < 71 μm; and average starting porosity 43%) as compared with results from simulations conducted at similar conditions and with the results of the analytical model by Pluymakers and Spiers (2014) (model “D3”). The starting porosity for the simulations was 41%.

The validation of the DEM model was done by comparing the macroscopic compaction rates produced by the model with those measured in the benchmark experiments, in the absence of grain boundary healing. The implementation of the pressure solution contact law was considered to be valid if the strain rates produced by the model at a given porosity would be in the same range as the strain rates measured in the experiments at identical values of porosity. To take into account the uncertainty of the kinetic parameters as reported by Spiers et al. (1990), the simulations were scaled to a minimum, a maximum, and a mean estimate of the parameter Z from equation (7b) (see Table 1 for the parameter values and Figure 2 for the results). The (scaled) kinetic parameter was adopted as reported by Spiers et al. (1990), and the model was not further “calibrated” by varying Z in the DEM simulations. Note that, although the value of Z was derived from laboratory experiments conducted in a similar way as reported here, there is no circularity in our approach: the analytical models that were compared with experiments to constrain Z were derived for a single contact and thus require upscaling to the size of the aggregate to simulate the evolution of particle contact area and stress with sample strain or porosity. In the DEM approach, this evolution is virtually unconstrained, as only the shape of the particle‐particle contact is determined by the DEM assumption of spherical particles. Hence, the characteristic evolution of strain rate with sample porosity emerges spontaneously from the DEM simulations, which is then directly compared with laboratory data.

The possible role that the particle stiffness may play is systematically investigated by repeating the above procedure for different DEM samples, with each simulation exhibiting a particle stiffness that is either identical to, 10 times lower, or 10 times higher than the value used in the reference simulation described above (simulation set A; see also supporting information Text S2). For each chosen particle stiffness, three unique samples (“sister samples”) were generated and deformed to get a measure of the numerical sample variability. At 32% porosity, the variability in strain rate between the sister samples lies within 1% of the average value (supporting information Figure S4). By comparison of the simulations exhibiting different stiffnesses, it is seen that even such a large variation in stiffness (2 orders of magnitude) has only a minor effect on the compaction rates after strains of a few percent, particularly when considering the uncertainty in (D C S) value and variability of the laboratory experiments. These results also show that the simulation outcomes are not very sensitive to the initial packing configuration, such as local porosity variations, and that the compaction behavior is mostly controlled by the IPS kinetic parameter Z .

After validation of the contact model by direct comparison with experiments, additional verification simulations were run to investigate the model dependency on stress and grain size. The aim of these numerical experiments was to compare the DEM model trends with theoretical predictions. For convenience, new samples were created for each simulation set following a somewhat different numerical procedure than was adopted for simulation set A. First, a total of 4,500 particles were generated at random locations in a rectangular simulation domain. The dimensions of the rectangular simulation box were scaled to the particle size, so that the lateral sides of the box had a length of 15 particle diameters. To eliminate possible boundary effects, periodic lateral boundaries were employed, that is, particles at a vertical side of the simulation box interacted with those at the opposite side. In the axial direction, regular (nonperiodic) boundaries were employed. To investigate the possible role of friction, two interparticle coefficients of friction were tested (μ =0.01 and μ =0.6). The DEM sample was precompacted by cyclically varying the load between 1 and 10 times the load that was used during the rest of the simulations. A total of five cycles were performed. In the simulations that were conducted to test the stress dependence of the model (from here on referred to as simulation set B), the axial stress was varied between simulations in the range of 0.1–1.0 MPa, while the particle diameter was set to unity for all particles. In the tests that were conducted to investigate the grain size dependence (simulation set C), the applied axial stress was 1.0 MPa and the particle radius was varied from simulation to simulation in the range of 0.1–1.0 units of length. Pressure solution was then activated, and the sample was compacted until a final porosity of 30% was reached. The value of the kinetic constant Z was set to 10−13 in simulation set B and to 10−14 in simulation set C. These values were chosen so as to maintain a reasonable computation time for the chosen particle sizes. Note that these particle sizes were arbitrarily chosen, and that therefore the spatial dimensions were not physically meaningful. However, for the compaction experiments, the inertia of the particles did not affect the observed compaction rates, and so the results are representative.

Testing Effects of Friction and Distributed Grain Size

Since shear stresses have been typically excluded from analytical expressions for IPS, it is worthwhile to systematically investigate the effect of friction on the macroscopic rate of compaction. To this end, a new set of simulations (simulation set D) was run following a similar procedure as was adopted for simulation set B. In these simulations, a sample consisting of 4,500 particles was generated with an interparticle friction coefficient of 0.01. After precompaction and load cycling, the aggregate was subjected to an axial stress of 1 MPa. Using this precompacted sample as a starting point, eight new simulations were initiated, each with a different interparticle friction coefficient of 0.01, 0.05, 0.10, 0.20, 0.30, 0.50, 0.70, and 1.00. Then, pressure solution was turned on and each simulation was run for 1,000 model seconds, and the resulting macroscopic compaction rates were compared. This procedure ensured identical starting microstructure (porosity, coordination number, etc.) for each simulation, regardless of the friction values chosen for “wet” compaction phase.

A final set of simulations was run to investigate the effect of a distributed particle size with DEM (simulation set E). Five samples were generated by rendering 10,000 particles in a box with periodic lateral boundaries. The radius of each particle was randomly sampled from a Gaussian distribution with unit mean radius and a prescribed coefficient of variation ν, defined as the standard deviation normalized by the mean value. The five samples were assigned a coefficient of variation of ν={0.0,0.1,0.2,0.3,0.4}. To prevent excessively small or large particles to be created, the Gaussian distribution was clipped at a distance ±2ν from the mean. A value of ν = 0.0 implies a single particle size, which was used as a reference. After the particles had settled by gravity, the sample was subjected to a load of 5.0 MPa and pressure solution was activated. The sample was then compacted until a final porosity of 20% was reached. The stress acting on each particle‐particle contact was outputted at 0.5% intervals of porosity, so that the relative contribution of each particle to the overall compaction rate could be analyzed. Note that the initial porosities varied between samples, owing to differences in particle size distribution. For an unbiased comparison between the various simulations, all macroscopic compaction data are reported as a function of porosity, rather than strain (as is common in many studies).

The details of the parameters and characteristics of all of the simulations are summarized in Table 2.

Table 2

The Parameters and Characteristics as Used in the Various Sets of Simulations

Parameter	A	B	C	D	E	Units
Contact stiffnessa	2.65 × 109	2.65 × 109	2.65 × 109	2.65 × 109	2.65 × 109	N m−1
Friction coefficient	0.6	0.01 or 0.6b	0.01 or 0.6b	0.01–1.0	0.6	‐
Mass density	2100	2100	2100	2100	2100	kg m−3
Radius (min/max)	(1/1.12)	0.5	0.1–1.0	1.0	variablec	m
Size distribution	uniform	monodisperse	monodisperse	monodisperse	Gaussian	‐
Lateral boundary conditionsd	rigid cylinder	periodic	periodic	periodic	periodic	‐
Applied stress	1.08	0.1–1	1	1	5	MPa
Wall stiffness	2.65 × 109	2.65 × 109	2.65 × 109	2.65 × 109	2.65 × 109	N m−1
Wall friction coefficient	0.1	0	0	0	0	‐

Note. See text for a description of the different simulation sets.

aIdentical values were used for both the normal and shear stiffnesses.

bTwo identical simulation sets were run with different friction values.

cParticle radii were sampled from a Gaussian distribution with a mean of 1 m and coefficient of variation ν of 0, 0.1, 0.2, 0.3, and 0.4.

dIn the axial direction, a constant stress boundary condition was applied by a rigid planar wall.

Results of Benchmark Simulations

Comparison With Experiments

To validate the DEM model, we compared the outcome of DEM compaction tests to the results of compaction experiments conducted on porous halite aggregates. Most experimental work on halite has focused on intermediate to long‐term compaction behavior (see e.g., Gratier et al., 2013, for an extensive summary); hence, little data are available for the compaction of aggregates of high porosity, which require a high temporal resolution in the first hour or so of a typical experiment. To serve as a benchmark for validation of the DEM model, a total of six 1 h compaction experiments were conducted on halite aggregates to obtain data with the required temporal resolution. The choice for short‐duration compaction tests was dually motivated: first of all, to simulate a short‐duration test in DEM, slower (scaled) kinetics can be employed to reduce the inertial effects that may affect the initially rapid stages of compaction (see section 3.2). Second, the relatively high porosities would ensure a minimal effect of grain boundary healing on the overall compaction behavior (Van Noort et al., 2008; Zubtsov et al., 2004). An extensive description of the experimental methods is given in the supporting information (Text S1; De Meer & Spiers, 1997; Gratier et al., 2013; Van Noort et al., 2008; Visser et al., 2012; Zhang et al., 2002; Zubtsov et al., 2004).

The strain rate evolution with porosity obtained in the compaction experiments is plotted in Figure 2. The final porosities obtained after approximately 1 h of compaction were similar, reaching values of 30.73 ± 2.02%. Samples with a higher initial porosity generally also showed a higher final porosity. The strain rates produced at a given porosity varied within a band of roughly half an order in magnitude, which we took as a range of the experimental uncertainty. This variability may have been caused by initial (local) variations in porosity, grain size, or packing geometry, which may have manifested more explicitly in the relatively small samples used in this study.

Since the experimental results reported by Spiers et al. (1990) provide a range of values for Z , one simulation was run in simulation set A and the resulting strain rates scaled back according to three different values: one for a low estimate of Z =0.503 × 10−26, one for a high value of 1.51 × 10−26, and one with an average value of 1.01 × 10−26m3 Pa−1 s−1. These estimates are based on the statistical uncertainties in the determination of (D C S)0 and ΔH reported by Spiers et al. (1990), and they reflect the variability of their experimental results. The strain rate versus porosity results of this simulation are directly compared to the experimental data in Figure 2. This shows that end‐member estimates of Z envelope the range of experimental strain rates for porosities in the range of 30–37%. At porosities higher than 37% the DEM model seems to severely overestimate the experimental rates of compaction. This can likely be explained by the fact that in the DEM simulations the initial area of contact between two particles, before initiation of pressure solution, is purely due to elastic compression and is very small owing to the spherical geometry of the particles. In the benchmark experiments, the shape of the grains, being nonspherical, leads to larger initial areas of contact. This results in unrealistically small particle contact areas and high contact stresses in the DEM model, and hence larger strain rates produced in the simulations during initial stages of compaction. Analytical models suffer the same discrepancy when the porosity function tends to infinity. At larger strains, this initial discrepancy vanishes and a much better agreement is found between the DEM simulations and the laboratory experiments. Order of magnitude estimates of the sample permeability reveal that expulsion of the pore fluid (i.e., viscous advection) is very unlikely to impede the rate of compaction at this range of porosities.

When we compare the DEM simulation results to an analytical model (model D3 in Pluymakers & Spiers, 2014) with the same kinetic parameters as used in the mean estimate, we find good agreement between the two. While the analytical model tends to predict slightly higher compaction rates, the evolution of the strain rate with porosity is similar. This suggests that the contact‐scale physics is implemented sufficiently well to capture the macroscopic response that is predicted by the analytical models, which are based on the same physical principles. Hence, we can conclude that the proposed DEM contact law is capable of reproducing the compaction behavior realistically and can be used to investigate more complex granular aggregates for which analytical expressions do not exist.

Stress Dependence

Let us now consider the dependence of the sample strain rate on the applied axial (effective) stress. From theory it is expected that (for low effective stress) the pressure solution creep rate is linearly proportional to the applied effective stress (see equation (1a); Lehner, 1995; Rutter, 1976). This is also observed in compaction experiments (e.g., Spiers et al., 1990), although sometimes a higher stress exponent is observed (see Gratier et al., 2013). It can be shown that equation (6a) produces a stress exponent of +1 for a cubic close‐packed system, in agreement with theoretical values. However, in randomly packed aggregates, each particle contact may exhibit a different orientation and magnitude of the locally acting normal stress, which continuously evolves due to growth of the contact area. The resulting macroscopic compaction rate then represents a bulk‐averaged outcome of all local contributions of contact stress to the strain rate. Because of the internal heterogeneity of stress, a proportional relationship between the stress imposed at the boundaries of the simulation domain and the overall strain rate is not guaranteed for the DEM simulations. A new set of simulations (simulation set B) was therefore run so as to test the stress dependence of the DEM model.

First, a series of simulations was run with varying axial stress (0.10, 0.18, 0.32, 0.56, and 1.0 MPa) and with a near‐zero coefficient of interparticle friction (μ =0.01). The DEM samples were compacted down to a porosity of 30%. At 30.5% porosity, the instantaneous strain rate was measured and plotted against the applied axial stress in Figure 3 (blue dots). The best fit line through the data has a slope of 1.004 ± 0.001, which indicates that , as was expected from theory.

Figure 3

DEM results showing the sensitivity of the overall strain rate to the applied stress, as measured at 30.5% porosity, for a low interparticle coefficient of friction of μ = 0.01 and for a high coefficient of friction of μ = 0.60. Although the simulation with the higher coefficient of friction compacts more slowly (lower strain rates), the stress dependence does not change measurably. The reference value of strain rate is taken at σ =1.0 MPa for μ = 0.60.

Next, a set of simulations were run with the same axial stresses but with a different interparticle friction coefficient of 0.6 (see Figure 3; red triangles). Although the enhanced friction coefficient reduced the rates of compaction by about half an order of magnitude, it did not alter the stress dependence in the model. For easy comparison, the strain rates reported in Figure 3 were all normalized to the corresponding strain rate of the simulation where σ =1.0 MPa and μ =0.6. This result suggests that any deviations of stress exponent from the theoretical value of 1, as observed in experiments, are unlikely to be related to interparticle friction.

Grain Size Dependence

Similar to the stress exponent, it is expected from theory and experiments that diffusion‐controlled IPS creep rates are inversely proportional to the cube of the grain size, that is, (see equation (1b)). With the same reasoning as was put forward for the stress exponent (see section 4.2), equation (6b) does not guarantee a grain size exponent of −3 for disordered aggregates. The procedure for testing the grain size dependence of the DEM model (simulation set C) was similar to that of the stress dependence simulations. A constant axial stress of 1 MPa was applied to each of the five samples with different tested uniform particle sizes (0.10, 0.18, 0.32, 0.56, and 1.0 units of length). The interparticle coefficient of friction was 0.01 in the first set of simulations and was 0.6 in the second set. The log of the strain rate, measured at 30.5% porosity, was normalized to the corresponding strain rate of the simulation, where r = 1 and μ =0.6, and was then plotted against the log of the particle radius in Figure 4. The slope of the best fit line was very close to −3 for both sets of simulations, which indicates that , as was expected from theory, and that the grain size exponent is independent of interparticle friction.

Figure 4

DEM results showing the sensitivity of the overall strain rate to the grain size, as measured at 30.5% porosity. Similar to the stress dependence simulations, two interparticle friction coefficients were tested (μ = 0.01 and μ = 0.60). The interparticle friction does not affect the grain size dependence, but it does decrease the strain rate by a factor of 3, going from a friction coefficient of 0.01 to 0.6. The reference value of strain rate is taken at r = 1.0 for μ = 0.60.

Potential Mechanisms for Strain Rate Reduction

Effect of Interparticle Friction

From the simulation results reported in Figures 3 and 4, it is apparent that an increase in interparticle friction effectively reduces the overall compaction rates, although it does not affect the stress and grain size dependence. The evolution of the porosity with time, as observed from simulation set D, and the strain rates measured at fixed porosities in the range of 28–34% for the different values of μ are plotted in Figure 5. These results demonstrate that the compaction rates are systematically reduced with increasing coefficient of friction, by up to a factor of 3 (half an order of magnitude). The reduction of strain rates with increasing coefficient of friction is also observed in the data plotted in Figures 3 and 4. Furthermore, the effect of friction seems to diminish only slightly with decreasing porosity in the range of investigated porosities. The order of magnitude difference in strain rate ( ) of μ =0.01 and that of μ =1.0 at 34% porosity is 0.510, and at 28% it is 0.472. This implies that the effect of interparticle friction is only weakly dependent on porosity in the range of porosities investigated.

Figure 5

The effect of interparticle friction on the compaction behavior of the DEM simulations. (a) Overview of porosity evolution for interparticle friction coefficients ranging from 0.01 to 1.00. (b) Overall strain rate as a function of the interparticle friction coefficient, as measured at various porosities in the range of 28–34%.

It should be noted that the kinetic parameter (D C S)0 for diffusion‐controlled IPS as determined by Spiers et al. (1990) was obtained by fitting an expression of a similar functional form as equation (1b). Since (1b) does not explicitly take into account any effects of interparticle friction, the value of (D C S)0 is convolved with any contributions of friction to the overall compaction rates. If the contribution of friction is assumed to be independent of porosity, then friction cannot directly explain the tendency of analytical models to overestimate experimental compaction rates, particularly at lower porosities.

Effects of a Distributed Grain Size

As discussed in section 2.2, Niemeijer et al. (2009) report enhanced compaction rates in granular aggregates with a wider grain size distribution. To gain some insights into the effects of a distributed grain size, a number of compaction simulations were run with various grain size distributions (simulation set E), characterized by the coefficient of variation ν. In Figure 6a the overall compaction rate is plotted against bulk sample porosity for each simulation. For porosities >32% the strain rates are clearly enhanced by the distributed grain size: for a given value of porosity, the strain rates increase with increasing width of the distribution. This is consistent with the results of Niemeijer et al. (2009), although the enhancement of the strain rate obtained in the simulations is less than what was observed in the experiments. At around 32% a crossover at which the overall strain rates decrease with increasing ν occurs in our simulations. The decrease in strain rate with increasing ν is a feature that is not observed in experiments but is predicted by the model of Ter Heege et al. (2004).

Figure 6

DEM results for creep tests performed on samples with a fixed mean particle size but varying values of the width of the grain size distribution, expressed in terms of the coefficient of variation, ν. (a) The evolution of macroscopic strain rate with porosity. During the initial stages of compaction, where the porosities are high (>32%), the strain rates are enhanced for higher values of ν, whereas during the final stages of compaction the compaction rates are diminished with increasing ν. (b) Evolution of contact parameter with porosity. A similar crossover is found in as in but at a lower value of porosity.

It would be instructive to first consider the microscopic state of stress, that is, the stress supported by individual particle contacts. To this end, we define the average particle stress as the total vertical stress of all contacts associated with a given particle, divided by its coordination number, that is where N is the coordination number of particle j, and is the is the vertical component of the normal stress (normal force F divided by contact area A ) exerted by particle i on j. Then, a probability density function P(x) can be defined that describes the distribution of over all the particles such that

Note that the probability density, unlike the absolute probability, is allowed to locally exceed 1. For convenience, we define , with being the average particle stress measured at a bulk sample porosity of 35%. Then, P(x) is calculated at intervals of 5% of absolute porosity as to obtain the evolution of the distribution of stress with increasing strain. The results of this analysis are displayed in Figure 7, which clearly show that with decreasing porosity (increasing strain) (1) the average stress supported by the particles decreases (curves shifting left over the x axis), and (2) the stress becomes more homogeneously distributed, as indicated by a significant local increase in P(x). The first observation can be directly related to the notion that the average contact area A grows with ongoing pressure solution, as discussed in section 3.1. In continuum descriptions of IPS (e.g., equation (1a)), this evolution of the contact area is represented by an arbitrarily defined porosity function, which likely does not capture the detailed evolution of stress as seen in Figure 7.

Figure 7

Probability density functions of interparticle stress measured at four different bulk porosities (see legend) and for four different particle size distributions (individual panels). With decreasing porosity, both the magnitude and heterogeneity of the particle stress decrease.

The second observation that the stress distribution becomes more homogeneous with increasing strain is readily explained by the negative feedback loop that is induced by equation (6a): when the stress on a given contact is higher than average, it will experience enhanced rates of pressure solution (high V ps), and so the stress on the contact will decrease more rapidly as A grows, in turn resulting in lower V ps. This also implies that force chains initially present in the system will broaden and dilute as compaction progresses, and that it is unlikely that such a system supports localization of strain.

With this additional knowledge of the internal distribution of stress, we can now attempt to relate this observed crossover in strain rates to the internal stress distribution. It is expected that the distribution of contact stresses within the sample is related to the distribution of the particle sizes (e.g., Marketos & Bolton, 2007; Torok et al., 2005). However, it is not directly evident from experiments how the average contact stress changes with the width of the particle size distribution, and how it relates to the overall compaction behavior. To evaluate this, we calculate the axial (vertical) component of the normal force (F ) acting on each particle‐particle contact, as well as the area of contact (A ), and record these quantities at regular intervals of 0.5% porosity during the simulations. The bulk‐averaged, axial contact stress is then obtained as

Most analytical expressions for IPS only consider a uniform (mean) grain size, and therefore, the contribution of the diffusion distance to the strain rate can be directly related to the grain size and total strain (or porosity). The predicted strain rates for a given grain size and porosity are then linearly proportional to the (applied) normal stress (see equation (1b)). However, it is nontrivial to achieve such relations for a nonuniform grain size, and therefore, we choose to explicitly include the evolution of the diffusion distance in our analysis by defining a new quantity, , following equation (6b):

The parameter describes the average properties of the particle contacts that drive pressure solution and is proportional to the macroscopic strain rate, that is, . The evolution of with porosity is plotted in Figure 6b and shows a similar crossover as the overall strain rate. This suggests that the change in compaction behavior in our DEM model is related to changes in the internal distribution of stress and contact area. However, the porosity at which this crossover of occurs is about 5% lower (in absolute value) than that of . This implies that for a given , decreases with increasing width of the particle size distribution. We were not able to elucidate this behavior with our current analysis, and it is possibly a consequence of the geometrical arrangement of particles within the sample. Yet it can be concluded that the evolution of the internal state of stress and contact area in a polydisperse aggregate may in part explain the change in strain rate dependence on ν, helping to reconcile the experimental observations of Niemeijer et al. (2009) to the analytical model of Ter Heege et al. (2004).

While there is a clear effect of grain size distribution on the overall rate of compaction, this effect is generally less than half an order of magnitude for the largest value of ν investigated. Hence, polydispersity on its own cannot directly explain the large discrepancies between analytical predictions and experimentally observed compaction rates. However, it should be noted that in this study only a Gaussian distribution was considered, and that other types of grain size distributions (e.g., power law or bimodal) may produce a different distribution of contact stresses and strain rates, for instance, when large particles are fully immersed in a matrix of small particles. Additionally, the shape of the grains and that of the grain‐grain contacts could strongly affect the magnitude and distribution of the contact stresses. We reserve further investigations of the various effect of grain size distributions and grain shapes for future studies.

Conclusions

In this study, we have proposed a new DEM contact law and constructed a DEM model accounting for friction and pressure solution at the grain contacts, using a pressure solution contact law extracted from previous macroscale expressions for IPS. The novelty of the implementation lies in that the DEM particle overlap is used to model dissolved material, of which the dissolution rate is directly linked to microscale physical models. By comparison of the DEM results with uniaxial compaction experiments conducted on granular halite with a narrow grain size distribution, we found that the DEM model is capable of realistically describing the macroscopic compaction behavior as seen in the experiments. The verified numerical model was then used to investigate the effects of intergrain friction and of nonuniform grain sizes on the creep rates of granular halite.

It was found that interparticle friction can reduce the compaction rates by up to half an order of magnitude. The amount of reduction of compaction rates was independent of porosity and absolute strain rate. Furthermore, interparticle friction did not seem to affect the dependency of strain rate on the grain size and stress, which were identical for simulations with an interparticle friction coefficient of 0.01 and 0.6. However, it is likely that the estimates of (D C S)0 (a DEM input parameter), as obtained from experiments, are already affected by any effects of friction; and hence, friction itself cannot explain the discrepancy between analytical model estimates and experimental compaction rates, which can exceed several orders of magnitude difference.

When considering the effect of a distributed grain size, it was found that with increasing degree of polydispersity (i.e., wider grain size distribution) the bulk‐average value of ω (i.e., contact stress over contact area) increased, promoting overall compaction rates. This was in agreement with the results of earlier experiments on halite with a distributed grain size by Niemeijer et al. (2009), who observed enhanced strain rates in more polydispersed samples. However, as the porosity decreased below 32%, decreased faster with increasing ν and the effect reversed, that is, with increasing degree of polydispersity the average contact stress and strain rate decreased. This is consistent with the analysis of Ter Heege et al. (2004), who derived a composite flow law accounting for a lognormally distributed grain size and predicted a reduction in strain rates with increasing polydispersity for bulk (low porosity) materials. However, this effect of the grain size distribution was small (less than half an order of magnitude) and is likely to play only a minor role in the overall compaction behavior in experiments, when a lognormal distribution is considered. In these simulations it was also observed that with increasing compaction, the contact stresses become more narrowly (i.e., more homogeneously) distributed over the aggregate by preferentially dissolving material at the more highly stressed particle contacts. This can have important implications for other phenomena that have been linked to IPS, such as shear deformation of granular fault gouges under conditions where fluid‐rock interactions are prevalent.

Another conclusion is that interparticle friction and a lognormally distributed grain size cannot explain the observed discrepancies between experimental compaction rates and predictions by analytical models. Hence, other mechanisms that may impede compaction by pressure solution need to be considered, particularly for understanding the long‐term evolution of large‐scale geodynamical systems, such as compacting sedimentary basins and hydrocarbon reservoirs. A possible candidate mechanism is the process of sealing of the grain contact by grain boundary healing due to growth of contact asperities (see Van Noort et al., 2008). It is left for future studies to investigate the role of this mechanism on the overall deformation rates. Nonetheless, we have demonstrated that the discrete element method can successfully be used to simulate pressure solution in granular aggregates, and that the model can be used to relate microscale processes to macroscopic compaction behavior.

Supporting Information S1

Click here for additional data file.