Physics of the Earth and Planetary Interiors 354 (2024) 107227

Contents lists available at ScienceDirect

Physics of the Earth and Planetary Interiors
journal homepage: www.elsevier.com/locate/pepi

Simulation of microtextural evolution in omphacite: Orderingtransformation kinetics as unexplored archives of slab eclogitization
Ryo Fukushima a, * , Tatsuki Tsujimori a, b , Nobuyoshi Miyajima caGraduate School of Science, Tohoku University, Sendai 980-8578, Japan
Center for Northeast Asian Studies, Tohoku University, Sendai 980-8576, Japanc
Bayerisches Geoinstitut, Universität Bayreuth, 95440 Bayreuth, Germanyb

A R T I C L E I N F O

A B S T R A C T

Keywords:
Omphacite
Ordering transformation
Antiphase domains (APDs)
EclogiteTransmission electron microscopy (TEM)
Phase-field approach

Earth’s subduction zone processes and surface environments are intricately governed by mass transfer phe­nomena at plate convergent boundaries. The determination of their rates and timings from high-pressure
metamorphic rocks (e.g., eclogite), or remnants of ancient convergent boundaries, remains an ongoing chal­lenge. Here, we proposed the potential and versatility of ordering transformation kinetics of omphacite, anessential mineral found in eclogite, as a dynamic recorder of the metamorphic history. Through macroscopic
phase-field simulation, we explored the growth of antiphase domains (APDs) in metastable disordered ompha­
cite, discussing the feasibility of constraining metamorphic reaction kinetics based on the size and morphology ofomphacite APDs in eclogitized oceanic crust. Our simulation corroborated that omphacite nucleating later during
the prograde metamorphism can exhibit an incompletely ordered state with sparsely distributed ordered do­mains, which suggests their usefulness in estimating the recrystallization timing of the omphacite. Additionally,
we confirmed that the APD formation dynamics are significantly influenced by the initial cation configuration of
the disordered matrix. This implies the APD morphology in natural omphacite under slab-surface conditions may
reflect their precipitation kinetics. These findings provide valuable insights into the microtextural evolution ofomphacite due to its ordering transformation, thereby enhancing our ability to interpret morphological features.

1. Introduction
Constraining the timings and rates of mass-transfer phenomena at
Earth’s plate convergent boundaries is of considerable importance tounderstand solid-earth processes and dynamic changes of surface envi­
ronments. In addition to geophysical observations of modern subduction
zones, traditional petrological approaches, combined with radiometricage dating of rocks that recrystallized at ancient plate boundaries, or
high-pressure metamorphic rocks such as eclogite, have provided someconstraints on such chemical geodynamics at convergent zones. How­ever, there is an insurmountable barrier to the application of thistraditional methodology, as the temporal resolution of the metamorphic
mineral ages remains at ~105–6 years (e.g., porphyroblastic garnet;
Pollington and Baxter, 2010; Gorce et al., 2021; Tual et al., 2022). This
limitation may hinder the detection of their nucleation–growth timings,
which are controlled by shorter-span fluid flushing or metasomatism (e.
g., Yokoyama et al., 2002; Viete et al., 2018; Broadwell et al., 2019).
* Corresponding author.
E-mail addresses: ryo.fukushima.p7@dc.tohoku.ac.jp(N. Miyajima).

(R.

Fukushima),

Hence, there is an evident need to more precisely resolve the reaction
kinetics and burial durations of the individual minerals.
To overcome this barrier, we revisited the ordering transformation
kinetics of omphacite, a primary constituent of eclogite. Over the pasthalf century, numerous studies have provided insights into the size,
morphology and formation processes of antiphase domains (APDs)
within omphacite—a microtexture formed by cation ordering (Champ­
ness, 1973; Phakey and Ghose, 1973; Carpenter, 1978, 1979a, 1979b,
1979c, 1982a, 1982b; Carpenter and Smith, 1981). The ordering, with aspace group change from C2/c to P2/n, is described by splitting of theindividual M-cation sites, M1 for (Mg, Fe2+, Al, Fe3+) and M2 for (Na,Ca), into M1/M11 and M2/M21, respectively. Due to the breakdown oftranslational symmetry, antiphase boundaries (APBs) appear to form the
APDs. Focusing on the aspects at which disordered omphacite firstmetastably appears, even at lower temperatures than the critical tem­
perature of ~865 ◦ C (Champness, 1973; Carpenter, 1978; Carpenter andOkay, 1978; Carpenter and Putnis, 1985) and at which APDs generally

tatsukix@tohoku.ac.jp(T.

Tsujimori),

nobuyoshi.miyajima@uni-bayreuth.dehttps://doi.org/10.1016/j.pepi.2024.107227
Received 25 December 2023; Received in revised form 9 May 2024; Accepted 4 July 2024
Available online 9 July 2024
0031-9201/© 2024 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/).

R. Fukushima et al.

Physics of the Earth and Planetary Interiors 354 (2024) 107227

coarsen due to thermal annealing, Carpenter (1982b) proposed anempirical geothermometer/geospeedometer using the mean size ofequiaxed APDs (δ), based on a positive correlation between APD sizes
and peak metamorphic temperatures in natural omphacites:()
Qδ8 = 8K0 exp −
t,
(1)
RT

textures, one cannot clearly define sharp APBs at least from thedark-field images, posing a question whether this assumption is valid in
all cases. Besides, almost fully ordered omphacites (Fig. 1c,d) and thosewith a plenty of columnar APDs (Fig. 1c) sometimes coexist with suchpossibly metastable, disordered omphacites, leading to a difficulty ininterpreting their formation histories. Strictly speaking, if incompletelyordered omphacites can remain in low-T eclogites, a simple applicationof the ripening theory without considering the incipient ordering stage
could be insufficient to estimate the annealing timescale.
Therefore, to establish a reasonable indicator of T–t histories basedon the omphacite microtexture, we should first understand how cation
ordering affects the final morphology of APDs, for a further calibration
with natural samples. It may prove challenging to reproduce this cation
ordering process with actual annealing experiments due to its sluggish
kinetics (Carpenter, 1982b) and, as such, numerical simulation is themost effective way to understand the whole picture of the microtexture
formation owing to cation ordering. To simulate the morphologicalchanges of ordered domains, phase-field modeling has been widely used
and recognized as a reliable approach in material science (e.g., Vai­thyanathan and Chen, 2000; Zhu et al., 2004). This method basicallyrequires the excess free energy function against one or several param­
eters to describe the local order-disorder state (and/or chemical com­positions), as well as the other excess energy term(s) derived from thegradient of each parameter field. In the case of omphacite, Carpenter
et al. (1990) proposed a tricritical Landau expansion for the excess freeenergy as a function of the order parameter, which can be directly
incorporated into the phase-field model. Thus, the application of phasefield modeling to the omphacite system is a reasonable choice.
In this study, we use a phase-field approach to investigate how cation
ordering affects the incipient APD formation in incompletely orderedomphacite. By leveraging a previously proposed expression of its excess
free energy, we replicate the morphology and size of ordered domains inomphacite mainly under low-T eclogite-facies conditions. This is aimed
at determining the feasibility of extracting relevant temporal informa­tion to constrain their nucleation timings from observable micro-tonanoscale textures in natural omphacites.

where K0 is a constant (6 × 1035 Å8 yr− 1), Q is the activation energy(3–4 × 105 J mol− 1), R is the gas constant, T is the annealing temper­
ature, and t is the annealing time. Several researchers have discussed thepeak metamorphic temperatures and metamorphic timescales of eclo­gites using this T–t–δ relationship (Lardeaux et al., 1986; Wu et al., 2000;
Brenker et al., 2003; Müller et al., 2011; Pollok et al., 2014; Xie et al.,2020).
However, this T–t–δ relationship is, obviously, not precise enough toestimate detailed T–t paths, because this was only roughly calibrated sothat mean equiaxed APD sizes in most of the observed samples were
within a range of the geological timescale (104–8 years). In addition, wemust reevaluate whether this relationship is truly applicable, especiallyfor those formed at a depth of <~80 km along the typical subductionzone thermal gradients (i.e., low-T eclogites: Tsujimori and Mattinson,2021). Although Eq. 1 was established under the assumption that APD
coarsening is solely driven by the reduction of the APB surface area, wehave confirmed the presence of omphacite patches with heteroge­
neously bright ordered clusters (Fig. 1a) and with isolated ordered do­mains (Fig. 1b,c) in low-T eclogites (Fukushima et al., 2021b). For these

2. Model setting
We adopted a macroscale, phenomenological phase-field model witha single order parameter for the omphacite M1 site, assuming that theordering on the M2 site is linked to that on the M1 site (e.g., Ross, 1988;
Carpenter et al., 1990). By solving a simple evolution equation with thepreviously proposed Landau expansion of the excess free energy in theisotropic system, we were able to simply examine how cation ordering
proceeds in initially disordered omphacite with the ideal chemicalcomposition in accordance with various input parameters. It is worth
noting that the anisotropy of the gradient energy is disregarded here dueto the rare occurrence of crystallographically oriented APBs in ompha­
cite (Carpenter, 1979c). We also dismissed the presence of any compo­
sitional variation and lamellae within omphacite, and this allowed us to
set the 4 principal parameters: (1) annealing temperature (T); (2) excessenergy for sharp APBs (equilibrium APB energy); (3) the initial scalar
field for the ‘almost’ disordered omphacite, which must be properly
defined with some fluctuation for ordering to proceed; and (4) mobility(M) that relates real time (t) and normalized time for calculations (t*) as
t* = MRTt. However, we further regarded the equilibrium APB energy,
whose value is unknown for omphacite so far, as a constant for
simplicity. This is because the possible temperature range for low-T
eclogite formation is limited to ~450–600 ◦ C, and because it might beaffected by the compositional variation in such a narrow temperature
range, as in the case of alloys (e.g., Yu et al., 1994; Gorbatov et al.,2016). Further model details and calculation procedures are described
below.

Fig. 1. Dark-field TEM images (g = 050) of natural omphacites from low-T
eclogites. a, omphacite included in a garnet porphyroblast in a lawsonite
eclogite from the South Motagua Mélange (SMM), Guatemala. Xaeg = 0.08; Xaug
= 0.48; Mg# = 0.71; FeT = 0.23. b, omphacite with isolated equiaxed domains
included in a garnet porphyroblast in an epidote-glaucophane eclogite from
Syros, Greece (enlarged image of Fig. 6e in Fukushima et al., 2021b). The mean
APD size (δ) and the ordered domain fraction (F) are ~100 nm and ~0.5,
respectively. Xaeg = 0.00; Xaug = 0.47, Mg# = 0.66; FeT = 0.19. c, anotheromphacite inclusion with both equiaxed (dashed square) and columnar (solid
square) APDs from the same garnet as b. Grt = garnet, Omp = omphacite.
Representative compositions of the areas 1–3 are: Xaeg = 0.09, 0.24, 0.04; Xaug
= 0.50, 0.48, 0.52; Mg# = 0.70, 0.62, 0.70; FeT = 0.25, 0.44, 0.21. d, typicalequiaxed APDs in a matrix-forming omphacite in epidote eclogite from Omi,
Japan. The mean APD size and the ordered domain fraction are ~140 nm and~0.8, respectively. Xaeg = 0.02; Xaug = 0.54, Mg# = 0.77; FeT = 0.16.
2

R. Fukushima et al.

Physics of the Earth and Planetary Interiors 354 (2024) 107227

2.1. Phase-field model

fully ordered antiphase domains (Cahn and Hilliard, 1958):√̅̅̅ ∫ ϕm √̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
2agbulk (ϕ, T) + gm dϕ.
s=v − ϕmOwing to the simplification above, we can employ a general method
for simulating the dynamics of the order-disorder phase transition withincipient ordering as in Ninomiya et al. (1990). Basically, our approachis the same as their method, except for the form of the excess free energyfunction. All the symbols are summarized in Table 1. We introduced asingle phase-field parameter, ϕ(r, t), whose value can range from − 1 to
+1. This parameter represents the degree of order for M1 sites (by its
absolute value) and the relative phase (by its positive/negative sign) at a
position r and time t. The total free energy of the system (Gtotal ) is given
by:
}
∫ {
1
a2
(2)
Gtotal [ϕ(r, t) ] =gbulk (ϕ, T) + |∇ϕ|2 dr,
2rv

2.2. Calculation procedures
We basically conducted simulations in a 2-dimensional (2D) spacewith a 2 × 2 μm2 area (400 × 400 grids) and employed periodic
boundary conditions. Additionally, we attempted 3-dimensional (3D)
simulations in smaller spaces (0.6 × 0.6 μm2 area with 120 × 120 grids)
to discuss the dimensional effects. The 3D results were compared withthe 2D results by generating a 2D slice for each time step.
As the initial configurations, zero-mean Gaussian noise with a vari­ance of σ2 was independently introduced for each grid, each with a sizeof 5 nm on a side (dx, dy, dz). It is important to note that the grid size
should be larger than the unit cell size to robustly define the macro­scopic order parameter and the bulk free energy for each grid. In addi­tion, we can assume that the approximate conservation length
corresponds to the interatomic distance between neighboring crystal­lographic sites (~0.5 nm) (Carpenter and Salje, 1989). This suggests thatthe initial representative degree of order in a given grid should be in­dependent of those in the adjacent grids. Therefore, to simplify matters,
we discussed the initial cation configuration in terms of the fluctuation
amplitude, based on the working hypothesis that the cation configura­tion of the metastable disordered omphacite can be essentially treated as
the assembly of small portions with independent noise.
For the discretization of Eq. 5, we used a forward difference at time t*
and a second-order central difference for the space derivatives. The
time-step size (dt) was fixed at 0.01 to ensure stable calculation results.
The default value for the APB energy (s = 0.50 J m− 2) was determined sothat simulation results at a temperature of 550 ◦ C replicated the realequiaxed APD texture shown in Fig. 1b, assuming it is incompletelyordered. By tentatively setting the square root of the initial fluctuation
variance as σ = 10− 4, we adjusted the s value to reproduce a mean APD
size (δ) of ~100 nm when the ordered phase fraction (F) is ~0.5 in the
2D model. The order of magnitude of this s value, namely 10− 1 J m− 2, isconsistent with those of anorthite (Carpenter, 1994). For the dark-fieldlike image conversion, we assumed a response curve of the CCD camera
with a certain saturation level (Supplementary Fig. S1), while ignoring
completely buried ordered domains in the specimen. It is worth notingthat we tentatively speculated a convex-upward function for slightly
ordered states (|ϕ| < 0.1) to replicate the contrast of the real dark-field
images. For |ϕ| ≥ 0.1, where the sharp superlattice reflection due to thelong-range order should be, at least, observed (Boffa-Ballaran et al.,1998), we used a convex-downward parabolic function, as implemented
in the previous study (e.g., Ninomiya et al., 1990).
To quantitatively evaluate the individual snapshots of the evolving
phase field, we calculated the mean degree of order (spatial average of|ϕ|), the mean APD size (δ), the ordered phase fraction (F), and thespatial wavelength after every 50 time steps. For some representative
cases, we also calculated histograms of |ϕ| to visualize the temporalchange of its frequency distribution (i.e., ‘power spectrum’ as set forth in
Salje and Wruck (1988)). Since the calculation of δ and F requires
binarization of the corresponding phase field, we set a tentative
threshold of |ϕ| = 0.4 to convert it to a binarized image. The δ value wasobtained by averaging all the ordered segment lengths along 50 nmspacing lines. For the spatial wavelength calculation, we utilized amethod for estimating the averaged grain size using the autocorrelo­gram of the sediment image (Buscombe et al., 2010). By calculating the
autocorrelation of the phase field (before binarization) for each lag
along one direction, we obtained its representative length scale.
All the simulations were implemented with a set of Python routines
called Om-COS (Fukushima, 2023). We ran the scripts on the Windows
11 operating system using Python 3.9, with an AMD Ryzen™ 7 5700 G
processor and 16 GB RAM.

where v is a constant equivalent to the molar volume of omphacite(6.3295 × 10− 5 m3 mol− 1; the average of jadeite and diopside values(Holland and Powell, 2011)), gbulk (ϕ, T) is the temperature-dependentbulk excess free energy of omphacite assuming a uniform degree oforder (in J mol− 1), and a is a parameter that determines the isotropic
gradient energy (in m J0.5 mol–0.5). The gbulk (ϕ, T) is derived from theLandau-type expansion for omphacite whose composition is jd50 aug50
by Carpenter et al. (1990):(3)

gbulk (ϕ, T) = 11.4(T − 1138)ϕ2 + 4317ϕ6 ,where T is the annealed temperature in K. At a constant temperature, theminimum value of the gbulk (ϕ, T) ( = − gm ) can be achieved when ϕ = ±√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
4T)
ϕm = ± 22.8(1138−, which correspond to the equilibrium states. As
25902

ϕ(r, t) should be approximated as a non-conserved field (Carpenter andSalje, 1989), we considered the Allen-Cahn equation (Allen and Cahn,
1979):∂ϕδGtotal= − Mv
,
∂t
δϕ

(4)

where M is the mobility in mol J− 1 yr− 1. Dividing both sides by MRT, we
obtained the evolution equation to be solved:}∂ϕ
1 {22.8(T − 1138)ϕ + 25902ϕ5 − a2 ∇2 ϕ ,
=−
RT
∂t *

(5)

where t * is the normalized time defined as t* = MRTt. Additionally, theisotropic gradient energy coefficient a was related to the APB energy (sin J m− 2) by considering the free energy of a flat interface between two

Table 1
List of the symbols regarding the simulation.
symboldescription

δ
R
T
t
M
t*
F

mean equiaxed APD sizegas constantannealing temperaturetime
mobility
normalized timeordered phase fraction
standard deviation of the initial
fluctuation
APB energy
phase-field parameter
total free energy of the system
molar volume of omphacite
isotropic gradient energy coefficient
excess free energy of omphacite
minimum of the gbulk
degree of order when gbulk = − gm
grid size
time step sizeσsϕ
Gtotal
v
a
gbulk
− gm
ϕm
dx,dy,dz
dt(6)

input value
8.3143 J K− 1 mol− 1
400, 500, 550, 600 ◦ C
MRTt
10− 2, 10− 4, 10− 6, 10− 8,
10− 100.50 J m− 2
See Eq. 2.
6.3295 × 10− 5 m3 mol− 1
Converted from s (Eq. 6.).
See Eq. 3.
(22.8 (1138 – T)/25902)0.25
5 × 10− 9 m
0.01

3

R. Fukushima et al.

Physics of the Earth and Planetary Interiors 354 (2024) 107227

3. Results

Even after the mode of the histogram (i.e., the most probable value)
switches into around the equilibrium value, both the mean orderparameter and the fraction of ordered domains continue increasing. As
inferred from a previous theoretical study of disordering kinetics (Salje
and Wruck, 1988), this could raise a difficulty in detecting any differ­ence in ‘transitional’ omphacites (e.g., t* = 10.5–12.5 in Fig. 3) onlywith conventional X-ray experiments, whose results would reflect themost probable state.
For the sensitivity analysis of each parameter, we summarized thetime evolution of the spatially averaged order parameter (Fig. 4a), andthe relationship between the mean APD size (δ) and the spatial fractionof the ordered phase (F) (Fig. 4b,c). Because the mobility (M) acts only asa scaling factor of time, in this study, we focus only on how the annealed
temperature and the initial cation configuration influence the resultant
texture. The degree of order was found to suddenly increase after acertain incubation period, and this depends on the magnitude of theinitial fluctuation (Fig. 4a). Specifically, a larger initial variance (or alarger σ value) resulted in the incubation period being shorter at aconstant temperature. Note that the temperature effect on the incuba­tion duration is hard to evaluate since the mobility term should be asensitive function of temperature. As to the δ–F relationship, both ahigher temperature condition and a smaller initial variance lead to a

3.1. Advent and growth of ordered domains
The representative results of our modeling with random initialconfiguration are presented in Fig. 2a. It is clearly observed thatpartially ordered domains with various sizes randomly form in thealmost disordered matrix and subsequently coarsen. At a later stage, theamplitude of the long-range order increases throughout the entire
calculated area to form APBs with the equilibrium degree of order. Thiscan be explained by the decreasing total contribution of the gradientenergy term to suppressing ordering due to the slightly-ordered domain
growth. By converting the order-parameter field into a corresponding
gray-scale image, we reproduced images similar to the real dark-field
images with abundant equiaxed APDs (Fig. 2b). It should be notedthat only a ~50 ◦ C difference of the annealed temperature leads todistinct textures in the ordered domain size; a slightly higher tempera­ture can result in significantly larger, although fewer, domains.
The overall ordering process can be also figured out by generating ahistogram for each representative time step (Fig. 3). The initial sharp
distribution at |ϕ| = 0 is broaden and extends to higher |ϕ| values, andfinally produces a sharp peak again at the equilibrium degree of order.

Fig. 2. Representative results of the 2D simulations of equiaxed APD formation (scale bar: 1 μm). a, the time evolution (t* = 2, 4, 6, 8, 10, 12) of the phase field at atemperature of 550 ◦ C. The positive/negative values of the phase field were plotted in reddish/bluish colors. For each snapshot, the field is normalized by the
maximum degree of order. b, dark-field-image like visualization of selected results of the 2D simulations for equiaxed APD formation. Simulation results at 550 ◦ C (t*= 11, 12) and 600 ◦ C (t* = 16) are presented (brighter regions: ordered; darker regions: disordered). All the calculations were performed with σ = 10− 4.
4

R. Fukushima et al.

Physics of the Earth and Planetary Interiors 354 (2024) 107227

Fig. 3. Temporal evolution of the frequency distribution of the order parameter for a representative 2D simulation of equiaxed APD formation (550 ◦ C, σ = 10− 4).
The inset in each histogram (t* = 5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 15) shows the simulated dark field image (brighter regions: ordered; darker regions:disordered), and the red vertical line indicates the spatially averaged order parameter. The bin width of the histogram is set to 0.005. Note that the maximum value ofthe y-axis (fraction of grids in each class) is different: 1.0 or 100% (t* = 5, 6.5, 7.5); 0.2 or 20% (t* = 8.5, 9.5, 10.5, 11.5, 12.5); and 0.5 or 50% (t* = 15). (Forinterpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)larger δ/F ratio (Fig. 4b,c), which indicates the presence of larger but
fewer ordered domains. A similar tendency was observed in the 3D
simulation results (Supplementary Fig. S2), except that the incubation
duration seems to be a little longer than that in the 2D case (Fig. 4a).
In addition, our quantification of their spatial wavelengths, which
are naturally expected to be twice of the representative ordered cluster
sizes, demonstrated that its growth would obey the parabolic rate law,regardless of the input parameters and the dimension (Fig. 4d). This isconsistent with previous studies of the cluster size evolution during
cation ordering in alloys (Hashimoto et al., 1978; Kawasaki et al., 1978).
Note that we interpreted the wavelength drop shown in Fig. 4d (e.g.,
depicted by the yellow marker) as being caused by the abrupt orderingand subsequent topological changes in the phase field (cf., Marek and
Demjen, 2017). We excluded the wavelength–time curves after suchsudden decreases for the purpose of this research.

Because their heterogeneous nucleation processes cannot be tech­nically simulated only with a simple phase-field approach, the columnar
APD formation cannot be perfectly modeled unless the individual APD
nucleation sites are successfully estimated. Here, we assumed that theinitial presence of small completely ordered domains (randomly inphase or antiphase to each other) were aligned along the bottom sideof the simulation area. For the remaining area, we set random initialconfigurations, similar to those set in the previous case (See Supple­
mentary Fig. S3 for the detailed model setup). This treatment serves tosimplify the heterogeneous nucleation of the ordered domains as if they
appeared immediately after the omphacite had recrystallized or
precipitated. This approach also disregards the internal strain field andsubsequent domain growth anisotropy, which would have influenced
the resultant microtexture.
Our simulation revealed that the nuclei grow perpendicularly to themodeled nucleation boundary (or surface), making their shapes
columnar, and finally interact with equiaxed ordered domains (Fig. 5a).
On the beginning, the length of the columnar APDs monotonously in­creases. Thereafter, equiaxed ordered domains appear in the almost
disordered matrix and each of them reaches the equilibrium degree oforder, similar to the results observed in the previous modeling. After thepinning of the columnar APDs to the equiaxed ones, their elongation
speed suddenly decreased, and the elongated domains subsequently
shrunk due to domain coarsening. This suggests there should be a
maximum elongation length for such columnar APDs that nucleated at
grain boundaries or dislocations. Under a constant temperature, asmaller initial variance leads to larger maximum length, probably due tothe longer incubation time for the equiaxed APD formation (Fig. 5b).
Note that, depending on the amplitude of the initial fluctuation, ripening

3.2. Interaction between equiaxed APDs and initial nuclei
Although our modeling with random initial configurations seems toreplicate equiaxed APD textures in natural omphacites, this approachstill cannot simulate microtextures with columnar APDs that terminate
at grain boundaries and dislocations, as shown in Fig. 1c. There are two
mechanisms proposed to explain such columnar APD growths: a het­
erogeneous nucleation process of ordered phases at a high supercoolingdegree (Champness, 1973; Carpenter, 1979c), and a rapid overgrowth ofordered omphacite grain (Fukushima et al., 2021b). However, the latter
would be rare, and has been reported only in omphacites included ingarnet porphyroblasts. Therefore, in this study, we attempted to model
how columnar APDs grow from the initial ordered nuclei.
5R. Fukushima et al.

Physics of the Earth and Planetary Interiors 354 (2024) 107227

Fig. 4. Results of the sensitivity analysis of the input parameters (σ and T) for the equiaxed APD formation simulation. Under each condition, the simulation was
performed 5 times to check the effects of the initial randomness. a, the time evolution of the spatially averaged order parameter with different σ values (10− 6, 10− 4,
10− 2) at 550 ◦ C. Two representative results of the 3D simulations (σ = 10− 4, 10− 2) are also presented. b, the δ–F relationship at different temperatures (500, 550,
600 ◦ C) with σ = 10− 4. c, the δ–F relationship with different σ values (10− 6, 10− 4, 10− 2) at 550 ◦ C. d, the log–log plot of the spatial wavelength versus the normalized
time with different conditions of (T, σ) = {(550 ◦ C, 10− 2), (550 ◦ C, 10− 4), (600 ◦ C, 10− 4)}. A 3D representative result with T = 400 ◦ C and σ = 10− 6 is also presented.
See details in the text for interpretation of the trends.

between the adjacent elongated domains can proceed before theirimpingement on the equiaxed APDs. The maximum length and theelongation speed of the columnar APDs with the corresponding pa­
rameters seems similar even in the 3D simulation, except that thinnercolumnar APDs diminish faster than in the corresponding 2D simulation(Supplementary Fig. S4).

difference in the incubation time can be compensated for by slightly
adjusting the values of the initial fluctuation variance, the following
discussion focuses on the incipient ordering process based on the 2D
simulation results.
4.2. Timescale for the transitional state

4. Discussion

Our simulation result of the equiaxed APD formation demonstrated
that the duration of ordering from the initially disordered omphacite
depends not only on the annealed temperature but also the initial cation
configuration. The presence of an incubation period, probably due to thecompetition between the gradient and bulk energy terms, may allow theoverall ordering rate to be slower than a simple assumption of uniform
ordering would lead us to expect. Nevertheless, our simulation at aconstant temperature suggests the timescale required for abrupt
ordering is constant (e.g., ~5 normalized time units at 550 ◦ C) regard­
less of the initial cation configuration (Fig. 4a).
Based on the observation that incompletely ordered omphacite withequiaxed APDs (δ ~100 nm, F ~0.5) remains in low-T eclogite from
Syros, Greece (Fig. 1b), we can roughly estimate the duration required
for the transition. Because the omphacite grain was included at the rim

4.1. Dimensional effects
If we consider the larger measurement uncertainty of each parameterin the 3D simulation due to the smaller system size, we did not observe asignificant difference between the 2D and 3D simulations in terms of the
δ–F relationship, the growth exponent of ordered domains, and the
maximum columnar APD length (or the elongation speed). On the other
hand, domain coarsening between columnar APDs is faster, and the
incubation time is slightly longer in the 3D simulation. We postulate that
this is due to the larger gradient energy contribution for the individual
ordered clusters or APDs in the 3D system. However, since our main
focus is not the later-stage domain coarsening, and because this small
6

R. Fukushima et al.

Physics of the Earth and Planetary Interiors 354 (2024) 107227

Fig. 5. Representative results of the 2D simulations of columnar APD formation (scale bar: 1 μm). a, the time evolution (t* = 0, 4, 8, 12, 16, 20) of the phase field at atemperature of 550 ◦ C and with σ = 10− 4. The maximum length of the columnar APDs is ~170 nm. b, selected results of the 2D simulations of columnar APD
formation at 550 ◦ C. Simulation results with σ = 10− 6 (t* = 16), σ = 10− 8 (t* = 20.5), and σ = 10− 10 (t* = 25) are presented. The maximum APD lengths are ~230
nm, ~300 nm, and ~370 nm, respectively.

of a garnet porphyroblast (Fukushima et al., 2021b), it should have beenannealed at the peak metamorphic temperature of ~550 ◦ C (Laurent
et al., 2018) for a shorter time than the timescale of the whole eclogitefacies metamorphism. According to the previous geochronologicalstudies of eclogite-facies rocks from the Cycladic Blueschist Unit (Lagoset al., 2007; Dragovic et al., 2012, 2015; Gorce et al., 2021; Tual et al.,2022), the age difference between the garnet-core formation and its
post-peak isothermal decompression would be at most a few million
years. Thus, we can, for example, envisage that the omphacite grain
retained the incomplete order despite being annealed at 550 ◦ C for 1
million years.
We found that the normalized time of t* ~11.5 explains the micro­
texture with δ ~100 nm and F ~0.5 when σ = 10− 4 (tentative, ~0.1% of
the equilibrium value). Assuming t* = 11.5 is equivalent to 106 years,
the normalized time unit equals ~9 × 104 years, and the mobility (M) at
550 ◦ C is calculated to be ~2 × 10− 9 mol J− 1 yr− 1. This subsequently
means that the ordered-phase fraction rapidly increases within ~4 ×
105 years (i.e., 5 normalized time units) at 550 ◦ C. This order-ofmagnitude estimate is consistent with the APD-size variation in a sin­
gle natural low-T eclogite (Fukushima et al., 2021b); such a limited but
not too-short timescale compared to the duration of isothermal
decompression could allow for various APD textures depending on thetiming of the individual omphacite crystallizations, especially in thelater-stage eclogite-facies metamorphism. Thus, ordered domains inincompletely ordered omphacite could have the potential to record
short-span events during slab eclogitization. Nevertheless, for a precise

estimation of the mobility, we must choose the most feasible values forthe model parameters, especially for the unknown initial fluctuation
variance.
4.3. Assessment of the formation mechanism of columnar APDs
In our simulation of the columnar APD formation, their maximum
lengths must be overestimated since fully ordered nuclei were intro­
duced in the initial phase field. However, our simulation with various
degrees of the initial fluctuation failed to reproduce longer columnar
APDs than 1 μm (Fig. 5b), which were indeed observed in the naturalomphacite specimens (Fig. 1c). Considering the incubation time mainly
depends on the initial fluctuation variance and not much on any otherparameters including the APB energy (Supplementary Fig. S5), thiscontradiction should be answered by the following two explanations
about the initial configuration of the columnar-APD-rich omphacite: (1)the fluctuation was too small to simulate with our present model; or (2)the initial presence of heterogeneous, columnar domains that werealready ordered to some extent (i.e., a non-random fluctuation in theinitial disordered matrix).
If the first interpretation was correct, a question arises as to why thefluctuation is so small compared to equiaxed-APD-rich omphacites thatescaped any heterogeneous nucleation of APDs along their grainboundaries/dislocations. Unfortunately, it is uncertain whether the
initial fluctuation variance can be different by tens of order of magni­
tude between omphacites with equiaxed APDs and those with a plenty of7

R. Fukushima et al.

Physics of the Earth and Planetary Interiors 354 (2024) 107227

columnar APDs. On the other hand, the second explanation suggests thatdetailed microtextural analysis of columnar-APD-rich omphacites might
be a key to estimating the initial cation configuration itself. For example,
if a certain metastable disordered omphacite was formed by a dis­solution–precipitation mechanism, such a local heterogeneity was likely
to be introduced when the omphacite nucleated; and the crystal struc­ture of parent minerals can remain in the newly-forming phases,
possibly a silica-rich amorphous phase, to finally form the partially or­dered omphacite (Hellmann et al., 2012; Konrad-Schmolke et al., 2018).
Such a succession of the precursor structure is also inferred from theprevious report of the irregular space group of omphacite replacing from
igneous augite (Carpenter, 1978).
In any case, considering the microtextural variations (i.e., equiaxed
and long columnar APDs) in natural omphacites included in the samegarnet porphyroblast, we can conclude that there should have beensome difference in the initial cation configuration among the individual
metastable ‘disordered’ omphacites, even if they formed under similar
P–T conditions. It should depend on its nucleation mechanism (e.g.,
direct nucleation from solute-bearing aqueous fluids, or recrystallizationof lower-grade minerals such as chlorite and amphibole), and possibly
the compositional heterogeneity of the precursors. This suggests an
intimate relationship between the omphacite precipitation kinetics andits micron-scale ordered texture.

expressions, such as omphacite solid-solution models based on sym­
metric formalism (Holland and Powell, 1996; Green et al., 2007).
Nevertheless, these models basically assume an identical value for theorder parameters of M1 (QM1) and M2 (QM2) sites, oversimplifying theempirical correlation of QM1 = 2QM2 (Carpenter et al., 1990; BoffaBallaran et al., 1998). Therefore, we need a more accurate expression forthe omphacite free energy, especially at low-T eclogite-facies conditions.
Furthermore, when applying this model to natural omphacites, it isessential to consider the impact of compositional deviations on the bulk
free energy (and potentially the gradient energy) as they are not negli­
gible when discussing variations in ordered domain size andmorphology. It is widely acknowledged that jadeite content influences
the final APD size (Carpenter and Smith, 1981; Müller et al., 2011), andthe addition of iron is known to suppress cation ordering (Cámara et al.,1998). Such deviated compositions lead to a reduction of the effective
critical temperature, lowering the equilibrium degree of order, andmaking APBs thicker (Cahn and Hilliard, 1958). Thus, one should becautious about discriminating incompletely ordered textures fromequilibrium ones based only on the ordered domain fractions, without
considering these effects.
For omphacite whose composition is extremely different from theideal one, it is easy to judge that its microtexture is likely to be atequilibrium even though its ordered phase fraction is small. In fact, thearea 2 omphacite in Fig. 1c is characterized by an extremely high ironcontent (See the caption and Supplementary Table S1), probably thereby
tiny APDs are sparsely distributed there. However, for omphacite whose
composition is relatively less deviated, the compositional criteria toassure incomplete ordering are still ambiguous. Further microtexturalinvestigations combined with both precise compositional measurements
and theoretical considerations are required.

4.4. Constraints on the initial fluctuation amplitude
Our numerical simulation has revealed that cation configurations in
the metastable disordered omphacite strongly affect the overall ordering
rate and resultant APD textures. Given that disordered omphacite innatural low-T eclogites can nucleate via solid–solid reactions, and thatits mineral reactants should be heterogeneously distributed on a
micrometer scale, its metastable crystal structure should have been
determined not only by temperature but also by bulk-rock chemicalcompositions, local chemical heterogeneity, and the precipitation ki­
netics of the omphacite grain. Considering these complex factors, itwould be promising if one could determine the modeled initial fluctu­
ation variance by analyzing the observed microtexture itself in a naturalomphacite. Our simulation result of the equiaxed APD formation sug­gests that a two-order-of-magnitude difference in fluctuation standard
deviation can be distinguished based on the observed domain size andfraction (Fig. 4c). Additionally, if there are any columnar APDs that
nucleated at grain boundaries or dislocations in the early stage, weshould be able to constrain the fluctuation amplitude from theirmaximum length (Fig. 5). In this sense, it may be possible to establish theinitial fluctuation amplitude based on the natural microtextures,
although the reliable calibration of all the other non-unique parameters,such as the APB energy, grid size, and ordered-phase threshold in thedark-field TEM images, would be still challenging.

4.6. Perspectives
Our simulation has revealed the versatility of incompletely orderedomphacite as a recorder of relatively short-span geodynamic histories
and the precipitation kinetics of metastable disordered omphacite. Our
results emphasize that the incipient ordering and later-stage domaincoarsening should be treated as different regimes. This suggests that thepreviously proposed equation of δ8 ∝ t is not necessarily applicable tothe low-T eclogite-facies omphacites. Even if such sluggish domain
growth is due to impurity partitioning into the APBs (Carpenter, 1982b),
its effects should be insignificant during the incipient ordering before
the development of sharp APBs. Therefore, the mean size of orderedclusters in natural low-T eclogite-facies omphacites might follow thegrowth rate law δ2 ∝ t before reaching a completely ordered state(Fig. 4d). This time-sensitive domain growth would allow us to propose
the potential of the microtexture in relatively young (i.e., laternucleating) omphacites or incompletely ordered omphacites as a geo­
speedometer. On the other hand, if certain metastable disorderedomphacites are strongly characterized by the heterogeneous cation
configurations and if later-stage domain coarsening is sluggish enough,the observable various APD morphologies may reflect different forma­tion mechanisms of individual omphacite grains, rather than their burial
durations.
Despite the invaluable insights of our simulations, we haveencountered the limitation in the applications of such a simple forward
modeling. We require reasonable calibration of the several parameters,
including the initial cation configuration, the gradient energy, and themobility values. To constrain the latter two variables, we might be able
to perform experiments to observe the ordering process by using ‘highly
fluctuated’ disordered omphacites. It has been reported that reorderingof high-T-annealed omphacites occurred at 600 ◦ C, probably due to theinitial presence of short-range order (Holland, 1983). This suggests thata relatively large initial fluctuation enables such rapid ordering, as one
can observe in a laboratory. Furthermore, we must consider that thelimited number of APD observations, due to their time-consuming

4.5. Reliability of the excess free energy expression
In addition to the difficulty in calibrating the model parameters andestimating the initial cation configurations, we must carefully consider
the reliability of the excess free energy expression of omphacite. The
most questionable aspect of the Landau type expression is that theequilibrium degree of order at 500–600 ◦ C (~0.69–0.75) is significantly
lower than that of M1 sites (QM1) in some omphacites from natural low-T
eclogites (~0.9). With the model we used here, such a highly orderedstate can be attained at the equilibrium temperature of <200 ◦ C.
Although Carpenter et al. (1990) explained this as ordering occurring
during cooling after peak metamorphism, there is no guarantee thatordering can proceed at such a low temperature. Indeed, a static energy
calculation based on empirical interatomic potentials demonstrated thecalculated equilibrium value of QM1 at 600 ◦ C is ~0.9 (Vinograd et al.,2007), suggesting a potential low-T inaccuracy in the model used in our
simulation. There are, of course, alternative choices for free energy
8Physics of the Earth and Planetary Interiors 354 (2024) 107227

R. Fukushima et al.

procedures, may generally bias the characterization of a given eclogite.
Both experiments and comprehensive analyses of natural omphacites
are required for the further development of microtextural interpretationin omphacite.
If the ambiguity of the initial cation configuration and the imaging
bias problem are solved, and if the compositional effect can be correctly
handled, APD textures in omphacite would be a promising tool toconstrain the kinetics of eclogite-facies slab dehydration. Specifically, its
combined use with geochronology for eclogite-forming minerals (e.g.,garnet, zircon) should open a new way for discussing short-span meta­morphic events. Even if one has not calibrated the actual mobility value
itself, the phase field model-based geospeedometer should enable its
physically meaningful interpolation into a shorter timescale, based on alonger timescale constrained by the complement radiometric dating. In
addition, one might be able to incorporate some effects of deformationand local strain, with the aid of TEM observations of various dislocations(e.g., Van Roermund and Lardeaux, 1991; Müller et al., 2011), latticepreferred-orientation analyses for matrix-forming omphacites (e.g.,
Brenker et al., 2002; Urlich and Mainprice, 2005), and spectroscopic
studies to measure residual stress in omphacite included in garnet(Baratelli et al., 2024). The establishment of the new indicator for theslab-eclogitization kinetics should help us understand the grain-scale
reaction process, even from natural samples whose protoliths weretotally replaced by the eclogitic mineral assemblages. Our investigation
here is therefore the first step to realize this goal, by providing acornerstone for disentangling its ordering transformation kinetics.

BGI at the University of Bayreuth, and the Joint Usage/Research Center
PRIUS, Ehime University, Japan (assisted by Hiroaki Ohfuji, Sayako
Inoué, and Takeshi Sakai), and was funded in part by grants from the
JSPS KAKENHI JP22KJ0304 to R.F. and JP21H01174 to T.T. The au­thors are also grateful to the International Joint Graduate Program in
Earth and Environmental Sciences (GP-EES), Tohoku University. The
Scios FIB and the Titan G2 TEM at BGI were financed by DFG Grant Nos.
INST 91/315-1 FUGG and INST 91/251-1 FUGG, respectively. We ex­press our deep appreciation to Mitsuhiro Toriumi and Tadao Nishiyama
for their valuable comments on the simulation. We thank two anony­mous reviewers for the constructive discussion, and Jie Li for her
thoughtful editorial handling.
Appendix A. FIB-TEM/STEM-EDS analysis of natural omphacite
specimens
We have selected natural omphacites from low-T eclogites in the
SMM (Guatemala: Tsujimori et al., 2006), Syros (Greece: Fukushima
et al., 2021b), and Omi (Japan) as representatives, as shown in Fig. 1.
The omphacites from the SMM (one grain) and Syros (two grains) arefound as inclusions in garnet porphyroblasts (Fukushima et al., 2021a),while the Omi omphacite (one grain) is a matrix-forming one. We pre­pared thin foils with a thickness of ~100–150 nm from each sample
using the focused ion beam (FIB) technique. Subsequently, we captured
dark-field TEM images of these foils using a transmission electron mi­
croscope. For the FIB fabrication process, we utilized the FEI Scios
DualBeam FIB-SEM system, which is equipped with a gallium ion-gun
and a thermal field emission electron gun, at Bayerisches Geoinstitut,
Universität Bayreuth (BGI), and Geodynamics Research Center, Ehime
University (GRC). The dark-field TEM images of the Syros omphacites
were obtained using an FEI Titan G2 80–200 S/TEM at BGI, while thoseof the other omphacites were obtained using a JEOL JEM-2010 TEM at
GRC.
Representative major element compositions of the individual sam­ples were obtained by performing energy-dispersive X-ray spectroscopy(EDS) with 4 silicon drift detectors attached to the FEI Titan G2 80–200
S/TEM. The EDS spectra data were processed by using the software
Esprit, with the Cliff-Lorimer method for the Z-number correction (Cliffand Lorimer, 1975) and the absorption correction based on experi­mentally calibrated k-factors of major elements including oxygen (vanCappellen and Doukhan, 1994). The ferric iron contents were estimatedwith Fe3+ = Na – Al, by speculating that Al is only located in the M1
sites. The site fractions of the other divalent cations were determined sothat (Ca + Mg + Fe2+)M2 = (Mg + Fe2+)M1 is maintained, assuming
equipartition of Mg/Fe2+ between the M1 and M2 sites. We thencalculated the following values to characterize each omphacite
composition: Xaeg = Fe3+, Xaug = (Mg + Fe2+)M1, Mg# = [Mg/(Mg +Fe2+)]M1, and FeT = (Fe2+ + Fe3+).

5. Conclusions
By revisiting microtextural formation processes due to the orderingtransformation in omphacite under low-T eclogite-facies conditions, wedemonstrate their potential for decoding the slab eclogitization kinetics.
With the macroscopic phase-field approach, we can illustrate micro­
textural evolution due to cation ordering in a metastable disorderedomphacite. Although it remains challenging to constrain the parametervalues, our simulation suggests the size and fraction of ordered domains
can significantly change within <~105 years. In addition, our simula­tion results are consistent with the previous observation of various
ordered-domain morphologies in natural omphacites under similar P–T
conditions; there is no surprise if a given initial cation configuration
determines the final microtexture, reflecting the nucleation–growth ki­
netics of the omphacite grain. This study proposes the new framework
for analyzing omphacite microtextures, and would enlarge the appli­
cability of their formation dynamics for petrogeological research of
convergent orogens.
CRediT authorship contribution statement
Ryo Fukushima: Writing – review & editing, Writing – originaldraft, Visualization, Software, Project administration, Methodology,
Investigation, Funding acquisition, Data curation, Conceptualization.
Tatsuki Tsujimori: Writing – review & editing, Supervision, Resources,
Funding acquisition, Conceptualization. Nobuyoshi Miyajima: Writing– review & editing, Investigation.

Appendix B. Supplementary data
Supplementary data to this article can be found online at https://doi.
org/10.1016/j.pepi.2024.107227.
References

Declaration of competing interest

Allen, S.M., Cahn, J.W., 1979. A microscopic theory for antiphase boundary motion andits application to antiphase domain coarsening. Acta Metall. 27 (6), 1085–1095.
Baratelli, L., Murri, M., Alvaro, M., Prencipe, M., Mihailova, B., Cámara, F., 2024. Raman
scattering of omphacite at high pressure: towards its possible application to elastic
geothermobarometry. Am. Mineral. https://doi.org/10.2138/am-2023-9140.
Boffa-Ballaran, T., Carpenter, M.A., Domeneghetti, M.C., Tazzoli, V., 1998. Structural
mechanisms of solid solution and cation ordering in augite-jadeite pyroxenes: I. A
macroscopic perspective. Am. Mineral. 83 (5–6), 419–433.
Brenker, F.E., Prior, D.J., Müller, W.F., 2002. Cation ordering in omphacite and effect on
deformation mechanism and lattice preferred orientation (LPO). J. Struct. Geol. 24(12), 1991–2005.

The authors declare no competing interests.
Data availability
Data will be made available on request.
Acknowledgments
This research received support from CNEAS at Tohoku University,
9Physics of the Earth and Planetary Interiors 354 (2024) 107227

R. Fukushima et al.
Brenker, F.E., Müller, W.F., Brey, G.P., 2003. Variation of antiphase domain size inomphacite: a tool to determine the temperature–time history of eclogites revisited.
Am. Mineral. 88 (8–9), 1300–1311.
Broadwell, K.S., Locatelli, M., Verlaguet, A., Agard, P., Caddick, M.J., 2019. Transientand periodic brittle deformation of eclogites during intermediate-depth subduction.
Earth Planet. Sci. Lett. 521, 91–102.
Buscombe, D., Rubin, D.M., Warrick, J.A., 2010. A universal approximation of grain size
from images of noncohesive sediment. Case Rep. Med. 115 (F2), F02015.
Cahn, J.W., Hilliard, J.E., 1958. Free energy of a nonuniform system. I. Interfacial freeenergy. J. Chem. Phys. 28 (2), 258.
Cámara, F., Nieto, F., Oberti, R., 1998. Effects of Fe2+ and Fe3+ contents on cation
ordering in omphacite. Eur. J. Mineral. 10 (5), 889–906.
Carpenter, M.A., 1978. Kinetic control of ordering and exsolution in omphacite. Contrib.
Mineral. Petrol. 67 (1), 17–24.
Carpenter, M.A., 1979a. Omphacites from Greece, Turkey, and Guatemala: composition
limits of cation ordering. Am. Mineral. 64 (1–2), 102–108.
Carpenter, M.A., 1979b. Experimental coarsening of antiphase domains in a silicate
mineral. Science 206 (4419), 681–683.
Carpenter, M.A., 1979c. Contrasting properties and behaviour of antiphase domains inpyroxenes. Phys. Chem. Miner. 5 (2), 119–131.
Carpenter, M.A., 1982a. Time-temperature-transformation (TTT) analysis of cation
disordering in omphacite. Contrib. Mineral. Petrol. 78 (4), 433–440.
Carpenter, M.A., 1982b. Omphacite microstructures as time-temperature indicators of
blueschist- and eclogite-facies metamorphism. Contrib. Mineral. Petrol. 78 (4),441–451.
Carpenter, M.A., 1994. Evolution and properties of antiphase boundaries in silicate
minerals. Ph. Transit. 48 (1–3), 189–199.
Carpenter, M.A., Okay, A., 1978. Topotactic replacement of augite by omphacite in ablueschist rock from north-west Turkey. Mineral. Mag. 42 (324), 435–438.
Carpenter, M.A., Putnis, A., 1985. Cation order and disorder during crystal growth: some
implications for natural mineral assemblages. In: Thompson, A.B., Rubie, D.C. (Eds.),
Metamorphic Reactions: Kinetics, Textures, and Deformation. Springer, New York,New York, pp. 1–26.
Carpenter, M.A., Salje, E., 1989. Time-dependent Landau theory for order/disorderprocesses in minerals. Mineral. Mag. 53 (372), 483–504.
Carpenter, M.A., Smith, D.C., 1981. Solid solution and cation ordering limits in hightemperature sodic pyroxenes from the Nybö eclogite pod, Norway. Mineral. Mag. 44(333), 37–44.
Carpenter, M.A., Domeneghetti, M.C., Tazzoli, V., 1990. Application of Landau theory tocation ordering in omphacite I: equilibrium behaviour. Eur. J. Mineral. 2 (1), 7–18.
Champness, P.E., 1973. Speculation on an order-disorder transformation in omphacite.
Am. Mineral. 58 (5–6), 540–542.
Cliff, G., Lorimer, G.W., 1975. The quantitative analysis of thin specimens. J. Microsc.
103 (2), 203–207.
Dragovic, B., Samanta, L.M., Baxter, E.F., Selverstone, J., 2012. Using garnet to constrain
the duration and rate of water-releasing metamorphic reactions during subduction:
an example from Sifnos, Greece. Chem. Geol. 314, 9–22.
Dragovic, B., Baxter, E.F., Caddick, M.J., 2015. Pulsed dehydration and garnet growthduring subduction revealed by zoned garnet geochronology and thermodynamic
modeling, Sifnos, Greece. Earth Planet. Sci. Lett. 413, 111–122.
Fukushima, R., 2023. Om-COS: cation ordering simulator for omphacite. Zenodo.
https://doi.org/10.5281/zenodo.7982097.
Fukushima, R., Tsujimori, T., Aoki, S., Aoki, K., 2021a. Trace-element zoning patterns inporphyroblastic garnets in low-T eclogites: Parameter optimization of the diffusionlimited REE-uptake model. Isl. Arc. 30, e12394.
Fukushima, R., Tsujimori, T., Miyajima, N., 2021b. Various antiphase domains in garnethosted omphacite in low-temperature eclogite: a FIB-TEM study on heterogeneousordering processes. Am. Mineral. 106 (10), 1596–1605.
Gorbatov, O.I., Lomaev, I.L., Gornostyrev, Y.N., Ruban, A.V., Furrer, D., Venkatesh, V.,
Novikov, D.L., Burlatsky, S.F., 2016. Effect of composition on antiphase boundary
energy in Ni3Al based alloys: ab initio calculations. Phys. Rev. B 93 (22), 224106.
Gorce, J.S., Caddick, M.J., Baxter, E.F., Dragovic, B., Schumacher, J.C., Bodnar, R.J.,
Kendall, J.F., 2021. Insight into the early exhumation of the Cycladic Blueschist unit,
Syros, Greece: combined application of zoned garnet geochronology,
thermodynamic modeling, and quartz elastic barometry. Geochem. Geophys.
Geosyst. 22 (8), e2021GC009716.
Green, E., Holland, T., Powell, R., 2007. An order-disorder model for omphacitic
pyroxenes in the system jadeite-diopside-hedenbergite-acmite, with applications to
eclogitic rocks. Am. Mineral. 92 (7), 1181–1189.
Hashimoto, T., Nishimura, K., Takeuchi, Y., 1978. Dynamics on transitional ordering
process in Cu3Au alloy from disordered state to ordered state. J. Phys. Soc. Jpn. 45(4), 1127–1135.
Hellmann, R., Wirth, R., Daval, D., Barnes, J.P., Penisson, J.M., Tisserand, D., Epicier, T.,
Florin, B., Hervig, R.L., 2012. Unifying natural and laboratory chemical weatheringwith interfacial dissolution–reprecipitation: a study based on the nanometer-scalechemistry of fluid–silicate interfaces. Chem. Geol. 294, 203–216.
Holland, T., 1983. The experimental determination of activities in disordered and shortrange ordered jadeitic pyroxenes. Contrib. Mineral. Petrol. 82 (2–3), 214–220.
Holland, T., Powell, R., 1996. Thermodynamics of order-disorder in minerals: II.
Symmetric formalism applied to solid solutions. Am. Mineral. 81 (11− 12),1425–1437.

Holland, T., Powell, R., 2011. An improved and extended internally consistentthermodynamic dataset for phases of petrological interest, involving a new equationof state for solids. J. Metamorph. Geol. 29 (3), 333–383.
Kawasaki, K., Yalabik, M.C., Gunton, J.D., 1978. Growth of fluctuations in quenched
time-dependent Ginzburg-Landau model systems. Phys. Rev. A 17 (1), 455–470.
Konrad-Schmolke, M., Halama, R., Wirth, R., Thomen, A., Klitscher, N., Morales, L.,
Schreiber, A., Wilke, F.D.H., 2018. Mineral dissolution and reprecipitation mediatedby an amorphous phase. Nat. Commun. 9 (1), 1637.
Lagos, M., Scherer, E.E., Tomaschek, F., Münker, C., Keiter, M., Berndt, J., Ballhaus, C.,2007. High precision Lu–Hf geochronology of Eocene eclogite-facies rocks from
Syros, Cyclades, Greece. Chem. Geol. 243 (1–2), 16–35.
Lardeaux, J.M., Caron, J.M., Nisio, P., Péquignot, G., Boudeulle, M., 1986.
Microstructural criteria for reliable thermometry in low-temperature eclogites.
Lithos 19 (3–4), 187–203.
Laurent, V., Lanari, P., Naïr, I., Augier, R., Lahfid, A., Jolivet, L., 2018. Exhumation of
eclogite and blueschist (Cyclades, Greece): pressure–temperature evolution
determined by thermobarometry and garnet equilibrium modelling. J. Metamorph.
Geol. 36 (6), 769–798.
Marek, J., Demjen, E., 2017. Is autocorrelation image analysis the proper method in
nanoparticle sizing? J. Nanopart. Res. 19, 208.
Müller, W.F., Xu, Z., Brenker, F.E., 2011. Transmission electron microscopy of omphacite
and other minerals in eclogites from the CCSD borehole, China: indications for theirdeformation and temperature history. Eur. J. Mineral. 23, 645–659.
Ninomiya, H., Eguchi, T., Kanemoto, H., 1990. Dynamics of pattern formation of
antiphase ordered domain structures in alloys. Ph. Transit. 28 (1–4), 125–131.
Phakey, P.P., Ghose, S., 1973. Direct observation of anti-phase domain structure inomphacite. Contrib. Mineral. Petrol. 39 (3), 239–245.
Pollington, A.D., Baxter, E.F., 2010. High resolution Sm–Nd garnet geochronology
reveals the uneven pace of tectonometamorphic processes. Earth Planet. Sci. Lett.
293 (1–2), 63–71.
Pollok, K., Heidelbach, F., John, T., Langenhorst, F., 2014. Spherulitic omphacite inpseudotachylytes: microstructures related to fast crystal growth from seismic melt at
eclogite-facies conditions. Chem. Erde 74 (3), 407–418.
Ross, C.R., 1988. Statistical mechanical modeling of the kinetics of order-disorder inomphacitic pyroxenes. Phys. Chem. Miner. 15, 274–282.
Salje, E., Wruck, B., 1988. Kinetic rate laws as derived from order parameter theory II:interpretation of experimental data by Laplace-transformation, the relaxation
spectrum, and kinetic gradient coupling between two order parameters. Phys. Chem.
Miner. 16, 140–147.
Tsujimori, T., Mattinson, C., 2021. Eclogites in different tectonic settings. In: Elias, S.,
Alderton, D. (Eds.), Encyclopedia of Geology, 2nd ed. Academic Press, London,
pp. 561–568.
Tsujimori, T., Sisson, V.B., Liou, J.G., Harlow, G.E., Sorensen, S.S., 2006. Petrologic
characterization of Guatemalan lawsonite eclogite: Eclogitization of subductedoceanic crust in a cold subduction zone. Geol. Soc. Am. Spec. 403, 147–168.
Tual, L., Smit, M.A., Cutts, J., Kooijman, E., Kielman-Schmitt, M., Majka, J., Foulds, I.,
2022. Rapid, paced metamorphism of blueschists (Syros, Greece) from laser-based
zoned Lu-Hf garnet chronology and LA-ICPMS trace element mapping. Chem. Geol.
607, 121003.
Urlich, S., Mainprice, D., 2005. Does cation ordering in omphacite influence
development of lattice-preferred orientation? J. Struct. Geol. 27 (3), 419–431.
Vaithyanathan, V., Chen, L.Q., 2000. Coarsening kinetics of δ′-Al3Li precipitates: phasefield simulation in 2D and 3D. Scr. Mater. 42 (10), 967–973.
van Cappellen, E., Doukhan, J.C., 1994. Quantitative transmission X-ray microanalysis ofionic compounds. Ultramicroscopy 53 (4), 343–349.
van Roermund, H.L.M., Lardeaux, J.M., 1991. Modification of antiphase domain sizes inomphacite by dislocation glide and creep mechanisms and its petrologicalconsequences. Mineral. Mag. 55 (380), 397–407.
Viete, D.R., Hacker, B.R., Allen, M.B., Seward, G.G.E., Tobin, M.J., Kelley, C.S.,
Cinque, G., Duckworth, A.R., 2018. Metamorphic records of multiple seismic cycles
during subduction. Sci. Adv. 4 (3), eaaq0234.
Vinograd, V.L., Gale, J.D., Winkler, B., 2007. Thermodynamics of mixing indiopside–jadeite, CaMgSi2O6–NaAlSi2O6, solid solution from static lattice energy
calculations. Phys. Chem. Miner. 34, 713–725.
Wu, X., Han, Y., Meng, D., Yang, W., Li, D., 2000. TEM and HRTEM study of
microstructures in omphacite from UHP eclogites at Shima, Dabie mountains, China.
Acta Geol. Sin. 74 (2), 154–162.
Xie, Z.J., Liu, X.W., Jin, Z.M., Li, Z.Y., 2020. Microstructures and phase transition inomphacite: constraints on the P-T path of Shuanghe eclogite (Dabie orogen). J. Earth
Sci. 31, 254–261.
Yokoyama, T., Nakamura, E., Kobayashi, K., Kuritani, T., 2002. Timing and trigger of arc
volcanism controlled by fluid flushing from subjecting slab. Proc. Jpn. Acad. Ser. B.
78 (7), 190–195.
Yu, H.F., Jones, I.P., Smallman, R.E., 1994. The effects of temperature, composition andstrain rate on the deformation microstructure of Ni3Al. Philos. Mag. A 70 (6),951–967.
Zhu, J.Z., Wang, T., Ardell, A.J., Zhou, S.H., Liu, Z.K., Chen, L.Q., 2004. Threedimensional phase-field simulations of coarsening kinetics of γ’ particles in binary
Ni–Al alloys. Acta Mater. 52 (9), 2837–2845.

10

Supplemental Material for:
Simulation of microtextural evolution in omphacite: Orderingtransformation kinetics as unexplored archives of slab eclogitization
Ryo Fukushima1, Tatsuki Tsujimori1,2, and Nobuyoshi Miyajima3
1

Graduate School of Science, Tohoku University, Sendai 980-8578, Japan

2

Center for Northeast Asian Studies, Tohoku University, Sendai 980-8576, Japan

3

Bayerisches Geoinstitut, Universität Bayreuth, 95440 Bayreuth, Germany

This file includes:Figs. S1–S5 and Table S1

1

Ymax = 4.05 + Y1

Calculated intensity

4
Y = 5(X – 0.1)2 + Y1
(X > 0.1)

321

Saturation level = 0.4Ymax ~ 1.67

Y = Y1 – 100Y1(X – 0.1)2(0 ≤ X ≤ 0.1)
Y1 = 0.03Ymax ~ 0.125

0
0.0

0.2

0.4

0.6

0.8

1.0

Absolute value of the order parameter
Fig. S1. The tentative response curve for the dark-field-like image conversion (black solid
curve). Note that the saturation level (red dotted line) was set to 40% of the maximum intensity(Ymax).

2

b

a

1.0
Ordered phase fraction (F)Ordered phase fraction (F)

1.0
0.8
0.6
0.4

σ = 10–4
Calculated in 3D
T = 550℃
T = 500℃

0.2

0.8
0.6
0.4

T = 550℃
Calculated in 3D
σ = 10–2
σ = 10–4

0.2
0.0

0.0
0

50
100
150
200
Mean APD size (), nm

250

0

50
100
150
200
Mean APD size (), nm

250

Fig. S2. Results of the sensitivity analysis of the input parameters (σ and T) for theequiaxed APD formation simulation in 3D systems. Each simulation was performed 5 times tocheck the effects of the initial randomness. a, the 𝛿–F relationship at different temperatures (500
and 550°C) with σ = 10–4. b, the 𝛿–F relationship with different σ values (10–4 and 10–2) at
550°C. Note that the 3D system size was smaller than that of 2D, which would have led to thelarge variations in the 𝛿–F values.

3

a

y

zero-Neumann boundary condition
5 nm × 400 grids (= 2 μm)

zero-mean Gaussian noise(with the variance of σ2)

periodic boundary condition
5 nm × 400 grids (= 2 μm)

periodic boundary conditionx

initially ordered nuclei (25 × 40 nm2)aligned with a spacing of 10 nm
zero-Neumann boundary conditionyczero-Neumann boundary condition
5 nm × 120 grids (= 0.6 μm)

periodic boundary conditionyperiodic boundary condition

zero-mean Gaussian noise(with the variance of σ2)

periodic boundary condition
5 nm × 120 grids (= 0.6 μm)

periodic boundary conditionxz

periodic boundary conditionb

initially ordered nuclei (25 × 40 × 40 nm3)aligned with a spacing of 10 nm
zero-Neumann boundary condition

periodic boundary condition

Fig. S3. Schematic figures showing the boundary conditions and initial configurations in
the columnar APD simulations. The white rectangles/squares depict initially completelyordered nuclei, whose positive/negative signs were randomly designated. a, the condition in the2D simulation. b, c, the condition in the 3D simulation on the z = 5 plane and x = 120 plane,
respectively. The red line indicates the z = 5 plane.

4–0.73t* = 4t* = 6t* = 8t* = 10t* = 12t* = 14

Order parameter+0.73

Fig. S4. A result of the 3D simulation of columnar APD formation (scale bar: 500 nm). The
time evolution (t* = 4, 6, 8, 10, 12, 14) of the phase field at a temperature of 550°C and with σ =10–4 is presented. The maximum length of the columnar APDs is ~170 nm, similar to the 2D
simulation.

5

a

b

0.7

t* = 11

Mean order parameter

0.6
0.5
0.4
T = 550℃
σ = 10–4
Calculated in 2D
s = 0.25 J m–2
s = 0.50 J m–2

0.3
0.2
0.1
0.0
0

5

10
15
Normalized time

20

25

Fig. S5. Representative simulation results with a smaller APB energy value (s = 0.25 J m–2,
half of the default) in the 2D system (T = 550°C; σ = 10–4). a, the time evolution of thespatially averaged order parameter (blue) in the equiaxed APD formation simulation. For
comparison, the result with the default APB energy is also presented (black). Each simulation
was performed 5 times to check the effects of the initial randomness. b, a snapshot of thecolumnar APD formation at a time of t* = 11 (scale bar: 1 µm). The maximum APD length is~100 nm, just ~60% of that in the default case, suggesting the difficulty in making longcolumnar APDs (e.g., 1 µm) with this order of magnitude (i.e., 10–1 J m–2) of the APB energy
value.

6

Table S1. Representative compositions of omphacites in Fig. 1 (O = 6).
Fig. 1aFig. 1bFig. 1c-1
Fig. 1c-2
Fig. 1c-3(SMM)(Syros)(Syros)(Syros)(Syros)Fig. 1d(Omi)

Si [T]

2.00

2.00

1.94

1.95

1.95

2.02

Al [M1]

0.44

0.54

0.47

0.33

0.50

0.42

Fe2+[M1]

0.14

0.16

0.15

0.18

0.16

0.12

Fe3+[M1]

0.08

0.00

0.09

0.24

0.04

0.02

Mg [M1]

0.34

0.31

0.35

0.30

0.36

0.42

Fe2+[M2]

0.01

0.03

0.01

0.02

0.01

0.02

Mg [M2]

0.02

0.05

0.03

0.03

0.03

0.06

Ca [M2]

0.45

0.39

0.47

0.43

0.48

0.46

Na [M2]

0.51

0.52

0.56

0.57

0.54

0.44

Total

4.00

3.99

4.06

4.05

4.05

3.98

7