Poroelastic Modelling of Gravitational Compaction

A dynamic in-silico model captures the kinetics of 1-d gravity driven instabilities, in gravity or centrifuge, of fluid-infiltrated poroelastic media in a partial differential equation (pde). The pde yields the porosity profile over height and time for the given initial and boundary conditions, during slow compaction in counter-current fluid drainage. Processes captured are amongst others sedimentation, creaming and subsidence. The most important limiting prerequisite is that the incompressible dispersed medium is sufficiently structured and/or concentrated that it compacts during slow drainage, without segregation in sizes or in components. For Unilever, modeling of gravitational instability of products is important to quantify or extrapolate long time behavior during shelf life or use centrifuge data to quickly predict long term shelf performance of products.


Introduction
A dynamic in-silico model captures the kinetics of 1-d gravity driven instabilities, in gravity or centrifuge, of fluid-infiltrated poroelastic media in a partial differential equation (pde). The pde yields the porosity profile over height and time for the given initial and boundary conditions, during slow isotropic isothermal compaction in counter-current fluid drainage.
Processes included are sedimentation, creaming, caking, synaeresis, (ultra) filtration or settling of sludge, pulp, and slurry, washing, sand and ice cream sculpting, filter pressing, injection moulding, consolidation and subsidence in civil engineering, in gas and oil recovery and in carbon dioxide sequestration, dynamics of biological tissue, etc. The dispersed medium can be anything from (even partially aerated) porous, fibrous, granular dispersions, suspensions, emulsions, fat crystal networks, clothes, soil, rock, or their mixes, and the fluid can be either oleic or aqueous. Typical products are pourable dressings, spreads, mayonnaise, jam, butter, margarine, ice cream, frozen drink, custard, mustard, cheese, ketchup, pudding, yogurt, rock, sand, soil, clothes, muscle tissue, shampoos, creams, soaps and gels, et cetera.
The most important limiting prerequisite is that the incompressible dispersed medium is sufficiently structured and/or concentrated that it compacts isotropi- For Unilever, modeling of gravitational instability of products is important to quantify or extrapolate long time behavior during shelf life or use centrifuge data to quickly predict long term shelf performance of products.
For now, the profile is created in a 1-d vessel or reservoir, e.g. vessel or an aquifer, of fixed volume. The medium has initially constant composition and porosity, without gravity or centrifuge effects, as if the medium is freshly mixed or in zero gravity. The boundary conditions (bc) are no flow through the boundaries.
After the compaction process starts, the vertical profile almost immediately converts into three zones separated by two fronts: a clear (depleted) supernatant fluid layer, which borders the initial mix centrally, which borders the caked (elastically or plastically compressed) layer at the other side of the vessel. On which side the fluid layer occurs, depends on the density difference between the fluid and dispersed particulates. In sedimentation, the dispersed solids are heavier, causing the fluid layer to build on top; in creaming, the reverse occurs. The initial dispersive mix in the middle moves at constant speed and constant porosity, subject to buoyancy, and decreases in length as the clear fluid layer and caked layer increase in size. The cleared fluid layer height increases first linearly in time, until the two fronts merge into one front that slows in time due to the elastic force that builds on compaction and counteracts the buoyancy (see Figure 1), until a no-flow balance is reached at equilibrium. As fluid displaces dispersed material in pure counter-current flow, at each level the volume rate of fluid is equal to the volume rate of the dispersed material, but in opposite directions.
Other boundary conditions will create different profiles, for example, during a clothes wash process, in synaeresis, or in filtering compaction, et cetera.

Assumptions
The continuous incompressible fluid (aqueous or oleic) is assumed Newtonian. The isotropic dispersed medium is concentrated or structured in such a way that all incompressible dispersed ingredients move together in the opposite direction of the fluid, while compacting during fluid drainage (counter-current flow). The dispersed medium can be any form, including dispersions, emulsions, clothes, soil, rock, or their mixes, and in any state (aerated, porous, fibrous, granular, et cetera). The fluid wets the dispersed medium. Flow can be driven by buoyancy, centrifuge, applied pressure, simple fall impact of a wetted cloth, et cetera. Temperature is assumed constant. Coarsening (Ostwald ripening of dispersed material and bubbles) is assumed absent. The flow rates are so slow that the fluid follows Stokes flow and Darcy's law applies for the flow through the deforming porous media. As we are considering the flow in a vessel/reservoir with a closed bottom and closed top (which can be a meniscus), the total sample volume is conserved. The vessel is assumed "wide", such that any potential side-wall drag is not interfering with the flow rates, and the pore averaged flow is therefore assumed purely 1-dimensional in a vertical upward or downward direction. The total volume below any level does not change: the Darcy speed of particulate downward in sedimentation is therefore at any level exactly equal to the Darcy rate of fluid upward (or both zero). The same applies in reversed mode (e.g. creaming). For the small air bubbles, we assume that the bubble volumes are preserved, which is a simplifying approximation. Volume changes of bubbles are small anyway in the relatively small pressure gradients. The total gas volume fraction must also be constant, as all other ingredients are assumed incompressible and total volume is preserved. Whether the dispersed phase acts elastically or plastically is irrelevant in slow 1-d sustained compression: both for snow and cotton we feel in compression tension the reaction force, but on release of tension snow does not relax back (e.g. is plastically deformed), while cotton expands back (e.g. is elastically deformed).

Method
The fluctuating pore-level Navier-Stokes equations for the fluid flow at sub-pore level are averaged to macroscopic flow using poroelastics. This leads to a Darcy type of fluid flow equation driven by buoyancy forces, slowed down by counteracting elastic forces when the deformable porous medium compacts. The flow equation, together with the volume balance, lead to the partial differential equation (pde), that needs to be solved numerically for the proper initial and boundary conditions. Scaling relations are used to express the pde in scaled dimensionless variables such that one solution actually covers a whole range of actual conditions, dimensions, timings, and variation in parameter values, which cover many processes in various formulations and products.
The pde uses four Pi numbers, based on dimensionless combinations of one value of the seven relevant physical parameters (permeability, elasticity, density, viscosity, height of the vessel, acceleration of free fall, and initial porosity), to characterize any porosity profile as a function of height and time; from the initial waiting and transient to ultimate equilibrium. Two Pi numbers are scaling of time and height, one is the initial porosity, and the last Pi number is the ratio of buoyancy and elastic force.

Flow Averaging in Multiphase Flow
The local balances of mass, volume and momentum of each phase, e.g. fluid phase and dispersed structured solid (fabric) phase, are spatially fluctuating functions on the sub-pore level. These equations are averaged to a larger volume, to damp out the local pore-level fluctuations, and find the equations for the mass balance and the porosity, thus vary over distance and time, and make this second order pde nonlinear and awkward to solve. Cleverer combinations of these two varying parameters can be derived to create better constants and make life easier. The permeability is a function of both porosity and surface area per unit volume. For sphere packings and less regular assemblies of granular media, Carman-Kozeny found a typical experimental correlation ( ) , where a is the grain surface area per unit of grain volume and α c an experimental correlation coefficient [6] [7]. We may assume that the specific surface area a of the dispersed phase does not change significantly during compaction, e.g. during compressional porosity changes. This is a good approximation even for the lower porosities reached in a compacting porous cake, hence we may write k 0 is the permeability at a reference porosity 0 φ . The elasticity of a fibrous network pack varies with compression. To describe the elasticity, we will use the Van Wyk's heuristic equation [8]. For a fibrous fabric, it relates the applied pressure P to the inverse cube of the volume V, the intrinsic uncompressed volume V 00 , corrected for the incompressible (high pressure solid limit) volume V s as follows [9] [10]: . The equation is independent of the fibre diameter or elasticity, but includes the Young's modulus E, the mass of the fibres m and the bulk density ρ at low pressure. Further, the equation comprises a dimensionless constant K c characterising the fibres, of order 0.01 that will vary with fibre orientation. Now, the overall elasticity of the gel network is defined by , thus using the modified Van Wyk's equation In vertical compression, the volume of particulates is conserved; . So, between two states in compaction 3 4 We selected the reference state as the porosity condition, where we can measure the value of porosity dependent physical parameters like k 0 and ε 0 . In practice, this is the initial condition. In the uncompressed state, the pressure is zero when no external forces, like gravity, are applied. In the initial condition, the pressure might not be zero, even without external forces. The material might tend to swell (if excess free fluid is available) or crimp in all directions, even without external force.
If the dispersed material swells initially in an exponentially decaying process, the initial fluid created by the gravitational compaction of the depleted layer is resorbed by the swelling, until the rate of swelling is lower than the compaction rate and the depletion layer starts to build. This leads to an effective extra waiting time before the observable depletion layer (supernatant layer) starts to build. If the dispersed material crimps initially, this can be interpreted, and modeled, as a synaeresis process with an evenly distributed non-zero pressure initially. We will for now assume that sufficient concentration of dispersed material is present initially to avoid crimp.

The Poroelastic PDE and Its Scaling
For our pde we may write where we introduced the height of the aquifer or vessel H as spatial scaling. If we now introduce a dimensionless scaling parameter ( ) , e.g. negative in sedimentation, a dimensionless time . If we define a dimensionless fluid flow rate Q′ as , the dimensionless volume balance becomes We also define a dimensionless height of the supernatant Θ as Θ = Ω/H.

Scaling Consistent with Buckingham Pi Theorem
The dimensionless pde is consistent with the Buckingham Pi theorem. Next to the dimensionless dependent variable φ , we have 7 independent parameters φ is dimensionless, it is a Pi number in itself. As z and H have the same dimension, their ratio is a Pi number [11]. Thus, we need only one consistent set of the 7 independent parameter values [ 0 φ , Δρ, g, H, k, ε, μ] to deter-mine the value of the 4 independent dimensionless Pi numbers which characterize the porosity profile over time, from initial waiting, during linear growth and slow down of the clear layer height, to final equilibrium. Alternatively, these parameters like permeability and elasticity, can be determined indirectly from the experimentally characterized profile, for instance, by measuring the height Ω of the clear supernatant layer as a function of time.
The Froude number F is the ratio of the inertia force over gravity force, for a system characterized by a density ρ, a length L, and a speed v: In the same way, the Cauchy number is the ratio of the inertia force over the elastic The F and C numbers are taken squared to linearize in v.
Hence, the ratio of the Cauchy number over the Froude number gives the ratio of the buoyancy over the elastic forces: . This is proportional to our ( ) The δ is therefore also a measure of the ratio of gravitational and elastic force. For 0 φ = 0.3 we expect a balance when δ = 4.2.

Initial and Boundary Conditions
The initial condition is (for now) constant porosity 0 φ at t = 0. In practice, any monotonic porosity profile, constant or increasing towards the depletion side, can be taken as the initial start profile. The pertinent boundary conditions (bc) are no flow over the boundaries of the vessel Q = 0 at ξ = 0 and ξ = 1 at all times, hence conservation of dispersed material in the vessel: , which means that in view of the volume conservation at that cake boundary. If the pde prevails up to the boundary, and φ can be differentiated twice, then at that no-flow boundary . This is certainly so at the cake side boundary. This is a nonlinear Robin type of boundary condition. At the depletion (supernatant) side, the porosity rises in the beginning from 0 . Thus, as long as the porosity function is not satu- rated ( φ < 1), the no-flow boundary condition at the depletion side is and thus ( ) . This is the nonlinear Robin type of boundary condition again. In the experiment, the time needed for rising porosity to unity at the depletion side is an effective part of the waiting time: no real visible depletion (supernatant) layer is formed during that period of time.
As long as there is no shock front formed ( 0 1 φ < < ), the porosity obeys the pde over the entire range, hence ( ) ( ) ( ) ( ) However, as the porosity reaches φ = 1 at an infinite thin depletion/supernatant layer, the boundary condition is a no-flow condition 0 Q′ = at the bottom flank of the depletion layer of thickness Θ, where φ = 1, while the flank forms a moving boundary (e.g. a type of Stefan problem), a sharp moving shock front, described by Rankine-Hugoniot shock conditions. In the depleted (supernatant) layer is φ = 1 constantly. The lower boundary of the depleted layer is an almost vertical moving front ending in a sharp nick towards φ = 1. We may interpret this moving boundary as a space shifting of a Dirichlet boundary condition φ = 1 (over a spatial range with an extend Θ from the top of the vessel).
If we write the pde as , we recognise the general shape of a nonlinear degenerate elliptic-parabolic pde, the Richards' equation. The left-hand-side (LHS) is the hyperbolic sedimentation equation for the non-diffusive/non-dispersive sedimentation process. The settling velocity or conservative flux convection coefficient is an increasing function of φ . The LHS leads therefore to a shock front (kinematic wave) propagating with a speed U( 0 φ ); the sharp front separates the homogeneous suspension of initial porosity 0 φ from the clear fluid (the depleted, supernatant layer). Hydrodynamic dispersion modifies the above picture and spreads out the front.
The ultimate balance of self-sharpening and dispersion depends on δ, but is difficult to predict a priori, since the pde is nonlinear. The nonlinear nature rules out analytical techniques and consequently the equation can only be solved numerically. See, for instance, the work of Concha and Bustos, and of Bürger, Karlsen, and Evje, on sedimentation and flocculation, around the turn of the century, that describe flocculated suspensions.
One important characteristic is that the effective dispersion coefficient  becomes infinitely small for high values of φ approaching un-ity. This means that the dispersive broadening becomes very small for the higher porosity values near the start of the cleared zone. We know from the boundary condition ( ) that the flank of the cleared zone is also very steep for high values of φ . This makes the flank sharp and steep, as observed experimentally. On the other end, near the base of the caked zone at ξ = 0, the porosities become small. This means that the caked porosity profile is reasonably flat according to the boundary condition, helped by the fact that the dispersion is large, mixing the porosities over a wide range of ξ (similar to a "random walk").
This flatness of the cake is also observed experimentally. Basically, sharpening and dispersion are opposing effects on the front, thus they can balance each other, leading to a stabilized dispersed profile that, after a transient time, will propagate with a constant velocity U, without changing its shape [12]. In the experiments, stable and sharp moving fronts were observed.
Another important characteristic of the porosity profile in gravitation occurs, when starting with a profile that is monotone, neither decreasing (in sediment) nor increasing (in creaming), it will retain that same monotonic characteristic up to equilibrium. Any time or spatial oscillations indicate numerical instabilities in the numerical procedure.
One solution of the pde delivers a full set of curves that cover variations in all parameters for a given initial porosity 0 φ and for a given ratio δ of buoyancy and elastic force. Repeating the calculations for other values of 0 φ and δ give all possible solutions that exist in nature for slow counter-current compaction of incompressible structured or concentrated particulates and small gas/air bubbles in an incompressible Newtonian fluid. For each scaled non-dimensional curve φ (T, ξ) determined by a combination ( 0 φ , δ), there is in real space and time a fan of curves φ (t, z) as determined by the proportional scaling of time dimension t = τT and of space dimension z = Hξ. Sedimentation (δ < 0) and creaming δ > 0) are mirror images in vertical direction around the φ = 0 φ value when acceleration g is a constant.

Centrifuge: Radial-Symmetric Flow in Axial Rotation Symmetry
We describe the flow in a centrifuge in time t and in cylindrical coordinates (r, θ, z). In a fast spinning centrifuge, the flow is radial symmetric in axial symmetry, hence flow and porosity are a function of time and of radial coordinate r. For flow in radial direction in radial symmetry In the gravitational field, the acceleration ˆz ge = g , is practically constant (g = −9.8125 m/s 2 ), but in a centrifuge the centripetal acceleration is a function of r: , where ω the cyclic frequency. The buoyancy force increases pressure with distance from the axis in radial flow for axial symmetry in the centrifuge. In this radial flow in axial symmetry, the volume balance becomes: . For a small rectangular vessel at larger distance from the axis , where r is an average value of the distance of the sample from the rotation axis. We define Z for the sedimentation of a small sample in the centrifuge as 0 Z r r = − . In sedimentation, r 0 is the outer radius of the sample and 0 2 r r H = + . In creaming, r 0 is the inner radius of the sample and 0 2 r r H = − and 0 r r Z = + , with Δρ being negative. So, for such a small vessel spinning at large distances from the rotation axis for the small sample in the centrifuge the same role as g grav , but often represents a much larger gravity in fast spinning centrifuges: where RGF a large positive number, or in such centrifuges: If we represent earth's gravity data as a function of z and compare them with centrifuge data for small samples far from the axis as a function of Z, in for instance a sedimentation experiment, the graphs would be (nearly) identical at the same imposed gravity forces. This also applies when scaled dimensionless, using the same four Pi parameters H, τ, 0 φ , and δ.

Equilibrium Profile
In the ultimate equilibrium, the fluid flow 0 Q′ = over the entire vessel, hence ( ) . Integration yields the equilibrium profile ( ) point. If we take the reference point at the point In radial symmetry in equilibrium ( ) ( ) At the axis r = 0, the porosity is (virtually) maintained at the initial value 0 φ φ = as there is still the initial pressure irrespective of rotation speed, and from the axial symmetry, hence const = 0, which leads to the same condition as in constant gravity, but now z is replaced by r: g. at equilibrium, using Van Wijck in axisymmetric counter-current flow of incompressible fluid and particulates in a closed vessel in the centrifuge: . The equilibrium profile is over a finite part of the vessel: the rest is depleted zone (in the centrifuge at the inner part of the vessel in sedimentation, at the outer part in creaming), such that the particulate volume initially present over the entire vessel is preserved and compacted over the shorter equilibrium profile in the caked compaction zone. In the strong gravity field of the centrifuge, there is significant compaction, even near the top of the cake, which seems counter intuitive, as there is no push from excess weight above. However, the buoyancy force is present, even at the top of the cake.

Weak and Strong Form
The strong form requires the function φ to be continuous and at least twice differentiable over its domain. The solution of our pde: may contain discontinuous jumps, which are weak, around which the function remains continuous, twice differentiable and is therefore strong. If we integrate once over space: This integral describes the same physics but needs to be only once differentiable. This is a weak formulation and is approximate.
finite difference formulation is also an approximate and weak solution of the full equation and boundary conditions. When the spacing is smaller, the solution gets stronger. The position of the shock is described by the displacement from the flow from both the upstream and downstream sides, e.g. Rankine-Hugoniot jump conditions, as imposed by the requirement of conservation of energy, momentum, mass and volume.

Initial Speed of Shock Front
The volume balance can be rewritten as a weak shock front speed . Over a shock front the volume balance leads to Rankine-Hugoniot condition d d For the initial linear displacement cleared front Θ between 0 φ and depleted layer with 1 with respectively ( ) Since the dimensionless clear zone length in sedimentation is we find for the initial linear growth

Initial Shape of Shock Front
We may now calculate the (stronger) shape of the stable linear displacing shock front that is expected to retain its shape as long as the intermediate layer and the supernatant layer are present, e.g. as long as shock's rate of displacement is only a function of u, e.g. φ (u): Hence Thus, our pde is then only a function of u: , where c 1 is an integration constant. But we know that for  As the wave travels with constant shape, we may shift to pass at T = 0 the value φ = 1 at ξ = 1 and follow the curve: These parameter choices will be shown later to represent a virtual "sedimentation", e.g. inverse of the actual experimentally observed creaming, of a Hellmann's mayonnaise sample structured with whole egg (WE 6.2%) in a centrifuge at 1000 The height of the clear supernatant layer Ω is easily monitored in the product and is often used to experimentally measure compaction over time in earth's gravity or centrifuge. In this subsequent slowdown from elastic counter forces, we expect a displacement characterized by ~t Ω from the similarity with the diffusion and heat transport equations. We may anticipate this from a simplified model.

Simplified Model in Earth's Gravity and/or at Small Density Differences in Limit δ → 0
For small δ, like in earth's gravity and or at small density differences, our pde is . For small changes of φ , we have approx- The major effect in the small earth's gravity is therefore that a value of φ will displace as square root of time z t , as observed in many other physical phenomena, like chemical diffusion, heat transport, et cetera. For the slow 1-dimensional poroelastic drainage in the earth's gravity field there is a functional dependence or scaling with a dimensionless group 0 0 This is an important scaling number for counter-current compaction drainage processes in the earth's gravity in poroelastic media and is the equivalent of the Fourier number for heat transport. A typical time for approach to equilibrium would be a displacement over height H, e.g. during the drainage process. Hence, the permeability reduction and the elasticity enhancement are assumed to be compensated by their multiplication. We take a fluid-filled poroelastic medium (like a gel, a porous or fibrous network, representing an emulsion and/or a particulate) that is compacted by drainage from one side (like pressing with an open filter that allows easy passage of fluid and does not allow any passage of the particulates). We may then solve the pde analytically.
In the present case, we assume that: the initial condition is It is therefore assumed that the material drains from one side in a one-dimensional way and that the material is so thick, or times are so short, that the drainage disturbance does not reach the back of the fabric, within the time span of observation. Standard mathematics leads to a solution ( ) 0 The amount of fluid flowing is ( ) . Liquid flows in the direction of lowest pressure and highest porosity. To make liquid flow out of the porous medium through z = 0, the pressure on the outside must be lower than The liquid continues to flow with time at a diminishing rate. This is the case because it was assumed that the disturbance has not (yet) reached the other side of the medium. In our case, the production of fluid at z = 0 leads to a compaction shortening (e.g. depleted layer) of length where Ω the compaction length. This simple model applies in processes like filtration compaction, machine and hand wash, et cetera. In this simplified gravity compaction model, the compaction shortening Ω is proportional to t during the entire time span that the compaction disturbance has not reached the other side.

Height of Supernatant during Creaming of a Fibrous Nata de Coco Gel Network
The supernatant height Ω(t) = h(t) of an aqueous fibrous system of micro fibrillated bacterial cellulose (BC, Nata de Coco) and Carboxy Methyl Cellulose (CMC) at ratio 14/1, with captured soy bean oil (SB) droplets was followed in earth's gravity as a function of time [13], see Figure 3.
The volume of the sample is 20 ml, height ca. 42 mm. The experimental data of the initial rate as function of BC and Oil concentration are shown in Table 2; The creaming data for 0.06 wt% BC/CMC with 1 wt% oil is reproduced in Figure 4 and Figure 5 (experiments in triplicate): The cake compacts proportionally to the square root of time: The linear slopes are given in Table 3.

Creaming of Mayonnaise
Dressings with different formulations and varying processing were tested in a centrifuge to accelerate gravitational instabilities (sedimentation and creaming). These data are needed to predict the shelf life and to test structurants that allow the development of cheaper, lower carbon footprint, clean label, or better sustainable dressings. The LumiFuge centrifuge accelerates the study of gravitational instability behaviour of over a year into just 8 hours, due to a factor 1100 magnified gravitational force used in the centrifuge compared to earth's gravity.
Such data is often difficult to obtain in real time as the mayonnaises/dressings might be attacked from bacteria and/or fungi, might become chemically instable, or might change consistency on a smaller time scale, not directly related to gravitational instability during shelf life, but obscuring the results. As the centrifuge data for creaming of the dressings is rather complex, including changing rate over time, changing rate with speed of revolution and with composition, it is difficult to pinpoint centrifuge conditions (time/speed) that can be directly related to the magnitude of the actual slow creaming in the earth's gravitational field during shelf life. Therefore, a more elaborate analysis is warranted.
Several experiments focussed on the creaming of a Hellmann's Real mayonnaise formulation WE 6.2%, both in earth's gravity and in the LUMiFuge at 1000, 2000, 3000, and 4000 rpm; physical characterisation as follows.

Initial Porosity of WE 6.2%
The 67.5 wt% bean oil in the WE 6.2wt% formulation has a density of 0.9192 at 20˚C. Fresh whole egg consists typically of 74 wt% water, 11.8 wt% lipids and 12.8 wt% protein, and 1.4 wt% dissolved material (carbohydrate and minerals). A fresh egg has therefore 75.4 wt% water phase and 24.6 wt% solids. Some 34% of egg weight is the yolk. A fresh egg has a density of 1.033 [14]. Therefore, 100 g of a WE 6.2wt% formulation consists of 6.2 g WE 6.2%, and contains 1.5 g eggs solids phase of density 1.132 g/cm 3 and 4.7 g egg water of density 1.00. The rest of the WE 6.2% dressing is water plus dissolved solids = 100 − 67.5 − 6.2 = 26.3 wt%, with density 1.00. If we assume 100 g formulation, we have 67.5 g bean oil and 1.5 g as egg dispersed material and 26.3 + 4.7 = 31 g water phase. The volume of the dispersed solid is 67.5/0.9192 + 1.5/1.132 = 73.43 + 1.325 = 74.76 cm 3 solids and 31/1.00cm 3 water phase, thus in total 105.76 cm 3 . The initial volume percentage of dispersed material is therefore 74.76/105.76 = 70.7 vol%, hence 1 − 0 φ = 0.707, and the volume percentage of water phase is 29.3, hence the original homogeneously mixed WE 6.2% formulation has an initial porosity of 0 φ = 0.293.

Densities of WE 6.2%
The density of the aqueous fluid phase is around 1000 kg/m 3 , and the bean oil is 919.2 kg/m 3 . The density difference is therefore ( ) 3 1 0.9192 1000 80.8 kg m ρ ∆ = − = . We neglect the influence of the 1.5 wt% egg solids on the density difference, because the contribution is small and the egg white and yolk solids may redistribute between water and oil. Open Journal of Physical Chemistry

Viscosity of Fluid in WE 6.2%
The viscosity of the aqueous fluid phase is close to water, e.g. μ = 1 mPas.

Gravitation and LumiFuge
In our LUMiFuge centrifuge experiments, RGF can be 6-4000, vessel height H typically 22 mm, and axis of gyration typically 131 mm at the vessel's base, or on average 120 mm r = . We may assume that the LumiFuge's centrifugal force , then RGF = 1207 for WE 6.2% at 3000 rpm.

LumiFuge
The LUMiFuge (Figure 6   meniscus, at the top surface of the sample. A typical collection of results for samples running at 4000 rpm for 24 h is given in Figure 10. Here we see the clear layer at the bottom of the sample. In some cases, here in the samples in the middle, we see a thin oil rim at the top of the sample. As these samples have been running for a long time (24 h), the clear layer has practically reached the maximum thickness. The clear basal layer indicates that even the tiniest particles are not only captured in the network but also move concurrently within it. When the water phase is not clear, the particles move more independently, Figure 9. A sample of the LUMiFuge optical transmission versus position for Hellmann's Real WE 6.2 wt% at 3000 rpm during 2 hours. Creaming is also observed in the earth's gravitational field with new samples of Hellmann's Real in test tubes shown in Figure 11.

Experimental Data LumiFuge WE 6.2% at 3000 rpm
In these centrifuge experiments, the measured creaming height is sometimes indicated by h, but it actually always the height of the depleted layer Ω (in creaming) of our previous modeling. The creaming data of WE 6.2% at 3000 rpm from Figure 9 is shown in Figure 12.   The two traces of the creaming height as function of the square root of time for WE 6.2% at 3000 rpm are given in Figure 13 and Figure 14 (duplicate experiments, axes reversed).
The graph shows the linear part (right-hand side) and the reduction of speed when the creaming zone reaches the top of the sample (left-hand side). The linear part allows us to calculate the creaming speed.

Experimental Data LumiFuge WE 6.2% at Several Speeds
The experimental LumiFuge data for WE 6.2% at 1000, 2000, 3000, and 4000 rpm are shown in Figure 15: We see that all curves have essentially the same slope There is a scatter in the waiting time. probably due to differences in spin up timing of centrifuge. The square root linear parts of ~RGF grav g T ξ ρ ∆ are Figure 13. The creaming height Ω of WE 6.2% as function of square root of time for 3000 rpm followed during 2 hours (experiments in duplicate).  Table 4 and averaged in Table 5.
Our data are thus averaged:

Experimental Data WE 6.2% in Earth's Gravity Field
We may compare these data with the visual observation of the creaming height Figure 15. Log/log plot of WE 6.2% centrifuge data. Ω in a test tube as function of storage time (in earth's gravitational field), see Table 6 and Figure 16. The first data points are tricky as the tube had a V-shaped base (first few mm), so the early rates are difficult to read.

Extrapolation Centrifuge Data of WE 6.2%
to Reveal Shelf Life Stability Combining the gravity and centrifuge rates versus square root of time in one graph, for WE 6.2%, we obtain Figure 17: This graph shows that we cannot use a linear extrapolation to predict the creaming in the earth's gravity field. But if we plot the data on a log/log plot, see Figure 18, we find surprisingly: When RGF increases, we would expect the experiment to go faster as we exert a larger force ω 2 R. It is therefore not unlikely that Ω 2 /t is some function of RGF, e.g. even appearing to be directly proportional to RGF. We expect similar responses for other structurants than whole egg at 6.2 wt%, with a different proportionality constant. This equation already allows us to predict the (square root of time) rate in earth's gravity (RGF = 1) from the rates measured in centrifuge by a simple Table 6. Creaming of WE 6.2% in earth's gravity field. WE 6.2% week 1 week 2 week 3 week 4 week 5 week 6 week 7 week 8 Height Ω/ml 0.5 1 1.    log/log extrapolation, e.g. from a linear fit on a log/log plot of Ω 2 /t versus RGF.
For WE 6.2% at 3000 rpm, the linear part in square root of time quantified to d d 268.96 h t = and 247.12 mm month , on average 258 mm month in these two repeat experiments. We use the average slope at 3000 rpm of 258 mm month t Ω = , e.g.  mD and μ = 1 mPas, we estimate ε 0 = 3870 N/m 2 . Thus, for creaming of WE 6.2% at 3000 rpm, δ = 261.
Likewise we find for 1000 rpm: Figure 20. The predicted equilibrium curve for WE 6.2% at 1000 rpm.  NB. In sedimentation, Δρ is positive and thus δ is just the negative value, and hence the porosity profile is the vertical mirror image, when the acceleration (buoyancy) can be considered constant (earth's gravity field) or approximated to be constant (small sample in centrifuge). We may conclude that for the mayonnaise centrifuge experiments with small H and Δρ and thus, small value of δ, the diffuse square root dependence dominates quickly in the response.

Waiting Time
In all experiments, we see some time has passed before the cleared layer becomes visible, e.g. the waiting time. Part of that is the short time needed to rise the porosity to unity in the top of the vessel. In our experiments another part is the spin up time of the centrifuge. If a structured sample is slowly shaken just before the experiment to even out the porosity (or set up vertically after first slowly rolling horizontally), it might take some time before the firmed structure is developed, which might contribute to the waiting time.
The waiting time can be prolonged in the experiment by adding hygroscopic ingredients, like starch, that adsorb the initial free water, by swelling the porous matrix. This adsorption technique is often applied to extend shelf life in mayonnaises. The swelling suppresses the build-up of a depleted layer until the diminishing rate of swelling is equal to the depletion rate, from which moment on, the depleted layer starts to build, quickly accelerating to a constant rate, as swelling is decaying exponentially in time. The swelling and prolonged waiting time is not in the model yet.

Subsidence in the Groningen Gasfield, The Netherlands
When we look at, for instance, subsidence in the earth's gravity field, like Groningen's gasfield [17], as shown in Figure 21: We see that the subsidence response Ω/t even after several years of constant gas production was still in its initial linear constant slope transient. Nowadays, the gas production has declined, so a new, smaller linear trend slope of constant Ω/t must have started, that will transform eventually in a dispersive trend Ω 2 /t before the equilibrium subsidence depth is approached after many years. An independent determined value of the average elasticity of the rock formation Open Journal of Physical Chemistry might be used to predict this later time Ω 2 /t approach to equilibrium in the future and thus predict the decay of subsidence in the far future.

Numerical Simulation
In our numerical simulations, we describe the centrifuge WE 6.2% creaming experiments in the equivalent "sedimentation" mode by simulating with the oppo- , e.g. a front displacing with a constant rate and supernatant thickness Θ growing linearly. When the height of the intermediate layer has become zero, the supernatant grows gradually slower in height Θ until the equilibrium porosity profile is reached, with a balance between gravity and elastic forces. Typical values for δ in "sedimentation" of WE 6% mayonnaise are the negative of those given in Table 7  We see that at T = 0.01 (end of first time-step), the cake has already reached the normal exit boundary at ξ = 1, indicating that in a normal vessel length ξ = 1 the intermediate layer has eroded at T = 0.01 for WE 6.2% at 1000 rpm and the cake had already reached the thin supernatant layer. Hence, T = 0.01 is near the end of the linear growth period of the supernatant.
The subsequent rate of cake formation Θ slows, initially proportional to square root of time, and must probably be generated in material coordinates as a weak solution with a Stefan boundary condition. We tried a strong solution, but that is not stable in the supernatant layer. By mere luck we generated a solution that was reasonably stable for longer time. This example, calculated with COMSOL5.3a for δ = −30, 0 φ = 0.3, approximately equivalent to WE 6.2% at 1000 rpm is shown in Figure 24. Ten time-steps of 0.005 to dimensionless time T = 0.05 remain reasonably stable: Here we see a build-up of the supernatant layer and a diminishing of the supernatant depth growth rate at later time (T > 0.01). The supernatant porosities are fluctuating ( φ ≠ 1), and not reaching φ = 1, with a steep slope and a sharp nick to steady φ = 1 in the supernatant, indicating numerical instability. The later time displacements, and other simulation attempts, were not accurate due to the numerical errors in the supernatant range.
During slow down, improving the accuracy of the simulations might confirm the expected Θ 2 ∝ T, e.g. Ω 2 ∝ t behaviour. Because the displacement profile up to equilibrium has some shape similarity, we might try a self-similar transformation to predict larger time behaviour. This is not yet done.

Transformation to Material Coordinates
A trick might be used to avoid the difficult Stefan moving boundary problem by stretching the caked layer to remove the supernatant layer from the domain of integration. We have a dispersed material volume balance such that any two material-bound spatial coordinates at a small distance dξ 0 compact as , with boundary conditions (bc) at χ = 0 ( ) and φ = 1 at χ = 1, and initial condition T = 0, φ = 0 φ for 0 < χ < 1. Expressed in coordinate χ, the zone with dispersed material is virtually stretched, such that the compaction is compensated to keep the total length of that zone constant. By rescaling Notably, at equilibrium , e.g. integrated . At the entrance φ = 0 φ ∞ for u = 0 or which can be solved for 0 φ ∞ . For 0 φ = 0.3 and δ = −30 we find 0 φ ∞ = 0.179, which is indeed close to the predicted equilibrium curve (based on an estimated value of Θ equil ). We may now calculate Θ equil accurately by using , thus, at equilibrium: 3, δ = −30, we find ξ flank = 0.70, e.g. Θ equil = 0.30. Hence, we expect the experimentally estimated (long time) "equilibrium" thickness of about 5.5 mm to grow to 6.6 mm at full equilibrium (H was 22 mm). The earlier estimate Θ equil = 0.25 was probably too short. There is no need to estimate that equilibrium supernatant length from the experiment, because it is determined by the 4 independent Pi parameters and can thus, together with the equilibrium intercept 0 φ ∞ , be directly calculated from those parameters.
Using solution profile φ (u, T) at any given T, we may re-tabulate the proper curve φ (ξ, T) by the transformation . It is important to realize that we were indeed able to calculate all characteristics of the compaction curves from waiting, transient time, up to equilibrium, with only the 4 independent dimensionless groups (t/τ, z/H, δ, and 0 φ ) as given by the Pi theorem. The parameters themselves were calculated using only one value of the seven physical parameters [ 0 φ , Δρ, g, H, k, ε, μ] that appear in determining those 4 dimensionless groups. The only two difficult physical parameters to determine independently are permeability k and porosity ε. However, these can be determined using the experimental data, from the early linear slope of Θ vs t and later linear slope of Θ vs t respectively. The simulation data in material height coordinates u for sedimentation WE 6.2% at several rpm's are given in Figure 25. These simulations are reasonably accurate but are also time consuming (several hours on a 64 bit W10 AMD FX-8150 Eight-Core 8 Gb ram 3.6 GHz bench top computer). An overview of several representative simulations (COMSOL5.3a): Apparently, the flank and slope become steeper and flatter for increasing δ , with both the value and trend in the porosity of the intercept at u = 0 agreeing with the experimental observations. As the vertical differences in subsequent time-steps rapidly diminish, at an order of half of the difference in the previous time-step, the equilibrium curve is expected to be lower, approximately less than or equal to the shift during the last recorded time-step, hence already very close to the last and lowest recorded curve. As noticed before, the proper curves φ (ξ, T) can be found by the transformation ( )   This gives an indirect way to resolve φ (ξ,T) without stepping into the intricate Stefan moving boundary problem, and still find a rather strong solution. The physics of the moving boundary problem keeps the jump sharp even in these rather strong solutions, as was indeed observed in many experiments. NB. The initial flat profile was modified by a small upward curve near the vessel's top, to suppress oscillations and blow-up at the start.

Bead and Spring Model
The poroelastic counter-current compaction can be seen as a porous elastic weak "Bead and Spring" model, subjected to a gravity field, with the coordinate system (e.g. coordinate z) attached to the dispersed medium (the bead and spring system) of lateral area A, as shown in Figure 27.
Rotate graph 90 degrees counter-clockwise for actual position in sedimentation, 90 degrees clockwise for creaming. The porous medium in the vessel is assumed to consist of N cells filled with a (repeating bead and spring) elastic porous medium. Initially, all cell "boundaries" are equidistant at distance Δz 0 = H/N and porosity 0 φ . The left side is the bottom wall attached to the first cell. The right-hand side N is initially a bead next to the other vessel wall, without being attached to the wall (or simply ending in a water-filled medium to the right of that point N). The clearing between the right-hand side wall and the right-hand side bead is the thickness Ω of the depletion layer (initially zero). When gravity is applied, the spring chain will fall (displace to the left on the graph) compressing the elastic springs (near the bottom wall) and pushing the fluid    consistent up to a reasonable time, as shown in Figure 28: The experimental data are scaled dimensionless by τ = 45,477 s and H = 0.022 m. In comparing the simulation with the experimental data at 3000 rpm, it appears that most misfit is due to the uncertainty of the zero time: a small waiting time is expected. The observed order of magnitude of 100 seconds seems rather large but not unrealistic considering the preparation and design of the centrifuge experiment where the centrifuge has to gain speed. The fit can be improved by adapting the roughly guessed simulation parameters τ, δ, and incorporating a waiting time, but this was not pursued further, since the reproducibility of the experiments is low (see the experimental section, Figure 15).
The detailed simulated profile for 3000 rpm is shown for δ = −261, N = 10 in Figure 29: The elementary bead and spring model is fast, but accuracy is uncertain, especially at longer time.

Discussion
The 1-dimensional consolidation theory in soil is often based on Terzaghi's theory that assumes that the soil is fully saturated with water, that water and soil are incompressible, the effective stress is the sum of the total stress in the rock and the pore pressure, strains in soil are small, Darcy's Law applies to hydraulic flow, permeability and volume compressibility are constant, and there is a unique relation between porosity and stress, independent of time. Our theory can cope with nonlinear changing of permeability and elasticity during the consolidation processes over larger, non-linear strains and stresses, and correctly incorporates the buoyancy force. be described by the same type of equations and relations, at least when the process is slow and reversible (elastic).

Impact and Outlook
There is a wide range of measures to stabilize dispersions/emulsions based on experience in the field that can now be understood and quantified according to our model: make fluid viscous and turn fluid into a gel during storage; add structurants to liquids (starch, alginate, micro-fibrillated cellulose, gums like LBG and guar); add a hygroscopic material (starch) to absorb the first free liquid; increase viscosity, decrease elasticity, permeability and density difference. We may now measure stability quickly in fast centrifuge and translate results straightforward into slow changes in shelf life in earth's gravity.
There are a large collection of models and literature on sedimentation, creaming and synaeresis in porous, granular, or fibrous media, for granulates, for suspensions, for dispersions, or for emulsions. The solutions presented there often apply to a smaller range and are usually based on semi-empirical interaction parameters for particle-particle interactions between specific particles. Our model uses a trick to circumvent the specific disperse interaction and focus on the fluid flow in poroelastic media, using smart scaling relations that make the theory more universal and applicable in many different formulations. Phenomenological models for thickening during sedimentation of flocculated suspensions are generally based on the Kynch model, where the compaction reaches a fixed lower porosity limit in the cake to avoid further compaction [19] [20] [21] [22]. In our model, that condition is reached by the local dynamic balance of buoyancy and elastic forces in the compaction drainage of porous media, using widely applicable smooth permeability and elasticity functions that at the same time determine the range of applicability, together with the simplifying assumptions mentioned in the introduction ( §1.1). The counter-current flow without segregation requires initial porosities below say typically 45%, or small values of parameter δ. Discrete droplet aeration asks for an overrun below 50% -150%, depending on level of air stabilizers. In soils, the residual gas saturations should be below typically 35% to trap the gas as discrete bubbles.
Extending the model by additional physics, like starting with a particular monotone profile, adding compressibility of fluid or solids, anisotropy, solving the full centrifuge equations in cylindrical symmetry, seems straightforward. Even replacing the Carman-Kozeny or Van Wijk equations respectively by the permeability and elasticity characterisations using other dedicated equations as a function of porosity can easily be done. The reverse process of drainage compaction, e.g. imbibition of fluid infiltrating into swelling porous media, is also important (in reservoir and civil engineering, like capturing rainfall), but is not modeled here, as plastic or elastic responses require different solutions. The reversible elastic mode of dynamic counter-current imbibition will resemble the reverse of our counter-current drainage compaction solutions. The extension to Open Journal of Physical Chemistry more than one flowing (connected) fluid phase (aqueous, and/or oleic, and/or gaseous) requires extra physics (capillarity, hysteresis, relative permeability) and additional models (Van Genuchten, Corey, Brooks). In our model, additional fluids are only present as disconnected dispersed bubbles (overrun) or droplets (emulsions) that are trapped in the dispersed phase and flow co-current with that dispersed phase, and counter-current with the draining fluid.