Fluid-driven metamorphism of the continental crust governed by nanoscale fluid flow.

The transport of fluids through the Earth's crust controls the redistribution of elements to form mineral and hydrocarbon deposits, the release and sequestration of greenhouse gases, and facilitates metamorphic reactions that influence lithospheric rheology. In permeable systems with a well-connected porosity, fluid transport is largely driven by fluid pressure gradients. In less permeable rocks, deformation may induce permeability by creating interconnected heterogeneities, but without these perturbations, mass transport is limited along grain boundaries or relies on transformation processes that self-generate transient fluid pathways. The latter can facilitate large-scale fluid and mass transport in nominally impermeable rocks without large-scale fluid transport pathways. Here, we show that pervasive, fluid-driven metamorphism of crustal igneous rocks is directly coupled to the production of nanoscale porosity. Using multi-dimensional nano-imaging and molecular dynamics simulations, we demonstrate that in feldspar, the most abundant mineral family in the Earth's crust, electrokinetic transport through reaction-induced nanopores (10-100 nm) can potentially be significant. This suggests that metamorphic fluid flow and fluid-mediated mineral transformation reactions can be considerably influenced by nanofluidic transport phenomena.

Metamorphic transformation processes and associated chemical differentiation of the solid Earth often require mass transfer over length scales that demand effective advective, fluid-mediated transport. It is therefore essential to understand the processes that generate pathways for fluid transport in the Earth’s lithosphere as well as alternate means of fluid motion. However, the intrinsic rock permeability within most of the lithosphere may well be too low (<10-17 m2)1 to transmit sufficient amounts of fluids to explain large-scale fluid-induced metamorphism and metasomatism. Thus, deformational processes such as grain-scale dilatancy, micro-cracking and hydraulic fracturing are often invoked to explain the formation of interconnected fluid pathways2. Mass transfer in such scenarios is commonly assumed to occur in response to gradients in the hydraulic potential1–3. In the absence of a continuous macroporous network, fluid connectivity along grain boundaries is controlled by surface energy reduction at grain-fluid interfaces4, 5. Mass transport at the scale of narrow grain boundaries (<1 nm) is limited to diffusive transport6. Substantial evidence, however, demonstrates that networks of fluid pathways (pores and fractures) in nominally impermeable rocks can emerge solely from coupled reaction-transport processes7–12, which indicates that reaction dynamics may have a more prominent role in spurring mass transfer than recrystallization driven by the minimization of surface energy differences.

At reaction fronts within complex porous networks, pressure gradient-driven fluid flow may become less effective in pores with nanoscale diameters (1 to 100s nm)11, 13–17 and the flow cannot be well described by the simple continuum fluid dynamics models applicable to permeable rocks with larger pore size. In this setting novel modes of fluid transport can arise from the combination of small pore sizes and electrostatic properties of the material interface18, 19. As the solid surface-to-volume ratio is substantially increased in nanoscale pores, surface-charge-governed transport phenomena18, which are ancillary to advective mass transfer in macroporous systems, begin to dominate. An example of this phenomenon is osmotically induced electric currents generated by the creation of salinity gradients through 40-nm pore channels with high surface energy. These currents have been shown to exceed charge transfer induced via fluid pressure gradients by many orders of magnitude20. This demonstrates that transport phenomena arising at the nanoscale can exert a physical influence on the ion transport properties of systems that are much larger in spatial extent.

The aim of the present work is to examine the role of ‘nanofluidic’ transport phenomena (nanofluidics refers to the study of the behaviour of fluids that are confined to structures of nanometer dimensions), which have initiated scientific and technological innovations20, 21–23, for mineral transformation processes. Herein, we investigate the development of nanoscale pore structures during pervasive (>60 km2) fluid-induced feldspar replacement reactions in the Larvik batholith24 situated in the Oslo rift area, Norway (Fig. 1a). Combining the pore structure information with molecular dynamics simulations we investigate the impact of different driving forces for fluid and mass transfer through feldspar nanopores (<1 µm).

Figure 1

Massive fluid-induced feldspar replacement in the Larvik batholith, Norway.

(a) Geological map showing the extent of feldspar replacement across 60 km2. Inset shows the geographic location of the study area. (b) Partly replaced larvikite. Unaltered feldspars are greyish to dark blue in b to d, whereas the replaced, secondary feldspars appear pink-coloured in b and ochre/greenish-coloured in c and d. (c) – (d) Magnified views showing the partial feldspar replacement. As the replacement goes to completion the original shape remains (pseudomorphism; dashed areas in (d)). Scale bars: b: 30 cm; c-d: 2 cm.

Reactive fluid flow through feldspathic igneous rocks

The Larvik batholith comprises a suite of plutons emplaced during intracontinental rifting in the early Permian (292-298 Ma)24. These igneous (syenitic to monzonitic) rocks predominately contain a single ternary feldspar (CaAl2Si2O8 – NaAlSi3O8 – KAlSi3O8) (Fig. 1b-d). Magma crystallization resulted in an igneous rock structure characterized by mm- to cm-sized, interlocking feldspar crystals with a strong blue iridescence (Fig. 1b-d). Permeabilities of igneous rocks similar to the ones reported here are between 10-17 to 10-21 m2 (ref1, 25) and can thus be regarded as nearly impermeable. Chemical modification of this larvikite leads to a prominent colour change from blue to ochre-greenish and red, making it a conspicuous visual reaction tracer. This discolouration proceeds along grain boundaries (Fig. 1c) and can result in the complete shape-preserving replacement (pseudomorphism) of an original feldspar crystal (Fig. 1d). Chemical analysis shows that the secondary, ochre/greenish, feldspar is an intergrowth of the feldspar end-members albite (NaAlSi3O8) and K-feldspar (KAlSi3O8) (Supplementary material 1). The initiation of feldspar replacement within the Larvik batholith is directly linked to the infiltration of externally-derived fluids along fractures and faults that are out of equilibrium with the host rock (Supplementary material 2). However, pervasive hydrothermal alteration occurred without significant deformation and affected an area of at least 60 km2 (Fig. 1). Fluids required to drive the hydrothermal alteration (300-500 °C, 1-3 kbar) were likely sourced from the intracontinental rift region. Comparison of whole-rock compositions (Supplementary material 3) suggest limited Na and K transport across scales much larger than the feldspar grain size (1-2 mm). Transport of Ca, released during the feldspar replacement, occurred across the scale of individual grains to form epidote (Ca2Al2(Fe3+,Al)(SiO4)(Si2O7)O(OH)) and other Ca-bearing phases. The absence of significant metasomatic mass transport, that is the co-infiltration of Na and K with externally-derived fluids or the massive removal of elements, makes the hydrothermal alteration of the Larvik batholith an ideal field example to extract the intrinsic characteristics of disequilibrium fluid-mediated feldspar replacement reactions. As there are a variety of possible feldspar replacement reaction pathways that can result in the observed replacement of albite (forward reaction) and K-feldspar (backward reaction) we only highlight the simplest reaction form Whether nearly isochemical or metasomatic (i.e., influx of Na and K) these fluid-mediated feldspar replacement reactions require the transport of aqueous species that can occur at the scale of individual feldspar grains (tens to hundreds of µm)26–28 to the scale of an entire plutonic complex (hundreds to thousands of m)12, 29, 30.

Reaction-induced feldspar nanopores

Electron microscopy analysis of partially replaced feldspar grains show that the original and secondary feldspar domains are divided by chemically and structurally sharp reaction interfaces (Fig. 2 a-d, Supplementary material 4). Nanopores (<1 µm) developed exclusively within the replaced, secondary feldspar (Fig. 2 e-g). Nanocrystalline Fe-(hydro)oxide precipitated within some of these pores resulting in the observed discolouration12 (Supplementary material 4 and 5). In three dimensions (Fig. 3 a) the pore network is defined by spatially aligned nanopore channels directed towards and ending at the reaction interface. Pore channel orientation is consistent with a preferred alignment parallel to the [100] feldspar direction28. Analysis by nanotomography (Fig. 3 b) in combination with atomic-scale imaging (Fig. 2 e-g) suggests that the average pore channel diameter is about 100 nm with a minimum diameter of 10-20 nm. In its present form, individual pore channels (x-axis in Figure 3) have been transformed into disconnected cavities with high aspect ratios. Disconnection of the pore space is likely due to a Rayleigh-type instability induced by surface energy minimization31. The above described microstructures are consistent with experimental26–28 and other natural12, 29 investigations of isochemical and metasomatic, fluid-mediated feldspar replacement reactions. The replacement took place via an interface-coupled dissolution-precipitation mechanism in which the supply of fluid that is out of equilibrium with the original ternary feldspar and the creation of porosity is essential for the reaction to proceed11. Porosity generation can be a result of a negative solid volume change between parent and product mineral and/or by differences in mineral solubility11, 26. Additionally, porosity generation is enhanced by mass transport of reaction products away from the reaction site to, e.g., form secondary minerals such as epidote and muscovite (KAl2(Si3,Al)O10(OH)2). The dynamic and transient nature of the pore network during fluid-mediated replacement reactions has been demonstrated using an ionic salt (KBr-KCl) analogue system11, 32. In this case, interface-orientated, dead-end pore channels also develop ahead of a more complex pore network, facilitating reaction progress through the continuous supply of fluid to the reaction front.

Figure 2

Micro- and nano-structures of hydrothermal feldspar alterations and mineral replacement interfaces.

(a) Backscattered electron image of a partly replaced ternary feldspar. (b)-(c) Show sharp reaction fronts dividing pore-free, original feldspar and pore-bearing, secondary feldspar. The replacement direction is given by arrows. (d) Bright-field transmission electron microscopy (TEM) image showing the replacement interface and the crystallographic continuity of the original feldspar twinning. (e) Pore distribution in secondary feldspar. (f) and (g) high resolution TEM images of two feldspar nanopores. Scale bars: a: 1 mm; b: 50 µm; c: 20 µm; d: 1 µm; e: 1 µm; f-g: 50 nm.

Figure 3

Nanopore channel network in replaced, secondary feldspar.

(a) Three-dimensional visualization of pore network using focused ion beam-assisted nanotomography (total volume: 14.22 µm × 12.65 µm × 11.92 µm). (b) Representative probability density of the pore channel diameter in the secondary feldspar. Pore channel analysis of the nanotomography volume shows that the average pore channel diameter is about 100 nm with a minimum diameter of 10-20 nm based on high-resolution transmission electron microscopy (Fig. 2 e-g).

Fluid and mass transport through feldspar nanopores

Transport through a reactive porous medium with dynamically evolving nanoporosity will likely not be limited only by gradients in fluid pressure ∇P. Fluid and mass transport may also be driven by other potential gradients, including dissolved ion concentrations ∇C and electric field −∇Φ (Fig. 4 a-d). Although each of these gradients produces an discrete force on a fluid and its dissolved constituents, which drives transport, all of these forces are inherently linked to each other. Whereas electro-osmotic fluid and mass transport is a well-studied phenomenon that is active in pores18, 19, flow as a result of concentration gradients in the fluid is commonly only associated with ion-selective osmosis through a semipermeable membrane. Here the H2O molecules can move only from low to high dissolved ion concentrations. However, in a fully permeable nanopore with charged wall surfaces, H2O molecule movement can occur as a result of a gradient in chemical potential of the dissolved ions (k T ∇ ln C, where k is the Boltzmann constant and T the temperature) between the pore ends in the direction opposite to the ‘usual’ direction of H2O molecule movement in a classical osmotic process20, 33. This diffusio-osmotic flow is produced by the electrostatic attraction of dissolved ions to the pore wall, which causes a locally high dissolved ion concentration (diffuse double layer)18,19,20. As a result, an osmotic pressure gradient is generated at the pore wall in the same direction as the gradient in dissolved ion concentration. This osmotic pressure gradient can drive fluid and mass transport with rates on the order of 1-10 µm per second34. Thus, to develop an atomistic understanding of fluid-mediated feldspar replacement reactions it is necessary to evaluate the impact of different driving forces (∇P, ∇C, −∇Φ) on the bulk fluid flux through feldspar nanopores. Here we use non-equilibrium molecular dynamics (NEMD) to simulate fluid and mass transport through feldspar nanopores with widths of 30, 50 and 100 nm, respectively. The term ‘non-equilibrium’ is used because fluid and mass transport occurs as a result of external fields applied to the molecules35, 36. The external fields were chosen so that the system response occurs in a linear fashion (Methods). In addition, we used equilibrium MD simulations to validate the choice of the external fields (Supplementary material 7). This approach allows us to determine individual transport properties on a molecular scale as a result of ∇P, ∇C and −∇Φ with minimal assumptions. All MD simulations were carried out with the open-source LAMMPS simulation package37 in a canonical (NVT) ensemble. The pore fluid was simulated as H2O containing K+ and Cl- corresponding to a solution with a concentration of 0.1 M at 500K and 20 MPa. Pore walls were modelled as orthoclase feldspar, representing one of the nano-porous replacement products (Fig. 2). Using orthoclase also allowed us to test the validity of our feldspar-H2O interface model against previous MD simulations38 and synchrotron X-ray reflectivity experiments39 that have characterized the orthoclase-H2O interface. For simplicity, we do not consider feldspar reactivity and atomic-scale surface roughness, thus neglecting hydrodynamic slip effects40. Further details about the feldspar surface generation and MD simulations as well as their boundary conditions can be found in Methods and Supplementary material 6-8.

Figure 4

Molecular dynamics simulations of hydrodynamic and electrokinetic fluid flow through a feldspar nanopore.

(a) Simulation box snapshot with 3D periodic boundary conditions: teal, Cl; blue, K; red, O; white, H; yellow, Si; magenta, Al. (b) Velocity profiles of H2O, K+ and Cl- induced by ∇P, k∇ 1n C and −∇Φ. (c) Comparison of volumetric fluxes achievable via and through a cylindrical feldspar pore (diameter 100 nm; length 1 µm). (d) Volumetric flux as a function of pore diameter with given boundary conditions of , Analytical predictions (An. pred.) are in agreement with NEMD.

Figure 4 b shows the velocity profiles of H2O and its dissolved ions (K+, Cl-), obtained by NEMD, through a 30-nm feldspar slit pore generated by gradients in ∇P, k T ∇ ln C, −∇Φ. Both, hydrodynamic and diffusio-osmotic velocity profiles are very similar. Electro-osmotic velocities show a separation of the dissolved ions and strong velocity changes very close (<5 nm) to the charged feldspar surface. It should be noted that the velocities produced by the simulation (Figure 4 b) are much larger than those expected in natural systems but are required to remove data noise caused by thermal oscillations of the H2O molecules. However, as the system responded linearly, the obtained values can be scaled to continuum descriptions. Our NEMD results are in agreement with analytical predictions34, 41 of transfer coefficients that couple fluxes of fluid and dissolved ions through a charged pore under the effect of ∇P, ∇ ln C, −∇Φ (Supplementary material 8). On the basis of the Nernst-Planck equation we can define the overall volumetric bulk fluid flux J [m3s-1] through a pore channel as a combination of the individual hydrodynamic diffusio-osmotic and electro-osmotic fluxes as follows Here K is the hydrodynamic [m4 Pa-1 s-1], K the diffusio-osmotic [m Pa-1 s-1], and K the electro-osmotic [m4 V-1 s-1] transport coefficient derived from the MD simulations (Supplementary material 6-7). Figure 4 c shows J through a cylindrical 100 nm feldspar pore with a length of 1 µm as a function of its individual contributions from gradients in ∇P, ∇ ln C and −∇Φ. Our simulation reveals that, e.g., a pressure drop of 0.1 MPa will produce the same volumetric bulk fluid flux J = 1.5⋅10-15m-3 s-1 as a concentration difference in dissolved ions (C2/C1) of only 18% between the pore ends. This strong contribution of diffusio-osmotic flow to J is found across the entire range of pore channel diameters from 5 to 100 nm (Fig. 4d). At pore diameters <50 nm the contribution of electro-osmotic flow to the overall volumetric flux becomes even more important (Fig. 4 d). The pressure drop of 0.1 MPa has been chosen as an example as it is well within the laminar regime of fluid flow (Reynolds number = 10-6) through the cylindrical feldspar nanopore.

The NEMD simulations above show that diffusio-osmosis can have a considerable effect on bulk fluid fluxes parallel to the flow direction through a feldspar nanopore and that neglecting electrokinetic flow can result in a large deviation from the actual fluid and mass flux in a reaction-induced pore network. A simple calculation, based on our MD simulations, shows that a 0.1 M KCl fluid within a 10x10x100 nm pore contains 1125 K atoms at 500 K and 20 MPa. Thus only ~200 atoms, equivalent to C2/C1 = 18%, must be released or consumed during the dissolution-precipitation-mediated feldspar replacement reaction to induce a diffusio-osmotic flux that is equivalent to a 0.1 MPa pressure difference (Fig. 4 c). An orthoclase pore wall surface with a diameter of 100 nm exposes 80,000 K atoms. Assuming a dissolution rate of 5-50 atoms per second42 a concentration gradient across a pore, connecting dissolving and precipitating mineral domains, will be established nearly instantaneously. Moreover, diffusio-osmotic bulk fluid transport during mineral replacement reactions will likely be aided by the exposure of negatively charged feldspar surfaces and transient dipole moments such as the one normal to the (010) feldspar surface39. However, our simulations with charge-neutral surfaces show that diffusio-osmosis is not solely dependent on charged surfaces (Supplementary Fig. S8), which may diminish due to environmental changes (e.g., pH) and surface re-organization, consistent with previous theoretical predictions35.

Reaction-induced pore networks will consist of both interconnected and dead-end pores. By using MD simulations of a single nanopore open at both ends we have established a general concept of the impact of different driving forces to the overall fluid flux through a feldspar nanopore. However, in dead-end pore geometries orientated towards the reacting fluid-mineral interface flow through ∇P is impossible and the driving forces associated with ∇C will fluctuate as a result of continuous dissolution and precipitation. In such a setting the solute gradient ∇C along the pore, simulated here (Fig. 4), can be viewed to be similar to a high concentration of dissolved ions close to the reaction interface due to mineral dissolution and a low concentration away from the mineral interface driving flow. Indeed, recent experiments on the KBr-KCl salt analogue system suggest that convective motion in dead-end pore channels occurs as a result of self-induced diffusio-osmotic transport during mineral replacement reactions43. By tracking quantum dots (20 nm in diameter) released from the reacting crystal, the authors of the study suggested that the salt replacement rate is limited by ionic transport in and out of pores terminating at the reaction interface. Within these pores ionic gradients occur between the bulk fluid and the reaction front facilitating convective transport as a result of a self-generated electric field that occurs to maintain local charge balance within the fluid. As silicate dissolution (e.g., feldspars) is considerably slower than salt dissolution, interfacial reaction kinetics will contribute to the overall control of diffusio-osmotically driven fluid transport in reaction-induced pores. To advance our understanding of these transport processes and evaluate their contribution to the kinetics of fluid-driven metamorphism and metasomatism future studies need to couple transport coefficients derived at a molecular level to ab-initio simulations of reaction rates. These can be used in continuum-scale transport models that account for the crystallographic evolution of the pore network, encompassing dead-end geometries close to the reacting mineral interface and interconnected porosity as well as grain boundaries in its wake.

Consequences for reactive fluid transport in the solid Earth

Our investigations show that feldspar replacement reactions induced by fluid infiltration provide a record of fluid-rock interactions that can affect large volumes of the continental crust. These pervasive replacement reactions are mediated by the development of nanoporosity (<1 µm). Interconnected porosity allowing fluid to penetrate to the reaction front emerges and evolves dynamically as a result of dissolution and precipitation and could thus be independent of tectonic forcing and deformation. Determination of fluid transport coefficients based on molecular dynamics simulations suggest that fluid and mass transport through reaction-induced feldspar nanopores can be self-generated through the emergence of concentration gradients and may not rely on fluid pressure gradients. Only by linking the feedback of different driving forces (equation 2) with molecular dissolution and precipitation processes and the dynamic evolution of porosity can we begin to evaluate the impact of electrokinetic transport processes on metamorphic reactions. Recent experiments and theoretical considerations43, 44 indicate that diffusio-osmotic flow may enable convective mass transport in and out of pore structures that terminate at the reaction interface independent of pressure and temperature gradients. With increasing recognition of the universality of dissolution-precipitation-mediated mineral transformations8, 11, 45, 46 electrokinetic transport may have a profound effect on systems where reactive fluids interact with the solid Earth. This may hold not only for metamorphic environments but also for engineered subsurface systems47, 48 where (reactive) fluid transport through low-permeability rocks is critical.

Methods

Electron microprobe analysis

Chemical analysis of the original and hydrothermally altered feldspar grains was conducted in a JEOL JXA-8600 superprobe using an acceleration voltage of 15 kV. Standardization was carried out using a selection of natural minerals and synthetic phases.

Focused ion beam scanning electron and transmission electron microscopy

Electron-transparent thin foils were prepared for (scanning) transmission electron microscopy ((S)TEM) by using a FEI Nova Nanolab focused ion beam - scanning electron microscope (FIB-SEM). The FIB-SEM was also used to acquire a slice-and-view series for 3D volume reconstructions. Slice imaging was carried out in backscattered electron mode at 2 kV and 0.84 nA with a voxel size of 8.33×8.33×20 nm³. All FIB-SEM nanotomography volumes were reconstructed and analysed using FEI Avizo 9. Pore channel diameters were obtained by using the cross-correlation diameter obtained via FEI Avizo 9. At the given microscope conditions, we determined a lower boundary limit of the detectable pore size of 50 nm for the analysed nanotomography volumes (Fig. 3B). Electron-transparent FIB foils were investigated in a FEI Talos F200X (S)TEM equipped with four energy-dispersive X-ray detectors (Super-X EDX). The FEI Talos F200X TEM information limit is 0.12 nm. All FIB-SEM and TEM analyses were carried out at the Microscopy Square, Utrecht University.

Feldspar surface model generation

We constructed a model for a slab of orthoclase (KAlSi3O8) with the following experimental unit cell parameters49: a=8.554 Å, b=12.97 Å, c=7.207 Å, α=90°, β=116.01°, γ=90°. A (001) orthoclase slab was generated and cut along the cleavage plane β (for more details see Kerisit et al.38). Hydrogen atoms were added to form hydroxyl groups on both slab surfaces. The surface charge of -0.25 e/nm2 was achieved through tetrahedral substitution of Si by Al. The slab consisted of 9×9 unit cells with a thickness of ~25 Å and was divided into two parts located at the top and bottom of the simulation box forming a section of a slit-like pore (Fig. 4A). The pore space was filled by water, K+ and Cl- corresponding to a 0.1 M solution at 500K and 200 bar. All parameters (Lennard-Jones, Buckingham and Coulomb potential) describing the fluid-fluid and fluid-solid interactions are taken from Kerisit et al.38. The positions of all orthoclase atoms were fixed.

Molecular dynamics simulations

Molecular dynamics (MD) simulations (including non-equilibrium dynamics) were carried out with the open-source LAMMPS37 simulation package in a canonical (NVT) ensemble. A constant temperature of 500 K was maintained with a Nosé-Hoover thermostat. The equations of motion were integrated using a 1 fs time step and the SHAKE algorithm. Following Kerisit et al.38 we described the short-range potential parameters and ionic charges by using the extended simple point charge (SPC/E) model50 for water and a modified CLAYFF force field51 for the feldspar. Long-range electrostatic interactions were calculated using the particle-particle particle-mesh (PPPM) method with a tolerance of 10-5. To investigate the fluid transport properties we simulated, using non-equilibrium MD, a pressure gradient, chemical potential gradient of the solute and electric field parallel to the surface by applying constant external forces to the fluid atoms36. The absolute values for the individual external forces Fx were chosen to result in a linear response of the system, i.e. Here v is the fluid velocity in the slab center. In addition, we used equilibrium MD simulations to validate the selection (Supplementary material 7). MD simulations were performed for 3D pore slabs in a 2D geometry (slit pores) with the following dimensional (x × y × z) boundary conditions; the x and y pore slab dimensions were set to 10 nm, whereas z varied between 30 to 100 nm. The (x,y) plane are solid feldspar walls with periodic boundary conditions. The driving forces applied across the pores are as follows. For pressure-driven flow we used forces of 1.25×10-5, 7.5×10-6 and 4×10-7 kcal/mol/A. The chemical potential gradient of the solute was achieved by applying 0.01, 0.03 and 3×10-4 kcal/mol/A to the ions. Electro-osmotic flow was simulated by applying an electric field of 0.03 V/A to all pores. The velocities of the particles are collected every ps for 10 ns, after an equilibration period of 25 ns. The estimated volume flow rates for cylindrical pores were obtained by integrating the 2D flux profiles J, obtained via MD, of the slit pores by where L is the pore length and z the direction perpendicular to the pore surface.