# Reactive Transport Modeling in Porous Fractured Media: Application to Weathering Processes

Sylvain Favier, Yoram Teitler, Fabrice Golfier, Michel Cathelineau

► To cite this version:

Sylvain Favier, Yoram Teitler, Fabrice Golfier, Michel Cathelineau. Reactive Transport Modeling in Porous Fractured Media: Application to Weathering Processes. Water Resources Research, 2024, 60 (2), pp.1-19. 10.1029/2023WR035283. hal-04715736

HAL Id: hal-04715736

https://hal.science/hal-04715736v1

Submitted on 1 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Water Resources Research®

## RESEARCH ARTICLE

10.1029/2023WR035283

**Key Points:**

* Dual Porosity modeling of a fractured rock mass recreates weathering heterogeneities

* Weathering heterogeneities originate from the physico-chemical differentiation of the fractures and the matrix

* Discrete Fracture Matrix modeling reveals the impact of the geometry and permeability of fractures on weathering heterogeneities

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

S. Favier,
favier.syl@wanadoo.fr

**Citation:**

Favier, S., Teitler, Y., Golfier, F., & Cathelineau, M. (2024). Reactive transport modeling in porous fractured media: Application to weathering processes.
*Water Resources Research*, *60*, e2023WR035283. https://doi.org/10.1029/2023WR035283

Received 5 MAY 2023
Accepted 26 JAN 2024

**Author Contributions:**

**Conceptualization:** Sylvain Favier

**Funding acquisition:** Michel Cathelineau

**Investigation:** Sylvain Favier

**Methodology:** Sylvain Favier,
Fabrice Golfier

**Project administration:** Fabrice Golfier,
Michel Cathelineau

**Resources:** Yoram Teitler, Fabrice Golfier

**Supervision:** Yoram Teitler,
Fabrice Golfier, Michel Cathelineau

**Validation:** Sylvain Favier,
Yoram Teitler, Fabrice Golfier,
Michel Cathelineau

**Visualization:** Sylvain Favier

**Writing – original draft:** Sylvain Favier

© 2024. The Authors.

This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.

<figure><img src="image_2.png" /></figure>

## Reactive Transport Modeling in Porous Fractured Media: Application to Weathering Processes

Sylvain Favier¹ ID, Yoram Teitler¹, Fabrice Golfier¹ ID, and Michel Cathelineau¹ ID

¹Université de Lorraine, CNRS, GeoRessources, Nancy, France

**Abstract** Complex permeability distribution is frequently encountered in weathered rocks with inherited heterogeneities, such as networks of discontinuities that impact fluid flow and water-rock interactions. Reactive transport modeling is a powerful tool for analyzing fractured media's chemical weathering dynamics. In this study, two different ways to model the fracture-matrix interaction were explored. A 1D dual-porosity simulation of reactive transport is designed to study the global impact of fractures on the weathering progression front, compared with a Single Porosity (S.P.) (unfractured) reference model. Then a 2D discrete model using a Discrete Fracture Matrix (DFM) representation is built to highlight the importance of medium geometry and fracture connectivity on weathering. The relevance of such models was evaluated in the case of weathering processes affecting ultramafic units in New Caledonia. Comparing the average behavior of the dual and S.P. models reveals a similar trend in both global mineralogical evolutions. Yet differences occur when confronting the weathering rates. By extending the weathering over a greater depth, the fractures favor a smoother transition between the saprolite horizon and the oxide ore. Furthermore, a partition of the newly formed phyllosilicates occurs in fractures and matrix: Ni-phyllosilicates preferentially precipitate in fractures, while Mg-phyllosilicates remain trapped within the matrix due to differences in residence time and liquid-to-solid ratio. The DFM modeling highlights the effect of the connectivity of a complex fracture network on the dissolution-precipitation processes. Depending on the permeability and geometry, fluid flow is modified, favoring the formation of weathering heterogeneities typically observed on the field.

## 1. Introduction

Weathering is defined as the chemical degradation of mineral phases under surface or subsurface conditions (Viers & Oliva, 2018). It is interconnected with several processes that shape the Earth's surface and subsurface constitution and function (Beinlich et al., 2020; Frings & Buss, 2019). Earth's surface evolution and climate regulation are thus controlled by silicate weathering and their ability to sequestrate CO₂ (Beerling et al., 2020; Goudie & Viles, 2012; West et al., 2005). Weathering is also responsible for the formation of economic supergene deposits resulting from the concentration of highly demanded metals (Ni-Co-Sc, Fe, Cu, REE, PGE, U, Zn) enriched by successive processes of dissolution, transport, and re-precipitation (Cathelineau et al., 2016; Dill et al., 2010; González-Álvarez et al., 2016; Sillitoe, 2005; Teitler et al., 2019; Valeton, 1994). Yet, deciphering the different processes controlling weathering remains a significant challenge (West et al., 2005). Notably, the impact of chemical and physical heterogeneities on the weathering front progression is still yet not fully understood (Maher & Navarre-Sitchler, 2019). A pure “per descensum” model based on the rainwater infiltration explains the downward progression of the weathering front but fails to reproduce the heterogeneities observed in the field, such as boulder formation (Chigira & Hirata, 2021; Md Dan et al., 2015; Ollier, 2010). Thus, recent studies have highlighted fractures' critical impact in controlling dissolution and redistribution processes by creating preferential fluid flow pathways (Borg, 2009; Cathelineau et al., 2017; Favier et al., 2022; Steefel & Hu, 2022). However, the differences in nature and intensity of chemical interaction along a fracture network and within the surrounding rock matrix remain poorly addressed.

Due to the complexity of mechanisms involved in weathering, reactive transport models (RTMs) are a powerful tool to understand, quantify and predict the coupling effect of chemical and physical processes at different scales and across large temporal scales on weathering (Goddéris et al., 2019; Steefel et al., 2015). They shed light on the weathering front progression from the sample scale (Navarre-Sitchler et al., 2011) to the outcrop scale by successfully recreating the geochemical and mineralogical evolution in surface environments (Domènech et al., 2017; Heidari et al., 2017; Myagkiy et al., 2017). Furthermore, they provide reliable estimates of weathering
**Writing – review & editing:**
Sylvain Favier, Yoram Teitler,
Fabrice Golfier, Michel Cathelineau

rates using simple conceptual fractured or unfractured rock models, which are upscaled to derive macroscopic laws at the hillslope scale (Brantley et al., 2017; Lebedeva & Brantley, 2017; Reis & Brantley, 2017).

Most RTMs are 1D models restricted to the vertical progression of weathering front, but recent studies were extended to 2D and 3D configurations. These models mainly focus on the topographic impact on dissolution processes by considering lateral flows (Lebedeva & Brantley, 2017; Myagkiy et al., 2019). Some occasionally study the occurrence of a single fracture on the weathering patterns (Deng & Spycher, 2019; Myagkiy et al., 2019; Trinchero et al., 2021; Zhang et al., 2022). Yet, these models remain rare and do not address the issues of fluid mixing and flow channelization driven by complex discrete fracture network (DFN) geometries because of the challenges of explicitly modeling fracture-matrix interactions (Berre et al., 2019). However, recent advances in cutting-edge modeling techniques and the ease of access to computational power have opened new possibilities to RTM for fractured porous media applications. Different approaches, involving variable representations of both fractures and matrix, are thus available, as illustrated in Figure 1. The first solution uses an equivalent continuum porous medium model (ECPM) to consider the fracture network's upscaled contribution to the model's hydrological parameters. Such a formulation successfully recreates the upscaled hydrological heterogeneities resulting from the DFN and, therefore, provides the possibility to study their influence as preferential flow pathways or barriers (Pandey & Rajaram, 2016; Steefel & Hu, 2022; Trinchero et al., 2017; Viswanathan et al., 2022). Nevertheless, these representations are limited by mesh accuracy and cannot claim to capture precisely the mineral reactions at the fracture walls (Maher & Navarre-Sitchler, 2019).

In contrast to ECPM, the DFN models offer an accurate description of fluid flow and reactive transfer processes but only within the fractures (Hyman et al., 2015; Molson et al., 2012; Sherman et al., 2021). The Discrete Fracture Matrix (DFM) or hybrid models divide the domain into regions for which an explicit fracture formulation is coupled with a continuum representation of the matrix (Viswanathan et al., 2022), adding more resolution to the description, albeit at the cost of heavier numeric simulations. As an alternative, the Multi-Continua models (e.g., mobile-immobile, dual-porosity), also referred to as multi-rate mass transfer models, represent a compromise in terms of computational efficiency and ease of use (Davit et al., 2019; Golfier et al., 2007; Haggerty & Gor elick, 1995). The underlying idea is the introduction of multiple subdomains corresponding to the different characteristic times of processes at stake with a mass exchange term between regions. These mass transfer functions, however, cannot fully capture the complex flow behavior in fractured media, and the application to reactive transfer processes is also delicate (de Dreuzy et al., 2013; Orgogozo et al., 2013).

This paper investigates the impact of fracture networks acting as chemical and hydrological heterogeneities on the weathering front migration through different models of increasing complexity. First, a general approach with a 1D dual-porosity model is developed to highlight the impact of fractures on weathering patterns by comparing the results with an unfractured porous domain (Single Porosity (S.P.) model). In essence, this model is an upscaled representation of fractured reservoirs so that a 1D description is theoretically sufficient to capture fracture contribution provided that the transfer function which drives mass exchanges between fractures and matrix is correctly determined. For this purpose, a fine calibration of mass transfer coefficient is conducted to maintain a quantitative link between the underlying geological model (i.e., the periodic fracture network used with the DFM) and the effective parameters involved in this transport model. Then, a discrete approach is presented, where fractures are explicitly defined as 1D interfaces embedded in a 2D porous rock matrix. The effect of the fracture network connectivity on the redistribution of metals is thus studied, and potential hotspots of interest are iden tified. The results are discussed in light of field observations from Ni-laterite deposits in New Caledonia. Indeed, the local mining activity provides a rare opportunity to access up to 50 m-thick laterite profiles extending from the bedrock (fractured, olivine-rich peridotite) to the oxide zone.

1.1. Weathering Processes of Peridotites in New Caledonia

In New Caledonia (Figure 2a), the Peridotite Nappe underwent, after its obduction, deep weathering under a hot, wet tropical climate, resulting in laterite formation (Maurizot et al., 2020; Sevin et al., 2014). The development of Ni-rich laterites (Figure 2b) results from (a) the solubilization of Mg, Si, and Ni through olivine dissolution, (b) the trapping and concentration of Ni in the saprolite horizon as neo-formed phyllosilicates and Fe-oxy hydroxides, and (c) the formation of an oxide-rich horizon (mostly Ni-bearing goethite) after phyllosilicate dissolution (Butt & Cluzel, 2013; Cathelineau et al., 2017; Dublet et al., 2012; Freyssinet et al., 2005; Golightly, 2010; Haldar, 2017; Santoro et al., 2022; Teitler et al., 2019).
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Different modeling strategies depending on the accuracy of the fractures' representation and the matrix of a porous fractured medium. Discrete fracture network (DFM), equivalent continuous porous medium, DFM, Multi-Continua (the Multi-Continua). Adapted from Viswanathan et al. (2022).</figcaption></figure>

Recent studies highlight the structural control of the weathering front pro-
gression (Favier et al., 2022; Iseppi et al., 2018; Jeanpert et al., 2019) as either a
drain or a permeability barrier depending on the fractures coating and their
orientation to the fluid flow. The impermeable fractures are serpentine (liz-
ardite) filled fractures. Furthermore, outcrops are typically structured by sub-
orthogonal fracture sets (Figure 3b), which cross-cut each other, forming
dihedral geometries (Iseppi et al., 2018). By controlling the fluid flow zone,
these geometrical patterns allow the formation of weathering heterogeneities
such as the boulders in saprolite (Myagkiy et al., 2019). Such organization is
common for weathered soils (Figure 3a) and is also found in granitic or basaltic
environments (Chigira & Hirata, 2021; Md Dan et al., 2015; Ollier, 2010; Sak
et al., 2010). Boulder typically display a concentric weathering zonation: the
outer part (zone 2 in Figure 3c) is more weathered (zone 1 in Figure 3c) and
presents higher concentrations in Ni than the core (Figure 3d).

## 2. Modeling Strategy

2.1. Conceptual Models of Weathering in Fractured Media

The formation of residual laterite soil results from the complete dissolution of olivine, which is the predominant phase in the bedrock. The medium is assumed to be fully water-saturated by the downward migration of rainwater, which progressively dissolves olivine. Released silica, magnesium, and nickel partly reprecipitate as newly formed minerals: Ni- and Mg-phyllosilicates and goethite, constituting the saprolite horizon. These phyllosilicates are ultimately dissolved into goethite to form the oxide horizon. While precipitation of quartz have been reported in such geochemical systems, recent studies have highlighted the impact of evaporation to reach the saturation with respect to quartz (Cathelineau et al., 2023). Furthermore, preliminary studies regarding its precipitation in fully saturated media have highlighted its absence in the system (Myagkiy et al., 2017). Hence quartz was not included in the geochemical model. In this study, the porosity and permeability of the medium are assumed to be constant. A similar assumption has been used in the previous 1D single and D.P. simulations carried out by Myagkiy et al. (2017). Porosity, and hence permeability, evolve in complex ways both from dissolution-precipitation processes but also through compaction. Only considering porosity changes induced by mineral reactions could underestimate pore velocity, and coupling reactive transport with an elastoplastic constitutive model of soil is a highly challenging task. Furthermore, fractures are assumed to be much more permeable than the matrix, even in the saprolite and oxide horizon. Indeed,

<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) Geological map of New Caledonia (b) Typical weathering profile in New Caledonia (after Myagkiy et al., 2017).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a) Typical weathering profile organization (after Ollier, 2010) (b) Typical fracture network in New Caledonian peridotites. (c) Weathering zonation in a typical boulder. (d) Ni distribution inside the boulder.</figcaption></figure>

the permeability of the saprolite horizon, which is the most permeable layer, is estimated at around $10^{-13}$ m² at the pit scale, while fracture permeability ranges from $10^{-12}$ m² to $10^{-10}$ m² (Jeanpert, 2017). The main part of the flow is likely to be channeled in the fracture network. Therefore, not accounting for the porosity and permeability changes, should not significantly impact the outcomes of the simulation in a first-order approach.

Three different models are designed to capture the impact of a fracture network on weathering processes: (a) a 1D unfractured or S.P. model, hereafter named as S.P. Model, which serves as a comparison reference with previous studies (Domènech et al., 2017; Myagkiy et al., 2017), (b) a 1D D.P. model, named D.P. Model, to study the global behavior of fractured rock mass, and (c) a DFM model, which explicitly includes the geometry of fracture network in the modeling.

## 2.2. Model Descriptions

### 2.2.1. Single Porosity Model

The S. P. model proposes a 1D discretization of advective-diffusive mass transport along a porous column coupled with chemical reactions (Parkhurst et al., 2013a). The equation governing mass transport for a reactive chemical species is solved with PhreeqC and written as (Parkhurst & Appelo, 1999):

$$ \frac{\partial q}{\partial t} + \frac{\partial c}{\partial t} + v \frac{\partial c}{\partial x} - D_e \frac{\partial^2 c}{\partial x^2} = 0 \qquad (1) $$

Where c is the concentration in the solution (mol/L), t is the time (s), v is the pore water flow velocity (m/s), x is the distance along the column (m), $D_e$ is the effective diffusion coefficient (m²/s), and q is the concentration of solid phase calculated through the chemical reactions (mol/L). Note that this change in concentration in solid phase $\frac{\partial q}{\partial t}$ is calculated for each species separately from the transport part, as the result of all equilibrium and kinetically controlled reaction rates. The geochemical model is described hereafter in Section 2.3 and is the same for all the different transport models.

### 2.2.2. Dual Porosity Model

In contrast, the D.P. model divides the fractured porous medium into mobile regions of highly transmissive fractures and immobile regions corresponding to the rock matrix and the secondary fracture network. Such duality is easy to translate into fluid flow and tracer transport models. However, when dealing with reactive transport
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Dual Porosity (D.P.) formulation of the fractures. (a) General fractured model characterized by a fracture spacing *s*. (b) Definition of the advective area defined by the fracture characterized by its hydraulic aperture *a*<sub>f</sub> and a damaged zone characterized by its penetration length *L*<sub>pene</sub>, and the diffusive area defined as the rest of the matrix. (c) Discretization of the D.P. formulation into interconnected cells. The advective or mobile cells were both advection and diffusion occur, are connected to the diffusive or immobile cells by diffusion. Diffusive exchange between mobile and immobile cells is formulated in terms of a mixing process expressed as a first-order rate expression.</figcaption></figure>

modeling, the contrast in properties is not so obvious. Indeed, the drastically different reactive surfaces between fractures and rock matrices lead to reaction times per region behaving oppositely to the characteristic transport times. Furthermore, only a small part of the rock matrix is conventionally assumed to participate significantly in the reaction, due to diffusion limitations. In addition, the presence of damaged zones in the vicinity of fractures can enhance this contrast as they are more transmissive than the matrix (Laurent et al., 2019; Neretnieks, 1980; Wogelius et al., 2020). Such zones may be considered by sub-gridding the rock matrix, as done in the Structured INteracting Continua approach proposed by Babey et al. (2016). Hereafter, rather than introducing a third sub-domain in our multi-continua model, we merged the damaged zone area with the fracture region and used a D.P. approach. In other words, we consider two zones: (a) an advective zone, corresponding to the fractures and a portion of the matrix in their vicinity, where advection drives the mass transport, and (b) a diffusion zone, corresponding to the matrix where the diffusion is the dominant phenomenon (Figure 4). The advective zone is therefore characterized by the fracture hydraulic aperture *a*<sub>f</sub> (m) and a penetration length *L*<sub>pene</sub> (m) corresponding to the advection-dominated portion of the matrix. Hence, an advective zone length *d*<sub>f</sub> (m) is defined as (Figure 4):

$$d_f = a_f + 2L_{\text{pene}} \qquad (2)$$

Both zones are connected by an exchange factor (α) which controls the mass transfer between advective and diffusive zones. Thus, the weathering solution first interacts with the advective zone (fractures) before reacting only by diffusion with the matrix. The RTM is therefore modified as follows (Laurent et al., 2019; Parkhurst & Appelo, 1999):

$$\varphi^{\text{adv}} \varepsilon_f \frac{\partial q^{\text{adv}}}{\partial t} + \varphi^{\text{adv}} \varepsilon_f \frac{\partial c^{\text{adv}}}{\partial t} - \varphi^{\text{adv}} \varepsilon_f D_e \frac{\partial^2 c^{\text{adv}}}{\partial x^2} + \varphi^{\text{adv}} \varepsilon_f v \frac{\partial c^{\text{adv}}}{\partial x} = -\alpha(c^{\text{adv}} - c^{\text{dif}}) \qquad (3)$$

$$\varphi^{\text{dif}} \varepsilon_m \left( 1 + \frac{\partial q^{\text{dif}}}{\partial c^{\text{dif}}} \right) \frac{\partial c^{\text{dif}}}{\partial t} = \alpha(c^{\text{adv}} - c^{\text{dif}}) \qquad (4)$$
Where $\phi^{adv}$ and $\phi^{dif}$ are respectively the volume fraction of the advective and the diffusive zones: $\phi^{adv} = \frac{d_f}{s}$ and $\phi^{dif} = \frac{s-d_f}{s}$, with s (m) the fracture spacing (Figure 4). $\epsilon_f$ and $\epsilon_m$ are the advective zone porosity and the matrix porosity respectively, with $\epsilon_f = \frac{a_f + 2L_{pene}\epsilon_m}{d_f}$.

$c^{adv}$, $c^{dif}$, correspond to the averaged fluid concentrations (mol/L) in the advective or the diffusive zones, while $q^{adv}$, $q^{dif}$ represent the concentrations in the solid phase (mol/L). This D. P. model was implemented in PhreeqC (Brusseau et al., 1989; Parkhurst & Appelo, 2013). In addition to the above-mentioned intrinsic averaged quantities in the advective and diffusive regions (hereafter referred to as D.P fractures and D.P. matrix quantities, respectively), we will define weighted average quantities on the whole domain (hereafter referred to as D.P. model quantities) to facilitate the comparison with S.P. model.

The D.P. model is parametrized by the volume fraction of the advective area $\phi^{adv}$ and the exchange factor $\alpha$. To ensure that the global behavior of the D.P. model is similar to the DFM model, the exchange factor is calibrated by fitting the D. P. model solution to the breakthrough curve predicted by the DFM model for a tracer transport simulation. The resulting $\alpha$ value is then applied to all the different solute species. Since $\alpha$ controls the mass transfer between fractures and matrix, it is expected that the results should be sensitive to its variations. The impact of the exchange factor on the simulation is discussed further in supporting information Text S9 in Supporting Information S2 and Figures S4 and S5 in Supporting Information S1. Note that the choice of $L_{pene}$ does not influence the results but at which spatial resolution we want to investigate the results. A large penetration length will focus on the behavior of the matrix, while a short penetration length will better describe the fractures-matrix interface. In our study, $L_{pene}$ was chosen to be in the range of the width of mineralized fractures observed on core samples.

A more detailed description of the D.P. model is discussed in supporting information Text S1 in Supporting Information S2.

### 2.2.3. Discrete Fracture Matrix Model

A DFM approach is applied to explicitly represent the fractures as lower dimensional interfaces embedded in the porous rock matrix. Since fractures are 1D interior boundaries within a 2D domain, a continuity boundary condition provides that the pressure and the mass flux are continuous, with the calculation of an out-of-plane flux as a source term in the mass and momentum conservation equations (COMSOL Multiphysics, 2022).

The mass transport equations in the matrix and fractures are written as follows:

$$ \text{Matrix: } \frac{\partial \epsilon_m q_m}{\partial t} + \frac{\partial \epsilon_m c_m}{\partial t} - \nabla \cdot (D_e \nabla c_m) + \mathbf{u}_m \cdot \nabla c_m = 0 \quad (5) $$

$$ \text{Fractures: } \frac{\partial \epsilon_f q_f}{\partial t} + \frac{\partial \epsilon_f c_f}{\partial t} - \nabla \cdot (D_e \nabla c_f) + \mathbf{u}_f \cdot \nabla c_f = \frac{J_{f,o}}{d_f} \quad (6) $$

where $c_m$ and $c_f$ correspond to the concentrations respectively in the matrix and fractures. $q_f$ and $q_m$ correspond to the concentrations in the solid phases respectively in the fractures and the matrix (mol/m³). $D_e$ corresponds to the effective diffusivity (m²/s), $\mathbf{u}_m$ and $\mathbf{u}_f$ are the fluid velocity (m/s). The mass coupling between the fractures and matrix is achieved through an out-of-plane flux $J_{f,o}$ (mol/m²/s). Such term is calculated as a function of the fluxes at the hanging wall and footwall of the fracture (COMSOL Multiphysics, 2022).

The fluid flow in both the matrix and the fractures are modeled using Darcy's law:

$$ \text{Matrix: } \left\{ \begin{array}{l} \nabla \cdot (\mathbf{u}_m) = 0 \\ \mathbf{u}_m = \frac{k_m}{\mu} \nabla p_m \end{array} \right. \quad (7) $$
<table><caption>**Table 1**

Hydrogeological Parameters</caption><thead><tr><th>Parameters</th><th>Symbol</th><th>Value</th><th>Unit</th><th>Reference</th></tr></thead><tbody><tr><td>Dynamic viscosity</td><td>μ</td><td>1e-3</td><td>Pa.s</td><td></td></tr><tr><td>Temperature</td><td>T</td><td>25</td><td>°C</td><td></td></tr><tr><td>Density</td><td>ρ<sub>w</sub></td><td>1,000</td><td>kg/m<sup>3</sup></td><td></td></tr><tr><td>Porosity: flow barrier</td><td>ε<sub>i</sub></td><td>1%</td><td>%</td><td>(Favier et al., 2022)</td></tr><tr><td>Porosity: matrix</td><td>ε<sub>m</sub></td><td>1%</td><td>%</td><td>(Favier et al., 2022)</td></tr><tr><td>Porosity: fracture</td><td>ε<sub>f</sub></td><td>3.5%</td><td>%</td><td>This study</td></tr><tr><td>Permeability: flow barrier</td><td>k<sub>i</sub></td><td>10<sup>-17</sup></td><td>m<sup>2</sup></td><td>(Jeanpert, 2017)</td></tr><tr><td>Permeability: matrix</td><td>k<sub>m</sub></td><td>2 × 10<sup>-15</sup></td><td>m<sup>2</sup></td><td>(Jeanpert, 2017)</td></tr><tr><td>Hydraulic aperture fracture</td><td>a<sub>f</sub></td><td>2.55 × 10<sup>-5</sup></td><td>m</td><td>This study</td></tr><tr><td>Permeability: fracture</td><td>k<sub>f</sub></td><td>5.4 × 10<sup>-11</sup></td><td>m<sup>2</sup></td><td>This study</td></tr><tr><td>Fracture thickness</td><td>d<sub>f</sub></td><td>1</td><td>mm</td><td>This study</td></tr><tr><td>Fracture spacing</td><td>s</td><td>30</td><td>cm</td><td>(Jeanpert, 2017)</td></tr><tr><td>Penetration length</td><td>L<sub>pene</sub></td><td>4.87 × 10<sup>-4</sup></td><td>m</td><td>This study</td></tr><tr><td>Rainwater rate</td><td>p0</td><td>4.76 × 10<sup>-8</sup></td><td>m/s</td><td>(Myagkiy et al., 2017)</td></tr><tr><td>Effective diffusivity</td><td>D<sub>e</sub></td><td>4.16 × 10<sup>-10</sup></td><td>m<sup>2</sup>/s</td><td>(Myagkiy et al., 2017)</td></tr><tr><td>PhreeqC fracture - matrix coupling term</td><td>α</td><td>1.1 × 10<sup>-6</sup></td><td>s<sup>-1</sup></td><td>This study</td></tr><tr><td>Volume fraction of the advective zone</td><td>φ<sup>adv</sup></td><td>0.33%</td><td>%</td><td>This study</td></tr><tr><td>Volume fraction of the diffusive zone</td><td>φ<sup>dif</sup></td><td>99.67%</td><td>%</td><td>This study</td></tr><tr><td>Molar weight of olivine</td><td>M<sub>0</sub></td><td>145</td><td>g/mol</td><td>This study</td></tr><tr><td>Density of peridotite</td><td>ρ<sub>p</sub></td><td>2,700</td><td>kg/m<sup>3</sup></td><td>(Kelemen et al., 2004)</td></tr><tr><td>Molar weight of sepiolite</td><td>M<sub>s</sub></td><td>647.83</td><td>g/mol</td><td>(Blanc et al., 2012)</td></tr><tr><td>Density of sepiolite</td><td>ρ<sub>s</sub></td><td>2,270</td><td>kg/m<sup>3</sup></td><td>(Blanc et al., 2012)</td></tr></tbody></table>

$$ \text{Fractures : } \left\{ \begin{array}{l} \nabla \cdot (\mathbf{u}_f) = \frac{Q_{f,o}}{d_f} \\ \mathbf{u}_f = \frac{k_f}{\mu} \nabla p_f \end{array} \right. \qquad (8) $$

where $k_m$ and $k_f$ are the permeabilities of the matrix and fractures (m$^2$), $p_m$ and $p_f$ (Pa) are the pressure respectively in the matrix and the fractures and $Q_{f,o}$ (m/s) is the out-of-plane momentum flux between fractures and matrix. The fracture permeability is defined after the cubic law as: $k_f = \frac{a_f^2}{12}$.

A more detailed description of the DFM model is discussed in supporting information Text S2 in Supporting Information S2. The DFM model was solved within the ICP (interface Comsol-PhreeqC) platform, which couples the finite element Multiphysics software COMSOL Multiphysics® v5.6 for fluid flow and transport with the geochemical code PhreeqC (Nardi et al., 2014).

## 2.3. Hydrogeological and Geochemical Settings

All the data required for the hydrogeological parametrization are taken from the literature and are detailed in supporting information Text S3 in Supporting Information S2 and Table 1.

The geochemical system studies the formation of a laterite profile resulting from the dissolution-precipitation process of a bedrock composed of olivine and unreactive serpentine by the downward progression of rain-water. Indeed, the dissolution rate of serpentine is estimated to be 1 to three order of magnitude lower depending on pH than the dissolution rate of olivine (Daval et al., 2013; Krevor & Lackner, 2011; Myagkiy et al., 2017). Serpentinization degree of the bedrock classically varies between 40 wt.% and 90 wt.% in New Caledonia
<table><caption>**Table 2**

Effective Surface Area ($A_0$), Volume of Solution per Rock Mass (V)</caption><thead><tr><th></th><th>A<sub>0</sub> (m<sup>2</sup>/kg)</th><th>V (L/kg)</th><th>a<sub>w</sub> = A<sub>0</sub>/V (m<sup>2</sup>/L)</th></tr></thead><tbody><tr><td>Matrix</td><td>2.81 × 10<sup>1</sup></td><td>3.73 × 10<sup>-3</sup></td><td>7.54 × 10<sup>3</sup></td></tr><tr><td>Fractures</td><td>2.90 × 10<sup>-1</sup></td><td>4.35 × 10<sup>-5</sup></td><td>6.67 × 10<sup>3</sup></td></tr></tbody></table>

(Iseppi, 2018). Hence, a composition of 88 wt.% of unreactive serpentine and 12 wt.% of olivine was chosen for the bedrock. Olivine composition was based on quantitative, electron probe micro analyses of unweathered olivine: Mg<sub>1.865</sub>Fe<sub>0.127</sub>Ni<sub>0.008</sub>SiO<sub>4</sub> (Teitler et al., 2019). The dissolution of olivine is assumed to be kinetically controlled during the simulation while the precipitation and dissolution of newly formed minerals are supposed to be driven by thermodynamic equilibrium (Myagkiy et al., 2017). These neo-formed minerals consist of Ni-phyllosilicates (Pimelite and Falcondoite), Mg-phyllosilicates (Kerolite, Sepiolite) and Ni-Goethite. After their precipitation, the neo-formed phyllosilicates are eventually dissolved by the incoming rainwater and only Ni-Goethite remains in the system, trapping all the Ni either in its lattice or through sorption (Myagkiy et al., 2017).

The incoming solution consists of rainwater in equilibrium with the atmospheric CO<sub>2</sub> (pH = 5.6). The chemical composition of rainwater is taken from (Myagkiy et al., 2017). The precipitation rate was assumed to be 1,500 mm/year, also in agreement with Myagkiy et al. (2017).

The specific dissolution rate r (mol/m<sup>2</sup>/s) and overall reaction rate R<sub>d</sub> (mol/L/s) for olivine were defined following Myagkiy et al. (2017) and Pokrovsky and Schott (2000). Under far from equilibrium conditions (see the calculations of olivine saturation index SI in supporting information Figure S6 in Supporting Information S1 which supports this assumption) and, for pH < 9, olivine dissolution kinetics can be approximated as follows:

$$ \log(r) = -0.5pH - 6.64 \qquad (9) $$

$$ R_d = r * \left( \frac{A_0}{V} \right) * \left( \frac{n}{n_0} \right)^{\frac{2}{3}} \qquad (10) $$

where A<sub>0</sub> is the initial effective surface area of olivine (m<sup>2</sup>/kg<sub>solid</sub>), V is the volume of solution per rock mass unit (L/kg<sub>solid</sub>) (Table 2), n<sub>0</sub> and n represent the initial moles of solid and the moles of solid at a given time, respectively. More details on the calculation of A<sub>0</sub> is given in supporting information Text S4 in Supporting Information S2.

The thermodynamic database for surface complexation of inorganic ions on goethite developed by Dzombak and Morel (1990) was used in simulations to allow the sorption of both Ni and Mg (Myagkiy et al., 2017). The llnl.dat thermodynamic database (Johnson et al., 2000) edited for Ni-Mg phyllosilicates minerals (Myagkiy et al., 2017) was used to conduct the geochemical calculations at 25°C (Table 3).

## 2.4. Numerical Validation

The correctness of the different S.P., D.P., and DFM models was assessed as detailed in supporting information Text S5 in Supporting Information S2, and Figure S1 in Supporting Information S1.

# 3. Results and Discussion

## 3.1. Weathering Rates

The formation of a lateritic profile is characterized by the propagation of two different weathering fronts: (a) the primary silicates dissolution front in the bedrock and (b) the neo-formed silicates dissolution front in the saprolite horizon (Cathelineau et al., 2016; Favier et al., 2022). In our study, the primary silicates dissolution front is depicted by the olivine dissolution (Figures 5a, 5b, 5e, and 5f), while the neo-formed silicates dissolution front is represented by the Mg-phyllosilicates dissolution (Figures 5i and 5j). A total weathering rate of bedrock into oxide horizon (goethite) is estimated at 20 m/My from both S. P. and D.P. models, based on the final dissolution time of phyllosilicates.

<table><caption>**Table 3**

Equilibrium Phases</caption><thead><tr><th>Mineral</th><th>Formula</th><th>Log K</th><th>References</th></tr></thead><tbody><tr><td>Kerolite</td><td>Mg<sub>3</sub>Si<sub>4</sub>O<sub>10</sub>(OH)<sub>2</sub>:H<sub>2</sub>O</td><td>25.79</td><td>(Stoessell, 1988)</td></tr><tr><td>Pimelite</td><td>Ni<sub>3</sub>Si<sub>4</sub>O<sub>10</sub>(OH)<sub>2</sub>:H<sub>2</sub>O</td><td>11.46</td><td>(Myagkiy et al., 2019)</td></tr><tr><td>Sepiolite</td><td>Mg<sub>4</sub>Si<sub>6</sub>O<sub>15</sub>(OH)<sub>2</sub>:6H<sub>2</sub>O</td><td>30.44</td><td>(Stoessell, 1988)</td></tr><tr><td>Falcondoite</td><td>Ni<sub>4</sub>Si<sub>6</sub>O<sub>15</sub>(OH)<sub>2</sub>:6H<sub>2</sub>O</td><td>12.31</td><td>(Myagkiy et al., 2019)</td></tr><tr><td>Goethite_Ni</td><td>Fe<sub>0.98</sub>Ni<sub>0.03</sub>OOH</td><td>-3.99</td><td>(Myagkiy et al., 2019)</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Dissolution rates of olivine (in black) at 630 years in panels (a) Single Porosity (S.P.) Model (S.P. Model), (b) Dual Porosity (D.P.) Model (D.P. Model), (c) Porous matrix part of the D.P. Model (D.P. Matrix), (d) fractures in the D.P. Model (D.P. Fractures) and at 6,300 years in panels (e) S.P. Model, (f) D.P. Model, (g) D.P. Matrix, (h) D.P. Fractures, and dissolution rates of Mg-phyllosilicates (in blue) at 29,000 years in panel (i) S.P. Model, (j) D.P. Model, (k) D.P. Matrix, (l) D.P. Fractures.</figcaption></figure>

However, in New Caledonia, the reported weathering rates are assessed by considering the oxide horizon formation rather than the bedrock dissolution. They range from 10 m/My to 50 m/My under steady-state conditions and without erosion (Freyssinet & Farah, 2000; Freyssinet et al., 2005; Maurizot et al., 2020; Nahon & Tardy, 1992; Trescases, 1975). To convert our weathering rates into an oxide horizon formation rate, the compaction of the weathering profile must be accounted for (Maurizot et al., 2020), as well as the serpentinization degree of the peridotite impacting the initial quantity of olivine. Depending on the weathering heterogeneities, a compaction factor between 1 and 6 must be applied (Fritsch et al., 2014). By assuming a linear dependency of weathering time to the serpentinization degree, the global weathering rate is estimated in a range of 0.6 m/My (serpentinization: 40%; compaction factor: 6) to 20.5 m/My (serpentinization: 90%; compaction factor: 1). Given the different approximations made in the modeling (notably, not considering the porosity changes with weathering), such result is close to previous estimates.

Furthermore, the weathering is not uniform along the column, as already mentioned in the literature (Reis & Brantley, 2017). The olivine dissolution rate was calculated as a function of depth and at different times (supporting information Text S6 in Supporting Information S2), to provide a better insight into the weathering dynamics as illustrated in Figure 5. In the case of the D.P. model, we differentiate the total dissolution rate (calculated over the total column) from the dissolution rate per region (denoted as D.P. matrix or D.P. fracture). The dissolution rate follows an exponential decay trend at early stages (630 years, Figure 5) in agreement with Reis and Brantley (2017). The weathering is more localized at the top of the column for the S.P. model with higher weathering intensities than in the D.P. model.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Mineralogical profiles at 83,000 years of (a) the Single Porosity (S.P.) Model (S.P. Model), (b) the Dual Porosity (D.P.) Model (D.P. Model), (c) the matrix (D. P. Matrix), and the (d) fractures (D.P. Fractures) of the D.P. Model.</figcaption></figure>

In addition, the D.P. model displays a more intense weathering rate in the fractures (D.P. Fractures), since they are advection-dominated, than in the porous matrix (D.P. Matrix) under diffusion control (Figures 5c and 5d). Accounting for 99.67% of the column volume, the matrix behavior drives the general dynamics of the whole profile (supporting information S1). At later stages (6,300 years), the olivine dissolution rate becomes mostly uniform for both S.P. and D.P. models.

Once the olivine has been dissolved, the dissolution of the Mg-phyllosilicates initiates, showing two very different dynamics between the S.P. and the D.P. models. An exponential decay without any constant term, typical in weathering systems governed by local equilibrium (Ma et al., 2010; Reis & Brantley, 2017; Vanwalleghem et al., 2013), is observed for the D.P. model. On the contrary, the S.P. model only shows an extremely localized dissolution of the phyllosilicates, materialized by a peak in the dissolution rate (Figures 5i and 5j). The fractures behave more like the unfractured column predicted by the S.P. model but with a faster dynamic (Figure 5l).

The difference in dissolution rate profiles of olivine and Mg-phyllosilicates is explained by the kinetics driving olivine dissolution along the whole profile, while the exponential decay is caused by the equilibrium with the pore-water solution, which entirely controls the Mg-phyllosilicates dissolution.

Thus, depending on the competition between mass transfer and reaction kinetics, two limiting behaviors occur for the dissolution or precipitation rate of a given mineral species.

## 3.2. Geochemical Evolution

The comparison of global mineralogical behavior between the S.P. and D.P. models highlights similarities in the overall mineralogical composition of a lateritic profile. Three different mineralogical areas are thus observed in Figure 6, from bottom to top: (1) a lower zone dominated by Mg-phyllosilicates, where the pH is buffered at ca. 9, (2) a transition zone characterized by the significant occurrence of Ni-phyllosilicates, while the pH decreases from ca. 9 to 5.5, and (3) an upper zone, entirely composed of goethite, where the pH is buffered at 5.5 (Figures 6a and 6b). Such mineralogical zoning from a phyllosilicate-dominated horizon to an oxide-dominated horizon is
<figure><img src="image_9.png" /><figcaption>**Figure 7.** (a) Chemical path in Mg-SiO₂ activity diagram in fractures and matrix. (b) Chemical path in Ni-SiO₂ activity diagram in fractures and matrix. In colors, the main mineralogical species. Black: Olivine; Blue: Mg-phyllosilicates; Green: Ni-phyllosilicates; Red: Goethite.</figcaption></figure>

consistent with field observations and previous numerical modeling studies (Cathelineau et al., 2016; Domènech et al., 2017; Favier et al., 2022; Myagkiy et al., 2017).

The mineralogical sequence, that is, evolving from initial olivine (forsterite) into Mg-Ni-phyllosilicates—either the kerolite—pimelite or sepiolite—falcondoite series, and then into goethite (Butt & Cluzel, 2013; Cathelineau et al., 2016, 2017; Dublet et al., 2012) is depicted by the chemical path of the geochemical system in the $Ni^{2+}/(H^{+})^2$-SiO₂ and $Mg^{2+}/(H^{+})^2$-SiO₂ activity diagrams built with Geochemist's Workbench (Bethke, 2007) (Figure 7 a and b). While constrained by the model assumptions and the choice of phases authorized to precipitate, these activity diagrams give insights on why the Ni-phyllosilicates persist longer in the profile than Mg-phyllosilicates. The boundary between quartz and Mg-phyllosilicates is situated at a higher $Mg^{2+}/(H^{+})^2$ ratio than the corresponding boundary between quartz and Ni-phyllosilicates (Figures 7b and 7c). Apart from the olivine dissolution during which ($Mg^{2+}$) highly decreases, ($Mg^{2+}$) remains mainly sub-constant during the simulation because of incoming rainwater rich in ($Mg^{2+}$). The decrease in $Mg^{2+}/(H^{+})^2$ ratio is, therefore, mostly explained by the increase of $H^+$, that is, the pH lowering of the solution. Hence for the same pH variation, the solution becomes undersaturated with respect to the Mg-phyllosilicate while it remains saturated with respect to the Ni-phyllosilicates. The smaller domain of stability of Mg-phyllosilicates explains their earlier dissolution, while Ni-phyllosilicates remain saturated when the pH front progresses downwards.

If the mineralogical sequence is identical in both S.P. and D.P. models, spatial distribution differs depending on the model (Figures 6a and 6b). The transition zone (2) is thicker (approximately 50 cm-thick) in the D. P. model in relation with a gradual decrease of the pH, from ca. 9 to 5.5 in ca. 50 cm, while this transition is sharper (10 cm) for the S. P. model. Furthermore, Ni-phyllosilicates, which form within the transition zone, are locally more concentrated in the S.P. model since the transition zone is thinner.

These similarities between S.P. and D.P. models are also found in the bulk rock content profiles, where the behaviors of Mg, Fe, and Si are substantially similar (Figures 8b–8d). The migration of weathering fronts in both models is visible. A careful analysis of Figure 8 confirms the three distinct geochemical sections previously identified. A smoother transition between the oxide horizon (zone 3) and the saprolite horizon in the D.P. model (zone 2) is again observed. In addition, the Ni enrichment is localized in a 10 cm depth window for the S.P. Model, while it extends over ca. 50 cm for the D.P. model (Figure 8a). Changes in the Mg-phyllosilicates dissolution rate patterns along depth (Figures 6i and 6j) explain the discrepancy observed between a sharp transition in the S. P. model and a smooth transition in the D. P. model (Figures 7a and 7b). Indeed, the concentration of fluids along the fractures leads to the mitigation of the fluid flow and the spreading of weathering along depth in the matrix.

These differences are discussed below by analyzing the fractures and matrix behaviors predicted by the D.P. model.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Bulk rock content profiles at 83,000 years for the Single Porosity (S.P.) Model (S.P. Model), the Dual Porosity (D.P.) Model (D.P. Model), the matrix of the D. P. Model (D.P. matrix), and the fractures of the D.P. Model (D.P. fractures) of (a) Ni, (b) Mg, (c) Fe, (d) Si. Apart from Ni, the D.P. Model and the D.P. matrix curves are confounded for the other elements.</figcaption></figure>

### 3.3. Impact of Fractures on the Dual Porosity Model

The mineralogical and chemical profiles within the matrix of the D.P. model are similar to the S.P. model, although slightly delayed (Figures 6a, 6c, and 8). However, the changes induced by weathering within the fractures of the D.P. model are clearly different. The transition zone (2) is much wider than in the matrix (ca. 3 m) and much more enriched in Ni-phyllosilicates with content ranging from 40 wt.% to ca. 100 wt.% (Figure 6d). Furthermore, fractures are depleted in Mg-phyllosilicates in the transition zone. Thus, the fractured bedrock acts as filtering mechanism for the newly formed phyllosilicates. Ni-phyllosilicates preferentially precipitate in fractures, whereas Mg-phyllosilicates form in the matrix. Such mineralogical distribution leads to fractures characterized by high Ni content ranging from more than 10 wt.% to 30 wt.%, and almost no Mg (<1 wt.%), while the matrix is characterized by Ni content ranging from 1 wt.% to 6 wt.%. These results are consistent with field observations where fractures present high Ni concentrations ranging from 10 wt.% to 40 wt.% (Cathelineau et al., 2016, 2017), while Ni grades do not exceed 2 to 4 wt.% in porous matrix (Butt & Cluzel, 2013; Cathelineau et al., 2016; Teitler et al., 2019, 2022).

The study of the mineralogical evolution through time at a specific depth in the middle of the profile (Figures 9a and 9b) provides further insights into the mineral partitioning between fractures and matrix. Significant enrichment in Ni-phyllosilicates occurs only when the Mg-phyllosilicates are mostly dissolved (Figure 6). Otherwise, the presence of large quantities of Mg-phyllosilicates makes the absolute variations of Ni-phyllosilicates negligible. In fractures, the Mg-phyllosilicates are dissolved much faster (at ca. 20,000 years) than in the matrix, which leads to an earlier enrichment in Ni-phyllosilicates, and thus a longer transition phase (2) of about 90,000 years (Figure 9 a). In the matrix, the Mg-phyllosilicates are more slowly dissolved. Hence the enrichment in Ni-phyllosilicates is delayed until ca. 110,000 years (Figure 9 b). Therefore, a significant Ni-phyllosilicate enrichment is only observed for a short period of ca. 20,000 years, when the dissolution of Mg-phyllosilicates increases, followed by a significant precipitation of Ni-phyllosilicates. The detailed impact of
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Molar composition of (a) fractures and (b) matrix at *z* = 2 m in the Dual Porosity model. In blue: Mg-phyllosilicates, in green: Ni-phyllosilicates, and in red: goethite. The sudden end of data before the simulation end for Mg-phyllosilicates and Ni-phyllosilicates is explained by their total dissolution from the domain at that time. (c) Dependency of the final dissolution time of the Mg-phyllosilicate in the matrix and the fractures with respect to the saturation index of the dissolving solution.</figcaption></figure>

the molar evolution of the Mg-Ni-phyllosilicates on the relative enrichment in Ni-phyllosilicates is discussed in
supporting information Text S7 in Supporting Information S2 and Figure S2 in Supporting Information S1.

An in-depth analysis of Mg-phyllosilicates dissolution time reveals why they dissolve faster in fractures, therefore leading to higher residual enrichment of Ni-phyllosilicates (Figure 9c). Indeed, the final dissolution time of Mg-phyllosilicates as a function of the negative saturation index (undersaturated solution with respect to the Mg-phyllosilicates) calculated from the lixiviation of a Mg-phyllosilicate cell by varying undersaturated solutions, is larger in the matrix, whatever the saturation index. Such a difference comes from a lower liquid-solid ratio and a larger residence time in the matrix which reduces weathering.

Then, the mass transfer occurs between matrix and fractures and enhances the fractures enrichment in Ni-
phyllosilicates, explaining the differences observed in terms of Ni distribution (see supporting information
Text S8 in Supporting Information S2, Figure S3 in Supporting Information S1 for a detailed discussion of this
mass transfer).

**3.4. Mineralization Hotspots Revealed by DFM Model**

To reproduce typical bedrock structuration, which commonly presents two sub-orthogonal fracture sets with a mean fracture spacing of about 30 cm (Jeanpert et al., 2019), a 4.5 m × 4.5 m DFM reference case model is designed with two orthogonal fracture sets dipping respectively 20° and 110°. The spacing of the sub-horizontal fractures is ca. 40 cm, whereas it is ca. 70 cm for the vertical ones (Figure 10a).

The impact of fractures on the weathering process was studied with two other geometries or mineral infills: (a) the dihedron (Figure 10b) and (b) the serpentine bench (Figure 10c) where low-permeability layers, due to serpentine coating, act as natural flow barriers (Iseppi et al., 2018).

The mineralogical distribution of Mg-Ni-phyllosilicates and goethite clearly defines two geochemical domains separated by a transition area delineated by the fracture network where the Ni-phyllosilicates are concentrated (Figure 11). Such organization is consistent with the D.P. Model (Figures 6 and 8) and field observations (Favier et al., 2022).

Furthermore, boulders with a concentric weathering intensity are found in the transition zone and the lower part of the oxide zone. They are characterized by a rim enriched in Ni-phyllosilicates and a core mainly composed of Mg-phyllosilicates. They result from the diffusion of the weathering solution from the fractures (Myagkiy et al., 2019). Such an organization is extremely similar to that observed in outcrops (Figure 3b), where boulders
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Discrete Fracture Matrix (DFM) conceptual models and their pressure and fluid flow distribution. (a) DFM reference case, (b) Dihedron, (c) Serpentine bench.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Geochemical and mineralogical distribution in Discrete Fracture Matrix model at 127,000 years of (a) Ni and Ni-phyllosilicates, (b) Fe and Goethite, (c) Mg and Mg-phyllosilicates, (d) Si, and (e) pH.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Ni distribution and outcrop organization in New Caledonia similar to Discrete Fracture Matrix (DFM) models at 95,000 years for (a) the DFM reference model, (b) the dihedron, and (c) the flow barrier.</figcaption></figure>

are characterized by a concentric Ni enrichment mainly localized at their rim, while their core remains poorly enriched.

The addition of low permeability fractures materialized by (a) a dihedron and (b) a sub-horizontal layer, highlights the impact of such heterogeneities on the flux distribution and, therefore, on the weathering processes. By concentrating the flux to the lowest point, the dihedron favors a more intense weathering of the matrix in the vicinity of the fractures situated below (zone (a)) and conversely favors the formation of poorly weathered areas away from the outlet location (zone (b)) (Figure 12b). Such a simulation well reproduces the typical regolith geometry in the presence of a dihedron, where a weathered area (zone 3) is visible between two unweathered planes delimiting a poorly weathered region (zone 1) (Figure 12b). Likewise, the sub-horizontal flow barrier drives lateral weathering by inhibiting the role of the sub-vertical fractures (Figure 12c). The weathering only occurs at the hanging wall of the flow barrier (zone 3), while the footwall remains mostly unweathered (zone 1), which is again consistent with field observations (Figure 12c).

4. Conclusions

Permeability heterogeneities control weathering processes (Jung & Navarre-Sitchler, 2018), which either favor or hinder them by creating preferential flow paths or barriers (Goddéris et al., 2019; Navarre-Sitchler et al., 2015). The multiple porosity formulation of the bedrock weathering considers the impact of such heterogeneities by connecting two or more different hydro-chemical areas. In our study, the fractures are depicted as preferential flow paths concentrating the fluid flow, while the matrix is connected to the fractures through diffusivity.

The impact of fracture networks as chemical and hydrological heterogeneities on the weathering front progression in a fractured rock mass is investigated through different models of increasing complexity.

Comparing the global behavior of the D.P. model with the S.P. one sheds light on the role of fractures on weathering rates and mineralogical distribution along a vertical profile. By extending the weathering over a greater depth, the fractures favor (a) a smoother transition between the saprolite and the oxide horizons and (b) a
less localized Ni enrichment in the matrix compared to the S. P. Model. Fractures also participate in forming Ni heterogeneities observed in the field as shown in the D.P. model. Because of the differences in residence time and liquid-to-solid ratio between the advective and diffusive area of the fractured bedrock, Ni-phyllosilicates preferentially precipitate in fractures while Mg-phyllosilicates remain trapped in the matrix. Hence, the D.P. formulation reveals the role of fractures in the distribution of weathering heterogeneities.

Yet, such formulation cannot precisely model complex geometries and hydrological systems, such as multiple permeability heterogeneities and complex fracture networks. In contrast, the DFM approach more accurately addresses the complexity of weathering processes through a better description of the fracture connectivity effect on the dissolution-precipitation processes. Depending on its geometry and permeability, the fracture network controls the fluid flow, hence favoring the formation of geometrical weathering heterogeneities. The heterogeneous distribution of alteration rates and metal enrichment in the regolith as a function of the initial distribution of fractures explains the specific localization of Ni hotspots of interest for mining companies.

The D.P. and DFM models offer excellent opportunities to model complex systems displaying hydrogeochemical heterogeneities and, therefore, better model and explain the weathering processes. Yet they need to be carefully designed to describe the weathering systems correctly. Thus, the mass transfer coefficient, intrinsically related to the fracture network connectivity, is tedious to adjust. This difficulty may significantly impact the fluxes between advective and diffusive regions and hence the weathering of the matrix. A dense and highly connected fracture network, corresponding to a high mass transfer coefficient value $\alpha$, increases the global weathering by favoring the transfer of an unequilibrated solution to the matrix. In contrast, a low connectivity of fractures (i.e., a low $\alpha$) decreases the weathering by reducing the exchanges between the weathering solution and the porous matrix (see supporting information Text S9 in Supporting Information S2 and Figures S4 and S5 in Supporting Information S1 for a detailed discussion on the mass transfer coefficient).

Furthermore, if ICP theoretically allows a full coupling of the geochemical-hydrological systems, the current modeling is only partially coupled. Indeed, the porosity and permeability, which should be dependent on the rock's mineralogical nature, are set constant during the whole simulation. Such a limitation reduces the rock-water interactions by preventing the increase of porosity and permeability linked to the bedrock dissolution. Therefore, a fully coupled model should create higher weathering rates, perhaps closer to the reported ones.

## Data Availability Statement

Data were interpreted using the open-source software PhreeqC version 3 (Parkhurst et al., 2013b), and the commercial softwares COMSOL MULTIPHYSICS® version 5.6 (CONSOL Multiphysics®, 2020) and ICP version 2.1 (Nardi et al., 2022). Data sets and codes related to PhreeqC software (i.e., Single and Dual Porosity models) can be found at https://doi.org/10.24396/ORDAR-137 (Favier et al., 2024) an open-source online data repository hosted at ORDAR (OTELo Research Data Repository).

## References

Babey, T., de Dreuzy, J.-R., & Ginn, T. R. (2016). From conservative to reactive transport under diffusion-controlled conditions: Nonlinear fluid-rock reactivity from conservative transport. *Water Resources Research*, *52*(5), 3685-3700. https://doi.org/10.1002/2015WR018294

Beerling, D. J., Kantzas, E. P., Lomas, M. R., Wade, P., Eufrasio, R. M., Renforth, P., et al. (2020). Potential for large-scale CO2 removal via enhanced rock weathering with croplands. *Nature*, *583*(7815), 242-248. https://doi.org/10.1038/s41586-020-2448-9

Beinlich, A., John, T., Vrijmoed, J. C., Tominaga, M., Magna, T., & Podladchikov, Y. Y. (2020). Instantaneous rock transformations in the deep crust driven by reactive fluid flow. *Nature Geoscience*, *13*(4), 307-311. https://doi.org/10.1038/s41561-020-0554-9

Berre, I., Doster, F., & Keilegavlen, E. (2019). Flow in fractured porous media: A review of conceptual models and discretization approaches. *Transport in Porous Media*, *130*(1), 215-236. https://doi.org/10.1007/s11242-018-1171-6

Bethke, C. M. (2007). *Geochemical and biogeochemical reaction modeling* (2nd ed.). Cambridge University Press. https://doi.org/10.1017/CBO9780511619670

Blanc, P., Lassin, A., Piantone, P., Azaroual, M., Jacquemet, N., Fabbri, A., & Gaucher, E. C. (2012). Thermoderm: A geochemical database focused on low temperature water/rock interactions and waste materials. *Applied Geochemistry*, *27*(10), 2107-2116. https://doi.org/10.1016/j.apgeochem.2012.06.002

Borg, G. (2009). The influence of fault structures on the genesis of supergene zinc deposits. In S. R. Titley (Ed.), *Supergene environments, processes, and products*. Society of Economic Geologists. https://doi.org/10.5382/SP.14.11

Brantley, S. L., Lebedeva, M. I., Balashov, V. N., Singha, K., Sullivan, P. L., & Stinchcomb, G. (2017). Toward a conceptual model relating chemical reaction fronts to water flow paths in hills. *Geomorphology*, *277*, 100-117. https://doi.org/10.1016/j geomorph.2016.09.027

Brusseau, M. L., Jessup, R. E., & Rao, P. S. C. (1989). Modeling the transport of solutes influenced by multiprocess nonequilibrium. *Water Resources Research*, *25*(9), 1971-1988. https://doi.org/10.1029/WR025i009p01971

## Acknowledgments

This work has been funded and logistically supported by the CNRT research contract CSF N° 9PS2017-CNRT.
GEORESSOURCE/TRANSNUM«
Facteurs d'enrichissements et transferts de Ni, Co-Sc dans les saprolites de Nouvelle Calédonie: approche géométrique, minéralo-géochimique et numérique », and by the French National Research Agency through the national program "Investissements d'avenir" of the Labex Ressources 21 with the reference ANR-10-LABX-21-01/LABEX RESSOURCES21.
Butt, C. R. M., & Cluzel, D. (2013). Nickel laterite ore deposits: Weathered serpentinites. *Elements*, 9(2), 123–128. https://doi.org/10.2113/gselements.9.2.123

Cathelineau, M., Boiron, M.-C., Grimaud, J.-L., Favier, S., Teitler, Y., & Golfier, F. (2023). Pseudo-karst silicification related to late Ni reworking in New Caledonia. *Minerals*, 13(4), 518. https://doi.org/10.3390/min13040518

Cathelineau, M., Myagkiy, A., Quesnel, B., Boiron, M.-C., Gautier, P., Boulvais, P., et al. (2017). Multistage crack seal vein and hydrothermal Ni enrichment in serpentinized ultramafic rocks (Koniambo massif, New Caledonia). *Mineralium Deposita*, 52(7), 945–960. https://doi.org/10.1007/s00126-016-0695-3

Cathelineau, M., Quesnel, B., Gautier, P., Boulvais, P., Couteau, C., & Drouillet, M. (2016). Nickel dispersion and enrichment at the bottom of the regolith: Formation of pimelite target-like ores in rock block joints (Koniambo Ni deposit, New Caledonia). *Mineralium Deposita*, 51(2), 271–282. https://doi.org/10.1007/s00126-015-0607-y

Chigira, M., & Hirata, Y. (2021). The Kui boulder fields: Formation processes controlled by columnar joints of granodiorite. *Catena*, 207, 105683. https://doi.org/10.1016/j.catena.2021.105683

COMSOL Multiphysics. (2022). Battery design module user's guide. COMSOL Multiphysics® v 6.1. COMSOL AB (pp. 377–378).

COMSOL Multiphysics® (5.6). (2020). [Computer software]. COMSOL AB. www.comsol.com

Daval, D., Hellmann, R., Martinez, I., Gangloff, S., & Guyot, F. (2013). Lizardite serpentine dissolution kinetics as a function of pH and temperature, including effects of elevated pCO2. *Chemical Geology*, 351, 245–256. https://doi.org/10.1016/j.chemgeo.2013.05.020

Davit, Y., Golfier, F., Latché, J.-C., & Quintard, M. (2019). A domain decomposition approach to finite-epsilon homogenization of scalar transport in porous media. *SIAM Journal on Applied Mathematics*, 79(5), 1797–1822. https://doi.org/10.1137/17M1157775

de Dreuzy, J.-R., Rapaport, A., Babey, T., & Harmand, J. (2013). Influence of porosity structures on mixing-induced reactivity at chemical equilibrium in mobile/immobile Multi-Rate Mass Transfer (MRMT) and Multiple IInteracting Continua (MINC) models: Influence of immobile porosity structure on reactivity. *Water Resources Research*, 49(12), 8511–8530. https://doi.org/10.1002/2013WR013808

Deng, H., & Spycher, N. (2019). Modeling reactive transport processes in fractures. *Reviews in Mineralogy and Geochemistry*, 85(1), 49–74. https://doi.org/10.2138/rmg.2019.85.3

Dill, H. G., Gerdes, A., & Weber, B. (2010). Age and mineralogy of supergene uranium minerals—Tools to unravel geomorphological and palaeohydrological processes in granitic terrains (Bohemian Massif, SE Germany). *Geomorphology*, 117(1–2), 44–65. https://doi.org/10.1016/j.geomorph.2009.11.005

Domènech, C., Galí, S., Villanova-de-Benavent, C., Soler, J. M., & Proenza, J. A. (2017). Reactive transport model of the formation of oxide-type Ni-laterite profiles (Punta Gorda, Moa Bay, Cuba). *Mineralium Deposita*, 52(7), 993–1010. https://doi.org/10.1007/s00126-017-0713-0

Dublet, G., Juillot, F., Morin, G., Fritsch, E., Fandeur, D., Ona-Gnuema, G., & Brown, G. E. (2012). Ni speciation in a New Caledonian lateritic regolith: A quantitative X-ray absorption spectroscopy investigation. *Geochimica et Cosmochimica Acta*, 95, 119–133. https://doi.org/10.1016/j.gca.2012.07.030

Dzombak, D. A., & Morel, F. M. M. (1990). Surface complexation modeling: Hydrous ferric oxide.

Favier, S., Teitler, Y., Golfier, F., & Cathelineau, M. (2022). Multiscale physical–chemical analysis of the impact of fracture networks on weathering: Application to nickel redistribution in the formation of Ni-laterite ores, New Caledonia. *Ore Geology Reviews*, 147, 104971. https://doi.org/10.1016/j.oregeorev.2022.104971

Favier, S., Teitler, Y., Golfier, F., & Cathelineau, M. (2024). Data for the paper 'Reactive transport modelling in porous fractured media: Application to weathering processes' (Version 1) [dataset]. OTELO. https://doi.org/10.24396/ORDAR-137

Freyssinet, P., Butt, C. R. M., Morris, R. C., & Piantone, P. (2005). Ore-forming processes related to lateritic weathering. In J. W. Hedenquist, J. F. H. Thompson, R. J. Goldfarb, & J. P. Richards (Eds.), *One hundredth anniversary volume*. Society of Economic Geologists. https://doi.org/10.5382/AV100.21

Freyssinet, P., & Farah, A. S. (2000). Geochemical mass balance and weathering rates of ultramafic schists in Amazonia. *Chemical Geology*, 170(1–4), 133–151. https://doi.org/10.1016/S0009-2541(99)00245-4

Frings, P. J., & Buss, H. L. (2019). The central role of weathering in the geosciences. *Elements*, 15(4), 229–234. https://doi.org/10.2138/gselements.15.4.229

Fritsch, E., Juillot, F., Dublet, G., Fandeur, N., Fonteneau, L., Martin, É., et al. (2014). Analyse fine des minerais latéritiques: Approches pétrographique, minéralogique, géochimique et isotopique. *Les gisements de nickel latéritique de Nouvelle-Calédonie, III*, 156. [PDF]. https://doi.org/10.34928/CYJB-4085

Goddéris, Y., Schott, J., & Brantley, S. L. (2019). Reactive transport models of weathering. *Elements*, 15(2), 103–106. https://doi.org/10.2138/gselements.15.2.103

Golfier, F., Quintard, M., Cherblanc, F., Zinn, B. A., & Wood, B. D. (2007). Comparison of theory and experiment for solute transport in highly heterogeneous porous medium. *Advances in Water Resources*, 30(11), 2235–2261. https://doi.org/10.1016/j.advwatres.2007.05.004

Golightly, J. P. (2010). Progress in understanding the evolution of nickel laterites. In R. J. Goldfarb, E. E. Marsh, & T. Monecke (Eds.), *The challenge of finding new mineral resources: Global metallogeny, innovative exploration, and new discoveries*. Society of Economic Geologists. https://doi.org/10.5382/SP.15.2.07

González-Álvarez, I., Boni, M., & Anand, R. R. (2016). Mineral exploration in regolith-dominated terrains: Global considerations and challenges. *Ore Geology Reviews*, 73, 375–379. https://doi.org/10.1016/j.oregeorev.2015.11.017

Goudie, A. S., & Viles, H. A. (2012). Weathering and the global carbon cycle: Geomorphological perspectives. *Earth-Science Reviews*, 113(1–2), 59–71. https://doi.org/10.1016/j.earscirev.2012.03.005

Haggerty, R., & Gorelick, S. M. (1995). Multiple-rate mass transfer for modeling diffusion and surface reactions in media with pore-scale heterogeneity. *Water Resources Research*, 31(10), 2383–2400. https://doi.org/10.1029/95WR10583

Haldar, S. K. (2017). Deposits of Europe. In *Platinum-nickel-chromium deposits* (pp. 217–231). Elsevier. https://doi.org/10.1016/B978-0-12-802041-8.00008-0

Heidari, P., Li, L., Jin, L., Williams, J. Z., & Brantley, S. L. (2017). A reactive transport model for Marcellus shale weathering. *Geochimica et Cosmochimica Acta*, 217, 421–440. https://doi.org/10.1016/j.gca.2017.08.011

Hyman, J. D., Karra, S., Makedonska, N., Gable, C. W., Painter, S. L., & Viswanathan, H. S. (2015). dfnWorks: A discrete fracture network framework for modeling subsurface flow and transport. *Computers & Geosciences*, 84, 10–19. https://doi.org/10.1016/j.cageo.2015.08.001

Iseppi, M. (2018). Fracturation polyphasée et contrôles des gisements de nickel supergène de Nouvelle-Calédonie: Nouvelles méthodes d'exploration et modèles de gisements [PhD Thesis]. Retrieved from http://www.theses.fr/2018NCAI0001/document

Iseppi, M., Sevin, B., Cluzel, D., Maurizot, P., & Bayon, B. L. (2018). Supergene nickel ore deposits controlled by gravity-driven faulting and slope failure, peridotite nappe, New Caledonia. *Economic Geology*, 113(2), 531–544. https://doi.org/10.5382/econgeo.2018.4561

Jeanpert, J. (2017). Structural and hydrogeologic characterization of the peridotitic massifs of New Caledonia (Issue 2017LARE0047) [Thesis, Université de la Réunion]. Retrieved from https://tel.archives-ouvertes.fr/tel-01845503
Jeanpert, J., Iseppi, M., Adler, P. M., Genthon, P., Sevin, B., Thovert, J.-F., et al. (2019). Fracture controlled permeability of ultramafic basement aquifers. Inferences from the Koniambo massif, New Caledonia. *Engineering Geology*, **256**, 67–83. https://doi.org/10.1016/j.enggeo.2019.05.006

Johnson, J., Anderson, G., & Parkhurst, D. (2000). Database from thermo.com.v8.r6.230 prepared at Lawrence Livermore National Laboratory (Revision: 1.11).

Jung, H., & Navarre-Sitchler, A. (2018). Physical heterogeneity control on effective mineral dissolution rates. *Geochimica et Cosmochimica Acta*, **227**, 246–263. https://doi.org/10.1016/j.gca.2018.02.028

Kelemen, P. B., Kikawa, E., Miller, D. J., Abe, N., Bach, W., Carlson, R. L., et al. (Eds.). (2004). 209 Initial Reports, *Proceedings of the ocean drilling program* (Vol. 209). Ocean Drilling Program. https://doi.org/10.2973/odp.proc.ir.209.2004

Krevor, S. C. M., & Lackner, K. S. (2011). Enhancing serpentine dissolution kinetics for mineral carbon dioxide sequestration. *International Journal of Greenhouse Gas Control*, **5**(4), 1073–1080. https://doi.org/10.1016/j.igggc.2011.01.006

Laurent, G., Izart, C., Lechenard, B., Golfier, F., Marion, P., Collon, P., et al. (2019). Numerical modelling of column experiments to investigate in-situ bioleaching as an alternative mining technology. *Hydrometallurgy*, **188**, 272–290. https://doi.org/10.1016/j.hydromet.2019.07.002

Lebedeva, M. I., & Brantley, S. L. (2017). Weathering and erosion of fractured bedrock systems: Weathering and erosion of fractured bedrock systems. *Earth Surface Processes and Landforms*, **42**(13), 2090–2108. https://doi.org/10.1002/esp.4177

Ma, L., Chabaux, F., Pelt, E., Blaes, E., Jin, L., & Brantley, S. (2010). Regolith production rates calculated with uranium-series isotopes at Susquehanna/Shale Hills Critical Zone Observatory. *Earth and Planetary Science Letters*, **297**(1–2), 211–225. https://doi.org/10.1016/j.epsl.2010.06.022

Maher, K., & Navarre-Sitchler, A. (2019). 12. Reactive transport processes that drive chemical weathering: From making space for water to dismantling continents. In J. Druhan & C. Tournassat (Eds.), *Reactive transport in natural and engineered systems* (pp. 349–380). De Gruyter. https://doi.org/10.1515/9781501512001-013

Maurizot, P., Sevin, B., Lesimple, S., Bailly, L., Iseppi, M., & Robineau, B. (2020). Chapter 10: Mineral resources and prospectivity of the ultramafic rocks of New Caledonia (Vol. 51, pp. 247–277). Geological Society. https://doi.org/10.1144/M51-2016-17

Md Dan, M. F., Muhamad, E. T., Komoo, I., & Alel, M. N. A. (2015). Physical characteristics of boulders formed in the tropically weathered granite. *Jurnal Teknologi*, **72**(3). https://doi.org/10.11113/jt.v72.4015

Molson, J., Aubertin, M., & Bussière, B. (2012). Reactive transport modelling of acid mine drainage within discretely fractured porous media: Plume evolution from a surface source zone. *Environmental Modelling & Software*, **38**, 259–270. https://doi.org/10.1016/j.envsoft.2012.06.010

Myagkiy, A., Golfier, F., Truche, L., & Cathelineau, M. (2019). Reactive transport modeling applied to Ni laterite ore deposits in New Caledonia: Role of hydrodynamic factors and geological structures in Ni mineralization. *Geochemistry, Geophysics, Geosystems*, **20**(3), 1425–1440. https://doi.org/10.1029/2018GC007606

Myagkiy, A., Truche, L., Cathelineau, M., & Golfier, F. (2017). Revealing the conditions of Ni mineralization in the laterite profiles of New Caledonia: Insights from reactive geochemical transport modelling. *Chemical Geology*, **466**, 274–284. https://doi.org/10.1016/j.chemgeo.2017.06.018

Nahon, D., & Tardy, Y. (1992). The ferruginous laterites. In *Handbook of exploration geochemistry* (Vol. 4, pp. 41–55). Elsevier. https://doi.org/10.1016/B978-0-444-89095-5.50010-9

Nardi, A., Idiart, A., Trinchero, P., de Vries, L. M., & Molinero, J. (2014). Interface COMSOL-PHREEQC (iCP), an efficient numerical framework for the solution of coupled multiphysics and geochemistry. *Computers & Geosciences*, **69**, 10–21. https://doi.org/10.1016/j.cageo.2014.04.011

Nardi, A., Idiart, A., Trinchero, P., de Vries, L. M., & Molinero, J. (2022). ICP (2.1) [Computer software]. https://en.amphos21.com/producto/icp/#

Navarre-Sitchler, A., Brantley, S. L., & Rother, G. (2015). How porosity increases during incipient weathering of crystalline silicate rocks. *Reviews in Mineralogy and Geochemistry*, **80**(1), 331–354. https://doi.org/10.2138/rmg.2015.80.10

Navarre-Sitchler, A., Steefel, C. I., Sak, P. B., & Brantley, S. L. (2011). A reactive-transport model for weathering ring formation on basalt. *Geochimica et Cosmochimica Acta*, **75**(23), 7644–7667. https://doi.org/10.1016/j.gca.2011.09.033

Neretnieks, I. (1980). Diffusion in the rock matrix: An important factor in radionuclide retardation? *Journal of Geophysical Research*, **85**(B8), 4379–4397. https://doi.org/10.1029/JB085iB08p04379

Ollier, C. D. (2010). Very deep weathering and related landslides. *Geological Society, London, Engineering Geology Special Publications*, **23**(1), 5–14. https://doi.org/10.1144/EGSP23.2

Orgogozo, L., Golfier, F., Buès, M. A., Quintard, M., & Koné, T. (2013). A dual-porosity theory for solute transport in biofilm-coated porous media. *Advances in Water Resources*, **62**, 266–279. https://doi.org/10.1016/j.adwwatres.2013.09.011

Pandey, S., & Rajaram, H. (2016). Modeling the influence of preferential flow on the spatial variability and time-dependence of mineral weathering rates: Preferential flow, long-term weathering. *Water Resources Research*, **52**(12), 9344–9366. https://doi.org/10.1002/2016WR019026

Parkhurst, D. L., Appelo, C., & others (2013a). Description of input and examples for PHREEQC version 3—A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations. *US Geological Survey Techniques and Methods*, **6**(A43), 497.

Parkhurst, D. L., & Appelo, C. A. J. (1999). User's guide to PHREEQC (Version 2): A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations. https://doi.org/10.3133/wri994259

Parkhurst, D. L., & Appelo, C. A. J. (2013). Description of input and examples for PHREEQC version 3—a computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations (Techniques and Methods).

Parkhurst, D. L., Appelo, C. A. J., & others (2013b). PhreeqC (Version 3) [Computer software]. USGS. https://www.usgs.gov/software/phreeqc-version-3

Pokrovsky, O. S., & Schott, J. (2000). Kinetics and mechanism of forsterite dissolution at 25°C and pH from 1 to 12. *Geochimica et Cosmochimica Acta*, **64**(19), 3313–3325. https://doi.org/10.1016/S0016-7037(00)00434-8

Reis, F. D. A. A., & Brantley, S. L. (2017). Models of transport and reaction describing weathering of fractured rock with mobile and immobile water. *Journal of Geophysical Research: Earth Surface*, **122**(3), 735–757. https://doi.org/10.1002/2016JF004118

Sak, P. B., Navarre-Sitchler, A. K., Miller, C. E., Daniel, C. C., Gaillardet, J., Buss, H. L., et al. (2010). Controls on rind thickness on basaltic andesite clasts weathering in Guadeloupe. *Chemical Geology*, **276**(3–4), 129–143. https://doi.org/10.1016/j.chemgeo.2010.05.002

Santoro, L., Putzolu, F., Mondillo, N., Boni, M., & Herrington, R. (2022). Trace element geochemistry of iron-(oxy)-hydroxides in Ni(Co)-laterites: Review, new data and implications for ore forming processes. *Ore Geology Reviews*, **140**, 104501. https://doi.org/10.1016/j.oregeorev.2021.104501
Sevin, B., Cluzel, D., Maurizot, P., Ricordel-Prognon, C., Chaproniere, G., Folcher, N., & Quesnel, F. (2014). A drastic lower Miocene regolith evolution triggered by post obduction slab break-off and uplift in New Caledonia: Slab break-off in New Caledonia. *Tectonics*, 33(9), 1787–1801. https://doi.org/10.1002/2014TC003588

Sherman, T., Sole-Mari, G., Hyman, J., Sweeney, M. R., Vassallo, D., & Bolster, D. (2021). Characterizing reactive transport behavior in a three-dimensional discrete fracture network. *Transport in Porous Media*, 146(1–2), 307–327. https://doi.org/10.1007/s11242-021-01568-4

Sillitoe, R. H. (2005). Supergene oxidized and enriched porphyry copper and related deposits. In J. W. Hedenquist, J. F. H. Thompson, R. J. Goldfarb, & J. P. Richards (Eds.), *One hundredth anniversary volume*. Society of Economic Geologists. https://doi.org/10.5382/AV100.22

Steefel, C. I., Appelo, C. A. J., Arora, B., Jacques, D., Kalbacher, T., Kolditz, O., et al. (2015). Reactive transport codes for subsurface environmental simulation. *Computational Geosciences*, 19(3), 445–478. https://doi.org/10.1007/s10596-014-9443-x

Steefel, C. I., & Hu, M. (2022). Reactive transport modeling of mineral precipitation and carbon trapping in discrete fracture networks. *Water Resources Research*, 58(9). https://doi.org/10.1029/2022WR032321

Stoessell, R. K. (1988). 25°C and 1 atm dissolution experiments of sepiolite and kerolite. *Geochimica et Cosmochimica Acta*, 52(2), 365–374. https://doi.org/10.1016/0016-7037(88)90092-0

Teitler, Y., Cathelineau, M., Ulrich, M., Ambrosi, J. P., Munoz, M., & Sevin, B. (2019). Petrology and geochemistry of scandium in New Caledonian Ni-Co laterites. *Journal of Geochemical Exploration, 196*, 131–155. https://doi.org/10.1016/j.gexplo.2018.10.009

Teitler, Y., Favier, S., Ambrosi, J.-P., Sevin, B., Golfier, F., & Cathelineau, M. (2022). Evaluation of Sc concentrations in Ni-Co laterites using Al as a geochemical proxy. *Minerals*, 12(5), 615. https://doi.org/10.3390/min12050615

Trescases, J. (1975). L'élvolution gélochimique supégèlne des roches ultrabasiques en zone tropicale: Formation des gisements nickélifères de Nouvelle-Calédonie: Mémoires ORSTOM (Office de la Recherche Scientifique et Technique Outre-Mer). In *Mémoires ORSTOM (Office de la Recherche Scientifique et Technique Outre-Mer)*.

Trinchero, P., Iraola, A., Bruines, P., & Gylling, B. (2021). Water-mineral reactions in a translated single realistic fracture: Consequences for contaminant uptake by matrix diffusion. *Water Resources Research*, 57(10). https://doi.org/10.1029/2021WR030442

Trinchero, P., Puigdomenech, I., Molinero, J., Ebrahimi, H., Gylling, B., Svensson, U., et al. (2017). Continuum-based DFN-consistent numerical framework for the simulation of oxygen infiltration into fractured crystalline rocks. *Journal of Contaminant Hydrology*, 200, 60–69. https://doi.org/10.1016/j.jconhyd.2017.04.001

Valeton, I. (1994). Element concentration and formation of ore deposits by weathering. *Catena*, 21(2–3), 99–129. https://doi.org/10.1016/0341-8162(94)90006-X

Vanwalleghem, T., Stockmann, U., Minasny, B., & McBratney, A. B. (2013). A quantitative model for integrating landscape evolution and soil formation: Modelling landscape and soil formation. *Journal of Geophysical Research: Earth Surface*, 118(2), 331–347. https://doi.org/10.1029/2011JF002296

Viers, J., & Oliva, P. (2018). Chemical weathering. In W. M. White (Ed.), *Encyclopedia of geochemistry* (pp. 237–241). Springer International Publishing. https://doi.org/10.1007/978-3-319-39312-4_143

Viswanathan, H. S., Ajo-Franklin, J., Birkholzer, J. T., Carey, J. W., Guglielmi, Y., Hyman, J. D., et al. (2022). From fluid flow to coupled processes in fractured rock: Recent advances and new frontiers. *Reviews of Geophysics*, 60(1). https://doi.org/10.1029/2021RG000744

West, A., Galy, A., & Bickle, M. (2005). Tectonic and climatic controls on silicate weathering. *Earth and Planetary Science Letters*, 235(1–2), 211–228. https://doi.org/10.1016/j.epsl.2005.03.020

Wogelius, R. A., Milodowski, A. E., Field, L. P., Metcalfe, R., Lowe, T., van Veelen, A., et al. (2020). Mineral reaction kinetics constrain the length scale of rock matrix diffusion. *Scientific Reports*, 10(1), 8142. https://doi.org/10.1038/s41598-020-65113-x

Zhang, Q., Deng, H., Dong, Y., Molins, S., Li, X., & Steefel, C. (2022). Investigation of coupled processes in fractures and the bordering matrix via a micro-continuum reactive transport model. *Water Resources Research*, 58(2). https://doi.org/10.1029/2021WR030578