Tracking time with ricequakes in partially soaked brittle porous media

See allHide authors and affiliations

Science Advances  12 Oct 2018:
Vol. 4, no. 10, eaat6961
DOI: 10.1126/sciadv.aat6961


When brittle porous media interact with chemically active fluids, they may suddenly crumble. This has reportedly triggered the collapse of rockfill dams, sinkholes, and ice shelves. To study this problem, we use a surrogate experiment for the effect of fluid on rocks and ice involving a column of puffed rice partially soaked in a reservoir of liquid under constant pressure. We disclose localized crushing collapse in the unsaturated region that produces incremental global compaction and loud audible beats. These “ricequakes” repeat perpetually during the experiments and propagate upward through the material. The delay time between consecutive quakes grows linearly with time and is accompanied by creep motion. All those new observations can be explained using a simple chemomechanical model of capillary-driven crushing steps progressing through the micropores.


When it comes to pouring milk into a bowl of cereal (Fig. 1A), there is more to it than meets the eyes. Just like cereal, it is commonly known that crustal rocks and dry snow can crumble under pressure (13) or dissolve in fluid over time (46). But what happens when brittle porous media are both subjected to high pressures and soaked at the same time? Can hidden failure patterns develop in crustal rocks under experimentally inaccessible high pressures and long time scales (i.e., up to millions of years) (2)? To answer these questions, we study the infiltration of milk into cereals that are constantly pressured. As a result, we reveal “ricequakes,” the gradual jerky collapse of brittle porous media in chemically active fluids.

Fig. 1 The problem of partially soaked brittle porous media.

(A) Breakfast cereal. (B) Sketch of rockfill dam collapse under capillary action [inspired from (13)]. (C) Schematic of the experiment. (D) Stress over time since soaking the puffed rice revealing recurring stress drops. (E) Bursts of acoustic emission correlated with the stress drops. a.u., arbitrary units.

The substitution of subject materials by surrogate materials has been widely effective in science. Specifically, dry cornstarch, carbomer, and silicone have been used to demonstrate transient deformations in correspondingly brittle, semibrittle, and ductile crustal rocks (7), while cereal (puffed rice) has been used to reveal previously hidden propagating compaction patterns in dry porous media (810). Following these first observations with puffed rice, the compaction of dry snow also disclosed similar patterns (11)—this is likely because rocks, dry snow, and cereals all share the properties of brittle porous media. The testing of puffed cereal grains, however, is simpler as they are highly porous, brittle, and soft, which renders their compaction feasible within limited spaces and short time scales suitable for laboratory experiments. It is therefore likely that the ricequakes we find in partially soaked cereals may also develop in wet crustal rocks under experimentally inaccessible geological conditions.

This paper begins by detailing the experiments and the spectacular properties of ricequakes. We then explain the mechanisms controlling the ricequakes using a simple chemomechanical model, with which we discuss whether ricequakes may also be related to the collapse of rockfill dams (12, 13), sinkholes (14, 15), and ice shelves (4, 16). For example, (12, 13) reported compaction failures of rockfill dams caused by upward moisture migration through capillary suction in Fig. 1B. It is understood that wetting-induced collapse of such media depends on interfacial energy changes (17) and chemical dissolution (1820).


In the interest of further detecting and separating the microscopic mechanisms controlling the problem from global geometric effects, we simplified the geometry of the rockfill dam problem using a vertical cylinder filled with puffed rice overlaying on top of a specially designed granular filter (Fig. 1C). We first increased the load on top of the puffed rice pack using a servo-controlled loading machine programmed to achieve a designated stress. Once the stress was closely constant p* and after the deformations essentially ceased, we injected a set amount of fluid through a tube connected to the bottom of the granular filter, which soaked the lower part of the puffed rice pack and immediately consolidated the sample. Consequently, this has expectedly reduced the stress temporarily below its target constant value, to be regained after ≈250 s (Fig. 1D).

We also noticed that after fluid injection, the stress kept on dropping repetitively, while the delay between these drops prolonged over time. Along with the stress drops, the pack also emitted audible sounds, which resembled the clicking beats of a metronome slowing down over time (listen to audio file S1 and see the acoustic bursts in Fig. 1E). These beats were loud and required no amplification to be noted all around the laboratory.

To discern the principal mechanism controlling this phenomenon, we used a Hirox KH-8700 digital microscope and recorded the deformation in the small interfacial region between the fully soaked and dry parts of the material, as seen in Fig. 2A (and in the corresponding movie S1). To analyze the local vertical deformations across that region, we formed a spatiotemporal plot (Fig. 2B) of the fabric along a central vertical line in Fig. 2A (here focusing on the time after the sample had established the recurring ricequakes). As revealed by this figure, while the upper unsaturated pack significantly compacted over time, it collapsed almost as a rigid body motion into the lower fully saturated zone. This occurred while the interface between the fully soaked and unsaturated parts of the pack remained closely stationary. The downward motion of the upper part of the pack shows recurring incremental drops interluded by creeping compaction periods. These incremental collapses synchronize perfectly with the instances of stress drops and acoustic bursts, thus all associated to the so-called ricequakes. The delay times between consecutive quakes can thus be unequivocally measured, as plotted in Fig. 3A against time from fluid injection. As seen from this figure, the delay times increase mostly linearly with time.

Fig. 2 Deformation patterns.

(A) Macrophotography of an area surrounding the interface between the soaked and unsaturated parts. (B) Spatiotemporal plot of the material texture under the red line in (A), which shows recurring incremental compactions during ricequakes interluded by creep phases. The red dashed line in (B) highlights the position of the interface.

Fig. 3 Tracking time with ricequakes.

(A) Delays between consecutive ricequakes against time since injection, with the linear line predicted by Eq. 2 with typical micropore unit size lm = 0.6 mm. (B) Frequency distribution of micropore sizes extracted from several puffed rice cross sections, with a vertical line highlighting the typical micropore unit size lm = 0.6 mm. (C) Macrophotography of a puffed rice sliced section. (D) Overall deformation with time since injection during the experiment and based on the numerical and analytical models. (E) Magnified view of the deformation that highlights the ricequakes and the creep phase in between.

We note from the spatiotemporal plot in Fig. 2B that the displacement drops associated with the ricequakes are much smaller than the size of a single puffed rice grain. Therefore, to identify the microscopic feature potentially related to these drops, we took more than 15 slices of puffed rice grains and imaged them using a digital microscope, as shown in Fig. 3C. Using these images, we manually measured over 300 micropore unit sizes, as given by the histogram in Fig. 3B, showing pore unit sizes mostly under a millimeter. Therefore, a pack of puffed rice has two characteristic sizes: a macropore size defined by the voids between the grains (of the order of a typical centimeter grain) and a micropore unit size. This double porosity is common for geomaterials, such as aggregate packs of rockfill dams (17), whose hydromechanical collapses motivated this study. Since the displacement drops are much smaller than a typical grain, we assume that the micropores control the crushability of puffed rice. This is consistent with the idea that, under fluid activity, the puffed rice collapses into its gradually wetting micropores.

To explain the experimental observations, we developed a simple “crushing wave model” that builds on the leading role of brittle micropores (see Materials and methods and Fig. 4). In this model, a train of elasto-brittle carriages (see Fig. 4A) with initial spring equilibrium length given by the typical micropore unit size lm can crush under the action of partial fluid saturation. As is known from the physics of capillarity (21, 22) under gravity, the upward wetting of porous media is carried out by micropores that reside within the particles (as they undertake larger surface tension), and thus, lm also controls the upward capillary pull up of fluid through the capillary gravity length (23)Embedded Image(1)where g is the gravity acceleration, ρ is the fluid density, and γ is the surface tension.

This length is strongly connected to the saturation profile. For double porosity media, this profile is best captured by the “effective saturation” Se [or normalized water content (22, 24)], which reflects the dominant role of micropores within the particles, and the insignificant contribution of the interparticle macropores to capillary rise (17, 25). For the same reason, most of the saturation is contained within the grains, which is thus not visible from the outside. However, at the end of the test, after carefully removing the upper particles, we clearly identified internally wet particles in the centimeter above the visible wetting front level. The effective saturation profile is described generally for hydrophilic materials by Se(z) = exp(−z/LG), with z being the distance from the wetting front upward (see Fig. 4B). A usual way to establish this profile is by using water retention curves (22), but this requires equilibrium which cannot be maintained for puffed rice, as they quickly degrade under fluid. However, the use of the exponential profile for Se(z) is strongly supported by comparison with other common highly porous particles used in agriculture, horticulture, and the building industry (see supplementary text and fig. S1).

Fig. 4 Crushing wave model.

(A) Micropore units as train carriages. (B) Effective saturation profile Se(z) = exp(−z/LG). (C) Stress-displacement relationship for individual carriages. The red line illustrates the deformation path under fluid activity, and the red triangle indicates the compressed length after crushing.

As movie S1 and Fig. 2B suggest, the decay in effective saturation could be taken constant in time with respect to the initial position of the visible wetting front, as it is mostly stationary. Since under this partial effective saturation the solid skeleton chemically degrades, both its constrained stiffness modulus (initially K0; see Materials and methods) and critical crushing strength reduce over time (see Fig. 4C), which explains the experimentally observed creep phases. Furthermore, it is assumed that once the critical crushing strength reduces below the stress a micropore unit carries, it will sustain a wetting-induced collapse and the said unit will establish a new and reduced equilibrium spring length lcr given by Eq. 9. After each micropore unit collapses, the next micropore unit in line is targeted. These events add abrupt increments to the overall deformation pattern, as is evidenced by the experimental, numerical, and analytical lines in Fig. 3 (D and E). The numerical and analytical lines (derived in Materials and methods) are based on models with respectively fixed degrees of saturation in either the laboratory or material frames. According to our model, the delay between consecutive quakes Δtcr linearly grows over time from the point of fluid injection t as (see Eq. 13)Embedded Image(2)which can be evaluated against the experimental data in Fig. 3A. For the water used, we have γ = 73 × 10−3 kg s−2 and ρ = 1000 kg m−3. As for the micropore unit size, we use lm = 0.6 mm, being roughly the mode of the size distribution in Fig. 3B. In that sense, our model explains the delay time between ricequakes, without any fitting parameter. As we have all the parameters to calculate our theoretical λ, the fact that it matches the experimental λ is quite remarkable. Testing other values of λ would necessitate different lm, requiring a different material. Similarly, changing the fluid can slightly modify ρ and γ, but also strongly alter the chemistry. The only parameter whose value could be varied significantly without changing the solid and fluid is g, by using a centrifuge, which would be interesting to explore in the future.

As seen in Fig. 3 (D and E), our model also successfully predicts the overall deformation over time, while revealing the two consolidation mechanisms of ricequake and creep compaction described above. To establish this comparison, we used the puffed rice constrained modulus K0 = κp* with κ = 470 previously established using P-wave propagation velocity measurements (10) with piezoelectric sensors. The remaining property was the critical fluid activity acr = 0.96 that we determined as a free fitting parameter. The theoretical model can be further simplified by considering the delay times as arbitrarily small integration time steps, leading to the additional smooth line in Fig. 3 (D and E) (see Materials and methods). An effective creep rate is then derived (see Eq. 20) for the limit of sufficiently large times prior to complete collapse of all the micropores (see Eq. 21)Embedded Image(3)as depicted by the slope in Fig. 3E. Again, just like λ, the logarithmic creep rate Cl is also only a function of the gravity constant g and material constants (ρ, γ, lm, κ, and acr).

Although often small, dynamic effects are always experimentally inevitable and go beyond quasi-static models. To depict more clearly the dynamic stress signatures between all the events, we replot the stress against the log of time (Fig. 5A). As shown, after each drop, the stress returned to its target value much faster than the delay time between the events. This means that new events were likely not affected by the servo-controller recovery time, which strengthens the use of a quasi-static model to explain the delay time statistics. Nevertheless, it did probably affect the deformation slightly. We note that the stress drops increased from about 2 to 5% of the desired stress. This inevitable dynamic effect was translated into a growth in the incremental ricequake displacements over time, from around 50 to 80 μm (see Fig. 5B). As asserted by our model, the extent of creep increments between those quakes remained practically constant (see Fig. 5B), thanks to the short response time of the servo relative to the creep time between the events. These observations therefore reinforce our model assumptions yet motivate future extensions to include dynamic effects.

Fig. 5 Dynamic effects.

(A) Stress drops over logarithm of time, which highlights the fast stress recovery after each ricequake all throughout the test. (B) Incremental displacements, which highlights the constancy of the creep contribution between the events.


In conclusion, a new form of incremental collapse was detected in partially soaked brittle porous media using puffed rice, which produces acoustic bursts and stress drops. These ricequakes start shortly after the puffed rice is soaked, repeat perpetually, and the delay time between consecutive quakes prolongs linearly over time. This linear relationship is so robust that one may use the ricequakes for tracking time.

We showed that the growth rate of the delay times depends on the ratio between the micropore unit size lm and the capillary gravity length LG. Since lmLG, Eq. 2 suggests Embedded Image, and thus, the time between consecutive quakes scales quadratically with the micropore unit size. This means that delays between quakes would be hardly distinguishable for solids with much smaller micropores, such as soft coastal carbonates (26) and Antarctic ice sheets (4). However, this does not exclude the effect of the gradual collapse of micropores on the creep rate in these materials. The time signatures of the instabilities in our experiments are reminiscent of observations of tidal icequakes from the Whillans Ice Stream, Antarctica (27), which exhibit an astonishingly similar pattern of incremental displacements as in Fig. 3 (D and E). The slowdown in the waiting time between icequakes happens over many years and is thus much milder than the slowdown of ricequakes, which may be reasoned by the much smaller lm expected in the ice. Therefore, through upscaling, our model may also apply for geological pressures and durations in crustal rocks (18, 26, 2830) and ice sheets (4, 16, 27). That said, before extrapolating data from cereal to predict wetting-induced phenomena and creep rates in other materials, the validity of Eq. 3 should be tested by conducting corresponding experiments under different gravity g and a wider range of brittle porous solids and chemically active fluids.



In the following, the crushing wave model is developed. The model is structured as a train of elasto-brittle carriages from i = 1 to N akin to a vertical series of micropores extending from the wetting surface upward, with i = 1 being the first carriage initially positioned adjacent to the wetting surface and i = N being the last carriage, which touches a loading plate applying a pressure p* that is initially positioned at height H (Fig. 4A). Neglecting inertia effects and any friction between the tested material and the container boundaries, force equilibrium imposes the pressure p(i) = p* for any i. Irrespectively, the relation between the ith element pressure [p(i)] and its elastic strain [ε(i), positive in compression] is defined asEmbedded Image(4)where K(i) is the ith constrained elastic modulus, with K0 being its initial value before any interaction with fluid; l(i) is the compressed length of the ith carriage; and Embedded Image is its corresponding equilibrium length. In the following, the equilibrium length will have two possible values for two distinct behavioral phases of the material: prior and upon and after crushing (Fig. 4C).

Prior to crushing. The equilibrium length is given by the uncrushed characteristic micropore length lm, while the stiffness degrades with fluid activity Embedded ImageEmbedded Image(5)where the fluid activity is defined in terms of the time (t) integrated effect of the chemical reaction of the porous medium with the fluid, i.e., in terms of the effective degree of saturation Embedded Image of node (i) at time t, and the reaction rate α. The assumed linear decay of the stiffness with water activity agrees with previous experiments on crunchy material for nonnegligible degrees of saturation (31). The reaction rate can be expressed most generally as an Arrhenius’ relationship (32) Embedded Image [where A, Ea, R, and T(i) are the activation rate constant, activation energy, ideal gas constant, and absolute temperature of the ith carriage, respectively]. Under the constant room temperature (≈20°C) maintained in our experiments, α is therefore constant in both space and time and thusEmbedded Image(6)

Accordingly, for fixed degrees of saturation in time, one should expect a linear relationship between aS and Se, which is also consistent with experimental observations with cereals (33).

Upon and after crushing. The equilibrium length changes to a new smaller constant value lcr. The crushing of the ith carriage occurs when its fluid activity reaches a critical value Embedded Image. At this point, the constrained modulus attains a constant fully degraded valueEmbedded Image(7)where Embedded Image denotes the critical time for the crushing of the ith carriage. Under constant temperature, its value can be found from the following implicit relationEmbedded Image(8)

Accordingly, the onset of crushing of the ith carriage depends on how its effective saturation Embedded Image evolves with time. Experimentally, however, it is noticed that in between ricequakes, the capillary wetting front remains quite static, whereas immediately after a quake, its motion establishes a new equilibrium position fairly quickly (at the order of a second or less). Therefore, in the interest of providing the simplest explanation for the recurring ricequakes, it could be safely assumed that the saturation is fixed in space relative to the initial position of the wetting front, while the solid points consolidate, and thus, the effective saturation in the ith carriage is dependent on that carriage distance from the wetting front, which slightly shortens over time. In addition, note that the crushing length of the carriages lcr can be related to the other constants by combining Eqs. 4 and 8 [for p(i) = p*]Embedded Image(9)

Effective saturation profile

The equilibrated shape of wetting fronts in porous media partially soaked in fluid under gravity depends on many parameters (23), from the level of gravity to fluid viscosity, topology of the pores, and their connectivity in the solid medium. We take z = 0 as the lowest point of the wetting front, where the effective saturation is [Se(0) = 1]. From typical profiles reported in the literature, Se decays over distance z away from the lowest point of the front (34, 35). The typical length of decay is frequently related directly to the capillary gravity length (23)Embedded Image(10)where g is the gravity, ρ is the fluid density, and γ is the surface tension. A natural way to consider the saturation decay and take this length into account is to set the effective saturation profile as Se(z) = exp(−z/LG), assuming the puffed rice is a hydrophilic material (Fig. 4B). See supplementary text and fig. S1 for illustrations of the relevance of this profile for systems of highly porous particles. The carriage i whose top is at location z(i) therefore experiences effective saturation Embedded Image. Under the pressure level p* at the start of the wetting part of the experiments, we find the initial position of the center of the carriages Embedded Image. According to Valdes et al. (10), for the order of the p* ≈ 50 kPa in our experiments, K0 ≈ 15 MPa from shear wave measurements. This gives Embedded Image, and thus, the initial positions of the carriages can be calculated as Embedded Image. The initial effective saturation of carriage i is thereforeEmbedded Image(11)

Critical crushing times

Assuming small solid deformations, it is possible to resolve analytically the critical crushing times over time. Under the small deformation assumption, the initial effective saturation in a carriage i remains constant and attached to that carriage and, thus, does not depend on time. Therefore, using Eqs. 8 and 11, we findEmbedded Image(12)such that, with the help of Eq. 10, the time gap between two consecutive quakes grows linearly with time asEmbedded Image(13)


Assuming again small solid deformations, it is also possible to resolve analytically the overall displacement over time. At any given time t, the last carriage n experiencing crushing could be identified with the help of Eq. 12Embedded Image(14)while during that time t, the total deformation D(t) (with t = 0 denoting the point of soaking the medium) can be calculated by summing up the contribution from all the n crushed carriages and the Nn uncrushed onesEmbedded Image(15)

Since the deformation in each crushed carriages di is independent of time, the overall contribution to the deformation from all the n crushed carriages is given asEmbedded Image(16)

On the other hand, the length at time t of an uncrushed carriage i > n can be calculated from Eq. 4 as l(n) = lm(1 − p*/K(n)), where K(n) is given by Eq. 5 in terms of the integration of its corresponding effective saturation. The deformation from the time of soaking in that uncrushed carriage is then Embedded Image. Under the assumption of small deformation, which implies that each of the carriages carries a constant effective saturation that is only a function of the initial position through Eq. 11, we may therefore approximate the contribution from all the (Nn) uncrushed carriages to the overall deformationEmbedded Image(17)such that by combining Eqs. 9, 15, and 17, and since the initial height of the pack is simply H = Nlm (with Hlm), the total deformation at time t after the soaking event could be calculated asEmbedded Image(18)while n grows discontinuously according to Eq. 14, with n = 0 at time t = 0.

The analytic expression of the deformation over time in Eq. 18 is discontinuous due to its dependence on the typical size of the micropore units lm and the integer n marking the transition level between crushed and uncrushed carriages. The overall deformation curve is given by two-phase response modes, including a discrete crushing phase with n and D jumping distinctly, and a smooth creep phase. It is possible to smoothen the discrete part of the deformation by replacing the sum operation with an integral equation. This could be done by taking lm = dz as the infinitesimal integration width only within that discontinuous summation part of the equation. Out of this summation, we keep lm to be finite and take H = lmN as the sample height. In addition, with the help of Eq. 12, the transition integer n and transition level zn+1 can be replaced by the smooth continuous time Embedded Image as Embedded Image and Embedded Image (such that at t = 0, we have both n = 0 and z1 = 0, with Embedded Image as a nondimensional continuous time). Using this step, the approximated smooth analytic expression for the overall deformation D(t) takes the formEmbedded Image(19)

Creep rate

The advantage of the smooth solution in Eq. 19 is that it can be used to establish an effective creep rate as the change of deformation under constant stress p* over the change in the logarithm of timeEmbedded Image(20)

We then define the effective limiting creep rate Cl for when the crushing cycles are well established over sufficiently long times 1 ≪ αt (yet before complete annihilation of all micropores in the column, i.e., for Embedded Image)Embedded Image

Notice that while α does not affect the asymptotic creep rate in Eq. 21, it still controls the initial creep rate (when αt ~ 1), as could be seen in Fig. 3. Following our previous work (10), we found that the stiffness of puffed rice is linearly pressure-dependent with K0 = κp*, and thusEmbedded Image(21)

Numerical model

The small strain assumption, which implies constant fluid saturation per carriage, can be removed by integrating the model numerically. For that purpose, we started by calculating the equilibrium length of the ith carriage, which follows a biphase value dependent on the Heaviside function H(x) = 1 if x > 0, H(x) = 0 otherwiseEmbedded Image(22)where, at the start of the simulation, Embedded Image. In removing the small strain assumption, we must now recall that, unlike the small strain analytic model before, the effective saturation corresponding to the ith carriage Embedded Image generally varies with time. This is because the position of that carriage Embedded Image varies as a function of carriage lengths l(i), itself varying with time. Specifically, during the simulation, the ith carriage effective degree of saturation and length areEmbedded Image(23)with the initial condition for the lengths given by Embedded Image. The constrained elastic moduli K(i) and the lengths l(i) are updated using a fourth-order Runge-Kutta algorithm, with a time step Δt such as Δtacr/α the characteristic time scale of the system. Simulations are nondimensionalized [with lm the unit length, Δt the unit time, and p*/(lmΔt2) the unit mass], but results are presented with dimensions for easier comparison with experiments.


Supplementary material for this article is available at

Movie S1. Video of milk injection into a puffed rice column under constant pressure, focusing on the region near the saturated/unsaturated interface (the video is 15 times faster than real time).

Audio file S1. Two-minute recording of the sound emitted during the same experiment (in real time), starting 7 minutes after fluid injection.

Supplementary Text

Fig. S1. Description of previous experiments from the literature, which reinforces the use of exponential decaying effective degree of saturation for systems with highly porous particles.

References (3638)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank J. R. Valdès for introducing us to the world of cereals. Funding: This work was supported by the Australian Research Council through project DP160104310. Author contributions: Both authors conceived the ideas, carried out the experiments, developed the model, and wrote the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article