# Migration of dissolution front in a fracture network – Implications for weathering of fractured bedrock systems and boulder formation

Sylvain Favier, Yoram Teitler, Fabrice Golfier, Michel Cathelineau

► To cite this version:

Sylvain Favier, Yoram Teitler, Fabrice Golfier, Michel Cathelineau. Migration of dissolution front in a fracture network – Implications for weathering of fractured bedrock systems and boulder formation. Journal of Hydrology, 2024, 644, pp.132056. 10.1016/j.jhydrol.2024.132056. hal-05241926

HAL Id: hal-05241926

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

Submitted on 6 Sep 2025

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
Research papers

# Migration of dissolution front in a fracture network – Implications for weathering of fractured bedrock systems and boulder formation

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

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

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

Check for
updates

## ARTICLE INFO

### Keywords:

- Weathering regimes

- Fracture Networks

- Boulders

- Reactive Transport Modeling

- Scaling laws

## ABSTRACT

Fracture networks contribute to the differential weathering of fractured bedrock at multiple scales. The influence of competing fracture networks is discussed in light of hydrogeochemical conditions to predict the dissolution pattern resulting from spheroidal weathering. The investigation of the dimensionless Péclet (Pe) and Damköhler (Da) numbers on a simple case of a dual-scale fracture network and its application to three fracture network weathering reveals the influence of fluid flow and reaction rate conditions on the dissolution patterns.

The conceptual model was numerically validated and compared with observations on lateritic nickel deposits in New Caledonia. A fracture network will contribute to boulder development if the associated fluid flow regime reaches a Pe threshold of around 10. On the other hand, predominant diffusion will inhibit the fracture set in favour of a higher permeability fracture network. The Da drives, for its part, the width of the dissolution halo. Owing to the dependency of the fluid flow regime on the fracture network characteristics, scaling laws were established to estimate the size of the boulders at different scales by predicting the fracture contribution as a function of their Pe.

## 1. Introduction

Located at the base of the Critical Zone, the regolith consists of weathered material blanketing unweathered rock. Regolith formation results from progressive alteration of the bedrock by downward percolating supergene waters; the process is controlled by the permeability and heterogeneities such as fractures acting as flow barriers or drains (Favier et al., 2022; Iseppi et al., 2018; Munday et al., 2003). Weathering processes occurring in the subsurface are essential in the flux dynamics across the earth-atmosphere interface. They were thus described and modeled to understand better the resulting vertical distribution of chemical elements (Brantley and Lebedeva, 2011). However, the complex distribution of structural heterogeneities leads to an uneven dissolution and rough weathering front at the bedrock-saprolite interface, and hence an irregular surface and thickness of the regolith (Lebedeva and Brantley, 2017; Navarre-Sitchler and Brantley, 2007). As a result, in supergene Ni-Co deposits, such heterogeneous regolith distribution implies an uneven Ni-grade distribution, which is still not understood (Cathelineau et al., 2016; Myagkiy et al., 2019; Teitler et al., 2019, 2022). Furthermore, weathering also affects rock mechanical

properties, notably their hardness, with important implications for mining activities. Hence, an appropriate evaluation of spatial heterogeneities in the regolith is critical for mineral resource development and production in the mining industry.

Boulders, often associated with spheroidal, concentric or onion-skin weathering (Fletcher et al., 2006; Md Dan et al., 2015), are common weathering heterogeneities encountered at the bedrock-saprolite interface. Boulders typically consist of a weathered outer rim and a mostly unweathered core. In the most weathered cases, they are usually shelled by a thin saprolite layer surrounded by laterite (Md Dan et al., 2016; Ollier, 2010; Patino et al., 2003). This chemical dissolution process is usually accompanied by spheroidal fracturing that forms rindlets or "onion-skins" (Fletcher et al., 2006; Goodfellow and Hilley, 2022). In this work, however, we will only investigate weathering rind formation (i.e. the chemical dissolution) and not consider the weathering-induced fracturing (i.e., the hydromechanical coupling). The boulder development occurs on various rock types, including gneiss, schist, andesite, sandstone, granites, basalts, and peridotites. It is considered a common mechanism initiating the transformation of bedrock to saprolite (R. Fletcher et al., 2006) and resulting from differential weathering

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: favier.syl@wanadoo.fr (S. Favier), yoram.teitler@univ-lorraine.fr (Y. Teitler), fabrice.golfier@univ-lorraine.fr (F. Golfier), michel.cathelineau@univ-lorraine.fr (M. Cathelineau).</footnote>
(Sarracino et al., 1987). Weathering starts from a fracture network delimiting bedrock blocks. Then, the weathering front progressively develops inwards, creating boulders (Hirata et al., 2017; Navarre-Sitchler et al., 2011; Sarracino et al., 1987). At the scale of the fracture network, a focus has been made on the balance between chemical weathering and erosion rates to predict the block size distribution (Brantley et al., 2011; Fletcher and Brantley, 2010; Lebedeva and Brantley, 2017). The bedrock fracture spacing controls the initial boulder size (Fletcher and Brantley, 2010; Hausrath et al., 2011; Sklar et al., 2017), which also drives the rough weathering front progression of rock mass (Lebedeva and Brantley, 2017). A heterogeneous block size distribution is described depending on the balance between the resulting weathering and erosion rates (Brantley and Lebedeva, 2021).

Despite considerable progress on several fronts, the intricate couplings at stake with many physicochemical parameters (e.g., temperature, pH, ionic strength, saturation index, water content) and mechanisms (evolution of reactive surface area with rock alteration, mineral armouring due to precipitation of insoluble secondary compounds) influencing chemical weathering rates (Anbeek, 1993; Emmanuel, 2022; Jung and Navarre-Sitchler, 2018b) represent a challenging task. Such difficulty is notably illustrated by the discrepancy observed between laboratory and field-measured dissolution rates with several orders of magnitude difference (Navarre-Sitchler and Brantley, 2007; White and Brantley, 2003), which reflects our inability to predict weathering rates at large scale. One of the main reasons behind this gap is the function of which limiting reaction mechanism prevails (kinetics limited or transport limited reaction). Laboratory experiments typically take place under well-stirred conditions and situated far from equilibrium, while typical field conditions exhibit a decrease in weathering rate over time due to the accumulation of dissolution products as the alteration front moves downward (Maher, 2010; Navarre-Sitchler and Thyne, 2007). Such observation has to be balanced in the presence of high erosion rates, which limit the reaction front extent and, hence, the evolution toward chemical equilibrium (Dixon et al., 2012; Schoonejans et al., 2016). Thus, kinetic-limited weathering is typically observed in actively eroding mountains (West et al., 2005), and the limiting mechanism is a good indicator of tectonic or climatic control on chemical weathering (Mao et al., 2022; Riebe et al., 2017). As an illustrative example, the weathering processes of peridotites in New Caledonia, which have resulted in laterite formation, are characterized by olivine dissolution under a hot, wet tropical climate with high erosion rates (Cathelineau et al., 2017; Teitler et al., 2019). Reactive transport simulations conducted by Favier et al. (2024) highlight that permeability heterogeneities control weathering propagation and olivine dissolution is driven by kinetics (see the calculations of olivine saturation index SI in Supporting Information S1 of Favier et al. (2024), as large negative values indicate that reaction occurs far from equilibrium). We will focus hereafter only on the kinetics dissolution regime, and the weathering of laterite supergene deposits in New Caledonia will be used for application purposes.

Pioneering works on dissolution processes in porous media (Golfier et al., 2002; Ladd and Szymczak, 2021; Szymczak and Ladd, 2006) and later in fractures (Deng and Peters, 2019; Detwiler et al., 2003; Detwiler and Rajaram, 2007; Elkhoury et al., 2013; Oltéan et al., 2013; Trinchero et al., 2021; Yuan et al., 2019; Zhang et al., 2022), have shed light on the influence of dimensionless numbers, respectively Péclet and Damköhler numbers, to predict dissolution patterns. At the field scale (10–100 m), most weathering models have initially investigated the 1D vertical dissolution front (e.g., see Goddéris et al., 2019, for a general overview of these models). With the advances in modelling techniques of flow and transport in fractured mass rock (we refer the reader to the state-of-the-art review papers of MacQuarrie and Mayer, 2005, or Deng and Spycher, 2019) and computational power, however, recent numerical studies have attempted to explore this complex interplay between weathering processes and structural pathways in 2D and 3D fracture-matrix systems (e.g., Andrews et al., 2023; Andrews and Navarre-Sitchler, 2021; Favier

et al., 2024; Jung and Navarre-Sitchler, 2018b, 2018a). Some results have emphasized the role of topography and lateral flows on weathering patterns in hills (Lebedeva and Brantley, 2017; Myagkiy et al., 2019), the competition between weathering and erosion and the impact of water table depth (Lebedeva and Brantley, 2020) or fracture densities (Lebedeva and Brantley, 2023). Simulations in the two-dimensional networks of Andrews and Navarre-Sitchler (2021) indicate that fractures act as fast flow paths through the rock mass where mineral dissolution remains driven by kinetics, while dissolution in the porous matrix occurs close to equilibrium and is limited by diffusion. Andrews et al. (2023) have extended these findings in three-dimensional fractured media. Steefel and Hu (2022) have shown how the weathering processes could affect the connectivity of a natural intersecting fracture network. Finally, simulations of Pandey and Rajaram (2016) and Favier et al. (2024) on Discrete Fracture Matrix (DFM) models have exhibited corestone features or boulders similar to the ones observed on the field. Both results highlight the role of the original fracture network on boulder sizes.

However, relationships between the boulder size distribution, the geometrical organization of fracture networks, and their hydrological characteristics have not yet been fully addressed. Indeed, depending on weathering fluids percolating through the fracture network, the weathering front geometry may be affected. More specifically, fracture networks tend to display a multiscale nature (Ehlen, 2000, 2002), which impacts flow properties, notably the network permeability at different scales (Davy et al., 2006; de Dreuzy et al., 2012). When multiple fracture sets are present, each one may have a different permeability, constraining the fluid flow. Thus, competition may occur between the different fracture sets. Highly permeable fractures channel the central part of fluid flow so that fluid circulation along the secondary network with less permeable discontinuities may be inhibited. To our knowledge, such competition has been studied in the context of flow patterns for oil recovery but never in the presence of dissolution (Gong and Rossen, 2018). This fluid-flow partitioning between fractures and rock matrix is exacerbated by the fluid-rock interactions which may occur. Therefore, this competition may result in different boulder size distributions depending on which fracture network is promoted relative to another Fig. 1).

This paper studies the impact of a geometrical organization (fracture spacing) and hydrological characteristics (permeability) of a multiscale fracture network composed of sub-orthogonal fractures on chemical weathering front propagation. The hydrogeochemical conditions under which a fracture network contributes to weathering at the expense of another are explored to predict better the dissolution pattern resulting from spheroidal weathering. The hydrogeochemical settings are first investigated in the light of dimensionless Péclet (Huysmans and Das-sargues, 2005) and Damköhler (Dijk et al., 2002) numbers, which describe the fluid flow and the chemical reaction regime, respectively. A simple case of a dual-scale fracture network was designed to study the impact of their variations on dissolution patterns. Then, the findings were upscaled to multiscale fracture networks using scaling laws. An application to lateritic deposits in New Caledonia was carried out to predict the boulder size observed at different scales. Indeed, local mining activity provides a rare opportunity to access laterite profiles developed after peridotite in which the fracture networks have been highly documented (Cathelineau et al., 2017; Favier et al., 2022; Iseppi et al., 2018; Jeanpert et al., 2019; Quesnel et al., 2016). Such application highlights the potentialities to predict highly variable-grade distribution resulting from weathering heterogeneities in supergene deposits.

## 2. First insight from a simple two-scale fracture network

A first analysis is conducted on a general simple dual-scale fracture network to highlight the physical mechanisms at stake. Different hydrological settings and dissolution rates were explored to reveal the conditions under which a secondary fracture network (F2) will constrain
<figure><img src="image_3.png" /><figcaption>Fig. 1. a) Multiscale fracture networks shaping different boulder sizes at outcrop scale in new caledonia. the primary fracture network in white is f1 at observation scale $l_1$ with a fracture spacing $s_1$. b) Zoom of the yellow rectangle. In yellow is the secondary fracture network F2 at observation scale $l_2$ with a fracture spacing $s_2$. c) Typical boulder shape in the saprolite horizon. d) Conceptual multiscale fracture networks composed of a primary fracture network F1 and a secondary fracture network F2. e) Conceptual model of boulder formation resulting from the percolation of an acidic solution in the primary fracture network F1 (in black). In such a configuration, the secondary fracture network (in green), F2, does not contribute to the dissolution pattern. f) Conceptual model of boulder formation where both F1 and F2 contribute to the dissolution pattern. In brown, the resulting boulders. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

with the main network (F1), the boulder size, and conversely. Thus, the competition between two fracture sets, F1 and F2, with different permeabilities, and the impact on boulder size formed from rock matrix dissolution was studied numerically.

Variations of fracture permeability are first carried out to investigate the impact of changes in fluid flow regimes on dissolution patterns. Then, variations of the constant kinetics rate are realized for a given hydraulic configuration to study its influence on dissolution patterns.

## 2.1. Conceptual model

### 2.1.1. Flow equations

A simple geometry of orthogonal fractures was used. The first network, F1, consists of high-permeability fractures spaced 8 m apart. The second one, F2, consists of low-permeability fractures spaced 4 m apart. Their permeability, $k_f$ [m²], is constrained by their hydraulic aperture $a_f$ [m], according to the cubic law (Snow, 1969):

$$k_f = \frac{a_f^2}{12} \qquad (1)$$

A pressure condition is applied to the upper left and lower right corners of the model to achieve an even flow between the vertical and horizontal fractures (Fig. 2). A periodic condition was applied to the remaining sides of the model. 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 (Favier et al., 2024).

The Darcy law models the fluid flow in both fractures and matrix:

$$ \text{Matrix :} \quad \left\{ \begin{array}{l} \nabla \cdot (\mathbf{u}_m) = 0 \\ \mathbf{u}_m = -\frac{k_m}{\mu} \nabla p_m \end{array} \right. \qquad (2.a) $$

$$ \text{Fractures :} \quad \left\{ \begin{array}{l} \nabla \cdot (\mathbf{u}_f) = \frac{Q_{f,0}}{a_f} \\ \mathbf{u}_f = -\frac{k_f}{\mu} \nabla p_f \end{array} \right. \qquad (2.b) $$

where $k_i$ [m²], $u_i$ [m/s], and $p_i$ [Pa] are respectively the permeability, the fluid velocity and the pressure of either the matrix or the fractures, and $\mu$, the dynamic viscosity [Pa.s]. The permeability $k_m$ of the matrix has been set to a very low value to ensure that all the fluid flow is channeled inside fractures. Its value agrees with the permeability values of un-weathered peridotites in New Caledonia measured at a sample scale to avoid being impacted by fractures (Favier et al., 2022). $Q_{f,0}$ [m/s] is the out-of-plane momentum flux between fractures and matrix (COMSOL Multiphysics, 2022b).
<figure><img src="image_4.png" /><figcaption>Fig. 2. Geometry and boundary conditions of the model.</figcaption></figure>

The hydrogeological parameters are detailed in Table 1.

According to the cubic law, the fracture hydraulic aperture directly influences the fracture permeability and flow rate. Different F1 and F2 aperture values are thus tested to assess their influence on the competition between the two networks. Hydraulic aperture ranges from 10 µm to 700 µm for F1 and from 0.5 µm to 100 µm for F2 to capture the different fluid flow regimes.

### 2.1.2. Reactive transport equations

The conceptual reactive transport model is intended to be as general as possible to tackle the broadest range of situations. Hence, a fictive simple geochemical system is studied (Eq. (4)). Bedrock is assumed to be dissolved by an acidic species A in solution coming from the upper left corner of the model (Fig. 2). Fracture networks are considered open or partially filled with inert material so that only bedrock blocks can be dissolved. However, porosity changes due to rock alteration are neglected in this model since we consider that the mass fraction of insoluble compounds inert to weathering is significant in rock material. Thus, porosity and permeability are kept constant. This critical assumption should tend to reduce the onset and growth of dissolution-driven instabilities but is not expected to drastically affect the weathering regimes found hereafter or the main outcomes.

<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>10<sup>-3</sup></td><td>Pa.s</td><td></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>30%</td><td>%</td><td>(Myagkiy et al., 2019)</td></tr><tr><td>Permeability: matrix</td><td>k<sub>m</sub></td><td>10<sup>-20</sup></td><td>m<sup>2</sup></td><td>(Favier et al., 2022)</td></tr><tr><td>Hydraulic aperture fracture</td><td>a<sub>f</sub></td><td>variable</td><td>m</td><td></td></tr><tr><td>Permeability: fracture</td><td>k<sub>f</sub></td><td>variable</td><td>m<sup>2</sup></td><td></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></tbody></table>

Reactive mass transport of the species in solution is written as follows using a discrete fracture matrix approach, where 1D fractures are embedded in a 2D matrix (Fumagalli and Keilegavlen, 2019; Hyman et al., 2022; Odsæter et al., 2019; Favier et al., 2024):

$$ \text{Matrix:} \quad \left\{ \begin{aligned} \mathbf{J}_A^m &= (-D_e \nabla c_A^m) + \mathbf{u}_m c_A^m \\ \frac{\partial \epsilon_m c_A^m}{\partial t} + \nabla \cdot \mathbf{J}_A^m &= R_A \\ \frac{\partial q_A^m}{\partial t} &= R_A \end{aligned} \right. \qquad (3.a) $$

$$ \text{Fractures:} \quad \left\{ \begin{aligned} \mathbf{J}_A^f &= (-D_e \nabla c_A^f) + \mathbf{u}_f c_A^f \\ \frac{\partial \epsilon_f c_A^f}{\partial t} + \nabla \cdot \mathbf{J}_A^f &= \frac{J_{f,0}}{a_f} \end{aligned} \right. \qquad (3.b) $$

where $\mathbf{J}_A^m$ and $\mathbf{J}_A^f$ are the total flux in the matrix and the fractures, $\epsilon_m$ and $\epsilon_f$ are the matrix and fracture porosity. $q_A^m$ [mol/m$^3$] corresponds to the molar concentration of rock in the matrix, while $c_A^m$ and $c_A^f$ [mol/m$^3$] correspond to the acid concentrations in the liquid phase in the matrix and the fractures. $J_{f,0}$ [mol/m$^2$/s] is the out-of-plane flux from the fractures (Favier et al., 2024). This source term accounts for the coupling between fractures and matrix. Such a term is calculated as a function of the total molar fluxes at the upside and downside of the fracture (COMSOL Multiphysics, 2022a) so that mass flux continuity condition is ensured at the matrix-fracture interface.

As specified in the Introduction, we will focus on weathering processes far from equilibrium conditions, such as the weathering of peridotite in New Caledonia, which will be proposed as an illustrative example. The olivine dissolution has been studied extensively (Anbeek, 1993; Oelkers et al., 2018) and is used here as a proxy for various ultramafic rocks (Favier et al., 2022; Freyssinet and Farah, 2000; Trescases, 1975). For such conditions, the olivine dissolution rate is essentially constrained by pH variations (Oelkers et al., 2018; Wogelius and Walther, 1992), contrarily to transfer-limited reactions classically described from Transition State Theory. For the sake of simplicity, the overall dissolution rate $R_A$ [mol/m$^3$/s] of the bedrock by solution weathering is assumed to follow a first-order kinetics such as:

$$ R_A = -k_A c_A^m \qquad (4) $$

where $k_A$ is the dissolution rate constant [s$^{-1}$] and the acidic solution is limiting. The simplicity of the geochemical model allows us to better highlight the impact of the fracture network properties, i.e. spacing and permeability, on the weathering regimes. Furthermore, despite its simplicity, such a model is consistent in a first-order approach, with more realistic geochemical models of mafic-to-ultramafic bedrock weathering (Myagkiy et al., 2017; Pokrovsky and Schott, 2000; Røyne et al., 2008). In addition, and for the sake of comparison, all results and variables will be presented as dimensionless, as defined in Section 2.1.2. Results will be thus interpreted in terms of governing dimensionless numbers to decipher the contribution of different mechanisms and physical parameters to the weathering pattern (Detwiler and Rajaram, 2007; Golfier et al., 2002). A more complex rate law of dissolution kinetics could eventually be analyzed by modifying the Damköhler number definition.

Using COMSOL Multiphysics software, we numerically solved the 2D problem (Eqs. (2) and (3)) with a finite element method. The effect of the flow regime on dissolution patterns is first investigated using a constant dissolution rate $k_A = 10^{-7} s^{-1}$. Then, we study the impact of kinetics by testing variable dissolution kinetic rates ranging from $10^{-5} s^{-1}$ to $10^{-9} s^{-1}$ for two different flow conditions.

### 2.1.3. Dimensionless driving parameters

A dimensionless analysis exhibits four dimensionless numbers that
drive the propagation of the dissolution front. First, two Péclet numbers, Pe₁ and Pe₂, respectively associated with F1 and F2 fracture sets, define the flow regime (advection versus diffusion-dominated flow) in fractures. Since diffusion is the prevailing transport mechanism within the matrix, this Péclet number quantifies the amount of fluid that diffuses into the matrix relative to the one transported along the fractures. It is classically defined as follows (Huysmans and Dassargues, 2005):

$$Pe = \frac{t_{\text{diffusion}}}{t_{\text{advection}}} = \frac{L u_f}{D_e} = \frac{\sqrt{k_f} u_f}{D_e} = \frac{a_f}{2\sqrt{3}} \frac{u_f}{D_e} \quad (5)$$

where the square root of permeability is used as the characteristic length of fracture (Guéguen et al., 1996). This permeability is a function of fracture aperture, using Eq. (1). A similar definition is used for all the simulations described in the following sections.

In addition, the two last dimensionless numbers referred to as the Damköhler numbers, Da₁ and Da₂, respectively associated with F1 and F2 fracture sets, assess the relative importance of reaction kinetics and solute transport by diffusion within the matrix. It is expressed as:

$$Da_i = \frac{t_{\text{diffusion}}}{t_{\text{chemical}}} = \frac{\left(\frac{s_i}{2}\right)^2}{D_e} \cdot k_A \quad (6)$$

where sᵢ/2 [m] is the half-length of the porous matrix block surrounded by the fracture set i. Depending on the fracture set, the characteristic length sᵢ may vary, hence the corresponding matrix block size, as illustrated in Fig. 2. Thus, Da₁ and Da₂ are respectively defined with s₁ = 8 m for F1 and s₂ = 4 m for F2.

## 2.2. Results

Many simulations have been performed for different fracture apertures (i.e., different Pe₁ and Pe₂ values, see Section 2.2.1) or different reaction rates (different Da, see Section 2.2.2). All the other simulation parameters are listed in Table 1. Such modifications of flow regimes and reaction kinetics impact the distribution of weathered rock masses and, thus, result in different dissolution patterns (i.e. different concentration fields qₐ) as discussed below.

### 2.2.1. Impact of flow regimes

Fig. 3 shows the different dissolution patterns obtained numerically from simulations and plotted in a Pe₁ – Pe₂ diagram. As previously mentioned, the Péclet number represents the competition between diffusion and advection, so a high Pe corresponds to an advection-dominated configuration. Thus, at low Pe (Pe₁ < 1 and Pe₂ < 1), diffusion guides the rock dissolution (patterns 13 and 14). In contrast, patterns 10, 11 and 12 represent intermediate cases where the impact of advection, and thus of the F1 and F2 fracture networks on rock dissolution, begins to be visible. These patterns are observed for Pe between 1 and 10. Pattern 10 corresponds to a case where only F1 contributes to the weathering in addition to diffusion. Pattern 11 corresponds to an intermediate case where F2 starts to impact the dissolution, while for pattern 12, F2 drives dissolution. This gradation from F1 to F2 control is related to the increase of Pe₂. At a Pe₁ > 10 and Pe₂ < 1, the diffusion effect is no longer observed on the F1 network, so F1 solely drives the boulder shape. F2 begins to contribute to the boulder shape with F1 when Pe₂ is comprised between 1 and 10 (patterns 3, 6 and 7), and beyond Pe₂ = 10, both F1 and F2 drive the dissolution pattern when Pe₁

<figure><img src="image_5.png" /><figcaption>Fig. 3. Dissolution pattern of the rock depending on the flow regime in primary fracture network F1, and secondary fracture network F2. Regime I: Diffusion-dominated weathering regime. Regime II: Primary fracture-dominated weathering regime. Regime III: Secondary fracture-dominated weathering regime. Regime IV: Transition regime between the three main weathering regimes.</figcaption></figure>
> 10 (patterns 4 and 8).

Hence, to enhance the contribution of secondary network F2 at the expense of F1, a Pe higher than 10 is required. On the contrary, to constrain weathering patterns only by F1 at the expense of F2, a Pe₂ lower than 1 is required.

Three main weathering regimes can be thus identified based on this analysis:

**Diffusion-dominated weathering regime** (regime I, Pe₁ < 1 and Pe₂ < 1, representative pattern: 14): Diffusion mechanism dominates over transport by convection, and instabilities cannot develop. A uniform dissolution front is observed.

**Primary fracture-dominated weathering regime** (regime II, Pe₁ > 10 and Pe₂ < 1; representative pattern: 1): Mass transport in fractures F1 is driven by advection, but diffusion-dominated flow still occurs in fractures F2. A shrinking core mechanism reveals a boulder-shaped dissolution pattern delineated by the F1 fracture network (Favier et al., 2024; Pandey and Rajaram, 2016).

**Secondary fracture-dominated weathering regime** (regime III, Pe₁ > 10 and Pe₂ > 10; representative pattern: 4): Mass transport in both fracture sets is dominated by advection. Progression of the dissolution front results in a boulder-shaped dissolution pattern delineated by both the F1 and F2 fracture networks. The complete imprint of the original fracture network is preserved.

A transitional regime between the main weathering regimes is also identified (regime IV). Different representative patterns (6, 11, 13) are observed depending on the transition between the main regimes.

Note that boundaries between the different regimes have been drawn based on a qualitative observation of the dissolution patterns, and some uncertainties are associated with the transition regions. Fig. 3 highlights that the transition between the different dissolution regimes is only driven by the values of the two Péclet numbers and does not depend on the Damköhler number. A similar finding has been found for wormholing regimes in homogeneous porous media (Golfier et al., 2002). In conclusion, the changes in fracture permeability through increased hydraulic aperture drive the dissolution pattern and, thus, the boulder size.

### 2.2.2. Impact of reaction kinetics

Fig. 4 shows the different dissolution patterns obtained numerically by varying Da at constant Péclet number. Two representative values of Pe₂ are considered, 1 (regime II) and 10 (regime III), respectively, while Pe₁ is kept constant and equal at 500. In regime II, only the primary fracture network is expressed, and Da₁ is considered for analysis, while in regime III, both fractures are represented, and Da₂ is preferentially calculated. A similar Da value between the two regimes does not necessarily correspond to the same reaction rate value. Indeed, in Regime III (pattern 7), for instance, Da₁ ≈ 400 corresponds to a rate k_A = 5·10⁻⁸ s⁻¹, while the same Da₂ value is obtained in regime II for k_A = 10⁻⁸ s⁻¹. However, both patterns are remarkably similar, illustrating the competition between the characteristic times of diffusion and reaction kinetics. Consequently, we do not differentiate Da₁ and Da₂ in Fig. 4.

We observe that Da variations induced by changes in dissolution rate

mainly influence the size of the dissolution halo but not the weathering regime itself (Fig. 4), as explained in Li et al. (2014). A high Da reduces the size of the dissolution halo and favors the formation of well-delineated patterns, whereas a low Da favors a more diffuse dissolution of the whole rock. Indeed, for Da ≫ 1, the reaction rate is far greater than the diffusion rate. The reaction is thus assumed to occur instantaneously at the boulder rim, creating a sharp transition with the undis-solved boulder part. On the other hand, for Da close to 1, reaction and diffusion co-occur. Hence, the dissolution is more uniform within all the boulders. Such results are consistent with previous studies regarding the impact of the Da and the reaction rates on the width of the reaction front (Genthon et al., 2017; Genthon and Ormond, 2008; Li et al., 2014; Panga et al., 2005; Szymczak and Ladd, 2009).

## 3. Weathering front in a multiscale fracture network: application to boulder size prediction in New Caledonia

### 3.1. A theoretical weathering model in fractured porous media

The simple case study of the two-scale fracture network highlights the influence of the fluid flow regime on fracture-driven weathering. We observed that this dependency can be expressed solely through the Péclet number. More explicitly, for Pe > 10, the fracture set contributes to weathering. Assuming a similar relationship can be observed across different scales ranging between a lower and an upper cutoff scale, we may predict the large-scale chemical weathering features and boulder formation. Indeed, the fluid flow regime in fracture networks is notably driven by the fracture spacing and their hydraulic aperture. It is noteworthy that these parameters tend to follow scaling laws (Davy et al., 2006; de Dreuzy et al., 2012), classically expressed as power-law distributions, such as

$$a_i(l) = \alpha l^b \qquad (7)$$

where α is a constant, b is a power-law exponent and l [m] is the observation scale of interest.

Deriving a similar scaling relationship for fluid flow in fractures is usually more complex (Pyrak-Nolte and Nolte, 2016). However, based on Eq. (5), the Darcy law given in Eq. (2.b) and the cubic law given in Eq. (1), Peᵢ for a fracture set i can be expressed as a function of only two parameters, i.e., aᵢ and the average hydraulic gradient ||∇p||, drives the fluid flow within the fracture network.

$$Pe_i = f(a_i, ||\nabla p||) \qquad (8)$$

Combining Eqs. (7) and (8), we may thus write a scalar relationship for Peᵢ so that we have:

$$Pe_i(l) = \alpha(||\nabla p||)l^b \qquad (9)$$

Assuming the previously determined threshold for Peᵢ (Peᵢ = 10) applies to any scale between a lower and an upper cutoff scale, estimating the fracture network size beyond which fracture features (spacing) are

<figure><img src="image_6.png" /><figcaption>Fig. 4. Impact of Da on the dissolution pattern.</figcaption></figure>
preserved and, hence, the average boulder size at the outcrop scale should be possible.

To illustrate this approach, the boulder size is discussed in the olivine dissolution case study, which uses the example of New Caledonia lateritic deposits.

## 3.2. Application to laterite deposits in New Caledonia

### 3.2.1. Geological description

New Caledonia is an archipelago located in SW Pacific (Fig. 5). Ultramafic rocks account for about 30% of the exposed surface of the main “Grande Terre” island, as relics of the Peridotite Nappe obducted during the Eocene period (Cluzel et al. 2001). After obduction, the Peridotite Nappe underwent deep weathering under a hot and wet tropical climate, resulting in laterite formation (Sevin et al., 2014; Maurizot et al., 2020). The development of lateritic deposits is caused by the leaching of soluble elements (Mg, Ni, Si) from the peridotite bedrock. Dissolved elements reprecipitate deeper in the saprolite horizon as newly formed phyllosilicates and Fe-oxy-hydroxides. Finally, the newly formed phyllosilicates are also dissolved to enable the formation of an oxide-rich horizon (Butt and Cluzel, 2013; Cathelineau et al., 2017; Dublet et al., 2012; Freysinet et al., 2005; Golightly, 2010; Haldar, 2017; Santoro et al., 2022; Teitler et al., 2019). Weathering is driven by permeability heterogeneities and a dense fracture network (Favier et al., 2022, 2024; Iseppi et al., 2018; Jeanpert et al., 2019).

Due to the intense fracturing of bedrock at different scales (Favier et al., 2022; Iseppi et al., 2018), weathering heterogeneities such as boulders in the saprolite zone (Myagkiy et al., 2019) are observed from the deposit scale to the sample scale (Fig. 1, Fig. 6). Boulders are characterized by a concentric weathering zonation: the outer part is more weathered than the core (Myagkiy et al., 2019).

At each scale, they display an evenly spaced distribution so that a scaling law of their size can be stated by measuring the distance between the centre of neighbouring boulders or, when visible enough, between the delineating fractures (Fig. 6).

### 3.2.2. Scaling of fracture network

Fractured rock masses exhibit inherent spatial variability with different levels of complexity at various scales, which are complex to consider in flow or geomechanical modelling (Ma et al., 2022b, 2022a). This study uses a fractal analysis approach to address this issue and transfer flow properties from one scale to another. On the other hand, while the scaling law of fracture spacing can be derived from field measurements at different scales (Fig. 7.a), the distribution of fracture permeability is a priori unknown. It is possible, however, to estimate this relationship from field permeability measures of the fractured porous medium under certain assumptions (Fig. 7.b). The measures come from various sources. The fracture spacing was collected from Favier et al. (2022), and image interpretation was done by the authors (Fig. 6). The field permeability was derived from Favier et al. (2022) and the project HyperK (Jeanpert, 2017; Join et al., 2017) of the National Center for Technological Research of New Caledonia on nickel (CNRT).

Assuming geometrical considerations, the equivalent permeability of a fractured porous medium, $k_{eq}$ [m$^2$], can be expressed as a function of fracture spacing and fracture and matrix permeability, respectively $k_f$ and $k_m$:

$$k_{eq} = r_f k_f + r_m k_m \qquad (10)$$

where $r_f$ and $r_m$ are the contribution of the fractures and the matrix, respectively, to the apparent permeability $k_{eq}$, $r_f$ and $r_m$ are related to the fracture spacing $s$ [m] and the hydraulic aperture $a_f$ [m] (Snow, 1969). Hence, as a first order approximation, in the case of parallel fractures, $r_f$ and $r_m$ can be expressed as (Snow, 1969):

$$\left\{ \begin{array}{l} r_f = \frac{a_f}{s} \\ r_m = 1 - \frac{a_f}{s} \end{array} \right. \qquad (11)$$

Furthermore, assuming an extremely low permeability matrix, as is the case in New Caledonian peridotites (Favier et al., 2022; Jeanpert, 2017), and according to the cubic law, the equivalent permeability of porous medium can be simplified as follows:

<figure><img src="image_7.png" /><figcaption>Fig. 5. Simplified geological map of New Caledonia.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 6. Boulder size distribution at a) deposit, b) outcrop, and c) sample scales. d) Scaling law of boulder size.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7. Scaling laws of a) fracture spacing, b) permeability of the fractured porous medium, c) fracture hydraulic aperture, d) fracture permeability.</figcaption></figure>

$$k_{eq} = \frac{a_f a_f^2}{s 12} \qquad (12)$$

At this point, a scaling law of the fracture hydraulic aperture and their permeabilities can be easily derived from the equivalent porous permeability of the medium and the fracture spacing (Fig. 7.c and d).

Finally, as explained in section 3.1, a similar scaling relationship can be established for the Péclet number Pei. It indicates the Pe value of the largest fracture network Fi at a given scale.

Assuming that the average pressure gradient is known, the Pe value may be written either as a function of the observation scale l or the fracture spacing s and the following relationships are thus obtained for peridotites in New Caledonia:

$$Pe_i(l) = (1.45 \cdot 10^{-4} \|\nabla p\|) l^{2.38} \qquad (13)$$

$$Pe_i(s) = (1.65 \cdot 10^{-2} \|\nabla p\|) s^{2.38} \qquad (14)$$
Since the Péclet scaling law was derived from the fracture spacing law and the permeability scaling law, the lower and upper cutoff scale should be, in a first-order approach, in the range of those of the fracture spacing laws (i.e., 0.1 m and 100 m respectively), as they are the narrowest range.

### 3.3. Boulder size prediction on a multiscale fracture network

Before applying our theoretical weathering model to geological observations, we must assess the applicability at any scale of the Pe threshold value required for promoting the boulder delineation. For this purpose, we extend our model to a three-scale fracture network and simulate the formation of the boulders in the presence of three fracture sets, F1, F2, and F3, yielding different aperture and spacing values. Fracture spacings and permeabilities (i.e. apertures) are derived from the previously described scaling laws (Fig. 7). Note that this dependency implies that fracture spacing, aperture and permeability are linked. In other words, a densely spaced fracture network will display a high permeability. Different weathering patterns have been simulated over a domain of 20 × 15 m (patterns 1–4, as identified in Fig. 8), and results illustrate the competition between the fracture networks.

The simulation is conducted on COMSOL Multiphysics, with the same hydrogeochemical settings as the dual-scale fracture network case

study. The only difference with the numerical conditions illustrated in Fig. 2 lies in the fracture network geometry and the 0.7 m/m hydraulic gradient corresponding to the median hydraulic gradient in the Koniambo massif (Jeanpert, 2017). The gradient is applied at the upper boundary of the model, while the lower boundary has an outflow condition. A periodic condition is used on the other sides of the domain.

Results support previous results as the contribution of any fracture is controlled by the value of Pe (Fig. 8). In pattern 1, a weathered bedrock is created where only F1 participates in boulder delineation. Permeabilities of F2 and F3, which have, in that case, the same spacing, are too small to reach a Péclet number of 10 (Fig. 8 and Fig. 9). Pattern 2 represents a fractured medium where both F1 and F2 participate in weathering. F1 boulders cluster F2 boulders, while F3 is not expressed. In pattern 3, all fracture networks are represented in agreement with their Pe above 10. Large boulder clusters are delineated by F1, inside which smaller clusters are observed, delimited by F2, which are themselves delineated into boulder clusters by F3. Pattern 4 is an extreme case where F1 is much larger than F2 and F3, which exhibit the same geometrical features. Hence, the large boulders shaped by F1 are then cut into smaller ones by F2 and F3, which belong to the same network. These patterns are also consistent with the spatial organisation of the boulders encountered on the field. For example, pattern 2 (Fig. 8) mimics the outcrop structure described in Fig. 1. A primary fracture

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Boulder formation and clustering as a function of the contribution of each fracture network F1, F2, and F3. A zoom on a part of the initial domain (6 × 6 m size) is done to illustrate the fracture sets F2 and F3.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Simulated fluid regime for the different fracture networks F1, F2, and F3 in the patterns 1 to 4. In the black line, the calculated scaling law for Peᵢ is obtained from Eq. (14).</figcaption></figure>

network F1 (white dashed lines in Fig. 1) clusters boulders delineated by the secondary fracture network F2 (yellow dashed lines in Fig. 1).

The different weathering conditions encountered for each simulated pattern and each fracture set are shown in Fig. 9 as a function of the Péclet number value and fracture spacing. All the configurations are consistently aligned with the scaling law, given in Eq. (14), also plotted as a reference. Fig. 9 shows that hydrodynamic conditions contributing to boulder formation are only above the critical Péclet value (Pe = 10).

Therefore, this three-scale case study reveals the applicability of the Péclet threshold method on a multiscale system and the spatial hierarchical distribution of boulders resulting from this competition. A clustering of boulders is indeed observed based on the permeability and spacing of the fracture network. Large fracture networks are characterized by high permeability. Thus, a high Pe delineates large boulders, divided by smaller fracture networks characterized by lower permeability and Pe, into clusters of smaller boulders. Such phenomenon is repeated successively for each fracture network whose Pe is above 10.

### 3.4. Estimation of the hydraulic gradient in the massifs of New Caledonia

The previous multiscale study validates using the Pe threshold to predict the expression of the fracture network and, hence, the boulder size in a multi-fractured system. This section attempts to confront our numerical predictions with in-situ measurements of boulder sizes in New Caledonia.

First, ranges of hydraulic gradient have been estimated for four different peridotite massifs) Goro, Koniambo, Nakety and Tiebgahi) of New Caledonia (Fig. 5), based on the work of Jeanpert (2017) and data communicated by the Societe des Mines de la Tontouta (SMT) mining company. For each massif, different hydraulic gradients are given to describe their variability: an upper and lower limit and a median value (Table 2).

The hydraulic gradient chosen for the Koniambo massif corresponds to a vertical gradient calculated on a vibrating wire piezometer (Jeanpert, 2017). The extensive gradient range can be explained by the seasonal rainfall variability in New Caledonia (Moron et al., 2016). The Goro massif is also characterized by a highly varying vertical gradient between 0.04 m/m and 1 m/m (Jeanpert, 2017). Such variability is explained by the existence of highly conductive structural discontinuities. In Tiebaghi massif, the vertical hydraulic gradient is reported at 0.04 m/m, while the apparent horizontal gradient varies between 0.2 m/m and 0.005 m/m (Jeanpert, 2017). The horizontal gradients in the Nakety massif were derived from the piezometer level located along the streamline. Their variability account for the seasonal variations.

<caption>**Table 2**</caption>

Upper, median and lower values of vertical and horizontal hydraulic gradients for the Koniambo, Goro, Tiebaghi and Nakety massifs. When available, horizontal gradients are marked with *.

<table><thead><tr><th></th><th></th><th>Hydraulic gradient [m/m]</th><th>Reference</th></tr></thead><tbody><tr><td rowspan="3">Koniambo</td><td>upper</td><td>1.33</td><td rowspan="3">(Jeanpert, 2017)</td></tr><tr><td>median</td><td>0.707</td></tr><tr><td>lower</td><td>0.086</td></tr><tr><td rowspan="3">Goro</td><td>upper</td><td>1</td><td rowspan="3">(Jeanpert, 2017)</td></tr><tr><td>median</td><td>0.1</td></tr><tr><td>lower</td><td>0.04</td></tr><tr><td rowspan="3">Tiebaghi</td><td>upper</td><td>0.2 *</td><td rowspan="3">(Jeanpert, 2017)</td></tr><tr><td>median</td><td>0.02</td></tr><tr><td>lower</td><td>0.005 *</td></tr><tr><td rowspan="2">Nakety</td><td>upper</td><td>0.3 *</td><td rowspan="2">SMT</td></tr><tr><td>median</td><td>0.06 *</td></tr></tbody></table>

### 3.5. Estimation of the minimum size of boulders in the massifs of New Caledonia

Since no information on the paleo-hydraulic gradient that prevailed during the formation of laterites is available, the current vertical hydraulic gradient value is used. Indeed, in the early stages of formation, weathering is assumed to be controlled by vertical infiltration before uplift and promotion of lateral flow through increased matrix permeability (Butt and Cluzel, 2013; Dewandel et al., 2017; Jeanpert, 2017). The current hydraulic gradient probably overestimates the paleo-gradient because of the permeability increase due to weathering. Hence, the evolution of Peᵢ as a function of fracture spacing is plotted for the minimum and the maximum vertical gradient measured in all massifs to account for uncertainties. As numerical simulations on discrete fracture networks show, the fracture will enhance weathering front propagation if the corresponding Péclet number is higher than 10. For Pe between 1 and 10, we should observe an intermediate pattern between a non-delineated network and the onset of boulders linked to the competition with diffusion (Fig. 10).

Therefore, such constraints on the Peᵢ value give a minimum boulder size expected for the upper and lower bounds of the gradient. Interestingly, the higher the hydraulic gradient, the smaller the minimum boulder size (Table 3).

These results are globally coherent with the field measurements (Fig. 6) regarding the different hypotheses and simplifications made. Boulders at the sample scale range from 5 to 10 cm and are poorly delimited (Fig. 6.c). Hence, the Pe should be close to 1, roughly in agreement with boulder sizes obtained at Pe = 1 for high hydraulic gradients such as in Koniambo. Boulder sizes at Pe = 10 tend to be larger than those measured on the field at the outcrop scale, where well-delimited boulders range from 15 to 30 cm (Fig. 6.b). As a higher

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Evolution of Peᵢ as a function of the fracture spacing for the maximum and minimum vertical hydraulic gradient. The different weathering regimes are plotted.</figcaption></figure>
<table><caption>**Table 3**
Minimum boulder size constrained by the lower and upper bounds of the vertical hydraulic gradient.</caption><thead><tr><th></th><th>hydraulic gradient<br>[m/m]</th><th>Boulder size for Pe<sub>i</sub><br>= 1 [m]</th><th>Boulder size for Pe<sub>i</sub> =<br>10 [m]</th></tr></thead><tbody><tr><td>minimum</td><td>0.02</td><td>0.46</td><td>1.20</td></tr><tr><td>maximum</td><td>1.33</td><td>0.10</td><td>0.28</td></tr></tbody></table>

hydraulic gradient leads to smaller boulders, we infer that boulders observed on the field should have formed under a higher hydraulic gradient than the value used in our prediction. Furthermore, it is possible to assume local variations in the hydraulic gradient, either from the topographic effect or the medium permeability (Dewandel et al., 2017), that cause the size of boulders to vary locally.

### 3.6. Sensitivity of the current approach

If the current approach allows deriving scaling laws to predict boulder size and organization, it is highly sensitive to the quality of the estimated scaling laws. The subjective choice of the permeability and fracture spacing observation scale can strongly affect the scaling laws. Furthermore, the scaling laws only apply to a limited range of scales that should be more precisely investigated. In addition, such an approach assumes that the permeability of porous fractured medium is only driven by fractures at all the investigated scales. Previous studies also highlight the difficulties of upscaling the flow properties from the fracture scale to the network scale due to the roughness of fractures and the network topology (de Dreuzy et al., 2012).

This method is also sensitive to the estimates of the hydraulic gradient. Indeed, as discussed above, the hydraulic gradient is highly variable in New Caledonia depending on the season and the presence of highly conductive discontinuities, and it is likely to have changed through geological times since the onset of weathering.

Regarding the geochemical aspect, the chemical model used for this study is reduced to a highly simplistic system of first-order dissolution for clarity. Yet, more complex systems with cycles of dissolution-reprecipitation may change the behaviour of the weathering front progression, as well as the porosity-permeability changes of the rock as a function of mineralogical changes.

## 4. Conclusions

This study investigates the impact of hydro-geochemical features on chemical weathering front propagation in fractured bedrock. Competition between multiple fracture networks results in different dissolution patterns and variable-size boulder formation. This work focuses on geological systems characterized by high erosion rates and kinetic-limited weathering, such as supergene deposits in New Caledonia. The geometry of the computational domain relies on idealized structures of fracture networks and simplified geochemistry. Still, we use this simplification to conduct a sensitivity analysis of physical parameters governing fluid flow and weathering processes. Numerical modeling accounting for two and three-fracture networks has shed light on the influence of these fluid flow and reaction rate conditions, assessed through the dimensionless numbers Pe and Da, on dissolution patterns. Despite its simplicity, the proposed conceptual model is numerically validated and confronted against in situ observations from lateritic deposits in New Caledonia. The primary outcomes can be summarized as follows:

A fracture network will contribute to delineating the boulder shape if the associated fluid flow regime reaches a Pe threshold value of around 10. On the contrary, predominant diffusion will inhibit this fracture set in favor of a higher permeability fracture network. The threshold value is reached when the network permeability is high enough. Therefore, fracture permeability, network connectivity, and spacing are important

parameters to consider when discussing the contribution of a given fracture network (de Dreuzy et al., 2012; Long and Witherspoon, 1985) to weathering and, thus, the shape of dissolution patterns.

In contrast, the Damköhler number controls the dissolution halo width in the boulder. For $Da \approx 1$, a uniform dissolution front is observed, while for $Da > 10^4$, a thin dissolution halo occurs. Such results are consistent with those obtained in the case of wormholes (Barker and Dipple, 2019; Cabeza et al., 2020), where an optimum dissolution regime can be predicted in a Pe-Da diagram (Golfier et al., 2002). An optimum $Da$ value, $Da_{opt}$, maximizing the weathering rate of porous fractured bedrock, should also exist in our simulations and will be investigated in further studies.

The dependency of the fluid flow regime on the fracture network characteristics (fracture permeability and spacing) can be expressed using scaling laws. We thus succeeded in estimating the size of the boulders at different scales by predicting the fracture contribution as a function of their Pe. Hence, any fracture network displaying a fracture spacing above the spacing threshold corresponding to a Pe of 10 can generate boulders. Those boulders are then clustered hierarchically according to the size of their developing fracture network. As a consequence, the boulder size distribution and spatial organization are not only controlled by topological effects and competition between weathering and erosion rates (Brantley et al., 2011; Fletcher and Brantley, 2010; Lebedeva and Brantley, 2013, 2017) but also by the competition between the inherited fracture networks.

Furthermore, such a conceptual model may have implications for mineral ore targeting. Indeed, the intensity of weathering controls Ni-grade distribution in lateritic deposits (Cathelineau et al., 2016; Myagkiy et al., 2017). Hence, by better predicting the boulder size distribution and spatial organization, such an approach may help to identify the specific localization of Ni-hotspots, which are of prime interest to mining companies.

### CRediT authorship contribution statement

**Sylvain Favier:** Investigation, Methodology, Software, Writing – original draft. **Yoram Teitler:** Investigation, Resources, Writing – review & editing. **Fabrice Golfier:** Conceptualization, Software, Supervision, Validation, Writing – review & editing. **Michel Cathelineau:** Funding acquisition, Methodology, Supervision, Writing – review & editing.

### Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Michel Cathelineau reports financial support was provided by National Center for Scientific and Technological Research (New Caledonia). Michel Cathelineau reports financial support was provided by French National Research Agency. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

This work has been funded and logistically supported by the CNRT research contract CSF N° 9PS2017-CNRT.GEORESSOURCE/TRANS-NUM « Facteurs d'enrichissements et transferts de Ni, Co-Sc dans les saprolites de Nouvelle Calédonie: approche géométrique, minéralogique, 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.
## References

Anbeek, C., 1993. The effect of natural weathering on dissolution rates. Geochim. Cosmochim. Acta 57 (21-22), 4963-4975. https://doi.org/10.1016/S0016-7037 (05)80002-X.

Andrews, E.M., Hyman, J.D., Sweeney, M.R., Karra, S., Moulton, J.D., Navarre-Sitchler, A., 2023. Fracture intensity impacts on reaction front propagation and mineral weathering in three-dimensional fractured media. Water Resour. Res. 59 (2), e2022WR032121. https://doi.org/10.1029/2022WR032121.

Andrews, E., Navarre-Sitchler, A., 2021. Temporal and spatial heterogeneity of mineral dissolution rates in fractured media. Geochim. Cosmochim. Acta 312, 124–138. https://doi.org/10.1016/j.gca.2021.08.008.

Barker, S. L. L., & Dipple, G. M. (2019). Exploring for Carbonate-Hosted Ore Deposits Using Carbon and Oxygen Isotopes. In S. Decrée & L. Robb (Eds.), Geophysical Monograph Series (1st ed., pp. 185–207). Wiley. https://doi.org/10.1002/ 9781119290544.ch8.

Brantley, S.L., Buss, H., Lebedeva, M., Fletcher, R.C., Ma, L., 2011. Investigating the complex interface where bedrock transforms to regolith. Appl. Geochem. 26, S12-S15. https://doi.org/10.1016/j.apgeochem.2011.03.017.

Brantley, S.L., Lebedeva, M., 2011. Learning to read the chemistry of regolith to understand the critical zone. Annu. Rev. Earth Planet. Sci. 39 (1), 387–416. https:// doi.org/10.1146/annurev-earth-040809-152321.

Brantley, S.L., Lebedeva, M.I., 2021. Relating land surface, water table, and weathering fronts with a conceptual valve model for headwater catchments. Hydrol. Process. 35 (2), e14010. https://doi.org/10.1002/hyp.14010.

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.

Cabeza, Y., Hidalgo, J.J., Carrera, J., 2020. Competition is the underlying mechanism controlling viscous fingering and wormhole growth. Geophys. Res. Lett. 47 (3). https://doi.org/10.1029/2019GL084795.

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). Miner. Deposita 51 (2), 271–282. https://doi.org/10.1007/s00126-015-0607-y.

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

Cluzel, D., Aitchison, J.C., Picard, C., 2001. Tectonic accretion and underplating of mafic terranes in the Late Eocene intraoceanic fore-arc of New Caledonia (Southwest Pacific): geodynamic implications. Tectonophysics 340 (1-2), 23-59. https://doi. org/10.1016/S0040-1951(01)00148-2.

COMSOL Multiphysics. (2022a). Battery Design Module User's Guide. COMSOL Multiphysics v 6.1. COMSOL AB, Stockholm, Sweden, 377-378.

COMSOL Multiphysics. (2022b). Subsurface Flow Module User's Guide. COMSOL Multiphysics v 6.1. COMSOL AB, Stockholm, Sweden, 105.

Davy, P., Bour, O., De Dreuzy, J.-R., Darcel, C., 2006. Flow in multiscale fractal fracture networks. Geol. Soc. Lond. Spec. Publ. 261 (1), 31–45. https://doi.org/10.1144/ GSLSP.2006.261.01.03.

de Dreuzy, J.-R., Méheyst, Y., Pichot, G., 2012. Influence of fracture scale heterogeneity on the flow properties of three-dimensional discrete fracture networks (DFN): 3D fracture network permeability. J. Geophys. Res. Solid Earth 117 (B11). https://doi. org/10.1029/2012JB009461.

Deng, H., Peters, C.A., 2019. Reactive transport simulation of fracture channelization and transmissivity evolution. Environ. Eng. Sci. 36 (1), 90–101. https://doi.org/ 10.1089/ees.2018.0244.

Deng, H., Spycher, N., 2019. Modeling reactive transport processes in fractures. Rev. Mineral. Geochem. 85 (1), 49-74. https://doi.org/10.2138/rmg.2019.85.3.

Detwiler, R.L., Glass, R.J., Bourcier, W.L., 2003. Experimental observations of fracture dissolution: the role of Peclet number on evolving aperture variability: fracture dissolution. Geophys. Res. Lett. 30 (12). https://doi.org/10.1029/2003GL017396.

Detwiler, R.L., Rajaram, H., 2007. Predicting dissolution patterns in variable aperture fractures: evaluation of an enhanced depth-averaged computational model: predicting dissolution patterns. Water Resour. Res. 43 (4). https://doi.org/10.1029/ 2006WR005147.

Dewandel, B., Jeanpert, J., Ladouche, B., Join, J.-L., Maréchal, J.-C., 2017. Inferring the heterogeneity, transmissivity and hydraulic conductivity of crystalline aquifers from a detailed water-table map. J. Hydrol. 550, 118–129. https://doi.org/10.1016/j. jhydrol.2017.03.075.

Dijk, P.E., Berkowitz, B., Yechieli, Y., 2002. Measurement and analysis of dissolution patterns in rock fractures: measurement and analysis of dissolution in rock fractures. Water Resour. Res. 38 (2), 5-1-5-12. https://doi.org/10.1029/2001WR000246.

Dixon, J.L., Hartshorn, A.S., Heimsath, A.M., DiBiase, R.A., Whipple, K.X., 2012. Chemical weathering response to tectonic forcing: A soils perspective from the San Gabriel Mountains, California. Earth Planet. Sci. Lett. 323-324, 40-49. https://doi. org/10.1016/j.epsl.2012.01.010.

Dublet, G., Juillot, F., Morin, G., Fritsch, E., Fandeur, D., Ona-Nguema, G., Brown, G.E., 2012. Ni speciation in a New Caledonian lateritic regolith: a quantitative X-ray absorption spectroscopy investigation. Geochim. Cosmochim. Acta 95, 119–133. https://doi.org/10.1016/j.gca.2012.07.030.

Ehlen, J., 2000. Fractal analysis of joint patterns in granite. Int. J. Rock Mech. Min. Sci. 37 (6), 909-922. https://doi.org/10.1016/S1365-1609(00)00027-7.

Ehlen, J., 2002. Some effects of weathering on joints in granitic rocks. Catena 49 (1-2), 91-109. https://doi.org/10.1016/S0341-8162(02)00019-X.

Elkhoury, J.E., Ameli, P., Detwiler, R.L., 2013. Dissolution and deformation in fractured carbonates caused by flow of CO2-rich brine under reservoir conditions. Int. J. Greenhouse Gas Control 16, S203-S215. https://doi.org/10.1016/j. ijggc.2013.02.023.

Emmanuel, S., 2022. Modeling the effect of mineral armoring on the rates of coupled dissolution-precipitation reactions: Implications for chemical weathering. Chem. Geol. 601, 120868. https://doi.org/10.1016/j.chemgeo.2022.120868.

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 Geol. Rev. 147, 104971. https://doi.org/10.1016/j.oregorev.2022.104971.

Favier, S., Teitler, Y., Golfier, F., Cathelineau, M., 2024. Reactive transport modeling in porous fractured media: application to weathering processes. Water Resour. Res. https://doi.org/10.1029/2023WR035283.

Fletcher, R.C., Brantley, S.L., 2010. Reduction of bedrock blocks as corestones in the weathering profile: Observations and model. Am. J. Sci. 310 (3), 131-164. https:// doi.org/10.2475/03.2010.01.

Fletcher, R., Buss, H., Brantley, S., 2006. A spheroidal weathering model coupling porewater chemistry to soil thicknesses during steady-state denudation. Earth Planet. Sci. Lett. 244 (1-2), 444-457. https://doi.org/10.1016/j.epsl.2006.01.055.

Freyssinet, Ph., 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, 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. Chem. Geol. 170 (1-4), 133-151. https://doi.org/ 10.1016/S0009-2541(99)00245-4.

Fumagalli, A., Keilegavlen, E., 2019. Dual virtual element methods for discrete fracture matrix models. Oil Gas Sci. Technol. 74, 41. https://doi.org/10.2516/ogst/2019008.

Genthon, P., Join, J.-L., Jeanpert, J., 2017. Differential weathering in ultramafic rocks of New Caledonia: the role of infiltration instability. J. Hydrol. 550, 268-278. https:// doi.org/10.1016/j.jhydrol.2017.04.059.

Genthon, P., Ormond, A., 2008. Infiltrative instability near topography with implication for the drainage of soluble rocks. Hydrol. Earth Syst. Sci. 12 (6), 1285-1293. https:// doi.org/10.5194/hess-12-1285-2008.

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., Zarcone, C., Bazin, B., Lenormand, R., Lasseux, D., Quintard, M., 2002. On the ability of a Darcy-scale model to capture wormhole formation during the dissolution of a porous medium. J. Fluid Mech. 457, 213-254. https://doi.org/10.1017/ S0022112002007735.

Golightly, J. P. (2010). Progress in Understanding the Evolution of Nickel Laterites. In R. J. Goldfarb, E. E. Marsh, & T. Monecke, 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.

Gong, J., Rossen, W.R., 2018. Characteristic fracture spacing in primary and secondary recovery for naturally fractured reservoirs. Fuel 223, 470-485. https://doi.org/ 10.1016/j.fuel.2018.02.046.

Goodfellow, B.W., Hilley, G.E., 2022. Climatic and lithological controls on the structure and thickness of granitic weathering zones. Earth Planet. Sci. Lett. 600, 117890. https://doi.org/10.1016/j.epsl.2022.117890.

Guéguen, Y., Gavrilenko, P., Le Ravalec, M., 1996. Scales of rock permeability. Surv. Geophys. 17 (3), 245-263. https://doi.org/10.1007/BF01904043.

Haldar, S.K., 2017. Deposits of Europe. In Platinum-Nickel-Chromium Deposits. Elsevier, pp. 217-231.

Hausrath, E.M., Navarre-Sitchler, A.K., Sak, P.B., Williams, J.Z., Brantley, S.L., 2011. Soil profiles as indicators of mineral weathering rates and organic interactions for a Pennsylvania diabase. Chem. Geol. 290 (3-4), 89-100. https://doi.org/10.1016/j. chemgeo.2011.08.014.

Hirata, Y., Chigira, M., Chen, Y., 2017. Spheroidal weathering of granite porphyry with well-developed columnar joints by oxidation, iron precipitation, and rindlet exfoliation: spheroidal weathering of granite porphyry. Earth Surf. Proc. Land. 42 (4), 657-669. https://doi.org/10.1002/esp.4008.

Huymans, M., Dassargues, A., 2005. Review of the use of Péclet numbers to determine the relative importance of advection and diffusion in low permeability environments. Hydrol. Geol. J. 13 (5-6), 895-904. https://doi.org/10.1007/s10040- 004-0387-4.

Hyman, J.D., Sweeney, M.R., Gable, C.W., Svyatsky, D., Lipnikov, K., David Moulton, J., 2022. Flow and transport in three-dimensional discrete fracture fracture matrix models using mimetic finite difference on a conforming multi-dimensional mesh. J. Comput. Phys. 466, 111396. https://doi.org/10.1016/j.jcp.2022.111396.

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. Econ. Geol. 113 (2), 531-544. https://doi.org/10.5382/ econgeo.2018.4561.

Jeanpert, J., Iseppi, M., Adler, P.M., Genthon, P., Sevin, B., Thovert, J.-F., Dewandel, B., Join, J.-L., 2019. Fracture controlled permeability of ultramafic basement aquifers. Inferences from the Koniambo massif New Caledonia. Eng. Geol. 256, 67-83. https://doi.org/10.1016/j.enggeo.2019.05.006.

Jeanpert, J. (2017). Structural and hydrogeologic characterization of the peridotitic massifs of New Caledonia (Issue 2017LARE0047) [Thesis, Université de la Réunion]. https:// tel.archives-ouvertes.fr/tel-01845503.

Join, J.-L., Jeanpert, J., Dewandel, B., Genthon, P., Bouichet, G., Sevin, B., Maurizot, P., Maréchal, J.-C., Ladouche, B., Adler, P. M., Morgenstern, U., & Taupin, J. D. (2017). Hydrogéologie minière—HYPERK: Hydrogéologie des Péridotites karstifiées de NouvelleCalédonie. Rapport final. Projet CNRT 'Nickel et son environnement' (p. 56).
Jung, H., Navarre-Sitchler, A., 2018a. Physical heterogeneity control on effective mineral dissolution rates. Geochim. Cosmochim. Acta 227, 246–263. https://doi.org/10.1016/j.gca.2018.02.028.

Jung, H., Navarre-Sitchler, A., 2018b. Scale effect on the time dependence of mineral dissolution rates in physically heterogeneous porous media. Geochim. Cosmochim. Acta 234, 70–83. https://doi.org/10.1016/j.gca.2018.05.009.

Ladd, A.J.C., Szymczak, P., 2021. Reactive flows in porous media: challenges in theoretical and numerical methods. Ann. Rev. Chem. Biomol. Eng. 12 (1), 543–571. https://doi.org/10.1146/annurev-chembioeng-092920-102703.

Lebedeva, M.I., Brantley, S.L., 2013. Exploring geochemical controls on weathering and erosion of convex hillslopes: Beyond the empirical regolith production function: hillslope evolution and regolith thickness. Earth Surf. Proc. Land. 38 (15), 1793–1807. https://doi.org/10.1002/esp.3424.

Lebedeva, M.I., Brantley, S.L., 2017. Weathering and erosion of fractured bedrock systems: weathering and erosion of fractured bedrock systems. Earth Surf. Proc. Land. 42 (13), 2090–2108. https://doi.org/10.1002/esp.4177.

Lebedeva, M.I., Brantley, S.L., 2020. Relating the depth of the water table to the depth of weathering. Earth Surf. Proc. Land. 45 (9), 2167–2178. https://doi.org/10.1002/esp.4873.

Lebedeva, M.I., Brantley, S.L., 2023. Using homogenized models to explore the effect of fracture densities on weathering. Am. J. Sci. 323. https://doi.org/10.2475/001c.68308.

Li, D.D., Jacobson, A.D., McInerney, D.J., 2014. A reactive-transport model for examining tectonic and climatic controls on chemical weathering and atmospheric CO2 consumption in granitic regolith. Chem. Geol. 365, 30–42. https://doi.org/10.1016/j.chemgeo.2013.11.028.

Long, J.C.S., Witherspoon, P.A., 1985. The relationship of the degree of interconnection to permeability in fracture networks. J. Geophys. Res. 90 (B4), 3087. https://doi.org/10.1029/JB090Ib04p03087.

Ma, G., Rezania, M., Nezhad, M.M., 2022a. Effects of spatial autocorrelation structure for friction angle on the runout distance in heterogeneous sand collapse. Transp. Geotech. 33, 100705. https://doi.org/10.1016/j.trgeo.2021.100705.

Ma, G., Rezania, M., Nezhad, M.M., 2022b. Stochastic assessment of landslide influence zone by material point method and generalized geotechnical random field theory. Int. J. Geomech. 22 (4), 04022002. https://doi.org/10.1061/(ASCE)GM.1943-5622.0002308.

MacQuarrie, K.T.B., Mayer, K.U., 2005. Reactive transport modeling in fractured rock: a state-of-the-science review. Earth Sci. Rev. 72 (3–4), 189–227. https://doi.org/10.1016/j.earscirev.2005.07.003.

Maher, K., 2010. The dependence of chemical weathering rates on fluid residence time. Earth Planet. Sci. Lett. 294 (1–2), 101–110. https://doi.org/10.1016/j.epsl.2010.03.010.

Mao, H., Cui, L., Zhang, Z., Xu, S., Liu, C., Zhao, Z., 2022. Influence of monsoon climate on chemical weathering of granitic regoliths. Global Biogeochem. Cycles 36 (5), e2022GB007362. https://doi.org/10.1029/2022GB007362.

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. *Geological Society, London, Memoirs*, 51(1), 247–277. 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. J. Teknologi 72 (3). https://doi.org/10.11113/jt.v72.4015.

Md Dan, M.F., Mohamad, E.T., Komoo, I., 2016. Characteristics of boulders formed in tropical weathered granite: a review. J. Teknologi 78 (8–6). https://doi.org/10.11113/jt.v78.9635.

Moron, V., Barbero, R., Robertson, A.W., 2016. Subseasonal-to-interannual variability of rainfall over New Caledonia (SW Pacific). Clim. Dyn. 46 (7–8), 2449–2468. https://doi.org/10.1007/s00382-015-2712-0.

Munday, T., Green, A., Meyers, J., Macnae, J., Rutherford, J., Cooper, M., 2003. The cawse nickeliferous laterite deposits, western australia – a case study on the application of airborne geophysics in targeting zones of supergene enrichment in a complex regolith setting. ASEG Extended Abstracts 2003 (2), 1–5. https://doi.org/10.1071/ASEG2003ab113.

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. Chem. Geol. 466, 274–284. https://doi.org/10.1016/j.chemgeo.2017.06.018.

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. Geochim. Geophys. Geosyst. 20 (3), 1425–1440. https://doi.org/10.1029/2018GC007606.

Navarre-Sitchler, A., Brantley, S., 2007. Basalt weathering across scales. Earth Planet. Sci. Lett. 261 (1–2), 321–334. https://doi.org/10.1016/j.epsl.2007.07.010.

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

Navarre-Sitchler, A., Thyne, G., 2007. Effects of carbon dioxide on mineral weathering rates at earth surface conditions. Chem. Geol. 243 (1–2), 53–63. https://doi.org/10.1016/j.chemgeo.2007.06.009.

Odšačer, L.H., Kvamsdal, T., Larson, M.G., 2019. A simple embedded discrete fracture-matrix model for a coupled flow and transport problem in porous media. Comput. Methods Appl. Mech. Eng. 343, 572–601. https://doi.org/10.1016/j.cma.2018.09.003.

Oelkers, E.H., Declercq, J., Saldi, G.D., Gislason, S.R., Schott, J., 2018. Olivine dissolution rates: a critical review. Chem. Geol. 500, 1–19. https://doi.org/10.1016/j.chemgeo.2018.10.008.

Ollier, C.D., 2010. Very deep weathering and related landslides. Geol. Soc. Lond. Eng. Geol. Spl. Publ. 23 (1), 5–14. https://doi.org/10.1144/EGSP23.2.

Oltéan, C., Golfier, F., Bués, M.A., 2013. Numerical and experimental investigation of buoyancy-driven dissolution in vertical fracture: buoyancy-driven dissolution in fracture. J. Geophys. Res. Solid Earth 118 (5), 2038–2048. https://doi.org/10.1002/jgrb.50188.

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 Resour. Res. 52 (12), 9344–9366. https://doi.org/10.1002/2016WR019026.

Panga, M.K.R., Ziauddin, M., Balakotaiah, V., 2005. Two-scale continuum model for simulation of wormholes in carbonate acidization. AICHE J. 51 (12), 3231–3248. https://doi.org/10.1002/aic.10574.

Patino, L.C., Velbel, M.A., Price, J.R., Wade, J.A., 2003. Trace element mobility during spheroidal weathering of basalts and andesites in Hawaii and Guatemala. Chem. Geol. 202 (3–4), 343–364. https://doi.org/10.1016/j.chemgeo.2003.01.002.

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

Pyrak-Nolte, L.J., Nolte, D.D., 2016. Approaching a universal scaling relationship between fracture stiffness and fluid flow. Nat. Commun. 7 (1), 10663. https://doi.org/10.1038/ncomms10663.

Quesnel, B., Gautier, P., Cathelineau, M., Boulvais, P., Couteau, C., Drouillet, M., 2016. The internal deformation of the Peridotite Nappe of New Caledonia: a structural study of serpentine-bearing faults and shear zones in the Koniambo Massif. J. Struct. Geol. 85, 51–67. https://doi.org/10.1016/j.jsg.2016.02.006.

Riebe, C.S., Hahm, W.J., Brantley, S.L., 2017. Controls on deep critical zone architecture: a historical review and four testable hypotheses. Earth Surf. Proc. Land. 42 (1), 128–156. https://doi.org/10.1002/esp.4052.

Rayne, A., Jamtveit, B., Mathiesen, J., Malthe-Sørensen, A., 2008. Controls on rock weathering rates by reaction-induced hierarchical fracturing. Earth Planet. Sci. Lett. 275 (3–4), 364–369. https://doi.org/10.1016/j.epsl.2008.08.035.

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 Geol. Rev. 140, 104501. https://doi.org/10.1016/j.oregeorev.2021.104501.

Sarracino, R.S., Prasad, G., Hoohlo, M., 1987. A mathematical model of spheroidal weathering. Math. Geol. 19 (4), 269–289. https://doi.org/10.1007/BF00897839.

Schoonejans, J., Vanacker, V., Opfergelt, S., Ameijeras-Marino, Y., Christl, M., 2016. Kinetically limited weathering at low denudation rates in semiarid climatic conditions. J. Geophys. Res. Earth 121 (2), 336–350. https://doi.org/10.1002/2015JP003626.

Sklar, L.S., Riebe, C.S., Marshall, J.A., Genetti, J., Leclere, S., Lukens, C.L., Merces, V., 2017. The problem of predicting the size distribution of sediment supplied by hillslopes to rivers. Geomorphology 277, 31–49. https://doi.org/10.1016/j.geomorph.2016.05.005.

Snow, D.T., 1969. Anisotropie permeability of fractured media. Water Resour. Res. 5 (6), 1273–1289. https://doi.org/10.1029/WR005i006p01273.

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

Szymczak, P., Ladd, A.J.C., 2006. A network model of channel competition in fracture dissolution. Geophys. Res. Lett. 33 (5), L05401. https://doi.org/10.1029/2005GL025334.

Szymczak, P., Ladd, A.J.C., 2009. Wormhole formation in dissolving fractures. J. Geophys. Res. Solid Earth 114 (B6), 2008JB006122. https://doi.org/10.1029/2008JB006122.

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. J. Geochem. Explor. 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.

Trescales, J., 1975. L'élution géochimique supégene des roches ultrabasiques en zone tropicale: formation des gisements nickeléfè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 Resour. Res. 57 (10). https://doi.org/10.1029/2021WR030442.

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

White, A.F., Brantley, S.L., 2003. The effect of time on the weathering of silicate minerals: why do weathering rates differ in the laboratory and field? Chem. Geol. 202 (3–4), 479–506. https://doi.org/10.1016/j.chemgeo.2003.03.001.

Wogelius, R.A., Walther, J.V., 1992. Olivine dissolution kinetics at near-surface conditions. Chem. Geol. 97 (1–2), 101–112. https://doi.org/10.1016/0009-2541(92)00138-U.

Yuan, Wei, Zhang, Qin, 2019. A numerical simulator for modeling the coupling processes of subsurface fluid flow and reactive transport processes in fractured carbonate rocks. Water 11 (10), 1957. https://doi.org/10.3390/w11101957.

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 Resour. Res. 58 (2). https://doi.org/10.1029/2021WR030578.