## Abstract

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.

## INTRODUCTION

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 (*1*–*3*) or dissolve in fluid over time (*4*–*6*). 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.

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 (*8*–*10*). 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 (*18*–*20*).

## RESULTS

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.

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 *l*_{m} 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, *l*_{m} also controls the upward capillary pull up of fluid through the capillary gravity length (*23*)(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” *S*_{e} [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 *S*_{e}(*z*) = exp(−*z*/*L*_{G}), 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 *S*_{e}(*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).

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 *K*_{0}; 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 *l*_{cr} 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 Δ*t*_{cr} linearly grows over time from the point of fluid injection *t* as (see Eq. 13)(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 *l*_{m} = 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 *l*_{m}, 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 *K*_{0} = κ*p** with κ = 470 previously established using *P*-wave propagation velocity measurements (*10*) with piezoelectric sensors. The remaining property was the critical fluid activity *a*_{cr} = 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)(3)as depicted by the slope in Fig. 3E. Again, just like λ, the logarithmic creep rate *C*_{l} is also only a function of the gravity constant *g* and material constants (ρ, γ, *l*_{m}, κ, and *a*_{cr}).

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.

## DISCUSSION

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 *l*_{m} and the capillary gravity length *L*_{G}. Since *l*_{m} ≪ *L*_{G}, Eq. 2 suggests , 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 *l*_{m} expected in the ice. Therefore, through upscaling, our model may also apply for geological pressures and durations in crustal rocks (*18*, *26*, *28*–*30*) 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.

## MATERIALS AND METHODS

### Model

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 *i*^{th} element pressure [*p*^{(i)}] and its elastic strain [ε^{(i)}, positive in compression] is defined as(4)where *K*^{(i)} is the *i*^{th} constrained elastic modulus, with *K*_{0} being its initial value before any interaction with fluid; *l*^{(i)} is the compressed length of the *i*^{th} carriage; and 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 *l*_{m}, while the stiffness degrades with fluid activity (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 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*) [where *A*, *E*_{a}, *R*, and *T*^{(i)} are the activation rate constant, activation energy, ideal gas constant, and absolute temperature of the *i*^{th} carriage, respectively]. Under the constant room temperature (≈20°C) maintained in our experiments, α is therefore constant in both space and time and thus(6)

Accordingly, for fixed degrees of saturation in time, one should expect a linear relationship between *a*_{S} and *S*_{e}, which is also consistent with experimental observations with cereals (*33*).

*Upon and after crushing*. The equilibrium length changes to a new smaller constant value *l*_{cr}. The crushing of the *i*^{th} carriage occurs when its fluid activity reaches a critical value . At this point, the constrained modulus attains a constant fully degraded value(7)where denotes the critical time for the crushing of the *i*^{th} carriage. Under constant temperature, its value can be found from the following implicit relation(8)

Accordingly, the onset of crushing of the *i*^{th} carriage depends on how its effective saturation 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 *i*^{th} 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 *l*_{cr} can be related to the other constants by combining Eqs. 4 and 8 [for *p*^{(i)} = *p**](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 [*S*_{e}(0) = 1]. From typical profiles reported in the literature, *S*_{e} 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*)(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 *S*_{e}(*z*) = exp(−*z*/*L*_{G}), 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 . 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 . According to Valdes *et al*. (*10*), for the order of the *p** ≈ 50 kPa in our experiments, *K*_{0} ≈ 15 MPa from shear wave measurements. This gives , and thus, the initial positions of the carriages can be calculated as . The initial effective saturation of carriage *i* is therefore(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 find(12)such that, with the help of Eq. 10, the time gap between two consecutive quakes grows linearly with time as(13)

### Deformation

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. 12(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 *N* − *n* uncrushed ones(15)

Since the deformation in each crushed carriages *d*_{i} is independent of time, the overall contribution to the deformation from all the *n* crushed carriages is given as(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)} = *l*_{m}(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 . 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 (*N* − *n*) uncrushed carriages to the overall deformation(17)such that by combining Eqs. 9, 15, and 17, and since the initial height of the pack is simply *H* = *Nl*_{m} (with *H* ≫ *l*_{m}), the total deformation at time *t* after the soaking event could be calculated as(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 *l*_{m} 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 *l*_{m} = *dz* as the infinitesimal integration width only within that discontinuous summation part of the equation. Out of this summation, we keep *l*_{m} to be finite and take *H* = *l*_{m}*N* as the sample height. In addition, with the help of Eq. 12, the transition integer *n* and transition level *z*_{n+1} can be replaced by the smooth continuous time as and (such that at *t* = 0, we have both *n* = 0 and *z*_{1} = 0, with as a nondimensional continuous time). Using this step, the approximated smooth analytic expression for the overall deformation *D*(*t*) takes the form(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 time(20)

We then define the effective limiting creep rate *C*_{l} 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 )

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 *K*_{0} = κ*p**, and thus(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 *i*^{th} carriage, which follows a biphase value dependent on the Heaviside function *H*(*x*) = 1 if *x* > 0, *H*(*x*) = 0 otherwise(22)where, at the start of the simulation, . In removing the small strain assumption, we must now recall that, unlike the small strain analytic model before, the effective saturation corresponding to the *i*^{th} carriage generally varies with time. This is because the position of that carriage varies as a function of carriage lengths *l*^{(i)}, itself varying with time. Specifically, during the simulation, the *i*^{th} carriage effective degree of saturation and length are(23)with the initial condition for the lengths given by . 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 Δ*t* ≪ *a*_{cr}/α the characteristic time scale of the system. Simulations are nondimensionalized [with *l*_{m} the unit length, Δ*t* the unit time, and *p**/(*l*_{m}Δ*t*^{2}) the unit mass], but results are presented with dimensions for easier comparison with experiments.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/4/10/eaat6961/DC1

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.

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.

## REFERENCES AND NOTES

**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.

- Copyright © 2018 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).