Fujita et al.
Progress in Earth and Planetary Science(2024) 11:40
https://doi.org/10.1186/s40645-024-00632-z

RESEARCH ARTICLE

Progress in Earth andPlanetary Science

Open Access

Imaging flow focusing and isolationof aqueous fluids in synthetic quartzite:
implications for permeability and retained fluid
fraction in deep‑seated rocks
Wakana Fujita1,5, Michihiko Nakamura1* , Kentaro Uesugi2, Philipp Eichheimer3,4, Marcel Thielmann3 andGregor J. Golabek3

Abstract
The microstructure of realistic fluid–rock systems evolves to minimize the overall interfacial energy, enabling local
variations in fluid geometry beyond ideal models. Consequently, the permeability–porosity relationship and fluid
distribution in these systems may deviate from theoretical expectations. Here, we aimed to better understandthe permeability development and fluid retention in deep-seated rocks at low fluid fractions by using a combined
approach of high-resolution synchrotron radiation X-ray computed microtomography imaging of synthesized rocks
and numerical permeability computation. We first synthesized quartzite using a piston-cylinder apparatus at different fluid fractions and wetting properties (wetting and non-wetting systems with dihedral angles of 52° and 61°–71°,
respectively) under conditions of efficient grain growth. Although all fluids should be connected along grain edges
and tubules in the homogeneous isotropic wetting fluid–rock system enabling segregation by gravitational compaction in natural settings, the fluid connectivity rapidly decreased to ~ 0 when the total fluid fraction decreased
to 0.030–0.037, as the non-ideality of quartzite, including the interfacial energy anisotropy (i.e., grain faceting),
became critical. In non-wetting systems, where the minimum energy fluid fraction based solely on the dihedral angle
is ~ 0.015–0.035, the isolated (disconnected) fractions was 0.048–0.062. A streamline computation in the non-wetting
system revealed that with decreasing total porosity, flow focusing into fewer channels maintained permeability,
allowing the effective segregation of the connected fluids. These results provide insight into how non-wetting fluids
segregate from rocks and exemplify the fraction of retained fluids in non-wetting systems. Thus, the findings suggest
a potential way for wetting system fluids to be transported into the deep Earth’s interior, and the amount of fluids
dragged down to the Earth’s interior could be higher than what was previously estimated.
Keywords Permeability, Quartzite, Digital rock physics, Dihedral angle, Microstructure development*Correspondence:
Michihiko Nakamura
michihiko.nakamura.e8@tohoku.ac.jp
Full list of author information is available at the end of the article© The Author(s) 2024. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to theoriginal author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images orother third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of thislicence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/.

Fujita et al. Progress in Earth and Planetary Science(2024) 11:40

1 Introduction
Segregation of geological fluids is an essential elementary
process that controls the chemical differentiation andphysical properties of the Earth’s interior. Melt segregation from a partially molten region controls the chemistry of magmas and the solid matrix in a two-phase flow
(e.g., McKenzie and Bickle 1988; Iwamori 1993; Faul
2001). The excess pore fluid pressure produced by thedehydration of a subducting slab is a key factor controlling seismic activity within the slab (Kita et al. 2006;
Peacock et al. 2011; Morishige and van Keken 2018). In
turn, the pore fluids left in the slab contribute to the total
amount of water downdragged into the Earth’s deep interior. Melt segregation is also important for the core formation and concomitant chemical evolution of the Earthand other planetary bodies (Stevenson 1990; Terasaki
et al. 2008; Bagdassarov et al. 2009; Watson and Roberts
2011; Solferino et al. 2020).
In deep-seated rocks below the brittle–ductile transition, where rocks do not sustain fractures, the rockmicrostructure develops to reduce the interfacial energy(Watson and Brenan 1987) through a process of “texturalequilibration”. The local textural equilibrium has been
generally premised by the establishment of the dihedral
angle (θ), the most elemental parameter that characterizes pore fluid geometry. In an ideal simplified fluid–rocksystem, with a uniform grain size and isotropic interfacial energy, θ is established at every solid–liquid interface with a constant mean curvature of the fluid–solidinterfaces (Beere 1975; Bulau et al. 1979; Waff 1980; Parkand Yoon 1985; Ghanbarzadeh et al. 2015; Rudge 2018).
When θ < 60°, all the interstitial fluids are connected, andgravity drives fluid segregation and compaction of ductile solid matrix on geological timescales (McKenzie andBickle 1988). When θ > 60°, to the contrary, fluids withless than the volume fraction at a percolation threshold
can be isolated and retained in the grain boundaries,
thereby escaping gravitational segregation. The lattercase has been explicitly studied to examine the possibility of water being downdragged into the Earth’s deep
interior along with the subducting slab (e.g., Watson andLupulescu 1993; Ono et al. 2002; Mibe et al. 2003; Liuet al. 2018), because the dihedral angles between aqueous
fluids and eclogite-constituent minerals in the slab dehydration condition are large. Here, unaddressed issues aremechanisms for the non-wetting, hard-to-connect fluids
segregation, and the resultant proportion of retained fluids in the rocks.
Furthermore, the geometries of melts and aqueous
fluids in the grain and interphase boundaries of realpolycrystalline rocks often deviate from ideal values.
Interfacial energy anisotropy results in the statistical diversity in the dihedral angles (Wark and Watson

Page 2 of 161998) and faceting of flat solid–fluid interfaces, causing
disruption of fluid networks and isolation of pore fluids (Waff and Faul 1992; Cmíral et al. 1998; Wark et al.
2003; Yoshino et al. 2006; Shangshang and Faul 2016;
Price et al. 2006; Huang et al. 2023). In addition, theformation of fluid pools surrounded by more than threeneighboring grains has been reported in olivine–basalt
melt, quartz–albite melt, olivine-FeS melt, and quartz–
H2O systems (Jurewicz and Watson 1984; Waff and Faul
1992; Cmíral et al. 1998; Shangshang and Faul 2016). In
these experiments, large fluid pools were formed adjacent to small and rounded grains, suggesting the merging of neighboring pore fluids via the consumption ofsmall grains. Similar fluid pool formation was also
observed in situ for the analogous fluid–rock system
of ethanol and norcamphor (Walte et al. 2003). Theseexperimental results indicate that the microstructureof realistic fluid–rock systems develops to minimize
the total interfacial energy of the system, allowing
deviation in the local fluid geometry from ideal models
with a single dihedral angle and a uniform liquid–solidinterfacial energy. Thus, how and to what extent geological fluids are retained in deep-seated rocks is not
fully constrained.
The permeabilities of geological bodies have beenpreviously estimated using hydrological models basedon geothermal and metamorphic data (e.g., Manning and Ingebritsen 1999; Ingebritsen and Manning2010), laboratory measurements (e.g., Wark and Watson 1998; Katayama et al. 2012; Ganzhorn et al. 2019),
and numerical simulations (e.g., von Bargen and Waff
1986; Zhu and Hirth 2003; Miller et al. 2014; Eberhard
et al. 2022). Several theoretical models have been proposed for the permeability–porosity ( k –ϕ ) relationship by assuming a single dihedral angle and uniform
liquid–solid interfacial energy in a system (hereafter,
the “reference model”) to establish a mechanically andthermodynamically stable topology of fluids in a texturally equilibrated monomineralic rock (Beere 1975;
Bulau et al. 1979; Waff 1980; Park and Yoon 1985;
Ghanbarzadeh et al. 2015; Rudge 2018). Here, porosity is often replaced by the volume fraction of fluidor melt depending on the context. This assumption
requires uniform space-filling tesselation of grains,such as a tetrakaidecahedral unit cell. The permeability
is primarily controlled by the minimum cross-sectional
area of grain-edge tubules (aperture size) Atb, which isa function of the grain diameter d , ϕ , and θ (Watson
1982; von Bargen and Waff 1986; Rudge 2018). An
increase in the fluid fraction increases Atb, and thus,
permeability. At a given ϕ , Atb is larger for smaller dihedral angles (von Bargen and Waff 1986). Therefore, theabsolute value of permeability ( kabs ) depends on theFujita et al. Progress in Earth and Planetary Science(2024) 11:40

Page 3 of 16

grain size of the rocks, even if the rocks are geometrically similar. The relationship between kabs and ϕ in textually equilibrated rocks can be expressed as:

kabs =

d2ϕn
,
C

(1)

where C is a geometrical factor and n is an exponent,
typically 2–3. The power-law exponent n depends on
porosity: n = 2 when porosity is low and fluid resides at
triple-junction tubules along grain edges, and n = 3 forhigher porosity conditions with fluid-wetting grain faces(Rudge 2018). A theoretical model that considers variations in the dihedral angle and grain-scale heterogeneity of the fluid distribution also requires n = 3 (Zhu andHirth 2003). Permeability decreases with increasing dihedral angle at the same porosity because a larger dihedral
angle decreased Atb by distributing more fluids at thegrain corners than grain-edge tubules and grain faces,
which are more effective for fluid interconnection. Notably, the dependence of permeability on the dihedral angle
is prominent at low porosities, indicating that permeability is more sensitive to fluid topology at low porosities
than at high porosities.
Laboratory measurements of permeability using synthesized, texturally equilibrated quartzite and calcite
have shown that permeability is reduced in moderately
faceted systems by segregating fluids into grain-scalepools, resulting in Atb decrease (Yoshino et al. 2006;
Price et al. 2006). For example, in a moderately faceted
quartz–water (brine) system, the values of n and C were
3 and 270, respectively (Wark and Watson 1998; Lianget al. 2001). Thus, the k –ϕ relationship showed that thepermeability decreased from those of the reference
models by a factor ranging 5–10 at the same porosity. These previous studies indirectly interpreted themeasured permeability by observing 2D cross sections.
Direct computation of permeability using the 3D tomographic images, a digital rock physics approach, would
bridge the gap and contribute to a better understanding of how pore microstructures control the permeability and isolated fluid fraction of deep-seated rocks. In
the present study, we used quartzite, a polycrystalline
quartz aggregate. This well-studied system works as asemi-analog of silicate rocks because quartz’s mineralogical and microstructural properties at elevated pressures and temperatures are well established. Althoughquartzite remains simpler than actual rocks, it exhibits
nearly “textually equilibrated” (interfacial energy-minimized) microstructures through normal grain growthin the piston-cylinder syntheses (e.g., Watson andBrenan 1987; Wark and Watson 1998; Nakamura andWatson 2001), thus allowing to investigate the microstructural evolution of natural rocks.

2 Methods
2.1 Quartzite synthesis

We utilized some of the quartzite samples, which weresynthesized to examine the capsule-scale chemical redistribution under high pressure and temperature in thecompanion paper (Fujita et al. 2022). Other details of
the piston-cylinder synthesis were reported by Fujita
et al. (2022). The starting material was a powder of finegrained (< 45 μm) quartz and amorphous silica mechanically mixed with a weight ratio of 3:1. The amorphoussilica, prepared from the sol–gel method, was adequate
for rapid achievement of steady grain growth (Nakamura and Watson 2001; Ohuchi and Nakamura 2007a,b). We added the fluid source to produce three different ­CO2 fractions (XCO2 in mol): XCO2 = 0 (pure water),
0.27–0.29, and 0.44–0.47. Brucite powder was used as thewater source for the pure water experiments. Althoughenstatite will be formed at the bottom of the capsule after
dehydrated brucite reacts with silica, it would not affectthe fluid distribution in quartzite (Fujita et al. 2022) orother microstructures of interest in this study. For the­CO2-bearing runs, oxalic acid dihydrate powder was used
as a fluid source. Two geometries were tested for mixing with ­SiO2 powder (i.e., initial mixing by grinding ina mortar or placing at the capsule bottom). However, nosystematic difference effect of the fluid distribution was
found, as reported by Fujita et al. (2022). The added volumetric fluid fraction ϕad was controlled between 0.019
and 0.101 at each experimental condition by changing
the amount of fluid sources. The fluid fraction calculation included the absorbed H­ 2O in the amorphous silica (~ 0.6 wt%). Hereafter, we refer to the C­ O2-free and­CO2-bearing systems as wetting and non-wetting systems, respectively. All the synthetic experiments wereconducted at 900 °C and 1 GPa. From the series of experiments with run durations of 24–382 h reported by thecompanion paper (Fujita et al. 2022), we selected and prepared 11 representative experimental products with run
durations of 192 h, added fluid fractions of 0.025–0.101
and XCO2 = 0, 0.28, and 0.39–0.44, for non-destructive
imaging of the interstitial pores with synchrotron radiation X-ray computed microtomography (SRμCT—e.g.,
Lindquist et al. 2000; Zhu et al. 2011; Miller et al. 2014;
Skemer et al. 2017; Thomson et al. 2018) (see Table 1).
2.2 Synchrotron X‑ray microtomography

The three-dimensional (3D) pore structure was obtainedby using SRμCT. Details of sample preparation for imaging are presented in Fujita et al. (2022). The SRμCT
was conducted at the BL20XU experimental hutch2 ofSPring-8 (Uesugi et al. 2012) with the 510 nm edge lengthof the cubic voxels in the obtained images. The details

0.101

0.094

0.094

0.032

0.032

0.057

0.057

0.032

0.032

0.032

0.032

0.054

0.054

0.054

0.054

0.031

0.031

0.028

0.028

0.028

0.028

0.028

0.028

0.025

0.025

0.025

QW13-13(#1)

QW13-11(#1)

QW13-11(#2)

QW14-11(#1)

QW14-11(#2)

QW14-21(#1)

QW14-21(#2)

QW15-11(#1)

QW15-11(#2)

QW15-12(#1)

QW15-12(#2)

QW15-13(#2)

QW15-13(#3)

QW15-14(#1)

QW15-14(#2)

QW15-21(#1)

QW15-21(#2)

QW15-22(#1)

QW15-22(#2)

QW15-22(#4)

QW15-22(#5)

QW15-22(#6)

QW15-22(#7)

QW15-23(#1)

QW15-23(#3)

QW15-23(#5)

0.014

0.037

0.048

0.048

0.062

0.014

0.013

0.031

0.074

0.062

0.056

0.036

0.032

0.059

0.053

0.042

0.044

0.056

0.055

0.064

0.069

0.053

0.048

0.079

0.091

0.042

Measuredporosity,
φmes

0

0.034

0.047

0.038

0.057

0

0

0

0.069

0.057

0.036

0

0

0.050

0.015

0

0

0.008

0.025

0.061

0.067

0.048

0.043

0.077

0.089

0.040

Connected
porosity,
φcon

0

0

0

0

0

0

0

0

0

0

0

0.39

0.39

0.39

0.39

0.28

0.28

0.28

0.28

0

0

0.28

0.28

0.44

0.44

0

Fluid
composition
XCO2 [mol/
mol]

–13

7.50 × ­10–13–12

1.01 × ­10–12

1.51 × ­10–12–13

5.26 × ­10–13–12

1.18 × ­10–12

1.02 × ­10–12–12

1.29 × ­10–13

5.38 × ­10–13

Permeability was not calculated

Proportion of the range of kx, ky, and kz to the mean permeability

Average number of channels in a slice of each subvolume (with φcon> 0) obtained from every 10 computational grid cells (Fig. 8)

§

&

n.d.

2.86 × ­10–13

4.70 × ­10–13

–

–

n.d.

n.d.

n.d.

1.85 × ­10–12

1.12 × ­10–12

9.24 × ­10–14

n.d.

n.d.

3.18 × ­10–13

7.73 × ­10–13

2.07 × ­10–12

1.07 × ­10–12

2.35 × ­10–13

3.44 × ­10–13

1.35 × ­10–13

9.15 × ­10–14

2.49 × ­10–13

1.08 × ­10–12

2.20 × ­10–12

5.66 × ­10–13

4.30 × ­10–13

2.32 × ­10–12

3.34 × ­10–12

2.12 × ­10–13

k ­[m ]

2

Permeability (mean)

129

109

64

9

115

123

16

38

58

25

21

25

35

[%]

§

Permeability
anisotropy

–

4.7

5.4

n.d.

n.d.

–

–

–

8.9

n.c

4.1

–

–

12.0

5.5

–

–

4.3

n.c.**

28.3

25.5

25.3

28.1

35

38.3

5.4

Mean numberof ­channels&,Nch(2024) 11:40**Not computed

The direction without connection from one side to the opposite side of the subvolume

n.d.

5.38 × ­10–13

1.31 × ­10–12
n.d.

–

–

–¶
–

n.d.

n.d.

n.d.

n.d.

n.d.

n.d.

1.52 × ­10

1.08 × ­10–12

3.62 × ­10–13
2.85 × ­10

2.50 × ­10–13
–12

n.d.

n.d.

n.d.

n.d.

3.41 × ­10–13

3.46 × ­10–13

n.d.

5.06 × ­10–14

2.17 × ­10–13

1.37 × ­10–13

n.d.

n.d.

9.15 × ­10–14

2.69 × ­10–13

1.06 × ­10–12

2.05 × ­10–12

4.21 × ­10–13

4.70 × ­10–13

2.56 × ­10–12

n.d.

n.d.

n.d.

n.d.

1.86 × ­10

4.58 × ­10–12

3.64 × ­10n.d.

kz ­[m ]

2

n.d.

n.d.

n.d.#

2.29 × ­10

2.70 × ­10

3.61 × ­10–13

2.08 × ­10–12–12

3.55 × ­10

1.75 × ­10–13
2.82 × ­10

2.49 × ­10–13–12ky ­[m ]2

kx ­[m ]2

Permeability (orthogonal directions)

¶

71.0

71.0

71.0

56.2

56.2

56.2

56.2

56.2

56.2

66.3

66.3

53.2

53.2

52.2

52.2

57.9

57.9

53.2

53.2

38.9

38.9

40.5

40.5

38.3

38.3

69.2

Grain
diameter[µm]#*The run information of the other samples is presented in Table 1 of Fujita et al. (2022)

Added fluid
fraction, φad

Subvolume no.*Table 1 Summary of experimental conditions and analytical results for subvolumes (stacks of SRμCT slices)

Fujita et al. Progress in Earth and Planetary Science
Page 4 of 16

Fujita et al. Progress in Earth and Planetary Science(2024) 11:40

of the beamline setting were described by Huang et al.
(2021).
2.3 Microstructure characterization

Cross sections of quenched samples that differed fromthose used for SRμCT were observed using backscatteredelectron (BSE) imaging using a field-emission type scanning electron microscope (FE–SEM; JSM-7100F, JEOL
Ltd., Akishima, Tokyo, Japan) equipped with an energy
dispersive X-ray spectrometer (EDS, JEOL JED-2300F).
The details of the sample preparation procedures arereported by Fujita et al. (2022). For the 3D characterization of the pore microstructure and numerical simulation of fluid flow, cubic subvolumes with a dimension
512 × 512 × 512 voxels (i.e., 261 × 261 × 261 µm3) were
cut out from the reconstructed images. For run QW1313, whose grain size was coarser than that of the others, a larger subvolume of 600 × 600 × 600 voxels (i.e.,
306 × 306 × 306 µm3) was cut out. The edge lengths ofeach subvolume were 4.4–6.7 times larger than theirmean grain size, which are too small to represent thecapsule-scale fluid distribution (Fujita et al. 2022) but canrepresent multi-grain-scale fluid distributions. To check
representativeness, we examined two to six subvolumes
from different positions in each sample. Only one subvolume was cut from QW13-13 because the CT-scannedarea was much smaller than those of the other samples.
A smoothing filter was applied using the ImageJ software(National Institutes of Health, USA; https://​imagej.​net/),and the subvolumes were converted into a binary format. Dilation cracks perpendicular to the pressurization
direction, formed upon quenching and pressure releasein piston-cylinder experiments, were found only in
QW15-13. For QW15-13, subvolumes were chosen fromareas without dilation cracks. The open grain boundaries
were broadly recognized especially along relatively large
grains, and hence, they were eliminated by dilating andsubsequently eroding the solid phase by 1-pixel layerusing the software package “Slice” (Nakano et al. 2006).
We measured the connected porosity (ϕcon ) in addition
to the total porosity (ϕmes) from the corrected 3D images
using “Slice”. Here, the ϕcon is a volume fraction of thepores connected from one side to the opposite side of the
cubic subvolume and thus should directly affect permeability. The relationship between ϕmes and ϕcon highlighted
the microstructural difference between the wetting andnon-wetting systems.
2.4 Numerical computation of flow velocity
and permeability

The velocity and flow patterns of a fluid can be describedby the Reynolds number Re, which is the ratio of inertial
forces to viscous forces and is defined as:Page 5 of 16

Re =

ρvdch
,
η

(2)

where ρ is the density, v is the velocity, η is the viscosity of
the fluid, and dch is the characteristic internal diameter ofthe fluid channel. For all simulations, we employed a fluid
viscosity of 1 Pa s as the computed permeability does notdepend on the viscosity value. Reynolds numbers typically range from ­10–9 to ­10–10 for magmas (Glazner 2014)
and from ­10–8 to ­10–5 for groundwater flow. Therefore,
owing to the narrow fluid channels of our samples, Re is
expected to be small (< 1); thus, the fluid flow should belaminar (Bear 1988). Under these conditions, the incompressible Navier–Stokes equations can be simplified to
the Stokes equations (ignoring gravity).

∂vi
= 0;
∂xi
∂vj
∂vi
∂
η
+
∂xj
∂xj
∂xi

(3)

−

∂P
= 0,
∂xi

(4)

where xi and xj are the spatial coordinates corresponding
to the Cartesian axes.
Permeability was computed using a finite-difference
Lithospheric and Mantle Evolution Model (LaMEM)(Kaus et al. 2016; Eichheimer et al. 2019). LaMEM
employs a staggered finite-difference discretization thatdefines the pressure P at the cell centers and the velocities v at the cell faces. The computational grid cell corresponded to an image voxel. A free-slip boundary
condition was employed on the sides of the subvolume.
The Eqs. (3) and (4) were non-dimensionalized, discretized, and then solved for a given pressure difference
using the multigrid solvers in the PETSc library (Balayet al. 2019) and the V-cycle geometric multiplicative
multigrid solver (Fedorenko 1964; Wesseling 1995). The
applied non-dimensional fluid pressure was 1 on a subvolume face at the inlet side perpendicular to the flow
direction, and 0 on a subvolume face at the outlet side.
A single-grid solver was used in the case of QW13-13,
which had a larger subvolume. The relative and absoluteconvergence tolerances were specified to be ­10−9 and­10−10, respectively, sufficiently low for calculation convergence (Eichheimer et al. 2019). From the computed
velocity field in the flow direction z, the volume-averaged
velocity component vm was calculated as follows (e.g.,
Osorno et al. 2015):
1|vz |dv,vm =(5)
Vf V fwhere V f denotes the volume of the fluid phase. Permeability values were then calculated using Darcy’s law:

Fujita et al. Progress in Earth and Planetary Sciencekabs =

ηvm L,
�P(2024) 11:40

Page 6 of 16(6)

where L is the height of the subvolume and P is thepressure difference. A more detailed explanation of thecomputational method is provided by Eichheimer et al.
(2019).
Because the absolute permeability is length-scaledependent, it is normalized by grain size to extract thegeometric features of the system. In this study, the computed permeabilities were normalized to a mean grain
diameter of 1 mm so the results can be compared withthose with those of Wark and Watson (1998) and Lianget al. (2001):k = koutput ×

10002
d2mean,(7)where k is the normalized permeability in ­m2, koutput is thecomputed absolute permeability for the subvolume withmean grain diameter of dmean in µm. The error bars ofϕmes were calculated as half the difference between coating and scraping the solid (quartz) surface by a one-pixel
layer. To investigate the anisotropy of the permeability, apressure difference was applied to the three orthogonal
directions (kx , ky , and kz) for all subvolumes. Here, z was
the vertical direction (the pressurized direction of thepiston-cylinder experiments), and x and y were the horizontal directions perpendicular to each other. The meanvalues were then obtained to represent the sample.
The flow velocity distributions of the four representative subvolumes (QW15-22#1, QW13-13#1, QW1311#1, and QW15-13#2) were investigated to understandthe mechanism which governs permeability. The streamlines of the four subvolumes were visualized using theopen-source visualization application, ParaView (https://​www.​parav​iew.​org/).

3 Results
3.1 Shapes and occurrences of pore fluids

The mean grain size, grain size distribution, and dihedral
angles measured in the polished sections of the synthesized quartzite are described in Fujita et al. (2022). The
grain size distribution indicated that near-steady-stategrain growth was achieved in both the ­CO2-free (wetting) and C­ O2-bearing (non-wetting) systems. The dihedral angles between the fluids and quartz (θ), given by
the median values excluding facetted pores, were 52°,61°, and 71° for XCO2 = 0, 0.28, and 0.44, respectively.
The fluid fractions that minimize the total interfacial energy of the fluid-bearing polycrystalline system(MEFF) are ~ 0.035 and 0.015 for the dihedral angles of
61°, and 71°, respectively (Park and Yoon 1985; Watson2015). Based on the BSE images of the polished sectionsof the run products (Fig. 1a, b) and the two-dimensional
SRμCT slices (Fig. 1c, d), we can classify the pore fluids inthe quartz aggregate into three types based on the number of surrounding grains: (1) fluid lenses between twograin boundaries, (2) grain-edge tubules surrounded bythree grains, and (3) fluid pools surrounded by more thanfour grains. The fluid lenses were occasionally distributed
as lenticular islands (Fig. 1a–d) and sometimes exhibited necklace microstructures (German 1985). Amongrun products with similar fractions of fluid remaining inquartzite, fluid lenses occurred more commonly in samples with high XCO2 than in the other two types of samples. Grain-edge tubules were commonly observed inboth wetting and non-wetting systems. However, dry triple junctions existed even in the wetting system, in which
the dihedral angle was < 60° (Fig. 1a, b), as reported byprevious studies (Yoshino et al. 2006; Price et al. 2006).
The fluid pools were often surrounded by several to multiple grains (Fig. 1c, d), forming irregular shapes. The
size of the fluid pools is generally larger than those of thefluid lenses and cross section of the grain-edge tubules,
sometimes comparable to or larger than the mean grain
diameter, and occasionally as large as 50–100 µm. The
plucking-free SRμCT slices demonstrated that the grainscale pools were often adjacent to small and rounded
quartz grains (Fig. 1c, d). The grain-scale pools appeared
to form within quartzites with relatively large added fluid
fraction ϕad . In all runs, faceted grain-edge tubules werecommonly found together with an equilibrium shape surrounded by three grains with a constant interface curvature, as assumed in the reference models. Faceting wasmore common in the ­CO2-bearing quartz aggregates.
3.2 Pore connectivity

The connected porosity (ϕcon ) roughly coincided withthe measured total porosity ϕmes at larger porosities,
while ϕcon became lower than ϕmes at lower porosities (ϕmes  0.07) and reached 0 (Fig. 2), being consistent with the prediction from the presence of pinch-off
threshold for the measured permeability and observations of 2D cross sections (Wark and Watson 1998). In
the wetting system, the pore interconnection was sharply
cut off at ϕmes between 0.031 and 0.037. At ϕmes ≥ 0.037,
most subvolumes except QW15-21#1 and QW15-22#7
had ϕcon /ϕmes > 0.9. In contrast, in the non-wetting system, ϕcon extensively varied from 0.008 to 0.050 against asmall range of ϕmes (= 0.048–0.059), which corresponded
to a rapid increase in ϕcon /ϕmes from 0.13 to 0.91. In the
subvolumes with ϕmes above this range (ϕmes ≥ 0.079),
ϕcon /ϕmes was above 0.98, whereas at ϕmes ≤ 0.044, ϕcon
was 0. Notably, QW15-11#2 had the lowest ϕcon /ϕmes(0.14), but still maintained pore interconnection.

Fujita et al. Progress in Earth and Planetary Science(2024) 11:40

Page 7 of 16

Fig. 1 Representative BSE, synchrotron X-ray microtomography, and CT images. Representative BSE images of polished cross sections (a, b; Table 1of Fujita et al 2022) and slice images of synchrotron X-ray microtomography images (c, d) of the run products. a QW5-13: XCO2 = 0, ϕad = 0.102, 382 h;
b QW9-12: XCO2 = 0.44, ϕad = 0.094, 192 h; c QW13-13: XCO2 = 0, ϕad = 0.101, 192 h; and d QW13-11: XCO2 = 0.44, ϕad = 0.094, 192 h. White arrows with “L”,
fluid lenses; “T”, triple-junction tubules; “P”, fluid pools; white dotted arrows, dry triple junctions. Polished sections inevitably suffered from grainplucking (black arrows). Many fluid lenses were observed as dimples on the quartz grain surfaces as in (a). The voxel size of CT images was 0.513µm3; fluid lens and triple-junction tubules that were smaller than 0.51 µm are not resolved in (c, d). Note that the CT images are free from grainplucking. The fluid pools were often adjacent to small and rounded quartz grains as indicated by the white arrows marked “R.” Scale bars represent
40 µm for all images

In the present study, the pinch-off thresholds were0.037 and 0.048 in the wetting and non-wetting systems,
respectively, which were significantly higher than the values in the reference model (i.e., 0 for θ < 60° and ~ 0.03 forθ = 70°; Rudge 2018). The strict threshold value depends
on the resolution of 3D images. However, the fine voxel
size of SRμCT for the present study (510 nm) ensuresthat the permeability significantly deceases at this fluid
fraction range.
The pore interconnection behavior, along with thegeometry of each pore fluid, were visualized throughvolume renderings. The largest interconnected pores ofthe subvolumes are depicted in the upper row of Fig. 3,
and the remaining pores are shown in the lower row. The
largest connected pore was one that connected from one
side of the subvolume to the opposite side (i.e., that penetrated the subvolume all the way through).

As observed in the 2D images (Fig. 1), the pores in the
non-wetting systems (Fig. 3d–f, j–l) were more rounded
and the fluid tubules along the grain edges were thicker
than those in the wetting system as expected from thelarger dihedral angle (Fig. 3a–c, g–i). In both systems,
fluid networks were established by interconnecting
the fluids at the grain corners and relatively large fluid
pools that often had more than four faces, whereas
small fluids were left unconnected. In the wetting system, the fluid networks became thinner and sharperwith decreasing total fluid fraction; however, the connectivity remained high as long as the interconnection
was maintained (Fig. 3a, b, g, h). In contrast, in the
non-wetting system, a portion of the pores [0.45 = 0.28/(0.28 + 0.34)] in Fig. 3e maintained the pore connection
from one side to the opposite side. The isolated fluid
fractions after the connectivity reached 0 were larger inFujita et al. Progress in Earth and Planetary Science(2024) 11:40

Fig. 2 Relationship between the measured total porosity (ϕmes)
and connected porosity (ϕcon). Error bars for both ϕmes and ϕcon were
calculated by shrinking and expanding the solid phase by half a pixel
layer. Blue and orange symbols represent subvolumes of the wetting(XCO2 = 0) and non-wetting (XCO2 = 0.28–0.44) systems, respectively.
Filled symbols represent subvolumes having pore interconnection
between the two facing sides. Open symbols show subvolumes
that have isolated pores without a flow pathway

the non-wetting system (Fig. 3c, i) than in the wetting
system (Fig. 3f, l).
3.3 Permeability

The anisotropy of computed permeability, i.e., the difference between the maximum and minimum directions,
was lowest (9% of the mean) for QW15-21#2 and highest (128%) for QW15-23#3 (Table 1). Here, the samples

Page 8 of 16

for which permeability was obtained in only one or two
directions were excluded. No systematic tendency was
found in the vertical and horizontal directions of the
subvolumes, indicating that the fluid distribution was
not affected by the pressurization direction duringhigh-pressure syntheses. The anisotropy showed a weaknegative correlation with ϕcon and was relatively high at
ϕcon < 0.040. Ten out of eleven samples had ϕcon > 0.040
and showed an anisotropy < 65%. The small difference inlog scale among kx , ky , and kz implied that the subvolume sizes (261 × 261 × 261 µm and 306 × 306 × 306 µm)
were reasonable and could represent the grain-scale fluid
distribution.
Figure 4 shows the computed permeability k against
the measured total porosity (ϕmes) and connected porosity (ϕcon ) for all subvolumes with pore connection acrossthe subvolumes. The computed permeabilities generally
increased with increasing ϕmes and ϕcon in both the wetting and non-wetting systems, except for subvolumes
with relatively low ϕmes in the non-wetting system. The
computed permeabilities were smaller by approximately
half an order of magnitude than those of the reference
ϕ − k model calculated for θ = 50°, 60°, and 70° (Rudge
2018). However, most subvolumes showed an excellent agreement with the ϕmes–k relationship obtainedby direct laboratory measurements of the quartz + ­H2O
(i.e., wetting) system (Wark and Watson 1998), except
for five subvolumes (QW15-13#3, QW15-11#1, QW1521#1, QW15-13#2, and QW15-11#2) whose permeability agreed with or were lower than the laboratory
measurements for the calcite + ­H2O (i.e., non-wetting)system (Wark and Watson 1998). These subvolumes

Fig. 3 Three-dimensional fluid distribution of representative subvolumes. The upper row (a–f) and the lower row (g–l) show the maximum
volume pore and other pores of each subvolume, respectively. Three subvolumes of the wetting system—QW14-21#1 (a, g), QW13-13#1 (b, h),and QW15-23#5 (c, i)—and non-wetting system—QW13-11#1 (d, j), QW15-11#1 (e, k), and QW15-14#2 (f, l)—are shownFujita et al. Progress in Earth and Planetary Science(2024) 11:40

Page 9 of 16

Fig. 4 Computed permeabilities. Computed permeabilities in log units as a function of a the measured porosity ϕmes, and b connectedporosity ϕcon. The permeabilities are mean values of kx, ky, and kz obtained from the flow simulations in three orthogonal directions (x, y, and z)
and normalized to a mean grain diameter of 1 mm. Error bars for both ϕmes and ϕcon were calculated by shrinking and expanding the solid
phase by half a pixel layer. Blue and orange symbols represent subvolumes of the wetting (XCO2 = 0) and non-wetting (XCO2 = 0.28–0.44) systems,
respectively. Gray dotted lines are from the reference model by Rudge (2018) for dihedral angles of 50°, 60°, and 70° from left to right. Gray dashedlines are the results of the direct permeability measurements of synthetic quartzite (left, wetting) and calcite (right, non-wetting), respectively (Warkand Watson 1998; Liang et al. 2001)

exhibited low pore connectivity except for QW15-13#3(Fig. 2). Thus, the variation in k in the non-wetting system stemmed from the variation in pore connectivity.
Notably, QW15-13#2 and QW14-11#2 had the sameϕmes (0.053) but different k , ranging from 1.35 × ­10–13 to5.66 × ­10–13 ­m2, owing to the differences in ϕcon of 0.015
and 0.048, respectively. QW15-13#3 lacked connected
pores in x directions, suggesting the microstructural heterogeneity was large for its relatively high ϕcon (0.050).
Among the subvolumes of the wetting system, the mostnoticeable gap in k at similar ϕmes was between QW1521#1 and QW15-23#1 (2.35 × ­10–13 ­m2 at ϕmes = 0.056 and7.73 × ­10−13m2 at ϕmes = 0.048, respectively).
In Fig. 4b, the variation in k at relatively low ϕmes
(Fig. 4a) disappears against ϕcon . The results mostly fall
onto the ϕcon − k relationship of quartz analog (Warkand Watson 1998) or are slightly higher at ϕcon > 0.03.
The two subvolumes with the lowest ϕcon (QW1513#2, ϕcon = 0.015 and QW15-11#2, ϕcon = 0.008) were
off to the higher side of the ϕcon − k relationship of thequartz + ­H2O system.

coincide with ϕcon in the wetting system, the derived nvalues were almost the same for ϕmes – k and ϕcon–k , at
3.0 and 3.1, respectively. This similarity is consistent withthe previous ϕ – k relation of the quartz (wetting) andcalcite (non-wetting) analogs (Wark and Watson 1998).
For the synthetic calcite, ϕ was replaced with ϕ − 0.015,
where 0.015 is the pinch-off threshold below which no
fluid interconnection is allowed. One subvolume of the
non-wetting system (QW15-11#2) with significantly low
connectivity (ϕcon /ϕmes < 0.15) was excluded from the fitting. The obtained n value was 3.4 for ϕmes–k . For ϕcon–kin the non-wetting system, the data did not fit well into
the single power-law relation of Eq. (1). Instead, it could
be fit by using two lines switching at around ϕcon = 0.04;
n was determined to be 0.76 and 3.0 for the low ϕcon(= 0.0075–0.028) and high ϕcon(= 0.042–0.09), respectively. The obtained power-law exponent n for the low
ϕcon range was extremely low, indicating that permeability was insensitive to ϕcon , but rather controlled by othergeometrical parameters such as aperture size.

3.4 Power‑law exponent n

The results of flow velocity visualization showed thatthe velocity distribution was not homogeneous andwas concentrated in some channels (Fig. 6). On average, these high-velocity channels were thicker thanthe other channels, suggesting that the fluid selectively flowed through wider pore spaces. Bottlenecked

The power-law exponent n in this study was derived forboth the ϕmes–k (Fig. 5a) and ϕcon–k (Fig. 5b) relationships. For the wetting system, n was calculated by applying the least-squares method to all subvolumes where
the permeability was obtained. As ϕmes values mostly

3.5 Fluid flow visualization

Fujita et al. Progress in Earth and Planetary Science(2024) 11:40

Page 10 of 16

Fig. 5 Plot for the power-law exponent n of Eq. (1) for the a ϕmes–k and b ϕcon–k relationships.The non-wetting system in (b) was fitted by using
two lines with n = 0.76 and 3.0 for the low ϕcon(= 0.0075–0.028) and high ϕcon(= 0.042–0.090), respectively

segments have a high velocity along the channels owing
to the flow concentration. The volumetric flow rate q
in a channel of radius R is given by integration of thevelocity of the fluid v , expressed as,R

q = ∫ 2π rvdr.
0(8)

Therefore, fluid flow occurs effectively in a channel witha large R, and v becomes highest at a point with a small r
in the channel. A remarkable difference in the comparedflow pattern was found in QW15-13#2 (Fig. 6d), where
the flow lines consisted of significantly fewer channels,
which is consistent with its low pore connectivity (Fig. 2,
ϕcon = 0.015).
The volumetric flow rate along the flow direction was
examined for the selected subvolumes where the permeability was computed. First, computational grid cells
were sliced into layers perpendicular to the applied pressure gradient. The computed local velocity of the fluid
was stored at the cell faces (see Fig. 1 of Eichheimer et al.
2019). In the solid phase, the velocity was 0. Second, allvelocities stored at the cell faces located in the ith crosssection were summed to calculate the volumetric flowrate in the ith cross section of the channel. The summedup volumetric flow rate in all cross section is equal to thetotal flow rate Q:qi = Q.
(9)iFig. 6 Streamlines of the fluid velocity using open-source application
ParaView for four representative subvolumes. a QW13-13#1, b
QW14-21#2, c QW14-11#1, and d QW15-13#2. The high and low
velocity channels are shown in red and blue, respectively. Highvelocity commonly occurs in thick channels, suggesting that fluids
selectively flow in wider pore spaces. Subvolume QW15-13#2 in (d)shows a remarkable difference in flow patterns compared with thoseof the other three subvolumes, as it consists of significantly fewer
channels

To ignore computational artifacts in the flow calculation, cross-sectional channels that accounted for < 1% of
Q were excluded. The values of qi /Q were obtained fromevery 10th slice of grid cells perpendicular to the pressure gradient (i.e., the distance between slices is 5.1 µm).
In Fig. 7, the values of qi /Q are shown along the flow
direction for the four subvolumes, whose streamlines
were also visualized. The values of qi were calculated forthree orthogonal directions (x, y, and z), except for theFujita et al. Progress in Earth and Planetary Science(2024) 11:40

Page 11 of 16

Fig. 7 Ratio of volumetric flow rate to the total flow rate. Ratio of volumetric flow rate in the i-th cross section of a channel qi to the total flow rate Q
along the flow direction showing the flow velocity distribution. a QW15-22#1, b QW13-13#1, c QW13-11#1, and d QW15-13#2

z-direction of QW13-13#1, in which the connected pores
were absent. In QW15-22#1, QW13-13#1, and QW1311#1, each channel occupied < 30% of the total flow rate(i.e., qi /Q < 0.3) from the inlet to the outlet. The valuesof qi /Q along the flow direction were the lowest (~ 0.1)
and most homogeneous in QW13-11#1. In contrast,qi /Q fluctuated significantly along the flow direction in
QW15-13#2, indicating flow focusing. The maximum
qi /Q reached 1 near the outlet of the flow in the x-, y-,and z-directions, showing that the flows were concentrated in one channel at these distances.
3.6 Number of channels

The permeability computation (Figs. 4 and 5), fluid flow
visualization (Fig. 6), and flow rate investigation (Fig. 7)
demonstrated that the heterogeneity and anisotropy
resulting from non-uniform grain size, interfacial energy
anisotropy of minerals affected fluid connectivity andpermeability especially when the fluid fraction is low.
To examine these effects in the present experiments, thenumber of channels in a slice of a subvolume Nch was
obtained from every 10 computational grid cells along
the flow directions for the subvolumes with ϕcon > 0. To
ignore the computational artifacts in the flow calculation,
only cross-sectional channels that accounted for > 1% of
the total volumetric flow rate Q were counted as Nch . The
mean counted Nch of the each subvolume was normalized to the mean grain diameter of all the subvolumes(55 μm). Assuming a packing of hexagonal grains with aside length of 27.5 μm, the number of triple junctions ina slice with a side of 261 μm is approximately 70. Figure 8
shows that the normalized Nch increased with increasing
ϕcon both in the wetting and non-wetting systems (from 4
to 28 and from 4 to 38, respectively).
3.7 Flow concentration

We evaluated the degree of flow concentration in asingle channel by calculating qimax /Q , where qimax is the
maximum volumetric flow rate qi in a slice of computational grid cells. In a similar manner to Nch counting,Fujita et al. Progress in Earth and Planetary Science(2024) 11:40

Page 12 of 16

Fig. 8 Number of cross-sectional channels (Nch) in each subvolume against connected porosity (ϕcon). a Wetting and b non-wetting systems. Barsand solid circles represent the range and mean of each subvolume, respectively

qimax was obtained from every 10th slice of the computational grid cells along flow directions with connected
pores from inlet to outlet sides of the subvolume( ϕcon > 0). The mean qimax for each subvolume is plotted in Fig. 9. In the wetting system, qimax /Q showed noor weak negative correlation with ϕcon , while qimax /Qdecreased clearly from 0.70 to 0.17 with increasing ϕcon in the non-wetting system. These results are
consistent with the tendency of increasing Nch withincreasing ϕcon , especially for the non-wetting system.
The two subvolumes of the non-wetting system withthe lowest ϕcon and highest qimax /Q , QW15-11#2 andQW15-13#2, consisted of fewer channels than thosein the other subvolumes, as shown for QW15-13#2 inFigs. 6d and 7d. These subvolumes typically represent
concentrated fluid flows rather than pervasive flows.
The other subvolumes of the non-wetting systemshowed a gradual shift to pervasive flow with increasing ϕcon.

4 Discussion
4.1 Reliability of permeability computation

The reliability of the permeability values computed using
the 3D thermomechanical code LaMEM has been examined in previous studies. The permeability of pipe flow
computed using LaMEM has been confirmed to converge with the analytically calculated permeabilities ofthe Hagen–Poiseuille flow with increasing computational
resolution (Eichheimer et al. 2019). Then, experimentallymeasured permeabilities of sintered glass bead media
with a wide range of porosities were compared with thepermeabilities computed by LaMEM using subvolumesof micro-computed tomographic images of the glass
bead media (Eichheimer et al. 2020). The computed values agreed to an order of magnitude with the measured
values, demonstrating that the numerical model works
well not only for a simple pipe flow but also for a porous
flow. As the quartzite pores in this study were finer thanthose of the glass bead media, the computed fluid velocity error was expected to be larger than that for the glass

Fig. 9 Flow rate concentration to a single channel qmax
/Q against connected porosity (ϕcon) of each subvolume. Bars and solid circles representithe range and mean of each subvolume, respectively

Fujita et al. Progress in Earth and Planetary Science(2024) 11:40

bead media. Nevertheless, the computed values fit wellwith the empirical permeability–porosity relationshipobtained from laboratory measurements (Wark and Watson 1998). This empirically indicates the applicability ofpermeability computation using the LaMEM for porous
media featuring small pore spaces.
4.2 Fluid channel radius

Yokoyama and Takeuchi (2009) demonstrated that thecharacteristic maximum pore throat radius strongly correlates with the permeability of porous volcanic pyroclasts because the flux of a single flow channel is limited
by its pore throat radius (i.e., the minimum aperture size)
along the flow, whereas the total flux through the sample is strongly controlled by the maximum throat radiusof the available channels in the sample. In the present
study, we used a simplified method to estimate the maximum pore throat radius characteristic of the subvolumes
with the SRμCT slices by gradually eroding the porepixel layers until the pore interconnection broke off. The
number of erosion layers at which the interconnections
between the two facing sides of subvolumes are finally lost(i.e.,, ϕcon = 0), defined as Rmin, approximates the characteristic maximum pore throat radius of the fluid tubules.
Figure 10 shows the relationship between the number of
erosion layers and the connected pore volume divided by
the total pore volume (Vcon /Vtotal). For example, if thepore interconnection is lost after two layers were eroded(Rmin = 2 voxels = 1.02 µm), the maximum pore throat area
corresponds to 5–16 voxels. The results shown in Fig. 10

Page 13 of 16

indicate that Rmin is not always higher for subvolumes withlarger connected porosities. For example, in the wetting
system, QW14-21#1 has no interconnection after three layers of erosion even featuring a relatively high ϕcon (0.067),
while QW15-21#2 with a lower ϕcon (0.057) still maintained the interconnection even after erosion of four layers.
This may explain why QW15-21#2 had higher permeability
than QW14-21#1. It is clear that the wetting system tends
to have a larger Rmin than that of the non-wetting system.
The maximum Rmin for most of the wetting system was
1.02–2.04 μm, whereas that for the non-wetting system
was 0.51–1.02 μm.
4.3 Small power‑law exponent

The power-law exponent n obtained from the k–ϕcon relationship of the non-wetting system was not determined
to be a single value but could only be fit by two straightlines with slopes increasing from 0.76 to 3.0 at aroundϕcon = 0.04 (Fig. 5b). A small n value indicates that permeability is insensitive to changes in porosity. Generally,
a power-law exponent n of < 2 reflects the failure of new
interconnections with increasing porosity (Zhu and Hirth2003), being consistent with the flow focusing in the small
ϕcon range. In a rock in which all channels are represented
by N pipes with radius r , the absolute permeability can be
expressed as follows:

kabs = N

πr 4
,
8A

(10)

Fig. 10 Degree of pore interconnection. Degree of pore interconnection expressed as the connected pore volume divided by the total porevolume (Vcon /Vtotal), for a wetting (XCO2 = 0) and b non-wetting (XCO2 = 0.28–0.44) systems, respectively. Only subvolumes having a pore connection
from one side to the opposite side are shown. The horizontal axes are the number of erosion layers of voxels removed from the pores’ surface.
Erosion was performed until the final pore lost interconnection. The minimum pore throat radius is given by the thickness corresponding
to the maximum number of erosion layers (the right dead end of the symbols) plus 0–1 layer (0.51 µm)

Fujita et al. Progress in Earth and Planetary Science(2024) 11:40

Page 14 of 16

where A represents the cross-sectional area of the rock
(Mavko et al. 2020), and the tortuosity is assumed to be 1
for simplicity. When the number of channels in the crosssection is N , it is related to ϕ as:

N π r 2 = Aϕ.

(11)

Combining these equations and Eq. (1), N as a functionof ϕ is solved as:

N (ϕ) =

AC.
8π d 2 ϕ n−2(12)

An n < 2 is needed for N to increase with increasing
porosity ϕ. Therefore, n = 0.76 (< 2) for the non-wetting
system is plausible. At ϕ ~ 0.05–0.06, where the variationin pore connectivity exists, permeability is controlled bypore connectivity rather than porosity.
4.4 Fluid isolation in the polycrystalline rocks and its
implications for the fluid transport to deep earth

The results of this study depicted the processes in whichinterstitial pore fluids are left isolated in the experimental polycrystalline rocks. In the reference models of bothwetting and non-wetting systems, all fluid tubules were
uniformly narrowed with decreasing fluid fraction until
the pinch-off threshold, where the pore interconnection was lost at every grain-edge tubule. Figures 2 and 3
demonstrate, however, that the grain-edge tubules werepartly cut off at a fluid fraction near the threshold, while
a limited number of fluid tubules across the subvolume
were maintained. This feature in the experimental systemsuggests that the aperture size of the grain-edge tubules
was not uniform. Relatively wider tubules maintained thepore interconnection, while the narrower tubules werepinched off. The normalized Nch increased with increasing ϕcon both in the wetting and non-wetting systems(Fig. 8) also supported a picture in which the fluids gradually wetted the new grain edges with increasing fluid
fraction, even when the mean dihedral angle was < 60°,
where all the grain edges were supposed to be wetted
in the reference model. The maximum Nch , the numberof channels in a slice of a subvolume, was smaller thanthat for homogeneous hexagonal packing (70) in all subvolumes even at the highest ϕcon . This supports the BSE
image observation (Fig. 1) that a considerable proportionof triple junctions were dry or narrower than the resolution of SRμCT.
This study experimentally simulated grain growth andgeological fluid segregation in deep-seated rocks, which
include dehydrating slabs, rocks under prograde metamorphism, and partially molten upper mantle. The rangesof fluid fractions in these settings are often lower than afew percent (e.g., Ganzhorn et al. 2019). The results ofour study confirm that the amount of unsegregated fluid
left in practical rocks with heterogeneity and interfacial
energy anisotropy is not determined solely by MEFF.
Even in a wetting system where all interstitial fluids are
expected to interconnect for gravitational segregationin an ideal isotropic case, a percolation threshold wasobserved for the total fluid fraction (0.037 in volume).
This result suggests that fluids can be transported into
the deep Earth’s interior in a wetting system, contrary to
the prediction made with a homogeneous isotropic system. In the non-wetting system, the heterogeneous cutoff of the pore interconnection during fluid segregation
resulted in an isolated fluid fraction of 0.048–0.062. Thisvalue indicates the amount of non-wetting fluid transported into the deep Earth’s interior in subduction zones,
although the specific amount depends on many parameters such as lithology, fluid composition, and fluid separation processes.

5 Conclusions
To better understand the segregation mechanisms andretention fraction of aqueous fluids in deep-seated rocks,
we applied a combined approach of 3D imaging using
SRμCT and LaMEM computations of permeability andflow velocity distribution in the synthetic quartzite as asemi-analog of polycrystalline silicate rocks. In the wetting system (dihedral angle θ = 52°) with ϕmes > 0.037,
we found that the calculated permeability–total porosity (ϕmes) relationship was mostly consistent with thoseobtained by laboratory measurements using quartzite(Wark and Watson 1998). Below this ϕmes (percolationthreshold), fluid interconnection was lost, mainly owing
to the interfacial energy anisotropy of quartzite, and thepore fluids were left isolated. In the non-wetting system(dihedral angle θ = 61–71°), grain-scale fluid pools wereeffectively connected to allow similar permeability to thatin the wetting system at a fluid fraction > 0.062. Below
this fluid fraction, the cutoff of heterogeneous fluid
tubules formed the isolated fluid, and the percolationthreshold was 0.048. Streamline visualization revealed
that the fluid flow just above this threshold was focused
into a few channels, establishing efficient channelized
fluid pathways. These results demonstrated that the isolated fluids after the interconnection was lost in the wetting system can reach ~ 0.037, although the specific value
depends on the anisotropy and heterogeneity of the system. In the non-wetting system, on the other hand, thefluid isolation process with decreasing fluid fraction was
visualized for the first time. As a result, the fraction ofretained fluid was presented to be 0.048, providing an
estimate for the fraction of downdragged fluid in a subducting slab. This fraction was determined through grain
growth of polycrystalline rocks, beyond the previous

Fujita et al. Progress in Earth and Planetary Science(2024) 11:40

assessment solely based on the dihedral angle. The water
circulation in the Earth is thus potentially influenced byfactors such as lithology and fluid separation processes insubducting slabs.
Abbreviations
EDS	Energy-dispersive X-ray spectrometer
FE–SEM	Field-emission type scanning electron microscopy
BSE	Backscattered electron
LaMEM	Lithospheric and Mantle Evolution Model
SRμCT	Synchrotron radiation X-ray computed microtomography
Acknowledgements
This study was conducted as part of the Ph.D. project of WF in the joint graduate program of Tohoku University and the University of Bayreuth. We thank NP
Walte and T Nishiyama for fruitful discussions and T Nakano for helpful advice
on image analyses. We also thank two anonymous reviewers for constructive
comments that helped to improve the manuscript. The synchrotron radiation
experiments were performed at SPring-8 with the approval of the Japan
Synchrotron Radiation Research Institute (JASRI; proposal nos. 2018A1471,
2018A1464, and 2019B1785).
Author contributions
WF and MN designed the study and prepared the manuscript. WF performed
the experiments, analyzed the experimental products, and carried out thenumerical computations. GJG, MT, and PE contributed to establishing andapplying the digital rock physics approaches, and KU contributed to the
SRμCT imaging.
Funding
This study was supported by the JSPS Japanese–German Graduate Externship,
International Joint Graduate Program in Earth and Environmental Sciences,
Tohoku University (GP-EES), JSPS KAKENHI (Grant Numbers 25610154,
16H06348, and 22H00162 to M.N.), and the Ministry of Education, Culture,
Sports, Science, and Technology (MEXT) of Japan through its “Integrated
Program for Next Generation Volcano Research and Human Resource Development” and “Earthquake and Volcano Hazards Observation and Research
Program” to M.N.
Availability of data and materials
All data generated or analyzed during this study are included in this publishedarticle.

Declarations
Competing interests
The authors declare that they have no competing interests.
Author details
1
Department of Earth Science, Graduate School of Science, Tohoku University,
Sendai 9808578, Japan. 2 SPring-8, Japan Synchrotron Radiation Research Institute, Hyogo 679‑5198, Japan. 3 Bayerisches Geoinstitut, University of Bayreuth,
95440 Bayreuth, Germany. 4 Present Address: Math2Market, 67655 Kaiserslautern, Germany. 5 Present Address: Nuclear Waste Management Organization of Japan, Tokyo 108‑0014, Japan.
Received: 5 January 2024 Accepted: 21 May 2024

References
Bagdassarov N, Solferino G, Golabek GJ, Schmidt MW (2009) Centrifuge
assisted percolation of Fe-S melts in partially molten peridotite: time
constraints for planetary core formation. Earth Planet Sci Lett 288:84–95.
https://​doi.​org/​10.​1016/j.​epsl.​2009.​09.​010
Bear J (1988) Dynamics of fluids in porous media. Dover Publications Inc., New
York, p 764

Page 15 of 16

Beere W (1975) A unifying theory of the stability of penetrating liquid phases
and sintering pores. Acta Metall 23:131–138. https://​doi.​org/​10.​1016/​
0001-​6160(75)​90078-4
Bulau JR, Waff HS, Tyburczy JA (1979) Mechanical and thermodynamic constraints on fluid distribution in partial melts. J Geophys Res Solid Earth
84:6102–6108. https://​doi.​org/​10.​1029/​JB084​iB11p​06102
Cmíral M, Gerald JDF, Faul UH, Green DH (1998) A close look at dihedral angles
and melt geometry in olivine-basalt aggregates: a TEM study. Contrib
Miner Petrol 130:336–345. https://​doi.​org/​10.​1007/​s0041​00050​369
Eberhard L, Thielmann M, Eichheimer P, Néri A, Uesugi K, Nakamura M, Golabek GJ, Frost DJ (2022) A new method for determining fluid flux at highpressures applied to the dehydration of serpentinites. Geochem Geophys
Geosyst 23:e2021GC010062. https://​doi.​org/​10.​1029/​2021G​C0100​62
Eichheimer P, Thielmann M, Popov A, Golabek GJ, Fujita W, Kottwitz MO, KausBJP (2019) Pore-scale permeability prediction for Newtonian and nonNewtonian fluids. Solid Earth 10:1717–1731. https://​doi.​org/​10.​5194/​
se-​10-​1717-​2019
Eichheimer P, Thielmann M, Fujita W, Golabek GJ, Nakamura M, Okumura S,
Nakatani T, Kottwitz MO (2020) Combined numerical and experimental
study of microstructure and permeability in porous granular media. Solid
Earth 11:1079–1095. https://​doi.​org/​10.​5194/​se-​11-​1079-​2020
Faul UH (2001) Melt retention and segregation beneath mid-ocean ridges.
Nature 410:920–923. https://​doi.​org/​10.​1038/​35073​556
Fedorenko RP (1964) The speed of convergence of one iterative process.
J Comput Math Math Phys 4:227–235. https://​doi.​org/​10.​1016/​0041-​
5553(64)​90253-8
Fujita W, Nakamura M, Uesugi K (2022) Chemical compaction and fluid segregation in piston-cylinder experiments. Chem Geol 614:121182. https://​doi.​org/​10.​1016/j.​chemg​eo.​2022.​121182
Ganzhorn AC, Pilorgé H, Reynard B (2019) Porosity of metamorphic rocks
and fluid migration within subduction interfaces. Earth Planet Sci Lett
522:107–117. https://​doi.​org/​10.​1016/j.​epsl.​2019.​06.​030
German RM (1985) Liquid phase sintering. Springer, New York. https://​doi.​org/​
10.​1007/​978-1-​4899-​3599-1
Ghanbarzadeh S, Hesse MA, Prodanović M (2015) A level set method formaterials with texturally equilibrated pores. J Comput Phys 297:480–494.
https://​doi.​org/​10.​1016/j.​jcp.​2015.​05.​023
Glazner AF (2014) Magmatic life at low Reynolds number. Geology 42:935–
938. https://​doi.​org/​10.​1130/​G36078.1
Huang Y, Guo H, Nakatani T, Uesugi K, Nakamura M, Keppler H (2021) Electricalconductivity in texturally equilibrated fluid-bearing forsterite aggregates
at 800°C and 1 GPa: Implications for the high electrical conductivity
anomalies in mantle wedges. J Geophys Res Solid Earth. https://​doi.​org/​
10.​1029/​2020J​B0213​43
Huang Y, Nakatani T, Sawa S, Wu G, Nakamura M, McCammon C (2023) Effectof faceting on olivine wetting properties. Am Mineral 108:2244–2259.
https://​doi.​org/​10.​2138/​am-​2022-​8808
Ingebritsen SE, Manning CE (2010) Permeability of the continental crust:
dynamic variations inferred from seismicity and metamorphism. Geofluids 10:193–205. https://​doi.​org/​10.​1111/j.​1468-​8123.​2010.​00278.x
Iwamori H (1993) Dynamic disequilibrium melting model with porous flow
and diffusion-controlled chemical equilibration. Earth Planet Sci Lett
114:301–313. https://​doi.​org/​10.​1016/​0012-​821X(93)​90032-5
Jurewicz SR, Watson EB (1984) Distribution of partial melt in a felsic system: theimportance of surface energy. Contrib Mineral Petrol 85:25–29. https://​doi.​org/​10.​1007/​BF003​80218
Katayama I, Terada T, Okazaki K, Tanikawa W (2012) Episodic tremor and slow
slip potentially linked to permeability contrasts at the Moho. Nat Geosci
5:731–734. https://​doi.​org/​10.​1038/​ngeo1​559
Kita S, Okada T, Nakajima J, Matsuzawa T, Hasegawa A (2006) Existence of aseismic belt in the upper plane of the double seismic zone extending
in the along-arc direction at depths of 70–100 km beneath NE Japan.
Geophys Res Lett 33:1–5. https://​doi.​org/​10.​1029/​2006G​L0282​39
Liang Y, Price JD, Wark DA, Watson EB (2001) Nonlinear pressure diffusion in a
porous medium: approximate solutions with applications to permeability measurements using transient pulse decay method. J Geophys Res
106:529–535. https://​doi.​org/​10.​1029/​2000J​B9003​44
Lindquist WB, Venkatarangan A, Dunsmuir J, Wong T (2000) Pore and throat
size distributions measured from synchrotron X-ray tomographic imagesof Fontainebleau sandstones. J Geophys Res 105:21509–21527. https://​doi.​org/​10.​1029/​2000J​B9002​08

Fujita et al. Progress in Earth and Planetary Science(2024) 11:40

Liu X, Matsukage KN, Li Y, Takahashi E, Suzuki T, Xiong X (2018) Aqueous fluid connectivity in subducting oceanic crust at the mantle transition zone conditions. J Geophys Res 123:6562–6573. https://​doi.​org/​10.​1029/​2018J​B0159​73
Manning CE, Ingebritsen SE (1999) Permeability of the continental crust:
implications of geothermal data and metamorphic systems. Rev Geophys
37:127–150. https://​doi.​org/​10.​1029/​1998R​G9000​02
Mavko G, Mukerji T, Dvorkin J (2020) The rock physics handbook, 3rd edn. Cambridge University Press, Cambridge. https://​doi.​org/​10.​1017/​97811​08333​016
McKenzie D, Bickle MJ (1988) The volume and composition of melt generated byextension of the lithosphere. J Petrol 29:625–679. https://​doi.​org/​10.​1093/​
petro​logy/​29.3.​625
Mibe K, Yoshino T, Ono S, Yasuda A, Fujii T (2003) Connectivity of aqueous fluidin eclogite and its implications for fluid migration in the Earth’s interior. J
Geophys Res 108:2295. https://​doi.​org/​10.​1029/​2002J​B0019​60
Miller KJ, Zhu W, Montési LGJ, Gaetani GA (2014) Experimental quantification of permeability of partially molten mantle rock. Earth Planet Sci Lett
388:273–282. https://​doi.​org/​10.​1016/j.​epsl.​2013.​12.​003
Morishige M, van Keken PE (2018) Fluid migration in a subducting viscoelastic
slab. Geochem Geophys Geosyst 19:337–355. https://​doi.​org/​10.​1002/​
2017G​C0072​36
Mu S (2016) Faul UH (2016) Grain boundary wetness of partially molten dunite.
Contrib Mineral Petrol 171:40. https://​doi.​org/​10.​1007/​s00410-​016-​1250-z
Nakamura M, Watson EB (2001) Experimental study of aqueous fluid infiltration
into quartzite: Implications for the kinetics of fluid redistribution and grain
growth driven by interfacial energy reduction. Geofluids 1:73–89. https://​doi.​org/​10.​1046/j.​1468-​8123.​2001.​00007.x
Ohuchi T, Nakamura M (2007a) Grain growth in forsterite–diopside systems. Phys
Earth Planet Interior 160:1–21. https://​doi.​org/​10.​1016/j.​pepi.​2006.​08.​003
Ohuchi T, Nakamura M (2007b) Grain growth in the system forsterite–diopsidewater. Phys Earth Planet Interior 161:281–304. https://​doi.​org/​10.​1016/j.​pepi.​
2007.​02.​009
Ono S, Mibe K, Yoshino T (2002) Aqueous fluid connectivity in pyrope aggregates:
water transport into the deep mantle by a subducted oceanic crust without
any hydrous minerals. Earth Planet Sci Lett 203:895–903. https://​doi.​org/​10.​1016/​S0012-​821X(02)​00920-2
Osorno M, Uribe D, Ruiz OE, Steeb H (2015) Finite difference calculations of
permeability in large domains in a wide porosity range. Arch Appl Mech
85:1043–1054. https://​doi.​org/​10.​1007/​s00419-​015-​1025-4
Park HH, Yoon DN (1985) Effect of dihedral angle on the morphology of grains ina matrix phase. Metall Trans A 16:923–928. https://​doi.​org/​10.​1007/​BF028​
14844
Peacock SM, Christensen NI, Bostock MG, Audet P (2011) High pore pressuresand porosity at 35 km depth in the Cascadia subduction zone. Geology
39:471–474. https://​doi.​org/​10.​1130/​G31649.1
Price JD, Wark DA, Watson EB, Smith AM (2006) Grain-scale permeabilities of
faceted polycrystalline aggregates. Geofluids 6:302–318. https://​doi.​org/​10.​1111/j.​1468-​8123.​2006.​00149.x
Rudge JF (2018) Textural equilibrium melt geometries around tetrakaidecahedralgrains. Proc Math Phys Eng Sci 474:20170639. https://​doi.​org/​10.​1098/​rspa.​
2017.​0639
Skemer P, Chaney MM, Emmerich AL, Miller KJ, Zhu W (2017) Network topologyof olivine–basalt partial melts. Geophys J Int 210:284–290. https://​doi.​org/​
10.​1093/​gji/​ggx160
Solferino GFD, Thomson P-R, Hier-Majumder S (2020) Pore network modeling ofcore forming melts in planetesimals. Front Earth Sci 8:339. https://​doi.​org/​
10.​3389/​feart.​2020.​00339
Stevenson DJ (1990) Fluid dynamics of core formation. In: Newsom HE, Jones JH(eds) Origin of the earth. Oxford University Press, Houston
Terasaki H, Frost DJ, Rubie DC, Langenhorst F (2008) Percolative core formation inplanetesimals. Earth Planet Sci Lett 273:132–137. https://​doi.​org/​10.​1016/j.​
epsl.​2008.​06.​019
Thomson PR, Aituar-Zhakupova AA, Hier-Majumder S (2018) Image segmentation and analysis of pore network geometry in two natural sandstones.
Front Earth Sci 6:1–14. https://​doi.​org/​10.​3389/​feart.​2018.​00058
von Bargen N, Waff HS (1986) Permeabilities, interfacial areas and curvatures ofpartially molten systems: Results of numerical computations of equilibrium
microstructures. J Geophys Res 91:9261–9276. https://​doi.​org/​10.​1029/​
JB091​iB09p​09261
Waff HS (1980) Effects of the gravitational field on liquid distribution in partial
melts within the upper mantle. J Geophys Res 85:1815–1825. https://​doi.​
org/​10.​1029/​JB085​iB04p​01815

Page 16 of 16

Waff HS, Faul UH (1992) Effects of crystalline anisotropy on fluid distribution in
ultramafic partial melts. J Geophys Res 97:9003–9014. https://​doi.​org/​10.​1029/​92JB0​0066
Walte NP, Bons PD, Passchier CW, Koehn D (2003) Disequilibrium melt distribution
during static recrystallization. Geology 31:1009–1012. https://​doi.​org/​10.​
1130/​G19815.1
Wark DA, Watson EB (1998) Grain-scale permeabilities of texturally equilibrated,monomineralic rocks. Earth Planet Sci Lett 164:591–605. https://​doi.​org/​10.​1016/​S0012-​821X(98)​00252-0
Watson EB (1982) Melt infiltration and magma evolution. Geology 10:236–
240. https://​doi.​org/​10.​1130/​0091-​7613(1982)​102.0.​CO;2
Watson EB (2015) Lithologic partitioning of fluids and melts. Am Mineral
84:1693–1710. https://​doi.​org/​10.​2138/​am-​1999-​11-​1201
Watson EB, Brenan JM (1987) Fluids in the lithosphere, 1. Experimentally determined wetting characteristics of C
­ O2H2O fluids and their implications for
fluid transport, host-rock physical properties, and fluid inclusion formation.
Earth Planet Sci Lett 85:497–515. https://​doi.​org/​10.​1016/​0012-​821X(87)​90144-0
Watson EB, Lupulescu A (1993) Aqueous fluid connectivity and chemical-transport in clinopyroxene-rich rocks. Earth Planet Sci Lett 117:279–294. https://​doi.​org/​10.​1016/​0012-​821X(93)​90133-T
Watson HC, Roberts JJ (2011) Connectivity of core forming melts: Experimental
constraints from electrical conductivity and X-ray tomography. Phys Earth
Planet Interior 186:172–182. https://​doi.​org/​10.​1016/j.​pepi.​2011.​03.​009
Yokoyama T, Takeuchi S (2009) Porosimetry of vesicular volcanic products by awater-expulsion method and the relationship of pore characteristics topermeability. J Geophys Res 114:B02201. https://​doi.​org/​10.​1029/​2008J​
B0057​58
Yoshino T, Price JD, Wark DA, Watson EB (2006) Effect of faceting on pore geometry in texturally equilibrated rocks: Implications for low permeability at
low porosity. Contrib Mineral Petrol 152:169–186. https://​doi.​org/​10.​1007/​
s00410-​006-​0099-y
Zhu W, Hirth G (2003) A network model for permeability in partially molten rocks.
Earth Planet Sci Lett 212:407–416. https://​doi.​org/​10.​1016/​S0012-​821X(03)​00264-4
Zhu W, Gaetani GA, Fusseis F, Montési LGJ, De Carlo F (2011) Microtomographyof partially molten rocks: three-dimensional melt distribution in mantle
peridotite. Science 332:88–91. https://​doi.​org/​10.​1126/​scien​ce.​12022​21
Balay S, Abhyankar S, Adams MF, Brown J, Brune P, Buschelman K, Dalcin L, Dener
A, Eijkhout V, Gropp WD, Karpeyev D, Kaushik D, Knepley MG, May DA,
McInnes LC, Mills RT, Munson T, Rupp K, Sanan P, Smith BF, Zampini S, Zhang
H, Zhang H (2019) PETSc.s https://​www.​mcs.​anl.​gov/​petsc. Accessed 12 Apr
2024
Kaus BJP, Popov AA, Baumann TS, Pusok AE, Bauville A, Fernandez N, Collignon
M (2016) Forward and inverse modelling of lithospheric deformation on
geological timescales. John von Neumann Institute for Computing (NIC) 48.
Forschungszentrum Jülich, Jülich
Nakano T, Tsuchiyama A, Uesugi K, Uesugi M, Shinohara K (2006) "Slice"-Softwares
for basic 3-D analysis-, Japan Synchrotron Radiation Research Institute(JASRI), http://​www-​bl20.​sprin​g8.​or.​jp/​slice. Accessed 12 April 2024
Uesugi K, Hoshino M, Takeuchi A, Suzuki Y, Yagi N (2012) Development of fast andhigh throughput tomography using CMOS image detector at SPring-8. Proc
SPIE 8506. https://​doi.​org/​10.​1117/​12.​929575
Wark DA, Williams CA, Watson EB, Price JD (2003) Reassessment of pore shapes
in microstructurally equilibrated rocks, with implications for permeability of
the upper mantle. J Geophys Res 108 (B1). https://​doi.​org/​10.​1029/​2001J​
B0015​75
Wesseling P (1995) Introduction to multigrid methods, Institute for Computer
Applications in Science and Engineering, NASA Langley Research Center,
Hampton, ICASE Report, No. 95-1. https://​apps.​dtic.​mil/​sti/​pdfs/​ADA29​4291.​
pdf. Accessed 12 April 2024

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims in publishedmaps and institutional affiliations.