D Attié

M Batkiewicz-Kwasniak

P Billoir

A Blanchet

A Blondel

S Bolognesi

D Calvet

M G Catanesi

M Cicerchia

G Cogo

P Colas

G Collazuol

A Delbart

J Dumarchez

S Emery-Schrenk

M Feltre

C Giganti

F Gramegna

M Grassi

M Guigue

P Hamacher-Baumann

S Hassani
email: samira.hassani@cea.fr

F Iacob

C Jesús-Valls

R Kurjata

M Lamoureux

M Lehuraux

A Longhin

T Lux

L Magaletti

T Marchi

A Maurel

L Mellet

M Mezzetto

L Munteanu

Q V Nguyen

Y Orain

M Pari

J.-M Parraud

C Pastore

A Pepato

E Pierre

B Popov

H Przybiliski

T Radermacher

E Radicioni

M Riallot

S Roth

A Rychter

L Scomparin

J Steinmann

S Suvorov
email: sergey.suvorov@lpnhe.in2p3.fr

J Swierblewski

D Terront

N Thamm

F Toussenel

V Valentino

G Vasseur

U Yevarouskaya

M Ziembicki

M Zito

Characterization of resistive Micromegas detectors for the upgrade of the T2K Near Detector Time Projection Chambers

Keywords:

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.

Introduction

T2K [START_REF] Abe | The T2K Experiment[END_REF] is a long-baseline neutrino oscillation experiment exploiting a muon neutrino beam produced by the JPARC accelerator complex in Japan.

The T2K experiment includes the secondary neutrino beamline, a set of near detectors (INGRID and ND280) and the far detector SuperKamiokande.

T2K provided the first evidence of non-zero mixing angle θ 13 [START_REF] Abe | Indication of Electron Neutrino Appearance from an Accelerator-produced Off-axis Muon Neutrino Beam[END_REF] and discovered the appearance of electron neutrinos in a muon neutrino beam [START_REF] Abe | Indication of Electron Neutrino Appearance from an Accelerator-produced Off-axis Muon Neutrino Beam[END_REF][START_REF] Abe | Evidence of Electron Neutrino Appearance in a Muon Neutrino Beam[END_REF][START_REF] Abe | Observation of Electron Neutrino Appearance in a Muon Neutrino Beam[END_REF]. Combining T2K data with precise θ 13 measurement from reactor experiments, T2K has recently reported hints of large Charge-Parity (CP) violation in the leptonic sector [START_REF] Abe | Constraint on the matter-antimatter symmetry-violating phase in neutrino oscillations[END_REF], excluding CP-conservation at about the 2σ Confidence Level.

The T2K collaboration is now preparing for the second phase of the experiment (T2K-II), starting in Fall 2022, which will exploit the upgrade of the beam from 500 kW to 750 kW. T2K-II will collect in total more than 10 22

Protons-On-Target (POT), including the 3.6 × 10 21 POT already collected, thus enabling 3σ sensitivity on CP-Violation, in case of maximally violated CP, as currently indicated by the T2K results [START_REF] Abe | Constraint on the matter-antimatter symmetry-violating phase in neutrino oscillations[END_REF]. In order to cope with such increased statistics, improved control of the relevant systematic uncertainties is needed. To this aim, an upgrade to ND280 is being constructed. ND280 is a magnetized multi-purpose detector that measures, before the oscillation, the neutrino differential rate versus the kinematics variables of the particles, mostly muons, that result from the neutrino interaction with nuclei in dense parts of the detector. This allows one to constrain the neutrino flux and the neutrino-nucleus interaction cross-section. The present ND280 consists of two main parts: an upstream π 0 detector (P0D) and a downstream tracker which includes two Fine Grained Scintillators (FGD) interleaved with three vertical Time Projection Chambers (TPCs). The P0D and the tracker are surrounded by an electromagnetic calorimeter (ECAL) and by the UA1 magnet providing a 0.2 T magnetic field. The magnet yoke is further instrumented with a Side Muon Range Detector (SMRD).

The role of ND280 in the T2K oscillation analysis is crucial, allowing one to constrain the uncertainty on the expected number of events at the far detector to 4-5%. ND280 measurements are performed on different targets (Carbon and Oxygen) and rely on the precise measurement of the muon momentum measured by the TPCs, with a momentum resolution of 10% at 1 GeV [START_REF] Abgrall | Time Projection Chambers for the T2K Near Detectors[END_REF]. Better resolution is not needed since the determination of the incoming neutrino energy from the outgoing muon momentum in ν µ charged-current interactions are limited by the uncertainty induced by the Fermi motion of the nucleon in the nucleus.

Another important result of ND280 is the measurement of the ν e contamination in the beam [START_REF] Abe | Measurement of the intrinsic electron neutrino component in the T2K neutrino beam with the ND280 detector[END_REF][START_REF] Abe | Measurement of the charged-current electron (anti-)neutrino inclusive cross-sections at the T2K off-axis near detector ND280[END_REF], which constitutes one of the two main backgrounds to the ν e appearance in the electron-like sample selected at Super-Kamiokande. This analysis, based on the Particle Identification (PID) capabilities of the TPCs and of ECAL, is possible thanks to the deposited energy resolution of 8% in the TPCs that allows sufficient e-µ separation between a few hundred MeV and ∼2 GeV.

The upgrade of ND280 consists in substituting the P0D with a new tracker, similar to the existent one but with a horizontal orientation parallel to the neutrino beam, and able to detect particles transverse to the beam, for which the acceptance of the current tracker is small. The new tracker includes a 3-dimensional scintillator target for neutrino interactions(Super-FGD) made of about 2 million 1 cm 3 cubes, read out by wavelength shifting fibers in the 3 directions [START_REF] Blondel | The SuperFGD prototype charged particle beam tests[END_REF]. On the top and the bottom of the Super-FGD, two High-Angle TPCs (HA-TPCs) will be installed. The new tracker system will be surrounded by six Time-of-Flight modules. This new detector configuration will allow one to improve the angular acceptance of ND280, being close to the full 4π phase-space accessible at SuperKamiokande. In addition, the better tracking performances of the super-FGD will allow improving the reconstruction of the hadronic part of the neutrino interactions, which will be exploited in combination with the muon kinematics. Also, more target mass for neutrino interactions yields more data statistics.

Each endplate of the new TPCs will be instrumented with eight Micromegas type charge readout modules. The Micromegas detector charge amplification structure uses a resistive anode to "spread" the charges over several pads to improve the space point resolution: a signal is induced in pads adjacent to the pad receiving directly electrons from the avalanche, and the information from those additional signals allows one to better measure the track position. This effect is different from the charge sharing between pads resulting from the diffusion of the primary electrons in the TPC gas. This innovative technique is combined with the bulk-Micromegas technology [START_REF] Abgrall | Time Projection Chambers for the T2K Near Detectors[END_REF] to compose the "Encapsulated Resistive Anode Micromegas" (ERAM) detector. Performances of a prototype of an ERAM detector exposed to a test beam at CERN were shown in [START_REF] Attié | Performances of a resistive Micromegas module for the Time Projection Chambers of the T2K Near Detector upgrade[END_REF]. This paper describes the performance of one prototype ERAM module, built with the same design that will be used for the HA-TPCs, exposed to an electron beam at DESY. With respect to our previous studies reported in [START_REF] Attié | Performances of a resistive Micromegas module for the Time Projection Chambers of the T2K Near Detector upgrade[END_REF], in this paper we characterize both spatial and dE/dx resolution as a function of the angle of the track projection reconstructed in the ERAM plane, with the side of the ERAM plane, corresponding to the sides of the ERAM pads.

The HA-TPC should match the performance of existing TPCs in terms of momentum and dE/dx resolution. In addition, the super-FGD will enable the reconstruction of low momentum protons and neutrons [START_REF] Abe | T2K ND[END_REF], thus allowing one to reconstruct the incoming neutrino energy more precisely and effectively correct for the Fermi motion: an improved track momentum resolution in the HA-TPC, even beyond the previous specifications, is therefore useful to improve the neutrino energy resolution and, ultimately, the precision on neutrino oscillation measurements.

Resistive Micromegas

The ERAM detector uses the bulk-Micromegas technology invented in 2004 by a collaboration between the CERN/EP-DT-EF PCB workshop and CEA-IRFU [START_REF] Giomataris | Micromegas in a bulk[END_REF]. It was developed and used for the construction of the 72 bulk-Micromegas modules (9 m 2 total area) which equip the three T2K/ND280 vertical TPCs. The woven micromesh of the Micromegas detector is sandwiched between two layers of the same insulating material (pyralux) and exposed to UV radiation at the location where the pyralux is kept on top of the pad-segmented anode PCB after chemical development.

The ERAM detector is a 128 µm amplification gap bulk-Micromegas using the standard SD45/18 304L woven micromesh built on top of a resistive anode PCB. When a charged particle crosses the gas volume of a TPC, it generates a cloud of ionized electrons. These electrons are drifted to the anode readout plane of the TPC under a uniform electric field. On the readout plane, an avalanche is generated by a high electric field in the Micromegas amplification region. The voltage that induced this field is notated as "DLC voltage" in the following of this paper. The DLC layer is powered at High Voltage on each of the 4 sides through the DLC-covered kapton foil. The 3D position of the track is then reconstructed from the arrival time and the position of the avalanche on the readout plane. In the case of the bulk-Micromegas, the electron avalanche in the amplification gap is much smaller comparing to the pad size (≈ 1 cm) and therefore the position resolution is often limited by the pad size. The resistive anode Micromegas, first introduced and extensively studied by the ILC-TPC collaboration [START_REF] Dixit | Position sensing from charge dispersion in micropattern gas detectors with a resistive anode[END_REF], provides a way to induce a signal on a larger number of pads allowing for a better reconstruction of the track position. A sketch of the bulk and the resistive Micromegas concepts are presented in Figure 1. The anode, segmented in pads, is covered by a foil of insulating material, which has a thin resistive layer on its top side. The ERAM detector uses a Diamond-Like Carbon (DLC) thin layer sputtered on a 50 µm thick APICAL (Kapton) insulator sheet. The avalanche is then naturally quenched because the potential difference locally drops in presence of a high charge density.

The resistive layer effect can be approximated with a 2-D RC network. The charge deposited by the avalanche induces a signal also on the adjacent pads.

Within the RC model the charge spread is following a Gaussian distribution.

For a point charge deposited at r = 0 and t = 0, the charge density as a function of radius r and time t is given by:

ρ(r, t) = RC 4πt e -r 2 RC 4t (1)
where R is the resistivity per unit area and C is the capacitance per unit area. For this structure, the capacitance C is defined by the distance between the resistive layer and the grounded pads (glue thickness plus APICAL foil).

These characteristics determines the width of the induced charge spread σ t at a given time. When measured, due to integration of the charge collected by a front-end electronics of shaping time t, this spread is of the order of

σ t = 2t RC .
The resistive anode provides mainly two advantages: by spreading the charge between neighbouring pads, it improves greatly the resolution with respect to the pitch/ √ 12 provided by a mere hodoscope, and it suppresses the formation of sparks and limits their intensity. The first one is the main advantage in the context of the T2K experiment. A further and novel improvement of this technique is a new High Voltage powering scheme, where the mesh is set to ground and the anode to a positive amplification voltage. The insulation of the resistive anode from the pads, hence from the electronics, ensures a safe operation by a capacitive coupling readout and thus allows us to remove the cumbersome anti-spark protection circuitry of electronic read-out boards necessary in the case of the standard bulk readout.

At the end of 2017, the first series of prototypes were produced and tested to assess the feasibility of large charge spreading by a low resistivity anode. The vertical TPCs PCB, with an active area of 36×34 cm 2 covered by 0.97×0.69 cm 2 pads, was adapted to build an ERAM structure with "onshelf" 2.5 MOhm/ DLC foils. This first prototype was tested with a particle beam at CERN in 2018 and its performances are summarized in [START_REF] Attié | Performances of a resistive Micromegas module for the Time Projection Chambers of the T2K Near Detector upgrade[END_REF].

In fall 2018, the global design of the ND280 upgrade detector was fixed.

The sub-detector envelopes were defined and the size of the ERAM modules fixed to be 420 × 340 mm channels AFTER ASIC per card). The AFTER chip, already used for the existing TPCs of T2K and for the test beam described in [START_REF] Attié | Performances of a resistive Micromegas module for the Time Projection Chambers of the T2K Near Detector upgrade[END_REF] allows tuning several parameters such as the gain, the shaping time and the sampling time.

This prototype was tested in a test beam at DESY in June 2019.

Experimental setup

The prototype has been exposed to an electron beam at the DESY II test beam facility [START_REF] Diener | The DESY II Test Beam Facility[END_REF]. DESY II provides electron beams of 1 -6 GeV/c at a rate of up to several kHz, depending on the chosen momentum. In the test beam area TB24/1, a large-bore superconducting solenoid, called PCMAG, provides a magnetic field of up to 1.25 T. The magnet is mounted on a movable platform, which allows the setup to be moved horizontally and vertically, perpendicular to the beam line, as well as rotated by ±45°in the horizontal plane. The platform can position the device under test with a precision of about 0.2 mm horizontally, 0.1 mm vertically, and within 0.1°in angle.

Inside the bore of the magnet, a rail system is installed on which test devices can be mounted at different positions within the magnet. The TPC prototype is supported on a sled, which can move in and out of the magnet and can be used to rotate the chamber around the magnetic field axis as illustrated in Figure 2.

Usually, the magnet is positioned perpendicular to the beam, so the magnetic field is also perpendicular to the beam, in a horizontal plane. The ERAM module is a vertical plane parallel to the beam, but perpendicular to the magnetic field. The walls of the magnet present about 20 % of a radiation length, so that an electron beam of 6 GeV/c easily penetrates the magnet and the device under test. A set of four consecutive scintillation counters, of which each has an area of approximately 2.5 cm × 2.5 cm, is mounted about 1.5 m in front of the magnet. The coincidence between them is used as a beam trigger. In addition, a second set of scintillation counters above and below the magnet provides a cosmic trigger for tests without beam.

The prototype has been placed inside the magnetic field of 0.2 T provided by the PCMAG magnet and has been operated for these measurements with a gas mixture of 95% argon, 3% tetrafluoromethane (CF 4 ), and 2% isobutane (iC 4 H 10 ) that is the same gas mixture as the one used in the existing TPCs of ND280. The gas quality was constantly monitored during the measurement. For the results reported in this paper, the oxygen contamination was around 30 ppm at a gas flow rate of 30 /h. The chamber was operated at atmospheric pressure. Ambient temperature and pressure were constantly monitored. The results presented in this paper were obtained with electrons with momenta varying from 0.5 to 5 GeV/c. We have carried out drift distance scans with seven points, spaced by 2 cm, at B = 0 T and B = 0.2 T, peaking time of 200 ns or 412 ns, and DLC voltage of 370 V and 360 V. In addition a scan of the DLC voltage, varied from 330 V to 400 V, was performed.

We have calibrated the Micromegas position with respect to the magnet frame coordinate system, the moving table position, and the drift velocities at different drift fields with a 4 GeV/c electron beam and a short peaking time of 116 ns. At the standard T2K drift field of 275 V/cm, we obtained a drift velocity of 7.68±0.03 cm/µs by drift distance scans with the accelerator beam. A Gas Monitoring Chamber (GMC), identical to the ones deployed at T2K's ND280 detector [START_REF] Abgrall | Time Projection Chambers for the T2K Near Detectors[END_REF], monitored the exhaust gas for the duration of the beam-time. The GMC measured a drift velocity of 7.81 ± 0.02 cm/µs at the T2K field, in agreement with the beam scans, see Figure 3. Under the electric field of 140 V/cm, which is associated with the minimum transverse diffusion, we found a drift velocity of 5.84 cm/µs. Despite the impact of the gas bottle changes on the gas properties, the drift velocity V drift under the electric field of E = 275V/cm varies less than 6‰.

Charge spreading characterisation

As explained before, the resistive Micromegas technology produces a spreading of the collected charge into neighbouring pads. The charge spreading phenomenon, which drives the waveform shape in the pads, is described in [START_REF] Dixit | Position sensing from charge dispersion in micropattern gas detectors with a resistive anode[END_REF]. The signal in the adjacent pads is smaller in amplitude, delayed and longer in time compared to the signal in the leading pad where the avalanche is located. The time delay of the waveform increases with the distance of the adjacent pad to the track. To study the phenomenon of interest, we define a "cluster" as a group of pads in the perpendicular direction to the track. A schematic view of such a cluster for horizontal tracks and waveforms in the leading and in adjacent pads are shown in Figure 4.

In the current work, the maximum of the waveform is used as an estimator of the charge in each pad, which corresponds to different times for the leading and adjacent pads. In Figure 5, the pad multiplicity per cluster and the ratio of the charge measured in the pad with the largest signal with respect to the sum of the charges of the pads in the cluster (q max /q cluster ) are shown. The other possible approaches e.g. fit of the waveform with an analytical function are under study.

Most of the clusters are formed by more than two pads and the pad with the largest signal contains typically 80% of the total charge of the cluster. The effect of the high voltage on the multiplicity is also clearly seen in Figure 5.

The cluster multiplicity increases with the high voltage, as the probability for the smaller signals in some pads to pass the threshold increases with the gain.

For tracks projections in the pad plane which are parallel to the sides of the pads, transverse spreading is defined precisely within the given column.

While for inclined tracks and large square pads, the separation between longitudinal (along the track) and transverse spreading is more complicated. In order to distinguish these two spreading topologies, we define more sophisticated cluster patterns shown in Figure 6. These clusters are repeated to pave the whole ERAM. In the case of square pads, these patterns are optimal for angles with tangents 0 (column), 1 (diagonal), 0.5 (2 by 1), 0.

Pattern recognition

In the test beam analysis, we focus on the studies of tracks crossing the whole detector as more complicated typologies (e.g. showering, curved low-energy tracks) are difficult to interpret. Hence simple reconstruction algorithms like DBSCAN [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF] are sufficient in our case.

We select a track if it is crossing the whole detector without breaks or splits. A split is defined as the case where there is more than one cluster in a given column. An event containing a split is thus a multiple track candidate and is rejected in our analysis. However, with the resistive spreading, two close parallel tracks may not be separated by a gap and thus misreconstructed as one single track. To reject such a topology, a cut on the pad multiplicity in each cluster was implemented. The cut value was optimised for each Micromegas voltage and electronics shaping time. Examples of the accepted and rejected events are shown in Figure 7 (a) and (b) respectively.

Spatial resolution

As described in section 1, the momentum resolution has to be better than 10% at 1 GeV/c but, to fully exploit the ND280 upgrade capabilities, even better momentum resolution would be desirable. The momentum resolution is directly connected to the spatial resolution obtained for one cluster through the Gluckstern formula [START_REF] Gluckstern | Uncertainties in track momentum and direction, due to multiple scattering and measurement errors[END_REF]. For tracks with 70 point measurements (using clusters), a maximum drift distance of 1 m, and a magnetic field of 0.2 T, a spatial resolution of ∼ 800µm would be sufficient to reach a momentum resolution of 10% at 1 GeV/c.

As we will show in this section, the resistive Micromegas technology allows one to significantly improve the spatial resolution with respect to the bulk Micromegas, even in presence of slightly larger pads thus allowing one to reduce by ∼ 30% the total number of electronic channels for the same active surface. The test beam data have been used to characterize the ERAM module performances for electrons with different angles with the pads sides.

The Pad Response Function method

The charge spread described in the section 5 results in charge detection in a few pads around the avalanche arrival point. The maximal charge of the waveform signal induced in the adjacent pads is related to the initial charge on the leading pad, in a very different way from the one that would result from a direct charge sharing on the pads, for example due to the diffusion in the TPC gas. That is why the barycentric method is not an optimized estimator of the position reconstruction. The dependence of the induced maximal signal in the adjacent pads versus the initial charge depends mostly on the charge spreading in the Micromegas module but also on other effects like the diffusion in the TPC gas, the electronics response, the integration of the charge over the pad surface area etc. The data allow to measure and calibrate the resulting dependence using the so-called Pad Response Function approach. The Pad Response Function (PRF) which characterizes the relation between observed charge ratios and track position with respect to the pad (Equation 2). This method improves the spatial resolution compared to the barycentric method for TPCs with the resistive anode [START_REF] Attie | Beam tests of Micromegas LC-TPC large prototype[END_REF]. The PRF is defined as:

P RF (x track -x pad ) = Q pad /Q cluster (2)
where x track is the reconstructed position of the track, x pad is the centre of the pad, Q pad is the charge collected on a given pad and Q cluster is the charge collected on the whole cluster. The definition of the cluster is the same as described in section 5: it's a group of pads where one receives a charge from the initial avalanche and the others detect the charge spread in the resistive foil.

To parametrize the PRF we used the empirical ratio of two symmetric 4 th order polynomials proposed in [START_REF] Boudjemline | Spatial resolution of a GEM readout TPC using the charge dispersion signal[END_REF]:

P RF (x, Γ, ∆, a, b) = A × 1 + a 2 x 2 + a 4 x 4 1 + b 2 x 2 + b 4 x 4 (3)
where the parameters a i and b i can be related to the more physical parameters: the full width at half maximum Γ, the base width ∆, and two scaling parameters a and b.

Spatial resolution estimation

This parametrized function of the PRF allows one to determine the parameters using an iterative method. To get a first guess of the PRF parameters, we use the track position reconstruction obtained from the barycentric method. After having defined the position of all the clusters, the global track is fit with a 2 degree polynomial function. The fit is based on many measurements using the different clusters along the track ( 34) thus it is considered as a true track position and the PRF scatter plot is filled (Figure 8 (a)). The scatter plot is profiled along Y axis to form a graph that is further fit with the analytical function from Equation 3 (Figure 8 (b)).

After having determined the PRF, the track position in each cluster is obtained with the following χ 2 minimization:

2 - 1 - 0 1 2 [cm]
χ 2 = pads Q pad /Q cluster -P RF (x track -x pad ) σ Q pad /Q cluster (4)
where σ Q is the uncertainty on the charge measurements. In our analysis, we assume that charge measurement probability follows a Poisson distribution, hence: σ Q pad /Q cluster = Q pad /Q cluster . We proceed through the iterative process of the PRF estimation until the track fit quality is not improving anymore. Typically this procedure converges after few iterations and for the results shown in this paper 10 iterations were used.

The spatial resolution is defined as the difference between the reconstructed position in a given cluster and the track global fit (residual). The particular cluster where the resolution is studied is excluded from the fit to prevent biases. The residuals distribution is fit with a Gaussian function whose standard deviation defines the spatial resolution. The DESY beamline allows changing the momentum of the electrons delivered to the test beam area. In this way, we studied spatial resolution as a function of the track kinematics. We find no significant changes for the position accuracy reconstruction for the tracks in the range between 1 and 5

GeV/c (Figure 10).

Spatial resolution dependence on the track inclination

Tracks inclined with respect to pad borders are expected to be reconstructed less precisely compared to horizontal ones. As an example, in the current ND280 TPCs, the resolution degrades as a function of the track angle from 600 µm to ∼ 1 mm [START_REF] Abgrall | Time Projection Chambers for the T2K Near Detectors[END_REF]. It is then particularly interesting to investigate In order to do this, the different cluster patterns described in section 5

were used for different track inclination. The χ 2 fit (Equation 4) is applied to each cluster to extract the track position. As for horizontal tracks, the positions in given clusters are fit together to form a global fit and the iterative analysis is applied: the barycentric estimation is used as a prior following with PRF calibration. The results of the inclined track spatial resolution estimations are shown in Figure 11.

As expected, the column clustering method leads to a severe spatial resolution degradation with the track slope, reaching a maximum of 1.4 mm for 45 • tracks. For these angles, the diagonal pattern provides a significant performance improvement. In the intermediate regions (20 By taking the best clustering algorithm we observe a spatial resolution better than 600 µm for all the angles. We understand the difference between horizontal tracks and inclined tracks as due to the larger effective pad size for diagonal clustering and to the rectangular shape of the pads, while diagonal clustering would work better for square pads. For a spatial resolution of 600 µm, for 70 point measurements and a magnetic field of 0.2 T we expect a momentum resolution of 6% at 1 GeV/c that scales linearly with the spatial resolution.

Bias measurements

As described above, we define the spatial resolution as a standard deviation of the difference between the reconstructed position in a given cluster and a global track fit. Meanwhile, the mean value of the residuals is also an important characteristic that shows the bias of our measurements. In particular, it is interesting to study the biases with respect to the track position in the pad. For that, we use the natural beam spread. The electron beam profile is nearly Gaussian with a standard deviation ≈ 1 cm. We sample the residuals with the reconstructed track position in the pad. Thus, we can study the resolution and biases in the different pad regions. In Figure 12, we show both the spatial resolution and bias per column. Individual PRFs were used for each column to analyze the behaviour in the different regions of the detector independently. The resolution undergoes some oscillations because if a track is close to one of the leading pad borders, the neighbouring pads see a larger signal and we thus have a more reliable input for the position reconstruction. The resolution and the bias of the track reconstruction over the position in the pad for various columns (represented in legend). The pad borders are represented with vertical lines. In this coordinate system, the beam is centred around Y = -0.5cm, but tracks are also measured on the neighbour pads due to the beam spread.

We generalize the bias study for the whole detector. Figure 13 represents the fluctuations of the resolution and biases in the given column. We conclude that for most of the detector the biases are under control and smaller compared to the spatial resolution. In the downstream part of the detector, we found larger biases that could be related to the non-uniformities in the resistivity of this ERAM detector that will be described in section 9.

Deposited energy resolution

One of the main goals of a TPC is to perform particle identification (PID) based on the measurement of the ionization produced by charged particles crossing the gas volume. The PID capabilities depend on the resolution in the ionization energy loss measurements.

In the case of T2K TPCs, the PID is mainly used to distinguish electrons (produced by ν e ) from muons (produced by ν µ ). In the momentum range studied by T2K, the amount of ionization between electrons and muons differs by ∼ 40%. Therefore, a resolution of less than 10% allows to efficiently distinguish these two particles. In general, the resolution depends on the number of independent ionization measurements (i.e. the number of clusters)

and on the amount of ionization in each cluster. For the existing TPCs, a resolution of 8% was obtained by combining measurements in 2 Micromegas detectors (72 independent measurements of ionization).

In this section, we will describe the performances observed with one single ERAM detector (36×32 pads).

The method used to estimate the energy loss of a given track is called the truncated mean method: the charges contained in each cluster of the track are sorted by increasing order and only a fraction of the lowest charged clusters is kept to compute the mean deposited energy per track. Such a method allows to reject clusters with a large amount of charge, due to fluctuations in the ionization processes, which would degrade the relative resolution on the mean value, and thus the power to separate different types of particles.

The dependence of the dE/dx resolution for the tracks parallel to the pad size on the truncation fraction is shown in Figure 14. The best resolution is obtained for values of truncation fraction between 50% and 80%. Therefore, for all the results presented in this section, a truncation factor of 70% is used. In the existing ND280 TPCs or in the results published in [START_REF] Attié | Performances of a resistive Micromegas module for the Time Projection Chambers of the T2K Near Detector upgrade[END_REF], the cluster charge is defined by summing the waveform maximum seen in each pad composing the cluster. This definition was shown to be a good estimator of the charge in the case of bulk Micromegas, or when the drift distance is large enough to allow the predominance of transverse diffusion over the charge spreading induced in the ERAM module.

However, in the present configuration, the limited size of the TPC implies that the transverse diffusion is small and that charge spreading dominates, except when the track is close to the pad border. Therefore, by summing the maximum of the waveform seen in each pad of the cluster, we are over counting the leading pad charge: first, it is seen in the leading pad before the spread, and then seen by the neighboring pads after the spread.

To remove this over counting effect, the charge per cluster can be defined in a different way. For each cluster, we build a summed waveform defined as the sum of the waveform amplitude seen in each pad of the cluster at all times. We then take as a cluster charge estimator the maximum of this summed waveform. We refer to this method as WF sum .

A comparison of the results of the two methods as a function of the drift distance is shown in Figure 15 for the two different peaking times and for data taken with 0.2 T magnetic field and with the beam parallel to the ERAM side. It can be seen that with both methods, the resolution is well below 10% for all the drift distances and, as expected, the WF sum method gives a better resolution comprised between 8.5% and 9.0%.

a b

Figure 15: dE/dx resolution with respect to the drift distance for the beam parallel to the pad side with a magnetic field of 0.2 T and peaking times of 200 ns (a) and 412 ns (b). Q sum method consists in summing the maximum of the waveform in each pad of the cluster while WF sum corresponds to the maximum of the summed waveforms in a cluster.

Deposited energy resolution for inclined tracks

As explained in section 5, in order to reconstruct inclined tracks, different clustering algorithms are used. In the case of the deposited energy resolution, the usage of such algorithms has two impacts: a larger number of clusters per track are reconstructed but the track will have different paths in different clusters.

In the column or row clusters (defined in Figure 6), tracks with the same angle with respect to the pad sides have the same dx in each cluster (neglecting the curvature induced by the magnetic field). This is not true for diagonal clusters in which the dx can vary between 0 and the diagonal of the pad (∼ 1.5 cm). The distribution of the charge as a function of dx for tracks inclined by 45 degrees with respect to the pad plane is shown in Figure 16.

It is clear from this figure that there is a dependence, although this dependence is not linear as one would expect from the simple consideration that the deposited energy should be proportional to the path.

The non-linearity is due to the fact that each cluster sees not only the direct charge due to the primary ionization, but it also sees some charge due to the spread on the resistive plane and to the transverse diffusion in the gas. In order to correct the dx in each cluster, we fit the charge in each slice of dx with a Landau function and we take the Most Probable Value (MPV).

The distribution of MPV as a function of dx is then parametrized with a third degree polynomial. In each cluster, the charge is corrected to take into account the real path length dx and then the truncated mean is computed.

The dE/dx distributions for diagonal clustering with and without the dx correction are shown in Figure 17. In the case of horizontal or vertical tracks, we observe similar behavior: the horizontal tracks being slightly better because of the larger pad size in the horizontal direction. Extrapolated to two ERAM modules, we obtain a deposited energy resolution of the order of 6% for all the angles. extract RC. Then a simultaneous fit of two neighbouring pads waveform use separate Q(t) functions, as the distance to the track can be different for the two pads, but we consider RC and the electronics response parameters as common in the fit. Since we are using horizontal tracks, we can only fit the y 0 -position of the track and we do not have any constraints on its x-position.

The track position y 0 is obtained with the PRF χ 2 minimization method (see equation ( 4)). Figure 20 shows an example of waveform fit results for the leading pad and its neighbours. Another method can be used to cross-check the RC values obtained with the analytical model. This alternative method studies two parameters related to the signal propagation: the time at which the signal in neighbour pads is maximal and the ratio of amplitudes of the neighbour pads and the leading pad. The time difference between the leading pad and one of its neighbour pads is found to be proportional to RC as can be seen in the following formula: ∆t 1 -∆t 2 = RC × L × y 0 , with L the pad length and y 0 the track position and ∆t 1,2 = t Leading Pad -t Neighbor the time difference between the leading pad and the neighbor pad.

The RC maps obtained using the analytical model and the time difference method described above are shown in Figure 21. The RC values are given in unit of ns/mm 2 . Both methods give results of the same order of magnitude although the second method is less precise than the fit procedure. A nonuniformity of RC up to 30% is observed using both methods.

This non-uniformity is confirmed with the observation of the charge col-lected in the second to leading pad. The charge observed in the leading pad is uniform across the detector while the measurements in the neighbours

593

Q second /Q leading demonstrate fluctuations, especially in the downstream de-594 tector region (Figure 22). The lower charge fraction is consistent with the 595 higher RC value. We measured the performance of the ERAM prototype with beam particles at DESY. We studied both spatial and dE/dx resolution as a function of the angle of the track with respect to the ERAM plane. We also characterized charge spreading and produced a RC map of the prototype. Spatial resolution better than 600 µm is obtained for all the angles using a dedicated clustering algorithm which is adapted to the track angle. Energy resolution better than 9% is obtained for all the angles. We expect ionization energy loss resolution to be better than 7% for tracks crossing two ERAMs. Such performances fully satisfy the requirements for the upgrade of the ND280 TPC.

Figure 1 :

 1 Figure 1: A sketch of the Standard Bulk Micromegas layout (left) and the Resistive anode Micromegas layout (right).

Figure 2 :

 2 Figure 2: View of the setup.

Figure 3 :

 3 Figure 3: The drift velocity under a constant field of 275 V/cm as a function of date (a) and as a function of electric field (b) with density corrections applied. The vertical lines on (a) correspond to the moments when the bottles of premixed gas were replaced. Transition regions at ramp-up and between bottle changes can be seen. The measurements at different fields in (b) correspond to a full day during the first bottle (14.06).

Figure 4 :

 4 Figure 4: A schematic view of a column cluster (left) and the waveforms of each pad composing this cluster (right) for 412 ns shaping time.

Figure 5 :

 5 Figure 5: The number of pads in a cluster versus the DLC high voltage (left) and the fraction of the cluster charge collected in the pad with largest signal at 360 V (right).

Figure 6 :

 6 Figure 6: The different cluster patterns that can be used depending on the track (red line) angle with the pads sides : (a) column, (b) diagonal, (c) 2 by 1, (d) 3 by 1. For each case, the colored pads correspond to one cluster, associated to the colored leading pad crossed by the track.

Figure 7 :

 7 Figure 7: Event displays of the (a) single track and (b) multi-track in the prototype.

Figure 8 :

 8 Figure 8: The Pad Response Function (PRF) obtained with (a) scatter plot and (b) results of its profile and fit with analytical function. The pad borders are represented by vertical lines.

7. 3 .

 3 Spatial resolution dependence on the drift distance, momentum, high voltage The spatial resolution was studied with different samples. The beam position was varied within the drift distance of the field cage keeping the tracks parallel to the Micromegas plane. For resistive Micromegas, the spatial resolution is expected to degrade slightly for a larger drift distance affected by the transverse and longitudinal diffusion. The observed dependence, for a DLC voltage of 360 V, is shown in Figure 9 (a). A resolution between 200 and 250 µm is observed for the whole drift length. In addition, the high voltage applied to the Micromegas mesh was varied to study the detector performance in different regimes. Higher voltage is expected to enhance the initial avalanche, thus increase the smaller charge spreading signal as well. Signals in the neighbour pads are then more likely to pass above the threshold and are less affected by statistical fluctuations. As track position reconstruction relies on the charge spreading measurements we expect better performance with higher voltage. The results are shown in Figure 9 (b).

Figure 9 :

 9 Figure 9: Spatial resolution with respect to (a) beam injection position and (b) Micromegas voltage (DLC voltage) for horizontal tracks parallel to the MM plane. Points represent the mean value over detector columns and errors represent the fluctuations (RMS).

Figure 10 :

 10 Figure 10: Spatial resolution as a function of the momentum of the electron beam.

Figure 11 :

 11 Figure 11: Spatial resolution with respect to track angle within the MM side obtained with different cluster definitions.

Figure 12 :

 12 Figure12: The resolution and the bias of the track reconstruction over the position in the pad for various columns (represented in legend). The pad borders are represented with vertical lines. In this coordinate system, the beam is centred around Y = -0.5cm, but tracks are also measured on the neighbour pads due to the beam spread.

Figure 13 :

 13 Figure 13: The spatial resolution and bias fluctuations observed for the different positions of the track in the pad. Dots represent the mean value in a given column, filled areas correspond to the RMS, and error bars represent minimum and maximum values.

Figure 14 :

 14 Figure 14: The dependence of the dE/dx resolution on the truncation factor α for the beam parallel to the pad sides with a magnetic field of 0.2 T and electronics peaking times of either 200 ns and 412 ns

Figure 16 :

 16 Figure 16: The distribution of the charge with respect to the dx of the track inclined by 45 degrees with respect to the pad plane. The points represent the most probable value from the fit of the vertical slice with Landau function.

Figure 17 :Figure 18 : 8 . 3 .Figure 19 :

 17188319 Figure 17: The dE/dx distribution for 40 deg inclined tracks at an electronics peaking time of 200 ns (a) and 412 ns (b) with and without the correction for dx.

Figure 20 :

 20 Figure 20: Example of waveform fit results for the leading pad (b) and its neighbors (a) and (c) in a given cluster.

Figure 21 :

 21 Figure 21: The RC (ns/mm2 ) map obtained using the fit from the analytical model (left) and the time difference method (right) described in the text.

Figure 22 :

 22 Figure 22: The distribution of the charge in the leading pad (left) and the charge fraction in the second to the leading pad (right) versus the column for horizontal beam tracks.

  2 with 32 × 36 rectangular pads of size 10.09 ×

	11.18 mm 2 . The ERAM module studied in this paper has a resistivity close
	to the required one of 200 kOhm/ using DLC foils stack on a 75 µm glue
	layer. The resistivity maps measured at critical steps of the manufacturing
	of the DLC foil shows a non-uniformity of 20%. The detector is readout
	with two analogue 576 channels AFTER based Front-End Cards (eight 72

  The asymmetry with respect to 45 • is caused by the rectangular pad shape 11.3×10.2 mm. Thus, the diagonal pattern is considered a better choice for 48 • tracks than for 45 • tracks. Hence, tracks inclined with 50 • are reconstructed more accurately compared to 40 • tracks. Similar behavior is observed for all the other patterns.

• , 30 • , 60

• , 70 • ), the best result is achieved with the more complex patterns: "2 by 1" and "3 by 1".

Acknowledgements

The measurements leading to these results have been performed at the Test Beam Facility at DESY Hamburg (Germany), a member of the Helmholtz Association. The authors would like to thank the technical team at the DESY II accelerator and test beam facility for the smooth operation of the test beam and the support during the test beam campaign. We acknowledge the support of CEA and CNRS/IN2P3, France; DFG, Germany; INFN, Italy; National Science Centre (NCN) and Ministry of Science and Higher Education (Grant No. DIR/WK/2017/05), Poland; MINECO and ERDF funds, Spain. In addition, participation of individual researchers and institutions has been further supported by H2020 Grant No. RISE-GA822070-JENNIFER2 2020, MSCA-COFUND-2016 No.754496, ANR-19-CE31-0001, RFBR grants #19-32-90100, the Ministry of Science and Higher Education of Russia (contract #075-15-2020-778) the Spanish Ministerio de Economía y Competitividad (SEIDI -MINECO) under Grants No. PID2019-107564GB-I00 and SEV-2016-0588. IFAE is partially funded by the CERCA program of the Generalitat de Catalunya.

RC map calculation

The quantity controlling the charge spreading over time is the product RC (section 2. To have a better understanding of our detector, we reconstruct the map of RC using horizontal tracks in data where a scan in y-direction was performed at a peaking time of 412 ns. This map is crucial to characterize our detector and its uniformity and is also needed for a detailed detector simulation of the electronics response.

To extract these RC values, we use an analytical model of the charge dispersion, adjusted to the waveforms measured in the pads. The induced charge on a rectangular pad below the resistive layer can be calculated by integrating the charge density function over the pad area [START_REF] Dixit | Position sensing from charge dispersion in micropattern gas detectors with a resistive anode[END_REF] :

with q e is the initial charge, (x 0 ,y 0 ) the track position, x high , x low , y high , y low the pad boundaries. In the denominator σ(t) = (2t/τ ) + ω 2 , the term τ = RC where R is the surface resistivity of the layer and C the capacitance determined by the spacing between the anode and readout planes. Finally, ω is associated to the transverse diffusion term.

To compare to data, the characteristics of the front-end charge preamplifiers need also to be included. Longitudinal diffusion increases the size of electron charge clusters in the drift direction. The longitudinal diffusion is neglected here since we have only 15 cm drift distance. The parameterization of the electronics shaping time effects I(t) is obtained from the simulation.

The convolution of I(t) and Q(t) results in the full theoretical model, which is compared to the data. This convolution is handled numerically. The fit is based on clusters of pads perpendicular to the track. Each cluster consists of a so-called leading pad collecting essentially the initial charge deposit, and in so-called neighbour pads sensitive mostly to the induced charge due to the resistive effect, as shown in Figure 4.

The fit procedure is as follows: we first fit the leading pad waveform with the electronics response function I(t), then we fit simultaneously the two neighbouring pads waveform with a convolution of I(t) with Q(t) in order to


Brahim K Benazzouz

Khac Hieu Ho

Phuoc The Nguyen

Hai Hoang

Guillaume Galliero
email: guillaume.galliero@univ-pau.fr

Diffusive Transport of Gases in Saturated Nanopores: Caprock Leakage from a Molecular Simulation Perspective

Keywords: Oversolubility, Molecular Diffusion, Nanopores, Gases, Molecular Simulations

In this work, molecular dynamics simulations have been performed to shed some light on possible gas leakages by molecular diffusion through reservoir caprocks. To do so, solubility and transport diffusivity of light hydrocarbon gases (methane and ethane) in water-saturated slit nanopores have been investigated at typical reservoir thermodynamics conditions.

Interestingly, it has been noticed that the gas solubility in the saturated nanopores is noticeably higher than in the bulk state. This "oversolubility" is mainly due to an adsorptiondriven effect for both gases. However, a non-negligible confinement-driven effect has also been noticed for ethane. As a consequence, the oversolubility linearly increases with the reciprocal of the pore size. Regarding the transport diffusivity, results have shown that it decreases with the decrease of the pore size. However, rather surprisingly, from the transport diffusivity point of view, the pore size has a weaker impact on ethane than on methane. This is probably related to the fact that ethane is preferentially localized at the depleted water region in the nanopore, whereas methane is mainly localized at the adsorbed water region.

Finally, these results have been used to estimate the influence of these nano-confinement effects on the diffusive gases leakages through a modelled water-saturated caprock over geological times. It has been found that the estimated leakage could be noticeably enhanced by the nano-confinement effects, emphasizing the importance of taking into account these effects on gas leakage estimates through caprock.

Introduction

Caprock is often viewed as an impermeable rock type trapping fluids in reservoirs as discussed in [START_REF] Grunau | A worldwide look at the cap-rock problem[END_REF][START_REF] Krooss | Light hydrocarbon diffusion in a caprock[END_REF][START_REF] Schlömer | Experimental characterisation of the hydrocarbon sealing efficiency of cap rocks[END_REF][START_REF] Selley | Elements of Petroleum Geology[END_REF]. When considering the gas trapping, the sealing efficiency of caprock arises from the very tight nature of the porous medium, composed mostly of nanopores saturated by water in which most gases hardly solubilize. However, recent studies have shown that the solubility of gases in water confined in nanoporous materials, such as clay, can be noticeably greater than that in bulk water [START_REF] Sun | An accurate model to predict the thermodynamic stability of methane hydrate and methane solubility in marine environments[END_REF][START_REF] Phan | Aqueous methane in slit-shaped silica nanopores: high solubility and traces of hydrates[END_REF][START_REF] Yuan | Experimental study and modelling of methane adsorption and diffusion in shale[END_REF][START_REF] Ho | Solubility of gases in water confined in nanoporous materials: ZSM-5[END_REF][START_REF] Li | Water distribution characteristic and effect on methane adsorption capacity in shale clay[END_REF][START_REF] Liu | The oversolubility of methane gas in nano-confined water in nanoporous silica materials[END_REF][START_REF] Grekov | Thermodynamic data of adsorption reval the entry of CH4 and CO2 in smectite clay interlayer[END_REF]. This "oversolubility" phenomenon, combined with molecular diffusion process, may induce a more important than scheduled leakage of gases (hydrocarbons, carbon dioxide…) through a reservoir caprock, when considering typical geological timescales, i.e. from thousands to millions of years which are the typical timescales occurring in geological engineering and in petroleum system [START_REF] Krooss | Light hydrocarbon diffusion in a caprock[END_REF]Krooss et al., 1992a;Krooss et al., 1992b;[START_REF] Krooss | Diffusive losses of hydrocarbons through cap rock. Experimental studies and theoretical considerations[END_REF][START_REF] Schlömer | Experimental characterisation of the hydrocarbon sealing efficiency of cap rocks[END_REF][START_REF] Song | Comprehensive review of caprock-sealing mechanisms for geologic carbon sequestration[END_REF]. Furthermore, considering the reservoir as an open system, i.e. subject to a gas leakage through the caprock, has a strong impact on the way the vertical fluid distribution should be modeled at the initial state in a reservoir as shown in [START_REF] Montel | Initial state of petroleum reservoirs: A comprehensive approach[END_REF]. Thus, to improve the pre-evaluation of hydrocarbon amounts and distribution in the reservoirs using up to date models [START_REF] Montel | Initial state of petroleum reservoirs: A comprehensive approach[END_REF][START_REF] Dembicki | Practical Petroleum Geochemistry for Exploration and Production[END_REF][START_REF] Baghooee | A new approach to thermal segregation in petroleum reservoirs: Algorithm and case studies[END_REF], or to better understand and determine sources and origins of groundwater contamination [START_REF] Nowamooz | Numerical investigation of methane and formation fluid leakage along the casing of a decommissioned shale gas well[END_REF], this overlooked possible leakage phenomenon may worth of investigations which is the purpose of this work.

Solubility and molecular diffusion of hydrocarbons gases in water-saturated nanopores, that are required to estimate the diffusion leakages through the caprock (Krooss et al., 1992a;Krooss et al., 1992b;[START_REF] Schlömer | Experimental characterisation of the hydrocarbon sealing efficiency of cap rocks[END_REF], have attracted interest and attention over the years [START_REF] Sun | An accurate model to predict the thermodynamic stability of methane hydrate and methane solubility in marine environments[END_REF][START_REF] Campos | A molecular dynamics study on natural gas solubility enhancement in water confined to small pores[END_REF][START_REF] Campos | A molecular dynamics study on natural gas solubility enhancement in water confined to small pores[END_REF][START_REF] Phan | Aqueous methane in slit-shaped silica nanopores: high solubility and traces of hydrates[END_REF][START_REF] Yuan | Experimental study and modelling of methane adsorption and diffusion in shale[END_REF][START_REF] Ho | Solubility of gases in water confined in nanoporous materials: ZSM-5[END_REF][START_REF] Li | Water distribution characteristic and effect on methane adsorption capacity in shale clay[END_REF][START_REF] Phan | Confined water determines transport properties of guest molecules in narrow pores[END_REF][START_REF] Bui | Transport mechanism of guest methane in water-filled nanopores[END_REF][START_REF] Gadikota | Hydrophobic solvation of gases (CO2, CH4, H2, noble gases) in clay interlayer nanopores[END_REF][START_REF] Grekov | Thermodynamic data of adsorption reval the entry of CH4 and CO2 in smectite clay interlayer[END_REF][START_REF] Liu | The oversolubility of methane gas in nano-confined water in nanoporous silica materials[END_REF]. Many works, including simulations and experiments, have found that these quantities can be significantly different from those in bulk liquid water. More precisely, in nanoporous structure the solubility of various gases in the confined water can be largely enhanced compared to what occurs in bulk water [START_REF] Sun | An accurate model to predict the thermodynamic stability of methane hydrate and methane solubility in marine environments[END_REF][START_REF] Campos | A molecular dynamics study on natural gas solubility enhancement in water confined to small pores[END_REF][START_REF] Ho | Solubility of gases in water confined in nanoporous materials: ZSM-5[END_REF], a phenomenon sometimes called "oversolubility", and the gas molecular diffusion can be highly reduced [START_REF] Phan | Confined water determines transport properties of guest molecules in narrow pores[END_REF][START_REF] Bui | Transport mechanism of guest methane in water-filled nanopores[END_REF][START_REF] Gadikota | Hydrophobic solvation of gases (CO2, CH4, H2, noble gases) in clay interlayer nanopores[END_REF]. These phenomena are attributed to direct and indirect effects [START_REF] Phan | Aqueous methane in slit-shaped silica nanopores: high solubility and traces of hydrates[END_REF][START_REF] Ho | Solubility of gases in water confined in nanoporous materials: ZSM-5[END_REF][START_REF] Li | Water distribution characteristic and effect on methane adsorption capacity in shale clay[END_REF][START_REF] Phan | Confined water determines transport properties of guest molecules in narrow pores[END_REF][START_REF] Bui | Transport mechanism of guest methane in water-filled nanopores[END_REF].

The direct effect results from the gas-wall interaction, which can lead to an adsorption of gases on the pore walls [START_REF] Ho | Solubility of gases in water confined in nanoporous materials: ZSM-5[END_REF][START_REF] Li | Water distribution characteristic and effect on methane adsorption capacity in shale clay[END_REF]. The indirect effect emerges from the spatial variation of local properties of the confined liquid water due to the water-solid surface interaction, which can have a noticeable impact on the solubility and molecular diffusion of the gases [START_REF] Phan | Aqueous methane in slit-shaped silica nanopores: high solubility and traces of hydrates[END_REF][START_REF] Phan | Confined water determines transport properties of guest molecules in narrow pores[END_REF][START_REF] Bui | Transport mechanism of guest methane in water-filled nanopores[END_REF]. In nanopores, where the surface-to-volume ratio is very high, these effects may be significant [START_REF] Karniadakis | Microflows and Nanoflows[END_REF][START_REF] Israelachvili | Intermolecular and Surface Forces[END_REF]. Thus, when considering nanopores, the pore size is among the major factors that influences both solubility and mass diffusion of gases [START_REF] Sun | An accurate model to predict the thermodynamic stability of methane hydrate and methane solubility in marine environments[END_REF][START_REF] Campos | A molecular dynamics study on natural gas solubility enhancement in water confined to small pores[END_REF]. The systematic evaluation of such an effect on idealized systems is the main goal of the present work.

Molecular dynamics (MD) simulation has shown to be an efficient numerical tool to investigate thermodynamic and transport properties of nanoscale systems and to deal with extreme thermodynamic conditions, e.g. high pressures as found in reservoirs [START_REF] Allen | Computer Simulations of Liquids[END_REF][START_REF] Smit | Understanding Molecular Simulation: From Algorithms to Applications[END_REF][START_REF] Hoang | Simultaneous Description of Equilibrium, Interfacial, and Transport Properties of Fluids Using a Mie Chain Coarse-Grained Force Field[END_REF]2019;2021). This tool has been used to investigate either the solubility or the molecular diffusion of various gases in solvents confined in nanopores, often at ambient conditions [START_REF] Campos | A molecular dynamics study on natural gas solubility enhancement in water confined to small pores[END_REF][START_REF] Phan | Aqueous methane in slit-shaped silica nanopores: high solubility and traces of hydrates[END_REF][START_REF] Phan | Confined water determines transport properties of guest molecules in narrow pores[END_REF][START_REF] Bui | Transport mechanism of guest methane in water-filled nanopores[END_REF][START_REF] Gadikota | Hydrophobic solvation of gases (CO2, CH4, H2, noble gases) in clay interlayer nanopores[END_REF], yielding results consistent with experiments [START_REF] Phan | Aqueous methane in slit-shaped silica nanopores: high solubility and traces of hydrates[END_REF][START_REF] Gadikota | Hydrophobic solvation of gases (CO2, CH4, H2, noble gases) in clay interlayer nanopores[END_REF]. Interestingly, this numerical tool can provide many information at the molecular scale, such as thermophysical and transport properties [START_REF] Allen | Computer Simulations of Liquids[END_REF][START_REF] Smit | Understanding Molecular Simulation: From Algorithms to Applications[END_REF]), so facilitating the interpretation of the results [START_REF] Phan | Aqueous methane in slit-shaped silica nanopores: high solubility and traces of hydrates[END_REF][START_REF] Phan | Confined water determines transport properties of guest molecules in narrow pores[END_REF][START_REF] Bui | Transport mechanism of guest methane in water-filled nanopores[END_REF]. Thus, in this work, MD simulation has been employed to investigate the effect of pore size on both solubilities and molecular diffusion of light hydrocarbons gases, methane and ethane, in water-saturated nano-pores, i.e. pore sizes typical of that encountered in a caprock, at typical reservoir conditions. This article is organized as follows. In section 2, some details on the MD simulations are described. Then, results obtained from the molecular dynamics simulations are presented and discussed in section 3. Finally, the conclusion that summarizes the main outcomes is provided in section 4.

Models and Simulation

Fluid and solid models

For the sake of simplicity, and because caprock complexity in terms of porous structure and fluid composition is not accessible to molecular simulations, all hydrocarbon gases (methane and ethane), water and solid molecules have been modelled using one of the simplest coarse-grained models, the Lennard-Jones one. Thus, all interactions between particles have been described by a classical Lennard-Jones (LJ) potential [START_REF] Lennard-Jones | The electronic structure of some polyenes and aromatic molecules I-The nature of the links by the method of molecular orbitals[END_REF] as:

= 4 - (1)
where, r is the distance between the two particles, ε is the potential depth, σ is the particle "diameter". This simple model has shown to simultaneously predict rather well thermodynamic and transport properties of small hydrocarbons gases, i.e. methane and ethane [START_REF] Galliero | Molecular dynamics comparative study of Lennard-Jones -6 and exponential -6 potentials: Application to real simple fluids (viscosity and pressure)[END_REF]. The LJ parameters of methane used are those coming from the study of [START_REF] Martin | Transferable potentials for phase equilibria. 1. United atom description of n-alkanes[END_REF], and for ethane, they were obtained from Galliero et al. [START_REF] Galliero | Molecular dynamics comparative study of Lennard-Jones -6 and exponential -6 potentials: Application to real simple fluids (viscosity and pressure)[END_REF]. They are given in Table 1.

When considering a LJ coarse-grained representation of species exhibiting polar or hydrogen bonding interaction, the effective LJ parameters should vary with the thermodynamic state to provide a good prediction for both thermodynamic and transport properties [START_REF] Müller | Molecular modeling of fluid-phase equilibria using an isotropic multipolar potential[END_REF][START_REF] Hadley | Coarse-Grained Molecular Models of Water: A review[END_REF]. Thus, the LJ parameters of water molecules at the two studied reservoir conditions (T=373.15 K, P=50 MPa and T=423.15 K, P=100 MPa) have been adjusted to simultaneously provide a good agreement with density and viscosity data taken from the NIST reference database [START_REF] Lemmon | Reference Fluid Thermodynamic and Transport Properties[END_REF]. The LJ equation of equation of [START_REF] Kolafa | The Lennard-Jones fluid: an accurate analytic and theoreticallybased equation of state[END_REF], and

the LJ viscosity correlation of [START_REF] Galliero | Molecular Dynamics Study of the Lennard-Jones Fluid Viscosity: Application to Real Fluids[END_REF] were employed to perform the adjustment of water LJ parameters. Results are listed in Table 1. The nanopore walls have been modelled to represent the macroscopic characteristics of typical silicate mineral surfaces from the surface energy point of view. They are made of identical spherical particles arranged on a faced centered cubic (FCC) lattice. The LJ diameter of solid particles and the lattice size have been chosen to quantitatively mimic the roughness of silicate mineral surface. The LJ energy of solid particles has been selected to yield reasonable Hamaker constant of silicate mineral surface of about = 1 × 10 Joule [START_REF] Iwamatsu | A molecular theory of solvation force oscillations in nonpolar liquids[END_REF]. The parameters of the solid walls are shown in Table 2. The LJ potential parameters of the interaction between unlike molecules were calculated by using modified Lorentz-Berthelot rules [START_REF] Schnabel | Unlike Lennard-Jones parameters for vapor-liquid equilibria[END_REF]:

! "# = $ "#% ! " ! # (2a) & "# = '& " + & # ) (2b)
where, i and j indicate the species, and $ "# is a pre-factor which was introduced to modulate the strength of energy between species [START_REF] Schnabel | Unlike Lennard-Jones parameters for vapor-liquid equilibria[END_REF]. Such a parameter is often used to better model cross interactions between asymmetric molecules [START_REF] Hamani | Excess volume, isothermal compressibility, isentropic compressibility and speed of sound of carbon dioxide+ n-heptane binary mixture under pressure up to 70 MPa. II. Molecular simulations[END_REF]. $ "# have been set to 1.0 for the interactions between the fluids (hydrocarbons gases and water) and solid walls. Between the gases and the water, $ "# have been adjusted to provide good solubilities of gases in bulk liquid water, as detailed in section 3.1.

Simulated configuration

To determine the properties of light hydrocarbon gases in the saturated nanopores, MD simulations in the grand canonical ensemble were performed by using a scheme similar to the one proposed in (Hoang and Galliero, 2012a), as sketched in Fig 1a . In terms of thermodynamics conditions, two reservoir conditions have been considered: (1) T=373.15 K and P=50 MPa, and (2) T=423.15 K and P=100 MPa. To address the effect of pore size on solubility and molecular diffusion of gases, MD simulations have been performed on systems with various pore sizes, W, ranging from 5 to 10 times the diameter of a gas molecule, i.e.

from about 1.90 nm to about 4.24 nm, as one would expect in an ideal caprock [START_REF] Schlömer | Experimental characterisation of the hydrocarbon sealing efficiency of cap rocks[END_REF][START_REF] Kampman | Observational evidence confirms modelling of the long-term integrity of CO2-reservoir caprocks[END_REF][START_REF] Espinoza | CO2 breakthrough -Caprock sealing efficiency and integrity for carbon geological storage[END_REF]. Such pore sizes correspond to permeability, k, of the micro-Darcy order of magnitude (k scales with W²) as those found in caprocks, even if defining correctly permeability in such nanopore from geometric argument only is questionable due to the strong variations of the confined fluid properties [START_REF] Hoang | Shear behavior of a confined thin film: Influence of the molecular dynamics scheme employed[END_REF].

In details, the solid walls are placed in the middle of the simulation box in the direction y, and at the top and the bottom in the direction x. Each solid wall is composed of 4×10×10 FCC cells in the x, y and z direction respectively. They are separated by a distance W in the x direction, defined as the distance between two innermost opposing atomic solid layers. The walls are of finite size in the x and y directions, whereas they are extended through the whole cell in the z direction. In the reservoirs that are located sufficiently away from the solid walls to avoid the effects of the finite solid wall size in the y direction (Hoang and Galliero, 2012a), the temperature and pressure of the fluids are kept constant.

The initial configuration to compute mass diffusion coefficient and solubility of gases in the nanopores were setup as follows. First, grand-canonical molecular dynamics (GCMD)

simulations were performed to equilibrate separately two systems: one consisting of water and solid molecules, and one containing gases (methane/ethane) and water, as illustrated in Figs.

1a and 1b, respectively. Then, the water molecules in the pore of the first system and the fluid (gases + water) molecules in the reservoirs of the second system were brought together to set up the initial configuration, see Fig. 1. The simulation systems contain at least 1000 water molecules and 5000 gas molecules.

Simulation details

During MD simulations, which were performed using an in-house code (Hoang and Galliero, 2012a;2013), the equations of motion of fluid particles were integrated using the velocity-Verlet algorithm with a timestep *+ ≈ 0.0023ps [START_REF] Verlet | Computer "Experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules[END_REF], whereas the solid particles were fixed at their sites on the FCC lattice. Periodic boundary conditions were applied in all the directions [START_REF] Allen | Computer Simulations of Liquids[END_REF][START_REF] Smit | Understanding Molecular Simulation: From Algorithms to Applications[END_REF]. The temperature and pressure in the reservoirs were maintained at reservoir conditions by using a

Berendsen thermostat and barostat, respectively [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF]. The sampling region, of a width equal to & 234 , where the properties of the confined fluids were computed, was chosen be located at the center of the pore in the y direction to avoid effects due to the finite size of the solid walls (Hoang and Galliero, 2012a;2013), see Fig. 1a.

After having built the initial configuration, as detailed previously, the GCMD simulations were performed in two steps. First, the systems were equilibrated on 2 × 10 time steps, during which data in the sample region were collected to analyze the molecular diffusion of hydrocarbons gases during the unsteady regime. Then, a run of 2 × 10 time steps was achieved to compute the gas solubility and all equilibrium properties of the confined fluids.

Results and Discussions

Determination of the $ "# between gases and water

Similarly to what is done when using a Soreide and Whitson Equation of State [START_REF] Soreide | Peng-Ronbinson predictions for hydrocarbons, CO2, N2, and H2S with pure water and NaCl brine[END_REF], the $ "# of the interactions between the light hydrocarbons gases (methane/ethane) and water molecules have been selected to reproduce the experiment solubility (hereafter, the solubility is defined as a molar fraction and is noted S) of the gases in water at the studied thermodynamic reservoir conditions. More precisely, the values of $ "# of methane-water interaction have been adjusted to yield 5 6 ,89:; = 0.39% and 5 6 ,89:; = 0.78% at (T=373.15 K, P=50 MPa) and (T=423.15 K, P=100 MPa), respectively [START_REF] Duan | A thermodynamic model for calculating methane solubility, density and gas phase composition of methane-bearing aqueous fluids from 273 to 523 K and from 1 to 2000 bar[END_REF]. For the ethane-water interaction, they were adjusted to provide 5 6 ,89:; = 0.13% at (T=373.15 K, P=50 MPa) and 5 6 ,89:; = 0.28% at (T=423.15 K, P=100 MPa) [START_REF] Mao | An accurate model for calculating C2H6 solubility in pure water and aqueous NaCl solutions[END_REF]. For this purpose, the gas solubility dependence on $ "# in water has been estimated at the studied thermodynamic reservoir conditions by performing GCMD simulations on systems containing gas and water as shown in Fig. 1b. Simulation results are illustrated in Fig.

2.

The results clearly show that the gas solubility monotonously increases with $ "# in all cases, allowing to estimate the optimal $ "# by interpolating the simulation data. Optimized $ "# values are provided in Table 3. It is worth noticing that the $ "# significantly varies with the thermodynamic conditions. This can be probably understood from the strong dependence of the LJ potential parameters of coarse-grained water molecules on the thermodynamic conditions, see Table 2. In the following, these $ "# values have been employed to estimate the solubility and the molecular diffusion of methane and ethane in the saturated nanopores.

Effect of the pore size on solubility

Methane

Figure 3 shows the solubility of methane in the saturated nanopores for the two studied reservoir conditions and for pore sizes ranging from 1.86 nm to 3.72 nm. As expected, the solubility of methane in the nanopores is noticeably larger than that in bulk water [START_REF] Campos | A molecular dynamics study on natural gas solubility enhancement in water confined to small pores[END_REF][START_REF] Ho | Solubility of gases in water confined in nanoporous materials: ZSM-5[END_REF], up to three times the bulk value for the smallest pore. To better understand this behavior, the local solubility profiles along the x direction have been computed. As shown in Fig. 4, the increase of solubility results from a local solubility which is significantly larger (up to about fifteen times) in the vicinity of the solid surfaces than in the center of the pore which exhibits solubility values similar to what found in bulk water, as expected.

Regarding the physical origins of the increase of solubility of gases in saturated nanopores, two possible molecular mechanisms have been proposed in the literature [START_REF] Ho | Gas uptake in solvents confined in mesopores: Adsorption versus enhanced solubility[END_REF]: (1) Adsorption-driven effect and (2) Confinement-driven effect. The adsorptiondriven effect corresponds to an adsorbed gas layer at a distance of about & 4B:@C D34 from the solid surface [START_REF] Ho | Gas uptake in solvents confined in mesopores: Adsorption versus enhanced solubility[END_REF]2015). The confinement-driven effect is related to a gas layer localized around the solvent depleted zone [START_REF] Ho | Gas uptake in solvents confined in mesopores: Adsorption versus enhanced solubility[END_REF]2015). In our systems, results illustrated in Fig. 5 indicate that the methane peak is mostly located at a distance around & 4B:@C D34 from the solid surface, indicating so that the adsorption-driven effect has a dominant effect on the increase of solubility of methane in the simulated systems.

To quantify the increase of solubility due to the confinement, the oversolubility, ∆5 234 , has been computed. It is defined as the difference between the solubilities of gas in nanopore and in bulk water by: ∆5 234 = 5 234,FBGH -5 234,89:;

(3) where 5 234,FBGH is the solubility of gas in the nanopore and 5 234,89:; the value in the bulk.

Results are shown in Fig. 3. It has been found that the oversolubility increases linearly with the inverse of the pore size. This feature is consistent with experiment results [START_REF] Sun | An accurate model to predict the thermodynamic stability of methane hydrate and methane solubility in marine environments[END_REF]Miachon et al., 2013) and is simply due to the adsorption-driven nature of this phenomenon and because the local increase of solubility does not depend much on the confinement, see Fig. 4. In addition, the oversolubility is the largest for the highest reservoir temperature. This effect is subtle and is probably related to temperature effects that are stronger than pressure effects on water in terms of structure, see Fig. 5.

Ethane

The solubility of ethane in the saturated nanopores for the two studied reservoir conditions and for pore sizes ranging from about 2.12 nm to 4.24 nm is displayed in Fig. 3.

Results clearly show that the trends for ethane are similar to those found for methane, keeping in mind that ethane is roughly three times less soluble than methane in bulk water.

More precisely, ethane exhibits a local solubility profile similar to that found for methane for both conditions. However, the local increase of solubility close to the solid walls is more pronounced than in methane, reaching up values 40 times larger than the bulk value.

Another interesting feature is that the oversolubility of ethane scales well with the inverse of the pore size even if, for the smallest pores, the oversolubility is nearly constant. This behavior is consistent with the adsorption-driven effect combined with a non-negligible contribution related to the confinement-driven effect for the smallest pores, as confirmed by the density profiles shown in Fig. 5.

Effect of the pore size on molecular diffusion

There exist many different formalisms to describe mass diffusion leading to different diffusion coefficients definitions [START_REF] Cussler | Diffusion -Mass Transfer in Fluid Systems[END_REF][START_REF] Wesselingh | Mass Transfer in Multicomponent Mixtures[END_REF]. This is even more complex in nanopores, even at the pore scale, as the transport coefficients may become space dependent, non-local and non-isotropic [START_REF] Botan | Hydrodynamics in Clay Nanopores[END_REF]Hoang and Galliero, 2012b).

To keep the problem simple, a 1D Fickian formalism has been chosen to describe the gas diffusion through the saturated nanopores. Such a formalism simply connects the mass flux of the gas to its compositional gradient. A corresponding diffusion coefficient, i.e. Fickian diffusion coefficient I 234 sometimes called transport diffusivity, has been estimated by using a transient non-equilibrium MD (T-NEMD) approach that is rather straightforward to implement in MD simulations [START_REF] Maginn | Transport Diffusivity of Methane in Silicalite from Equilibrium and Nonequilibrium Simulations[END_REF][START_REF] Arya | A critical comparison of equilibrium, nonequilibrium and boundary-driven molecular dynamics techniques for studying transport in microporous materials[END_REF], see Fig. 1. This method consists in setting up an initial gas concentration gradient in the simulation system, and then to monitor the relaxation of the system to the equilibrium state assuming a purely diffusive behavior, as described in section 2.2.

More precisely, to estimate I 234 , the evolution of the molar fraction of the gas in the sampling region J 234 + (averaged over 100 time-steps) has been fitted by a simple 1D diffusion equation as [START_REF] Witherspoon | Diffusion of Methane, Ethane, Propane, and n-Butane in Water from 25 to 43°[END_REF][START_REF] Cussler | Diffusion -Mass Transfer in Fluid Systems[END_REF]:

K LMN O K LMN P = 1 - Q R S ∑ UV S exp -X2Y + 1Z O X S R S ⁄ \ LMN Z U (4)
where, J 234 ] is the molar fraction value at equilibrium and L is the characteristic length of the pore over which diffusion occur. The latter has been taken equal to half of the pore length in the y direction plus one molecular diameter of gas to take into account a finite size effect due to the initial configuration setup.

As confirmed by the simulation results shown in Fig. 6, the evolution over time of J 234 can be well fitted by Eq. ( 4) for both gases, indicating that the T-NEMD approach allows to obtain a good estimation of I 234 , with an uncertainty of about 25%.

Methane

Figure 7 provides the variation in the transport diffusivity of methane relatively to the pore size. Results indicate that the pore size has a slight, but non-negligible, effect on the transport diffusivity. More precisely, I 6 tends to decrease with the decrease of the pore size, but this decrease is not very pronounced. This result is consistent with what has been found by many authors on similar systems [START_REF] Hoang | Molecular dynamics study on the effect of solution-wall interaction potential on the properties of solution in uniformly charged hydrophobic channel[END_REF]Hoang and Galliero, 2012a;[START_REF] Zaragoza | Molecular dynamics study of nanoconfined TIP4P/2005 water: How confinement and temperature affect diffusion and viscosity[END_REF][START_REF] Tran | Molecular dynamics simulations on aqueous solution confined in charged nanochannels: Asymmetric effect of surface charge[END_REF][START_REF] Zhang | Transport Property of Methane and Ethane in K-Illite Nanopores of Shale: Insights from Molecular Dynamic Simulations[END_REF], which is that the local molecular diffusion of a confined fluid in the plane parallel to the solid surface does not vary significantly with the pore width when dealing with not too strong fluid-solid interaction. Furthermore, the transport diffusivity of methane in the water-saturated nanopores has been compared with the self-diffusion coefficient of methane in bulk water which is assumed here to be equivalent to the Fickian diffusion coefficient due to the low solubility of gases in water. This diffusion coefficient was computed by performing equilibrium molecular dynamics (EMD) simulations in combination with the usual Einstein's relation [START_REF] Allen | Computer Simulations of Liquids[END_REF][START_REF] Smit | Understanding Molecular Simulation: From Algorithms to Applications[END_REF]:

D _ = lim b→] b 〈Xe @ t -e @ 0 Z 〉 (5)
where, D _ is the self-diffusion coefficient, t is the time, 〈… 〉 is the average over molecules and e @ t is the position of molecule i at the time t.

As expected, results indicate that the transport diffusivity of methane is smaller in the confined water than in the bulk water, an effect which is more pronounced at the highest temperature, see Fig. 7. Furthermore, I 6 tends to its bulk value for the largest pores.

However, the relative effect for the two thermodynamic states is somewhat unexpected as it seems to be more important at T=423 K and P=100 MPa which corresponds to a less structured confined water, see Fig. 5.

It is worth noticing that the transport diffusivity of methane obtained here ranges from 4.5×10 -9 to 10.0×10 -9 m 2 /s, which is typically one order of magnitude higher than the reported experiment results in porous medium [START_REF] Schlömer | Experimental characterisation of the hydrocarbon sealing efficiency of cap rocks[END_REF][START_REF] Schloemer | Molecular transport of methane, ethane and nitrogen and the influence of diffusion on the chemical and isotopic composition of natural gas accumulations[END_REF]. However, here I 6 has been estimated at the pore scale, whereas the transport diffusivity was experimentally deduced at the porous medium scale, i.e. including the effect of tortuosity and porosity (Krooss et al., 1992a;1992b;[START_REF] Krooss | Diffusion of methane and ethane through the reservoir cap rock; implications for the timing and duration of catagenesis: Discussion[END_REF][START_REF] Song | Comprehensive review of caprock-sealing mechanisms for geologic carbon sequestration[END_REF].

Ethane

Pore-size dependence of the transport diffusivity of ethane in the saturated nanopores is shown in Fig. 7. Results indicate that its dependence on pore size is similar to that for methane, showing a slight decrease when the pore size decreases and values tending to the bulk data for the largest pores. In addition, as expected, I 6 is slightly smaller than I 6 in particular at the highest temperature. The only unexpected result is that the effect of the confinement seems less pronounced, in relative, on ethane than on methane, whereas the ethane molecule is bigger than the methane molecule. This is probably related to the localization of the gases relatively to water in the nanopore. More precisely, the ethane molecules are preferably located around the depleted water region, whereas methane molecules are mostly located in the adsorbed water region, see Fig. 5.

Leakage of gases through silicate caprock

In the most favorable case from the caprock sealing point of view, the leakage from a gas reservoir is only due to molecular diffusion across the low permeability barrier formed by the caprock. This contribution is often assumed to be negligible in most cases [START_REF] Schloemer | Molecular transport of methane, ethane and nitrogen and the influence of diffusion on the chemical and isotopic composition of natural gas accumulations[END_REF], even if it could play a role in some specific cases [START_REF] Wang | Anomalies of natural gas compositions and carbon isotope ratios caused by gas diffusion-A case from the Donghe Sandstone reservoir in the Hadexun Oilfield, Tarim Basin, northwest China[END_REF]. Thus, it is interesting to revisit the quantification of diffusive gas reservoir leakage over geological times using the previous results obtained on solubility and diffusion of gases in nanopores. To do so, it has been assumed that the leakage through the caprock is simply governed by a 1D diffusion process. In this framework, the cumulative amount of gas leakage per unit of caprock surface over a time t can be estimated from (Krooss et al., 1992a;1992b;[START_REF] Schlömer | Experimental characterisation of the hydrocarbon sealing efficiency of cap rocks[END_REF]:

i 234 + = j \ LMN kll K LMN m nMo + + K LMN m nMo R S ∑ U S 1 -exp j- U S \ LMN kll R S nMo S +p Uq p (6)
where, r is the porosity, s 63t is the thickness of the caprock, and I uv wxx is the effective transport diffusivity that is defined as I 234 y ⁄ with y being the tortuosity [START_REF] Song | Comprehensive review of caprock-sealing mechanisms for geologic carbon sequestration[END_REF]. As an example, the model caprock considered in this work is assumed to possess a thickness of 100m, a porosity of 1% and a tortuosity of 50, and is composed of watersaturated nanopores having a pore size distribution peaked on 3nm [START_REF] Schlömer | Experimental characterisation of the hydrocarbon sealing efficiency of cap rocks[END_REF][START_REF] Kampman | Observational evidence confirms modelling of the long-term integrity of CO2-reservoir caprocks[END_REF][START_REF] Espinoza | CO2 breakthrough -Caprock sealing efficiency and integrity for carbon geological storage[END_REF].

Evolutions of cumulative amount of gases leakage through the model caprock are shown in Fig. 8. They have been computed from the 1D diffusion process Eq. ( 6) using the solubility and the transport diffusivity values of gases in the saturated nanopore and in the bulk water over 50 millions of years. Interestingly, results clearly indicate that the nanoconfinement effects are able to noticeably increase the leakage of gases compared to what would be deduced without them.

The increased leakage compared to what could be estimated without taking into account the effects of nano-confinement reflects the fact that the oversolubility is dominant over the decrease of diffusion coefficient. More precisely, on geological time, this leakage increase can be approximated, see Eq. ( 6), by the ratio of the product I 234 wxx J 234 under confinement relatively to its value in bulk water. From our data this means an increase in gas leakage of about 165% for methane and 265% for ethane. Even if evaluated on an oversimplified system, this result indicates that the estimation of the diffusive leakage of gases through caprock should take into account the influence of nano-confinement effects.

Conclusions

In this work, molecular dynamics simulations have been performed in grand canonical like ensemble to investigate the solubility and transport diffusivity of light hydrocarbon gases Very interestingly, it has been shown that the solubility of hydrocarbon gases in the water-saturated nanopores is significantly higher than that in bulk water (so-called oversolubility). This results from the local gas solubility which is strongly enhanced in the vicinity of the solid surface and can be up to 15 times (40 times) higher than the solubility of methane (ethane) in bulk water. Such a behavior is mainly induced by an adsorption-driven effect for methane, whereas there is an additional non-negligible contribution due to a confinement-driven effect for ethane. As a consequence, it has been shown that the oversolubility linearly increases with the reciprocal of the pore size. Regarding the transport diffusivity of hydrocarbon gases in the water-saturated nanopores which has been estimated by using a "transient" non-equilibrium molecular dynamic (T-NEMD) approach, it has been shown that the pore size has a weak, but nonnegligible impact on it. More precisely, the transport diffusivity decreases with the decrease of the pore size. However, this effect is more pronounced on methane than on ethane which is rather surprising as methane molecules are smaller than ethane molecules. This is probably due to the fact that methane is mostly located at the adsorbed water region in the nanopore, whereas ethane is mainly located in the depleted water region.

Finally, these results have been used to quantify the diffusive gas leakage through an idealized water-saturated caprock over geological time thanks to a 1D diffusion model.

Results have shown that the nano-confinement effects can noticeably enhance the leakage of gases by diffusion, because the effect of oversolubility of gases under confinement is dominant over the decrease of diffusivity transport due to pore size. This result, even if deduced from an oversimplified system compared to real cases, emphasizes the importance of taking into account nano-confinement effects on gas leakage estimates.

Figure 1 :

 1 Figure 1: (a) Equilibrated system of water with the nanopore. (b) Equilibrated system containing gas and water. (c) Initial configuration of the simulated system to determine diffusion and solubility of gases in the nanopore. Red color corresponds to the solid particle, green color to water molecules and blue color to gas molecules.

Figure 2 :

 2 Figure 2: Dependence of hydrocarbon gas solubility in bulk liquid water on $ "# (a): Methane. (b): Ethane. (Red color) Squares corresponds to T = 373.15K and P = 50MPa. (Blue color) and circles to T = 423.15K and P = 100MPa. Dashed-Dotted curves are exponential fitting functions. Solid horizontal and vertical lines correspond to the experimental solubility of gases in water and the optimized $ "# values, respectively.

Figure 3 :Figure 4 :

 34 Figure 3: Variation in solubility and over-solubility (inserts) of gas with the inverse of the pore size in the confined region. (a): Methane. (b): Ethane. (Red color) Squares corresponds to T = 373.15K and P = 50MPa. (Blue color) Circles to T = 423.15K and P = 100MPa. Solid lines are linear fitting functions.

Figure 5 :

 5 Figure 5: Density profiles of gas and water in the sampling region (a): Methane in the nanopore of W=5σC1. (b): Ethane in the nanopore of W=5σC2. (c): Methane in the nanopore of W=9σC1. (d): Ethane in the nanopore of W=9σC2. Solid lines correspond to T = 373.15K and P = 50MPa and Dashed-dotted lines correspond to T = 423.15K and P = 100MPa. (Green color) Vertical dashed lines correspond to a distance of one gas molecular size from the solid surface.

Figure 6 :

 6 Figure 6: Evolutions of J 234 + J 234 ] ⁄ in the sampling region for a pore size equal to 9 times the corresponding gas diameters at the reservoir thermodynamic condition T = 423.15K and P = 100MPa. (a): Methane. (b): Ethane.

Figure 7 :

 7 Figure 7: Dependence of transport diffusivity of hydrocarbon gas with the pore size. (a): Methane. (b): Ethane. Squares correspond to the thermodynamic condition T = 373.15K and P = 50MPa. Circles correspond to T = 423.15K and P = 100MPa. Lines correspond to the self-diffusion coefficients of gases in bulk water.

Figure 8 :

 8 Figure 8: Cumulative amount of hydrocarbon gases leaking through the caprock by the molecular diffusion at the reservoir thermodynamic condition T = 373.15K and P = 50MPa. (a) Methane. (b) Ethane. (Red color) Solid lines correspond to D 234 and C 234 of gases in the saturated nanopores (with taking into account the effect of confinement). (Blue color) Dashed-Dotted lines correspond to D 234 and C 234 of gases in bulk water (without taking into account the effects of confinement).

(

  methane and ethane) in water-saturated slit nanopores at reservoir thermodynamics conditions, (T=373.15K, P=50MPa) and (T=423.15K, P=100MPa). For simplicity, all fluid and solid molecules have been modelled by Lennard-Jones (LJ) particles. Fluid LJ parameters have been adjusted to yield a set of target properties, such as density, viscosity and solubility, in good agreement with the experiments. Solid parameters, including the LJ potential and structure, have been chosen to quantitatively represent the roughness and surface energy characteristics of silicate mineral surfaces.

Table 1 :

 1 Lennard-Jones parameters of fluid molecules.

	Fluid	σ (nm)	ε (kJ/mol)
	Water (T=373K)	0.2876	3.500
	Water (T=423K)	0.2805	3.200
	Methane	0.3723	1.245
	Ethane	0.4239	1.978

Table 2 :

 2 Lennard-Jones parameters and lattice size of solid wall.

	σ (nm)	ε (kJ/mol)	Lattice size (nm)
	0.3160	0.762	0.4469

Table 3 :

 3 Optimized $ @A between methane/ethane and water.

		Water (T=373K)	Water (T=423K)
	Methane	0.72	0.47
	Ethane	0.64	0.48

Acknowledgements

We gratefully acknowledge TOTAL S.A. for the post-doctoral grant awarded to two of us (BKB and HH) and for letting us publish these results. We also thank the Pau University and the MCIA for providing computational facilities.