Abstract
In far-from-equilibrium chemical systems, self-organizing diffusion-reaction processes can give rise to complex patterns. Such self-organizing patterns are commonly found in diverse rocks, emerging as periodic stripes, rings, and various polygons. While these patterns are well documented, the drivers of their diversity and the mechanisms behind pattern selection remain unclear. This study investigates how geometric heterogeneity influences the formation of Liesegang patterns in natural geological materials such as Zebra rocks, pyrite ores, and orbicular granites. Using numerical simulations based on phase-field modeling, we replicate various Liesegang pattern morphologies observed in nature, such as circular, triangular, and interacting bands, and analyze their dependence on initial geometry and boundary conditions. We demonstrate that the spatial distribution of reactive fluids and the shape of nucleation zones play critical roles in determining the final pattern morphology. Our results provide insight into the self-organization processes in geochemical systems and offer a predictive framework for understanding rhythmic mineral banding in rocks.
1 Introduction
Liesegang patterns—periodic precipitation structures formed via reaction–diffusion mechanisms—are widely observed in geological settings [1–5]. These patterns arise from spontaneous chemical self-organization under non-equilibrium conditions and are commonly found in sedimentary and igneous rocks. Understanding the factors that govern their formation is crucial for deciphering paleofluid flow dynamics and geochemical evolution.
This work explores the influence of geometric constraints on Liesegang pattern [6] development through numerical simulations. We focus on three representative geological examples: hematite banding in Zebra rocks, pyrite zoning in ore deposits, and orbicular granites with concentric mineral rings. By varying the initial configuration of fluid sources and domain boundaries, we reproduce key features of these natural patterns and identify deviations from classical Liesegang behavior.
Rhythmic patterns in rocks are prevalent, originating from self-organized processes within a non-equilibrium, nonlinear geochemical system [1–5]. These periodic patterns appear across various scales, ranging from microscopic to millimeter, and even up to kilometer scales. Such periodic patterns are frequently associated with the Liesegang phenomenon, which involves rhythmic precipitation in gels. The Liesegang pattern formation is understood through the interaction between the diffusion of dissolved reactants and the kinetics of precipitate formation, including nucleation, growth, and ripening [6]. Analogous patterns in rocks include parallel mineral bands, concentric rings, oscillatory zoning, concretions, orbicules, among other forms. Numerous rock types exhibit patterned textures with periodic banding, potentially involving a Liesegang-type scenario. Examples include banded iron oxide minerals [7, 8], Zebra rocks [9–11], orbicular granites [12–14], and periodic oscillatory zoning [15, 16]. A thorough review of studies elucidating the analogy between Liesegang bands and geochemical banding is presented in [1, 2, 17].
Two primary models have been proposed to explain the mechanism behind Liesegang band formation. The “pre-nucleation” model [18] hypothesizes that nucleation and growth of precipitates occur when the reactant activity product exceeds a certain threshold, such as the critical supersaturation for nucleation or the solubility product for growth. Conversely, the “post-nucleation” model suggests that pattern formation arises from the instability and subsequent break-up of a uniform precipitate sol, driven by Ostwald ripening [19] and the Lifshitz-Slyozov instability [20]. Among post-nucleation models, one introduced two decades ago [21, 22] is based on phase separation theory. This model assumes that pattern formation results from spontaneous phase separation (unmixing) of the mixture. The phase separation model, underlying by the Cahn-Hilliard dynamics, has since gained prominence in multiple disciplines due to its conceptual simplicity and capacity to simulate irregular patterns and transitions between stripes and spots [23–26]. Our recent research has extended the phase separation model to replicate hematite patterns observed in Zebra rock from northern Western Australia [10, 11].
However, these studies have predominantly focused on regular geometries, overlooking the impact of geometric heterogeneity—defined here as the diversity in the spatial distribution, number, and shape of nucleation (initialization) regions—on Liesegang patterning in geological materials. In our context, geometric heterogeneity encompasses both the arrangement (e.g., regular or random) and the geometric shapes (e.g., circular, triangular) of nucleation centers, which collectively influence the resulting pattern complexity and interaction. This approach better reflects the intricate geochemical patterns observed in nature. Recent studies have also highlighted the importance of domain geometry and boundary conditions in shaping Liesegang-type structures. For example, [27] employed a vertex-based finite volume method to explore post-nucleation Liesegang patterns in confined geometries, demonstrating how symmetry and reversibility are affected by initial spatial constraints. Similarly, [28] combined reaction-diffusion equations with a Cahn-Hilliard-type precipitation model to simulate geochemical banding in limestones, including fractal analysis of pattern interfaces—an approach closely aligned with our interest in linking numerical simulations to natural textures. While our work builds upon spinodal decomposition principles similar to those used in [28], we emphasize the role of heterogeneous initialization regions in controlling pattern emergence, extending previous findings to more complex geological morphologies.
In this study, we integrate the phase separation model with explicit forms of geometric heterogeneity to replicate several complex geological patterns in rocks. First, we design a circular initialization region to reproduce the formation of single rhythmic pyrite rings observed in the Xingqiao sediment-hosted massive sulfide deposits [16]. Next, we model triangular-like Liesegang banding found in Zebra rocks from the Western Australian East Kimberley region. We then simulate more intricate Liesegang rings in rocks with multiple and heterogeneous initialization cores, using orbicular granites as a case study to replicate Liesegang ring interactions with various regular distributions. Finally, we reexamine the formation of heterogeneous pyrite rings observed in Xiaoqiao deposits with randomly distributed nucleation regions. A schematic figure is provided to categorize and illustrate these types of geometric heterogeneity.
The paper is structured as follows. In Section 2, we first present the phase separation theory for Liesegang pattern formation, followed by a discussion of three types of Liesegang patterns in rocks, casting their formation mechanisms within the phase separation model framework. Then, we state the problem and detail the numerical setup. Section 3 explores the influence of geometric heterogeneity on the formation of single and complex Liesegang patterns in three selected rock types, containing both regular and randomly distributed initialization shapes. Finally, we conclude our study in Section 4, highlighting potential directions for future research.
2 Methods
2.1 Phase separation model
Consider a closed system comprising two distinct components, denoted as and , characterized by concentration functions and , respectively. The system is divided into two distinct regions: an interior subdomain where component is uniformly distributed at a constant concentration , while component is entirely absent, and an exterior subdomain where component is absent and component is uniformly distributed at a constant concentration . At time , the two components initiate diffusion and undergo an irreversible bimolecular chemical reaction, governed by the reaction kinetics.in which the constant determines the reaction rate as described within the framework of mean-field theory, and is the reaction product and its concentration denoted by . The reaction-diffusion equations governing the evolution of two reactants in Equation 1 readwhere, the notation is the Laplacian operator while and correspond to the effective diffusion coefficients of two reactants and , respectively.
The reaction product experiences segregation into two distinct phases: a low-concentration phase and a high-concentration phase, without the formation of intermediate complexes in the phase separation model [21]. This disparity in concentration between the phases drives the emergence of precipitation patterns. In the present investigation, the concentration of generated reaction product transitions into an unstable regime, referred to as the spinodal region, where it spontaneously separates into a low-concentration domain (devoid of precipitate) and a high-concentration domain (containing precipitate). This process is mathematically described by the Cahn-Hilliard equation [29, 30].where denotes the mobility of the product. Additionally, accounts for stochastic influences, such as heterogeneity arising from initial concentration distributions, thermal fluctuations within the system, or external contributions from the surrounding environment. The total free energy density , which governs the phase separation process, is defined by a generalized chemical potential [29]. This potential comprises two contributions: the bulk free energy density and the interfacial energy density , expressed as:where, represents the gradient parameter, which is associated with the interfacial thickness. The terms and denote the bulk free energy functional and its density, respectively. When chemical reactions induce spontaneous segregation of the material into solid and fluid phases, an effective non-equilibrium chemical potential must be introduced to address potential microscopic deviations from action-reaction symmetry. The phase-separation model employed in this study describes the formation of stationary patterns observed in geological rocks as time approaches infinity.
The bulk free energy functional of the system exhibits two minima corresponding to equilibrium states, which are associated with the low and high concentrations of precipitate production, denoted as and , respectively. For simplicity, a Landau-Ginzburg double-well potential is employed to describe the free energy functional with minima located at and , and a maximum at , as shown in Figure 1.in which and represent parameters specific to the system, with the ratio determining the locations of the minima in the bulk free energy . The conditions and are necessary to facilitate phase separation within the system.
FIGURE 1
Taking the derivative of the free energy of with respect to the concentration , we can derive the generalized chemical potential hence writes asWe rewrite Equation 2, 3 in a dimensionless form by defining the normalized concentration, time, and length scales asand shifting the normalized concentration as
Coupling Equations 2–8, we obtain the dimensionless final set of equations governing the pattern-forming process:where represents the dimensionless gradient parameter. Additionally, and are redefined as the normalized mobility and normalized noise parameter, respectively, associated with the reaction product.
2.2 Three Liesegang patterns in rocks
This section provides three periodic patterns in rocks that potentially arise from cyclic precipitation as mineral-rich water infiltrates a porous rock and reacts to form an insoluble product, namely pyrite bands, Zebra rocks, and orbicular granites, see Figure 2. This occurs through the coupling of the precipitation reaction with diffusion and shows similar features to the Liesegang patterns. Here, we cast their formation mechanisms into the phase separation model framework.
FIGURE 2
2.2.1 Pyrite bands
One notable class of rhythmic geochemical patterns is observed in pyrite, the most prevalent sulfide mineral near the Earth’s surface, which appears in various forms, from microscopic framboids to macroscopic structures (see Figures 2a,b). This paper examines the formation of pyrite rings within the Xingqiao sediment-hosted massive sulfide deposits. These pyrite patterns exhibit periodic fluctuations in concentration within a single crystal, from core to rim, characterized by two main zones: dark zones rich in siderite and bright zones dominated by pyrite. Despite extensive research on the formation mechanism of pyrite rings, the process remains debated. Self-organization processes, linked to the Liesegang ring phenomenon, offer a promising explanation for these patterns. Recently, Qiu et al. proposed a Liesegang model based on pre-nucleation to explain pyrite ring formation and replicate textures in the Xingqiao deposits [16]. We reexamine this possibility using our post-nucleation model to reproduce similar pyrite patterns in the same geological context.
Our phase separation model hypothesizes that periodic pyrite rings form through the reaction of two soluble reactants . This involves a high concentration of outer reactant , i.e. , in contact with a permeable medium initially containing a low concentration of inner reactant , i.e. . As diffuses into the porous medium and reacts with , a solid precipitates in the wake of the diffusion front. The specific chemical reaction for pyrite precipitation is as follows:
We note that while (mackinawite) may transiently form as an intermediate during the above chemical reaction in Equation 10, its rapid conversion under near-neutral pH and sustained supply renders it outside the rate-limiting regime in our system [31, 32]. Thus, for phase separation modeling, we focus on the overall precipitation dynamics under conditions typical of sedimentary environments, where pyrite is stable and remains short-lived.
2.2.2 Hematite patterns
Ferric oxide (hematite and goethite) concentrations can develop in various sedimentary rocks, forming distinct spatial patterns such as concentric bands, planar bands, or columns. Hematite concentric bands are commonly found in chert concretions and nodules within limestones, indicating true self-organization during diagenesis, as these concretions are of diagenetic origin. Although these hematitic patterns are geologically minor, understanding their formation could highlight the complex kinetic interactions between reaction and transport in rocks, providing significant geochemical insights.
We examine the Zebra rock pattern from the northern region of Western Australia, known for its distinctive and prominent rhythmic hematite concentrations. Among all hematite banding patterns, Zebra Rock stands out due to its clarity and richness, making it an excellent model for studying pattern-forming processes related to subsurface environmental evolution. Zebra rock features reddish-brown bands, rods, and elliptical spots on a white or light background, as depicted in Figure 2c. Our previous research connected Zebra Rock formation to the Liesegang phenomenon, where supersaturation, nucleation, and depletion interact to create banded patterns. This connection was drawn because certain Zebra rock patterns resemble Liesegang bands. However, our earlier study did not account for the impact of complex geometries on pattern formation. Therefore, we now investigate how heterogeneous geometries influence morphological patterns in Zebra rocks.
Hematite in Zebra rock is a late-stage transformation product of Fe-Oxyhydroxide, precipitated from ferrous ions in an acidic geothermal environment, as described by the following chemical reaction in Equation 11:
2.2.3 Orbicular granites
Igneous rocks such as granites and diorites, including more basic plutonic rocks, sometimes feature ovoid structures known as orbicules [12]. These orbicules, ranging from several centimeters to about 1 m in diameter, consist of concentric, alternating light and dark mineral layers. They are believed to form through rhythmic selective crystallization of minerals around a foreign rock nucleus. The light-colored layers primarily contain plagioclase crystals with some dark minerals like hornblende or biotite, while the dark-colored layers are mainly composed of biotite or hornblende with minor plagioclase, see Figures 2d,e. The mineral distribution within each shell is irregularly circular or semi-circular, not symmetrically centered. It is suggested that orbicules form through counter-diffusion of heat and silicate components under conditions of rapid growth, potentially relating to self-organization and the Liesegang phenomenon.
Similarities between Liesegang patterns and multi-shelled orbicular granites, such as the geometric progression of band positions, imply that a Liesegang-like mechanism may partly explain their formation. In [13], Wang and Merino proposed an isothermal transport-reaction model based on a pre-nucleation framework to study orbicular granites. Although this model has been applied to the orbicular granite system, they concluded that a post-nucleation mechanism is essential for shell formation. Additionally, [33] noted that post-nucleation or hybrid pre- and post-nucleation models can more accurately reproduce realistic patterns in igneous layering or orbicular granite. In this study, we adapt the reaction-diffusion model by Wang and Merino to a post-nucleation phase separation model. The generic pathway for mineral precipitation follows this form:in which is one of the species in the melt such as , , , , , , , and is the stoichiometric coefficient of in reaction .
As for the oscillatory plagioclase in the orbicular granite, the specific chemical reaction in Equation 12 becomes:
Taking the above chemical reaction into the framework of , the outer reactants are the metal irons in the melt, i.e. , , and , whereas the inner reactant is the silicate melt. The product is the plagioclase in the orbicular granite.
2.3 Numerical analysis
This section presents the numerical setup and finite element discrete formulations to apply the proposed phase separation model to simulate periodic patterns in stones and rocks.
2.3.1 Problem statement
We examine a two-dimensional initial-boundary value problem defined within a square domain as displayed in Figure 3, where the boundary of the domain is represented by and is associated with a unit outward normal vector . The domain size is 100 or 200 depending on the single pattern formation or multiple pattern interactions, which is sufficiently large to present different categories of geological patterns observed in rocks. Table 1 lists the parameters used for all simulations. We impose Neumann boundary conditions on four edges for all variables as.
FIGURE 3
TABLE 1
| Property | Symbol | Value |
|---|---|---|
| Diffusion coefficient | 1.0 | |
| Diffusion coefficient | 1.0 | |
| Reaction rate | 1.0 | |
| mobility | 1.0 | |
| Parameter | 1.0 | |
| Parameter | 0.2 | |
| Gradient parameter | 0.2 | |
| Noise effect parameter | 0 |
Material properties for liesegang band simulation.
Initially, the inner electrolyte , characterized by a higher concentration , is confined to specific geometric configurations within inner domains, while the outer electrolyte , with a comparatively lower concentration , occupies the remaining regions. In this study, diverse shapes and configurations are employed for the inner domains, e.g. circular domain in Figure 3. These geometries are selected to reflect configurations that are either observable in natural settings or amenable to experimental validation. There are two groups of initial geometries in this studies depending on area of reactants in the entire domain. First group contains a single geometry occupied by the reactant to study the single pattern forming process in Section 3.1, while second categories have several inner areas filled by the reactants to examine the complex pattern interaction in Section 3.2. It should be noted that the dimensions and positions of the inner geometries will be explicitly defined in the corresponding section where each specific case is analyzed. Overall, the system is assumed to be uniform and constant within the stable region, expressed asThe scaled initializations used are , , , and in our study.
2.3.2 Finite element discretization
The coupled nonlinear partial differential equation (Equation 9) is solved using the open-source, high-performance phase-field software PRISMS-PF [34], based on the finite element method (FEM) through the deal.II library [35]. This framework supports adaptive mesh refinement, massively parallel computations, and matrix-free finite element simulations—features particularly advantageous for resolving fine-scale pattern evolution under complex geometric constraints. In contrast, [27] employed a vertex-based finite volume method (vb-FVM) to simulate Liesegang-type patterns in arbitrary domains, emphasizing flexibility in handling unstructured meshes and geometric confinement. While both FEM and vb-FVM are suitable for simulating reaction–diffusion systems on non-uniform grids, FEM offers superior accuracy in capturing sharp interfaces and gradient-driven dynamics due to its variational formulation and higher-order basis functions. Additionally, the matrix-free implementation in PRISMS-PF enhances computational efficiency, enabling scalable simulations of complex, heterogeneous nucleation scenarios that are central to our study.
The weak formulation of the coupled Equation 9, incorporating the Neumann boundary Equation 14 and initial conditions Equation 15, is derived using the following trial and test functions in Equation 16:and integration by parts. By adopting the forward Euler method for time discretization, the corresponding weak forms become
The subsequent numerical analysis is conducted using an implementation of Equation 17 within the PRISMS-PF framework. The computational code is parallelized utilizing the Message Passing Interface (MPI) from the PETSc library (available at https://www.mcs.anl.gov/petsc). To enhance computational efficiency, a matrix-free finite element framework and an adaptive meshing technique are utilized, offering superior performance compared to conventional finite element methodologies.
3 Results and discussion
The numerical algorithm is developed and validated across various domains where the reactants ( and ) are initially separated. Parameters are selected to generate Liesegang band patterns, which are characteristic of such systems. This section investigates both single Liesegang pattern formation and complex Liesegang pattern interaction. In this investigation, the initial geometries are categorized into two groups based on the spatial distribution of reactants within the domain. The first group comprises configurations with a single region occupied by reactant , enabling the analysis of individual pattern formation processes, as discussed in Section 3.1. The second group involves multiple distinct regions filled with reactant , facilitating the examination of complex interactions between patterns, as detailed in Section 3.2.
3.1 Single pattern formation
3.1.1 Circular geometry
Figure 4a gives one sample of concentrically zoned pyrite aggregates (CZPA) of the pyritic ores in the Xingqiao deposit. This zoning pattern can be explained by Liesegang diffusion-reaction through the interaction between the diffusion of soluble reactants and the precipitation of solids in a permeable medium. Specifically, from hydrothermal fluids infiltrates early diagenetic sediments with high porosity and then reacts with in pore fluids, which are produced by bacterial sulfate reduction during organic matter degradation. Consequently, the outer solute , inner solute , and product can be generalized into the form . We design a numerical setup and use our proposed phase-separation model to simulate the forming processing in Figure 4b, where a high concentration of diffuses into a -bearing permeable medium.
FIGURE 4
Figure 4c illustrates the generated patterns from our post-nucleation phase-separation model. One can observe that our simulation can reproduce the oscillatory zoning observed in the CZPA of the pyritic ores in the Xingqiao deposit. Figure 4a depicts the location and band widths of the selected CZPA. We assume that the diffusion direction of comes from the center nucleation point to the outer space, indicating a reverse Liesegang phenomenon [16]. It is worth noting that [16] measured the location and widths from the opposite direction and obtained a positive Liesegang ringing. Owing to the limited geological evidence, we cannot determine the real diffusion direction the same as the subjective selection in [16]. Our findings provide evidence that the CZPA in the Xingqiao deposit has reverse Liesegang ringing features. Later, we shall investigate heterogeneous distribution of the pyrite rings in the CZPA with irregular initialization.
3.1.2 Triangle geometry
We subsequently compare the primary features captured by our simulations with those observed in a specific Zebra rock sample, as depicted in Figure 5a. One can observe a triangular-shaped Liesegang band pattern diffuses from the lower-left region throughout the entire sample. We postulate the presence of a fluid channel in the core region, which supplies Fe-bearing fluid flow, thereby generating the rhythmic hematite. Upon examination, two primary deviations from classical Liesegang patterns are noticeable. Firstly, the band thickness decreases from the core to the rim, and secondly, the distance between adjacent bands diminishes as the diffusion process progresses. To replicate these counterintuitive phenomena, we devised a numerical setup featuring a triangular-shaped initialization region for the Fe-bearing fluid, as illustrated in Figure 5b. In this configuration, the constant outer reactants diffuse from the center to the periphery, interact with the dolomite, and eventually form periodic hematite.
FIGURE 5
In Figure 5c, we present the spatial concentration distribution of the precipitate . The resultant pattern closely resembles that observed in the Zebra rock sample, with the band thickness decreasing as the Fe-bearing fluid flows from the center to the boundaries. This contrasts with the increasing band width observed in the wake of the Fe-bearing diffusion front. The primary reason for this discrepancy is that we initialize the constant Fe-bearing fluid in the triangular region rather than continuously supplying the outer reactant source, as in our previous setup [10, 26, 36]. In our simulation, the limited reactants give rise to a counterintuitive phenomenon that diverges from the traditional Liesegang pattern. Furthermore, the Zebra rock sample exhibits a decreased distance between periodic bands behind the diffusion front, a phenomenon known as reversed Liesegang patterns. Similarly, our simulation reveals a slight reduction in the distance between neighboring bands, in particular when the bands are away from the core. Overall, the simulated pattern closely approximates the field observations in the Zebra rock sample, despite the two primary deviations from the standard Liesegang pattern. It is important to note that while our phase separation model successfully reproduces the non-standard Liesegang pattern observed in the field, further geological investigations are required to refine our numerical simulations.
3.2 Multiple pattern formation and interaction
Apart from the single Liesegang ring formation, we examine our phase separation model for multiple Liesegang ring formation. Specifically, we first present the two and multiple orbicular granite formations in igneous systems. We then model the irregular distribution of the pyrite rings formation in the Xiaoqing deposit.
3.2.1 Multiple circular geometries
Many previous hypotheses regarding the pattern formation in orbicular granite have assumed that crystallization initiates at the core, with subsequent band formation propagating outward [2, 37]. The presence of orbicular granites also indicates an interaction between different connected orbicular rings. Figures 6b, 7b illustrate the interactions of two and nine orbicular granites within igneous systems, respectively. One difference is the appearance of morphological patterns; Figure 6b shows two oval Liesegang rings, whereas Figure 7 shows more circular Liesegang patterns. To distinguish the oval and circular morphological occurrence in orbicular granite, we conduct two numerical simulations. One simulation includes two oval cores with initialized metal melt, while another contains nine circular cores. The metal melt is capable of diffusing and reacting with the silicate melt, thereby forming rhythmic orbicular rings, as described by Equation 13. Detailed numerical geometries and boundary conditions are provided in Figures 6b, 7b, respectively.
FIGURE 6
FIGURE 7
The resultant pattern morphologies for the interaction of two orbicular granites are depicted in Figure 6c. It is observed that the thickness of the bands is greater near the core, which is consistent with observations of natural orbicular granites. Additionally, the distance between neighboring orbicular rings increases from the core to the rim. Notably, when the orbicular rings approach one another, the formation of one ring is impeded by the presence of another. This simulated interaction behavior qualitatively mirrors field observations of orbicular granites.
To further validate our phase separation model, we increase the number of initial cores from two to nine and change the initialized shape from ovals to circles, as shown in Figure 7b. A significant difference from the interaction of two orbicular rings is the occurrence of the Ostwald ripening phenomenon. Several spots between the connected orbicular rings can be observed in the field, as illustrated in Figure 7a. This phenomenon is also captured in our simulation, where only the outer rings interact with rings from other orbicular structures. Overall, while further quantitative analysis is warranted, the primary features observed in the field can be replicated in our numerical simulations of orbicular pattern formation.
3.2.2 Randomly distributed circular geometries
In Section 3.1.1, we simulate the formation of zoning rings using a single circular pattern. However, the CZPA typically exhibits more complex distributions within the depositional environment, which suggests a heterogeneous distribution of concentrations and spatial positions, as illustrated in Figure 8a. To further investigate this complexity, we have developed a more complex numerical model to simulate these complex zoning patterns, as shown in Figure 8b. In this model, we randomly distribute 15 seeds throughout the domain to replicate the heterogeneous distribution. The specific locations and sizes of these seeds are detailed in Table 2. The boundary conditions and geometric configurations are consistent with those used for multiple circular tests described in Section 3.2.1.
FIGURE 8
TABLE 2
| Region | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 30 | 30 | 40 | 60 | 60 | 60 | 90 | 90 | 90 | 120 | 120 | 150 | 150 | 180 | 180 | |
| 30 | 54 | 160 | 10 | 90 | 120 | 50 | 90 | 150 | 90 | 180 | 20 | 150 | 50 | 100 | |
| Radius | 5 | 7 | 9 | 10 | 15 | 13 | 9 | 15 | 7 | 6 | 10 | 11 | 8 | 12 | 6 |
The coordinates and sizes of the circular initialization regions of the Fe-bearing fluid. The circular region is ordered by increasing -coordinate, and in cases of identical -coordinates, by increasing -coordinate.
Figure 8c displays the simulated irregular distribution of pyrite rings. These irregular pyrite patterns closely resemble the field observations of pyrite distribution depicted in Figure 8a. A key question concerning the CZPA at Xinqiao is whether it formed from a core outward or from a shell inward. The presence of bright cores suggests the former, while dark cores imply the latter. Qiu and colleagues have proposed two different setups to replicate both the dark and bright cores observed in CZPA. One hypothesis involves the diffusion of into the sediment forming a pyrite shell, with diffusing outward to this shell, thus creating a lower concentration of pyrite in the center, resulting in dark cores. The other hypothesis posits a constant concentration of at the center of the domain, which reacts with to form CZPA with bright cores. Notably, our simulation can reproduce both dark and bright cores simultaneously without altering the numerical setup. In our model, we assume that the highly concentrated outer electrolyte typically diffuses outward from the cores and reacts with the , which initially was evenly distributed in the porous sediment.
4 Conclusion
This study employs a phase separation model to examine the rhythmic patterns in rocks, with a particular focus on oscillatory zoning in pyrite, hematite patterns in Zebra rocks, and orbicular granites. The simulation results elucidate the influence of geometric heterogeneity on the development of both simple and complex periodic patterns in these three rock types, considering both regular and randomly distributed initial conditions. The simulated patterns qualitatively replicate the Liesegang patterns observed in natural settings, thereby demonstrating the capacity of the phase separation model to reproduce periodic rock patterns. Nonetheless, further research is necessary, particularly through detailed geochemical analyses of these patterns and quantitative comparisons between the simulated patterns and those observed in the field. This study establishes a foundational framework for utilizing the phase separation model to analyze self-organized patterns in heterogeneous geological contexts.
Statements
Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
Author contributions
CL: Validation, Conceptualization, Visualization, Methodology, Formal Analysis, Investigation, Writing – review and editing, Writing – original draft, Software. KR-L: Supervision, Funding acquisition, Writing – review and editing. MH: Funding acquisition, Writing – review and editing, Supervision.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. We acknowledge the funding support from the Research Grant Council of Hong Kong (MH: ECS 27203720 and GRF 17206521) and the Australian Research Council (KR.L.: ARC DP200102517, DP170104550, LP170100233).
Acknowledgments
The computations were performed using the research facilities offered by Information Technology Services, the University of Hong Kong.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
References
1.
OrtolevaPJMerinoEMooreCChadamJ. Geochemical self-organization I; reaction-transport feedbacks and modeling approach. Am J Sci (1987) 287:979–1007. 10.2475/ajs.287.10.979
2.
L’HeureuxI. Self-organized rhythmic patterns in geochemical systems. Philosophical Trans R Soc A: Math Phys Eng Sci (2013) 371:20120356. 10.1098/rsta.2012.0356
3.
YatsudaYTsushimaKFangQNabikaH. Chemical model for pattern formation in rocks via periodic precipitation of iron oxide minerals. ACS Earth Space Chem (2023) 7:2042–9. 10.1021/acsearthspacechem.3c00165
4.
WalimbePItataniMKulkarniPLagziIKulkarniS. A quasi universal matalon–packter law for a periodically precipitating system of iron (II) hydroxide involving volumes and concentrations of the invading electrolyte. Langmuir (2023) 39:13420–9. 10.1021/acs.langmuir.3c00778
5.
CowardAJBruggerJWilsonSSlimACWilliamsTPillansBet alMineralogy and geochemistry of pattern formation in print stone from the Pilbara, Australia. Mineralogical Mag (2025) 1–33. 10.1180/mgm.2025.11
6.
LiesegangRE. Eine scheinbar chemische Fernwirkung. Annalen der Physik (1906) 324:395–406. 10.1002/andp.19063240210
7.
WangYChanMAMerinoE. Self-organized iron-oxide cementation geometry as an indicator of paleo-flows. Scientific Rep (2015) 5:10792–15. 10.1038/srep10792
8.
YoshidaHKatsutaNSironoSNishimotoSKawaharaHMetcalfeR. Concentric fe-oxyhydroxide bands in dacite cobbles: rates of buffering chemical reactions. Chem Geology (2020) 552:119786. 10.1016/j.chemgeo.2020.119786
9.
KawaharaHYoshidaHYamamotoKKatsutaNNishimotoSUmemuraAet alHydrothermal formation of Fe-oxide bands in zebra rocks from northern Western Australia. Chem Geology (2022) 590:120699. 10.1016/j.chemgeo.2021.120699
10.
LiuCCaloVRegenauer-LiebKHuM. Coefficients of reaction-diffusion processes derived from patterns in rocks. J Geophys Res Solid Earth (2023) 128:e2022JB026253. 10.1029/2022jb026253
11.
LiuCCaloVMRegenauer-LiebKHuM. Deriving flow velocity and initial concentration from Liesegang-like patterns. J Geophys Res Solid Earth (2024) 129:e2024JB029379. 10.1029/2024jb029379
12.
LevesonDJ. Orbicular rocks: a review. Geol Soc America Bull (1966) 77:409–26. 10.1130/0016-7606(1966)77[409:orar]2.0.co;2
13.
WangYMerinoE. Oscillatory magma crystallization by feedback between the concentrations of the reactant species and mineral growth rates. J Petrology (1993) 34:369–82. 10.1093/petrology/34.2.369
14.
GhaouiJL’HeureuxI. Model of layered pattern formation in binary igneous systems. Solid Earth Sci (2021) 6:80–94. 10.1016/j.sesci.2021.04.001
15.
OrtolevaPJ. Self-organization and nonlinear dynamics in sedimentary basins. Philosophical Trans R Soc Lond Ser A: Phys Eng Sci (1993) 344:171–9. 10.1098/rsta.1993.0085
16.
QiuWJZhouMFWilliams-JonesAE. Numerical simulation of the self-organizational origin of concentrically zoned aggregates of siderite and pyrite in sediment-hosted massive sulfide deposits. J Geophys Res Solid Earth (2024) 129:e2023JB028101. 10.1029/2023jb028101
17.
OrtolevaPChadamJMerinoESenA. Geochemical self-organization II; the reactive-infiltration instability. Am J Sci (1987) 287:1008–40. 10.2475/001c.60405
18.
OstwaldWLehrbuch der Allgemeinen Chemie, vol. 1 (1897).
19.
OstwaldWLehrbuch der allgemeinen Chemie, vol. 2 (1902).
20.
LifshitzIMSlyozovVV. The kinetics of precipitation from supersaturated solid solutions. J Phys Chem Sol (1961) 19:35–50. 10.1016/0022-3697(61)90054-3
21.
AntalTDrozMMagninJRáczZ. Formation of liesegang patterns: a spinodal decomposition scenario. Phys Rev Lett (1999) 83:2880–3. 10.1103/physrevlett.83.2880
22.
AntalTDrozMMagninJPekalskiARáczZ. Formation of liesegang patterns: simulations using a kinetic ising model. The J Chem Phys (2001) 114:3770–5. 10.1063/1.1342858
23.
HantzPBiróI. Phase separation in the wake of moving fronts. Phys Rev Lett (2006) 96:088305. 10.1103/physrevlett.96.088305
24.
ThomasSLagziIMolnárJFRáczZ. Probability of the emergence of helical precipitation patterns in the wake of reaction-diffusion fronts. Phys Rev Lett (2013) 110:078303. 10.1103/physrevlett.110.078303
25.
DayehMAmmarMAl-GhoulM. Transition from rings to spots in a precipitation reaction–diffusion system. RSC Adv (2014) 4:60034–8. 10.1039/c4ra11223g
26.
LiuCHuMRegenauer-LiebK. Liesegang patterns interpreted as a chemo-hydromechanical instability. In: Proceedings of 12th international workshop on bifurcation and degradation in geomechanics. Springer (2022). p. 59–66. 10.1007/978-3-031-22213-9_7
27.
Abi MansourAAl-GhoulM. Vertex-based finite volume simulation of liesegang patterns on structureless meshes. Phys Rev E (2014) 89:033303. 10.1103/physreve.89.033303
28.
Al-GhoulMSultanR. Simulation of geochemical banding: theoretical modeling and fractal structure in acidization-diffusion-precipitation dynamics. Phys Rev E (2019) 100:052214. 10.1103/physreve.100.052214
29.
CahnJWHilliardJE. Free energy of a nonuniform system. I. Interfacial free energy. The J Chem Phys (1958) 28:258–67. 10.1063/1.1744102
30.
CahnJW. On spinodal decomposition. Acta Metallurgica (1961) 9:795–801. 10.1016/0001-6160(61)90182-1
31.
RickardDLutherGW. Chemistry of iron sulfides. Chem Rev (2007) 107:514–62. 10.1021/cr0503658
32.
SchoonenMBarnesH. Reactions forming pyrite and marcasite from solution: I. Nucleation of FeS2 below 100°C. Geochimica et Cosmochimica Acta (1991) 55:1495–504. 10.1016/0016-7037(91)90122-l
33.
FeeneyRSchmidtSStrickholmPChadamJOrtolevaP. Periodic precipitation and coarsening waves: applications of the competitive particle growth modela. The J Chem Phys (1983) 78:1293–311. 10.1063/1.444867
34.
DeWittSRudrarajuSMontielDAndrewsWBThorntonK. PRISMS-PF: a general framework for phase-field modeling with a matrix-free finite element method. npj Comput Mater (2020) 6:29. 10.1038/s41524-020-0298-5
35.
ArndtDBangerthWBlaisBClevengerTCFehlingMGrayverAVet alThe deal.II library, version 9.2. J Numer Mathematics (2020) 28:131–46. 10.1515/jnma-2020-0043
36.
LiuCCaloVMRegenauer-LiebKHuM. Dendritic growth patterns in rocks: inverting the driving and triggering mechanisms. J Geophys Res Solid Earth (2023) 128:e2023JB027105. 10.1029/2023jb027105
37.
GrossePToselliAJRossiJN. Petrology and geochemistry of the orbicular granitoid of Sierra de Velasco (NW Argentina) and implications for the origin of orbicular rocks. Geol Mag (2010) 147:451–68. 10.1017/s0016756809990707
Summary
Keywords
liesegang pattern, self-organization, geometry heterogeneity, zebra rock, orbicular granite, pyrite rock
Citation
Liu C, Regenauer-Lieb K and Hu M (2025) Influence of geometry heterogeneity on Liesegang patterns in rocks. Front. Phys. 13:1628328. doi: 10.3389/fphy.2025.1628328
Received
14 May 2025
Accepted
16 July 2025
Published
11 August 2025
Volume
13 - 2025
Edited by
Rabih Sultan, American University of Beirut, Lebanon
Reviewed by
Mazen Al-Ghoul, American University of Beirut, Lebanon
Masaki Itatani, Hokkaido University, Japan
Mohd Hariri Arifin, National University of Malaysia, Malaysia
Updates
Copyright
© 2025 Liu, Regenauer-Lieb and Hu.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Chong Liu, chongliu@connect.hku.hk
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.