N Rougemaille

F Montaigne

B Canals

A Duluard

D Lacour

M Hehn

R Belkhou

O Fruchart

S El Moussaoui

A Bendounan

F Maccherozzi

Artificial kagomé arrays of nanomagnets: a frozen dipolar spin ice

Keywords: numbers: 75.10.Hk, 75.50.Lk, 75.70.Cn, 75.60.Jk

Magnetic frustration effects in artificial kagomé arrays of nanomagnets are investigated using X-PEEM microscopy and monte carlo simulations. Spin configurations of demagnetized networks reveal unambiguous signatures of long range, dipolar interaction between the nanomagnets. As soon as the system enters the spin ice manifold, the kagomé dipolar spin ice model captures the observed physics, while the short range kagomé spin ice model fails.

Complex architectures of nanostructures are nowadays routinely elaborated using bottom-up and/or microfabrication processes. This technological capability allows scientists to engineer materials with properties that do not exist in nature, such as negative refractive index [1] or tunable photonic bandgap [2], but also to elaborate model systems to explore fundamental issues in condensed matter physics. Frustrated spin models are one example for which theoretical predictions can be tested experimentally [START_REF] Gingras | Highly Frustrated Magnetism[END_REF]. In these models, spins are placed on an infinite network of a given geometry, and interact through a short or long range coupling of ferromagnetic or antiferromagnetic nature. Under certain conditions, magnetic frustration effects occur, giving rise to an infinite degeneracy of the ground state, i.e. a finite entropy at zero temperature. Several experimental studies on two-dimensional, artificial arrays of ferromagnetic nanomagnets have been reported recently [START_REF] Wang | [END_REF][5][6][7][8], bringing together the communities of magnetic frustration, statistical physics and nanomagnetism. Properties of these arrays seem to be well described by short range spin ice (SI) models: they show disorder, satisfy ice-like rules, and correlation coefficients fit with predictions [6].

In the kagomé SI model, spins interact through first neighbor exchange coupling only. Experimentally however, nanomagnets are coupled via the magnetostatic interaction, which extends beyond nearest neighbors. While the main interest for frustrated compounds arises from the massive degeneracy of their ground state, this degeneracy is fully lifted when long range, dipolar interactions are included in the model [9,10]. Understanding whether or not these long range interactions influence the local spin configuration in artificial arrays of nanomagnets is thus essential, especially because these networks are often considered as a playground to study magnetic frustration effects on a mesoscopic scale. Combining experimental and theoretical investigations, we show that the kagomé dipolar spin ice (DSI) model captures the physics we observe, where the kagomé SI model fails. Honeycomb arrays of nanomagnets have been fabricated from a Si//Ta(5nm)/Co(10nm)/Ru(2nm) thin film. An aluminium mask, defined by e-beam lithography and lift-off, has been used to protect the structures during ion beam etching. The Al mask has been subsequently removed by chemical etching. The nanomagnets have a typical size of 470 × 70 × 10 nm 3 , ensuring a high uniaxial shape anisotropy, so that they can be considered, locally, as Ising pseudo spins. The distance between nearest nodes in the network is set to 500 nm (Fig. 1a). The Co nanomagnets are thus physically disconnected, and are coupled through the magnetostatic interaction. The networks of nodes have kagomé topology (Fig. 1b) and are composed of 342 Co nanomagnets. The magnetic configuration of each element is resolved using the X-ray PhotoEmission Electron Microscope (X-PEEM, SOLEIL instrument) of the Nanospectroscopy beamline at the ELETTRA synchrotron facility (Trieste, Italy). Magnetic images are acquired at the Co L 3 edge at remanence, after demagnetizing the networks ex situ. The importance of the demagnetizing procedure to approach the ground state spin configurations has been highlighted in previous works [11,12]. To demagnetize our networks, the sample is rotated with a period of about 20 ms in a damped, in-plane sinusoidal magnetic field (1 s period, 0.11 Oe amplitude variation between two periods). From the X-PEEM images, the magnetic state of each individual nanomagnet is measured (Fig. 1a). Data presented here have been obtained from 9 nominally identical networks, which have been demagnetized 6 times, out of which 34 magnetic configurations have been determined. The kagomé ice rule, i.e. the preferential selection of 1in/2out or 2in/1out low-energy spin configurations (see Fig. 1a), is globally well obeyed: over all magnetic configurations, less than 0.6% of the vertices have a 3in or 3out high-energy spin configuration. These "forbidden" vertices are likely due to structural defects in the networks and often appear at the same locations. Residual magnetization after demagnetization is low, of the order of a few percents, consistent with refs. 11 and 12. From the analysis of our magnetic images, we deduce the spinspin correlation coefficients C ij = S i . S j between normalized S i and S j spins (indices are relative, identical to those used in previous literature [6,13], and defined in Fig. 1b). The experimental values, averaged over the 34 demagnetized networks, are reported in Table I for the 7 nearest neighbors. These values are compared to the coefficients expected from the SI model, described by the Hamiltonian H SI = -J αβ ij S i . S j that couples nearest neighbors only, through the J αβ coefficient [13,14]. Although differences are visible, as previously observed [6], experimental results and theoretical predictions from the SI model agree relatively well at first sight.

Considering the nanomagnets as Ising spins [15], we can derive from the C ij values the network energy

E = -(1/N ) i =j N ij J ij C ij .
In this expression, N is the total number of spins, N ij and J ij being respectively the number of interactions and coupling coefficients of type ij (Fig. 1b). Within the point-dipole approximation, the coupling coefficients J ij can be calculated: for two magnetic moments µ S i and µ S j located at r i and r j :

J ij = -D ( 1 | r ij | 3 -3 ( S i . r ij )( S j . r ij ) | r ij | 5 S i . S j
), with r ij = r i -r j , and D = (µ 0 /4π)µ 2 the strengh of the dipolar coupling. The corresponding J ij values are reported in Table I and compared to the values obtained using micromagnetic simulations [START_REF]OOMMF User's Guide, Version 1[END_REF]. Good agreement is found between the two approaches, except for first neighbor interactions, for which the minimal distance between the nanomagnets is small compared to their extension.

The system energy can be split in two terms: E = E short + E long where E short is the energy related to the interaction between first neighbors only (E short = -N αβ J αβ C αβ /N ), and E long the energy related to the interactions beyond the nearest neighbors. We will show that, once the system enters the spin ice manifold, C αβ is configuration-independent, i.e E short is constant. In the following, we thus focus our discussion on the E long term. The presence of a few "forbidden" vertices in our networks is then not critical as E long excludes energy contributions from nearest neighbors. The experimental E long distribution is shown in Fig. 2a, and compared to the distribution expected from the SI model (dotted line). Interestingly, this distribution is shifted towards negative energy. To understand where energy is gained in our networks, we sort the 34 measured configurations in 3 subsets corresponding to the lowest (red), intermediate (green) and highest (blue) energy configurations (Fig. 2a), and we average the C ij coefficients within these 3 subsets (Fig. 2b). Two main features are then observed: i) independently of the configuration energy, the C αν and C αη coefficients have higher values than expected (also see Table 1), and ii) the values of the C αγ and C αδ coefficients strongly vary with energy, contrary to what is observed for the 5 other measured coefficients.

As depicted in Fig. 1b, the C αν and C αη coefficients reflect correlations along a line in the network. We might think that these higher correlation coefficients originate from the demagnetization procedure that may favor ferromagnetically coupled nanomagnets along these lines. However, although the system has threefold symmetry, the three main directions of lines are not equivalent, and we observe anisotropy in the C αν and C αη values (one direction is different from the two others). We thus attribute the high C αν and C αη values to variations of coupling coefficients between first neighbors, presumably due to anisotropy in the fabrication process. Replacing the theoretical correlation coefficients along a line by the experimental ones, keeping all the other coefficients unchanged, leads to a larger shift of the distribution towards negative energies (continuous line in Fig. 2a). This new energy shift reproduces well our data. Thus, the total energy shift we observe experimentally has two main origins: an intrinsic shift associated with the ice model itself (E long is lower in the SI manifold than in the paramagnetic state in which all spin configurations are equiprobable [START_REF]We show here a phenomenon[END_REF]), and an extrinsic shift due to the lithography process. Both phenomena are related to first neighbor interactions, and do not involve long range coupling.

In contrast with the other correlation coefficients, the values of C αγ and C αδ strongly vary with E long (see Figs. 2b). For several configurations, C αγ and C αδ reach absolute values that can be surprisingly large compared to the values expected from the kagomé SI model. Magnetostatic coupling between the (α, γ) and (α, δ) spins favors antiparallel alignment of the corresponding nanomagnets, and thus negative values of their spin-spin correlators. Intuitively, this could explain why the measured C αγ and C αδ values are not in quantitative agreement with the values expected from the kagomé SI model.

To test whether this interpretation better describes our findings, monte carlo simulations of the kagomé DSI model have been performed, using the following Hamiltonian:

H DSI = -J 1 ij S i . S j -i = J ij S i . S j
, where J 1 is an additional isotropic nearest neighbor coupling. This J 1 term can be tuned so that the effective nearest neighbor interaction J 1 + J αβ = 5 J αβ , to account for the spatial extension of the nanomagnets (see Table 1). The simulations were done for L × L × 3 sites lattices (L ranging from 11 to 144) with periodic boundary conditions, using a single spin flip algorithm and a simulated annealing procedure from T /J αβ = 100 to T /J αβ = 0.1. In these simulations, 10 5 modified monte carlo steps (mmcs) are used for thermalization, where 1 mmcs corresponds to a set of local updates sufficiently long to achieve stochastic decorrelation, itself determined using the on the fly calculated spin-spin autocorrelation time. Measurements follow the thermalization process and are computed with 10 5 mmcs. Results of the simulation are shown in Fig. 3 where spin-spin correlators are plotted as a function of temperature. In the high temper- ature regime, the system is in a paramagnetic state. As temperature is lowered, the system enters the short range kagomé spin ice manifold (C αβ = 1/6), and the 7 measured correlation coefficients nearly reach the values expected from the SI model. We emphasize that in absence of long range interaction, these values would then remain temperature-independent (see inset of Fig. 3). The situation is different when dipolar interactions are included in the calculation. At very low temperatures, a transition to an ordered phase takes place [9,10]: all correlation coefficients vary abruptly (not shown here), and strongly deviate from the SI values. In the intermediate regime, i.e when 2 ≤ T /J αβ ≤ 5 (see Fig. 3), mean values of the correlation functions do not allow easy discrimination between the SI and DSI models. However, two correlation coefficients slowly decrease when lowering the temperature (see the two lowest curves). These two coefficients, C αγ and C αδ , are precisely those for which we observe a change experimentally, when E long decreases. This strikingly contrasts with predictions from the SI model and fits well our experimental data.

An alternative description of these spin models is useful. Identifying each spin S i to a magnetic dipole allows the mapping of any spin configuration on the kagomé lattice onto a classical magnetic charge configuration on the honeycomb lattice (Fig. 1a). The high temperature paramagnetic phase of the SI and DSI models then translates into a paramagnetic charge state. When temperature is lowered, both models enter the spin ice manifold: the kagomé ice rules are satisfied, and charge values are constrained to ±1 (instead of ±1 and ±3 at high temperature), therefore corresponding to an emergent Z 2 field. In the SI model, this Ising field remains weakly correlated and temperature independent. The DSI model however behaves differently: at lower temperatures, it leaves a first spin ice manifold (SI-I) and enters a more constrained phase (SI-II), in which spins still fluctuate but where charges crystallize in an antiferromagnetic pattern [10]. During the phase transition of the hidden Z 2 field, the stochastic dynamics is strongly modified, with exponentially suppressed single spin flips, leaving only loop-like moves to operate efficiently in the SI-II manifold. We have thus investigated the temperature dependance of the short range magnetic charge correlation Q i .Q i+1 , and compare monte carlo predictions from the SI and DSI models (Fig. 4a) to our measurements (Fig. 4b). In the highest temperature regime, where all charge states are accessible (Q=±1 and Q=±3), the unconstrained paramagnetic charge phase is associated with an asymptotic correlator Q i .Q i+1 = -1 (Note that a Q=±3 charge state corresponds to a monopole defect, see for example refs. 18 and 19). When the temperature is lowered, the SI model enters the spin ice manifold, with constrained Q = ±1 charges and a temperature independant short range charge correlation Q i .Q i+1 =-1/9. In the DSI model however, the charge correlator first increases to reach a maximum value around -0.15, and then steadily decreases with temperature. Note that there exists a temperature range within the spin ice manifold, in which the mean values of the correlator, as well as their associated standard deviations, are mutually exclusive for the two models (Fig. 4a). Analysis of our X-PEEM images reveals strong deviations from the -1/9 value of the SI model (Fig. 4b), clearly outside standard deviations of the correlator distributions, and agrees well with predictions from the DSI model. We emphasize that in all 34 networks, the spin configurations are charge neutral, which rules out a global extrinsic breaking of the Z 2 symmetry. Therefore, we attribute the deviations of the short range charge correlations to the proximity of the SI-II regime, in which the Z 2 and translation symmetries are broken. Unambiguously, our results reveal that dipolar interactions are at work in our arrays.

In summary, numerical computation and magnetic imaging show that the kagomé SI model fails to interpret our experiments, while the kagomé DSI model captures the physics we observe. The fact that artificial arrays of nanomagnets may look like short range kagomé spin ice systems at first sight is fortuitous, and related to the critical slowing down at low temperature of the single spin flip dynamics. Therefore, while artificial kagomé arrays of nanomagnets should not exhibit finite entropy at zero temperature [9,10], they are most likely described by a frozen dipolar spin ice regime with an emergent symmetry breaking of the hidden Z 2 field associated with the classical magnetic charges.

FIG. 1 :

 1 FIG. 1: (Color online) (a) 15 µm X-PEEM magnetic image of a demagnetized network. Direction of magnetization is deduced from the black and white contrast, as illustrated for two vertices that satisfy the ice rule, i.e. the local 1in/2out (magnetic charge Q=-1) or 2in/1out (magnetic charge Q=+1) spin configuration. (b) Sketch of the array indicating the relative indices used for the 7 nearest neighbors. The partly superimposed dots highlight the nodes of the kagomé lattice.

FIG. 2 :

 2 FIG. 2: (Color online) (a) Experimental histogram deduced from the 34 demagnetized networks. Color code is used to sort the configurations in 3 subsets. Curves show the normalized E long distributions from the SI model (dotted line) and after replacing the expected Cαν and Cαη coefficients by their measured values (continuous line). (b) Experimental correlation coefficients averaged within each energy subset.

FIG. 3 :

 3 FIG. 3: (Color online) Correlation coefficients Cαi derived from the DSI model as a function of temperature T . Greek index i = β, γ, . . . identifies each correlator. Inset : temperature dependance of correlators within the SI model.

  FIG. 4: (a) Theoretical short range magnetic charge correlators Qi.Qi+1 (bold line) and their standard deviations (hatched regions) calculated from the SI and DSI models as a function of temperature. (b) Experimental Qi.Qi+1 values measured for the 34 demagnetized networks as a function of E long . (c) Histogram of the experimental charge correlators compared to the nearly gaussian distribution computed in the SI model, within the spin ice manifold.

TABLE I :

 I Theoretical (SI model) and experimental spin-spin correlators, number of ij interactions in our networks, and coupling coefficients in arbitrary unit deduced from the pointdipole approximation and from micromagnetic calculations.

		αβ	αγ	αν	αδ	ατ	αη	αφ
	Cij (theo.) 0.167 -0.062 0.101 -0.075 0.012 0.019 0.023
	Cij (exp.) 0.164 -0.056 0.151 -0.063 0.013 0.056 0.021
	Nij	648	612	612	303 1152 576 273
	Jij (dip.)	1	-0.137 0.045 -0.036 0.014 0.037 0.014
	Jij (µmag.) 4.905 -0.130 0.044 -0.030 0.013 0.040 0.015


Céline Moriou

Adilson D Da Silva

Marcos Joel

Vianelli Prado

§,|| Clément Denhez

Oleksandr Plashkevych

Jyoti Chattopadhyaya

Dominique Guillaume

Pascale Clivio

C2'-F Stereoconfiguration as a puckering switch for base stacking at the dinucleotide level

ABSTRACT: Fluorine configuration at C2' of the bis(2'-fluorothymidine) dinucleotide is demonstrated to drive intramolecular base stacking. 2'-β Fconfiguration drastically reduces stacking compared to the 2'-α series. Hence, base stacking emerges as being tunable by the C2'-F stereoconfiguration through dramatic puckering variations scrutinized by NMR and natural bond orbital analysis. Accordingly, 2'-β F-isomer photoreactivity is significantly reduced compared to that of the 2'α-F-isomer.

C2'-F nucleosides are therapeutically important compounds used as anticancer (gemcitabine, clofarabine) or antiviral (clevudine, sofosbuvir) agents. 1 Several radiolabelled 18 F-C2' nucleosides are clinical candidates currently evaluated to be used in positron emission tomography as radiotracers to image cancer, autoimmune or inflammatory diseases. 2 DNA-or RNA-oligonucleotides (ONs) containing 2'-F nucleosides are also clinically investigated. [START_REF]F = 14 Hz, C1' Tβ F p-), 85.1 (d, J 4',F or J 4',P = 11 Hz, C4' Tβ F p-), 82.1 (C4' -pTβ F ), 77.6 (d, J 3',F = 26 Hz, C3' Tβ F p-), 75.1 (d, J 3',F = 27 Hz, C3' -pTβ F ), 66.5 (C5[END_REF] For the treatment of age-related macular degeneration, the aptamer pegaptanib is commercialized 4 whereas ARC1905 is in clinical trial. 3a RB006, an aptamerbased anticoagulant is in Phase II, 3a as is inclisiran, a siRNA presently evaluated as a cholesterol-lowering drug. 5 Actually, most of siRNAs in clinical trials today contain a F-atom at some of their C2'-position. [START_REF]F = 14 Hz, C1' Tβ F p-), 85.1 (d, J 4',F or J 4',P = 11 Hz, C4' Tβ F p-), 82.1 (C4' -pTβ F ), 77.6 (d, J 3',F = 26 Hz, C3' Tβ F p-), 75.1 (d, J 3',F = 27 Hz, C3' -pTβ F ), 66.5 (C5[END_REF]6 In parallel, damaged ONs incorporating a 2'-F sugar residue at their altered base site serve as bio-chemical tools to study the recognition and cleavage reaction mechanism of glycosylases involved in DNA repair. 7 Finally, 2'-F ONs are also used in structural biology to manipulate DNA i-motifs 8 and Gquadruplex folding. 9 Conveniently, secondary structure transitions can be monitored by 19 F NMR. 10 C2'-F nucleoside sugar conformation is overwhelmingly driven by its C2' configuration. In solution, 2'-α substitution leads to a predominant North (N, C3'-endo) conformation 11a,b whereas 2'-β substitution generally leads to a major South-East (SE, C1'-exo) conformation. 11c-e Such trend is also observed at the oligonucleotide (ON) level in a double strand context where the 2'-β sugar pucker oscillates between South and East (O4'-endo). 10b,12 Actually, 2'-α-F and 2'-β-F ONs are considered as RNA and DNA mimics, respectively. 13a Divergent observations have been reported regarding the influence of the C2' configuration of 2'-F ONs on duplex 12a,13 or G-quadruplex stability. 9b-d This may be due to the type and number of 2'-F nucleoside residues incorporated, as well as the length and sequence context. Consequently, it remains still difficult to anticipate the influence of the C2'-F stereoconfiguration on ON structure/stability. Probing the effect of the C2'-F stereochemistry on very simple models would certainly help to unravel the intrinsic basic and minimal conformational behavior of the sugar phosphate backbone scaffold, and nucleobase stacking.

Recently, we have demonstrated that the doubly 2'-α-labeled F analogue (1, Tα F pTα F ) of the dinucleo-tide of thymine (2, TpT) exhibits a major (>92%) N pucker at both ends. 14 Concomitantly, 1 was found much more sensitive to UV light than 2. 14 Such high photochemical reactivity was attributed to a high, and favorably driven by the N sugar pucker conformation, propensity to intramolecular base stack.

In order to evaluate the importance of the stereochemistry of the fluorinated atom on ON stacking at the dinucleotide level, we have synthesized Tβ F pTβ F (3), the C2' epimer of 1, and investigated the influence of its sugar conformation on its intramolecular stacking. For this purpose, analysis of 3 J (PSEUROT+J HF ) 15 and high-level ab initio DFT calculations have been performed. We have then studied the intramolecular base stacking behavior of 3 using CD and photochemistry.

Synthesis of 3 started from 2'-β-fluorothymidine (Tβ F , 4) accessible through a Vorbrüggen coupling step, according to the procedure of Wilds & Damha. 11c Our synthetic route closely followed the strategy reported for the synthesis of 1 14 using known intermediates 5 and 7 [START_REF] Kois | [END_REF] (Scheme 1). The conformational population distribution of each sugar of 3 has been estimated using the temperature dependent analysis by PSEUROT+J HF 15 (Table 1, S17-S26). Sugar moieties in 3 adopt predominantly SE conformation at the 5′-end (66%) and a near-equal distribution of SE (51%) and N (49%) conformations at the 3′-end. The pseudorotational torsion angle (P) and puckering amplitude (Ψ m ) of SE and N conformers are P SE =~130°, Ψ m =35°; and P N =-11°, Ψ m =34°; respectively which are consistent with those obtained as the results of gas phase optimizations of 3 (Tables S1,S5) calculated at different levels of theory and sizes of the basis set. Correlation between HF/6-31G** total energy and endocyclic torsions from geometry optimizations along C1'-C2'-C3'-C4' (5'-end) rotation pathway is shown Figure S19. NMR data have shown surprisingly weak dependence of 3 J HH and 3 J HF on the temperature (Table S2), suggesting an entropy controlled process. Interestingly, the temperature insensitivity of the sugar of Tβ F (4) has been reported earlier. 11b The SE conformation calculated for each sugar residue of 3 is consistent with the preferred sugar conformation depicted for 4 in solution. 11d Therefore, it can be affirmed that 1 and 3 present a drastically different preferred sugar conformation. Additionally, puckering preference of 1 and 3 has been studied by high-level ab initio DFT methods (IEFPCM/M06-2X/6-31G(d,p)) using 2'-F-thymidine 3'-methylphosphates 9 and 10 as models (Tables S4,S5).

Table 1. Experimental and PSEUROT+HF calculated NMR data of the sugar moieties of 3.

3 J Jexp Jcalc PN ° ΨN ° PS ° ΨS ° XN% RMSHH RMSd ΔG 298,a ΔH a ΔS b
Tβ F p-J 1′2′ 3.9 3.9 a : in kJ.mol -1 ; b : in J.mol -1 .K -1 ; c : Estimated values.

-
Low-energy calculated conformations of 9 and 10 have been found to be in accordance with NMRderived sugar conformational preferences of 1 and 3, respectively (Table S4). Orbital interactions, which are important factors to drive the sugar puckering, 15b have been examined by natural bond orbital (NBO) analysis (Tables S6,S7).

For 9 and 10, the nO4'→σ*C1'-N1 anomeric effect involving one lone pair orbital of O4' has been calculated to be the major contributor to the stability of the N conformers. To a lesser extent, the S conformers have been found to be stabilized by two other anomeric effects; nO4'→σ*C1'-N1 involving the other O4' lone pair and nN1→σ*C1'-O4'. Sugar pucker has also been found driven by antiperiplanar effects (APEs). 15b Most prevalent APEs (E>12 kJ/mol) found for the N and S conformers of 9 and 10 (Table 2) are highly consistent with those reported for conformers of 2'-F nucleosides. 11d,17 Importantly, the dihedral angles formed between X and Y of the fragment XCCY found in major APEs correspond to a trans geometry domain (Table 2). Such stereochemical arrangement induces a maximal σ/σ* orbital overlap geometry, stabilizing the conformer concerned. The correlation observed between these trans dihedral angles and stabilizing energy contributions (Table 2) supports our calculations. Interestingly, concomitantly and independently to our study, APEs of several nucleosides including 2'-deoxy-2'fluorouridine and 4 have recently been reported, validating our study. 18 In a 2',2'-difluorinated nucleoside, the existence of hydrogen bonds between the 2'αF atom and H4' (N conformer) and 2'βF atom and H6 (S conformer) has been suggested as conformational stabilizing factor. 17a In our case, corresponding H … F distances: 2.66 Å (N conformer of 9), and 2.75 Å (S conformer of 10) satisfied H-bond requirements. However, the participation of such electrostatic attractions was discarded for lack of directionality 19 (∠C-H … F=78.8° and 97.9° for 9 and 10, respectively). CD melting spectra of 3 have been recorded to get an insight on its intramolecular stacking status (Figure S1). Absence of significant signal variation upon heating and similarity of these spectra with that of 4 provide definitive evidence that 3 is predominantly unstacked in aqueous solution at room temperature. It is in sharp contrast with the stacking ability of 1 14 calculated to be at least 50% stacked by two different curve fitting methods (Table S3). 20 Unfortunately, mathematical treatment of melting CD data of 3 could not be further exploited due to insufficient amplitude variation. However, assuming that 2 and 3 share a similar stacked geometry, the stacked population of 3 could be estimated to be a half lower (ca 10%) than that of 2 21 by comparing their CD diff spectrum 22 (Figure S2). Because CD studies clearly demonstrated that 3 is nearly 90% unstacked, in sharp contrast with 1, 14 and since stacked thymine dinucleotides are considered to be prone to afford photoproducts (PPs) upon 254 nm exposure, 14,[23][24][25] we have investigated the photoreactivity of 3 to correlate it with the 2'-F configuration.

In a previous study, 14 we observed that UV exposure of 1 (Tα F pTα F ) leads to the cis-syn CPD and (6-4) PP with a quantum yield (Φ) calculated to be (2.7±0.3)x10 -2 and (0.13±0.02)x10 -2 , respectively. UV irradiation of 3 at 254 nm gave rise almost exclusively to 11, identified as the cis-syn CPD (Scheme 2), together with traces of 12 (not represented) tentatively assigned to a (6-4) PP from its UV absorption at 329 nm. However, this latter PP could not be isolated in preparative scale. Its Φ was estimated 24 to be (0.067±0.002)x10 -2 which is significantly lower than that of TpT (Φ (6-4) =(0.10±0.05)x10 -2 ) 26 and 1. 14 Scheme 2. Photochemical behavior of 3.

Photolysis and PP formation kinetics of 3 have been studied with respect to those of 2 (Figure 1). Results indicate that 2'-β-F substitution drastically lowers the rate of radiation-induced dinucleotide disappearance as well as the initial rate of CPD formation (Figure 1). Φ of 11 formation has been calculated to be (0.48 ± 0.08)x10 -2 . This formation efficiency is a half of that reported for TpT-derived CPD (Φ CPD =(1.1±0.05)x10 - 2 ), 26 in contrast to what had been observed for 1. 14 The low photoreactivity of 3 can be explained by its low stacking propensity, and consequently the limited formation of the conformation state necessary to allow PP formation. Interestingly, CD studies have shown that 3 was half stacked compared to 2. Table 2. Major APEs contributing to conformer stability in 9 and 10 from NBO analysis (corresponding dihedral angle in ° and energy in kJ/mol).* 18 (156, 17) σC2'-H2' → σ*C3'-O3' 18 (-167, 17) σC3'-C4' → σ*C2'-F 17a, 18 (-158, 18) σC2'-H2" → σ*C1'-N1 (-144, 12) σC1'-C2' → σ*C3'-O3' 17a,18 (154, 16) σC3'-C4' → σ*C2'-F 17a, 18 (155, 17) σC1'-H1' → σ*C2'-F 17,18 (164, 17) σC2'-H2" → σ*C1'-O4' 11d,17b (162, 12) * Bonds involved in APEs are colored and are highlighted in bold (donors) and dashed (acceptors). The origin of this quantitative stacking difference might come from the conformation of the sugar puckers as well as from steric hindrance introduced by the substituent. It is known that a 5'-end N sugar conformation favors intramolecular stacking and a 3'-end N sugar conformation strengthens this effect. 22,27 Assuming that 1-3 present a sufficiently similar stacking geometry, this might explain why the 2'-α-F diastereomer is much more stacked than the 2'-β-F diastereomer. That 3 is less stacked than 2, despite a higher N population, indicates that the 2'β-F substitution per se impedes the stacking phenomenom.

9 10 N conformer S conformer N conformer S conformer σC3'-H3' → σ*C2'-F 17a,18 (168, 18) σC1'-C2' → σ*C3'-O3' 17a,
In conclusion, we have evidenced the dramatic effect of the C2' F configuration on intramolecular base stacking. This different stacking behavior is likely driven by the specific puckering preference of each 2'-F diastereomers, consistent with the interplay of specific stereoelectronic effects. When embedded in a double stranded ON context, these 2'-F motifs have already been shown to have the ability to modify duplex stability. 12a More interestingly, incorporated in a ss ON context, where no base pairing and adjacent base stacking occur, 2'-F-dinucleotide residues might have the capacity to tune intramolecular stacking at their site of incorporation. Such behavior is of utmost importance and could be useful for manipulating the topology/folding of functional DNA-or RNA-containing loops such as ss DNA or RNA ONs, 28 or aptamers; loops being often an element of stabilization of the structure or a recognition element. 29 The latter is supported by the recent incorporation of 4 or 2'-F purine residues in loops of aptamers to provoke a preorganization of the loop for better ligand binding. 9f,30 EXPERIMENTAL SECTION General Remarks. Solvents and chemicals used for the reactions were purchased from commercial suppliers.

Dichloromethane, acetonitrile, N,Ndiisopropylethylamine (DIPEA), ethyl acetate and heptane were dried by distillation from calcium hydride. Pyridine and triethylamine (TEA) were dried by distillation from KOH and kept over KOH. TLC was performed on silica gel plates (Kieselgel 60, F 254 ) with detection by UV and visualization by spraying with a methanolic solution of sulfuric acid. Chromatography was performed on silica gel 60, particle size 35-70 µm, unless otherwise stated. Medium pressure reverse phase chromatography was performed on LiChroprep RP 18. 1 H NMR and 13 C NMR spectra were recorded on 300, 500 or 600 MHz spectrometers. Observed chemical shift (δ) values are given in ppm and coupling constants (J) in Hz. Following abbreviations are used: s (singlet), br s (broad singlet), d (doublet), t (triplet), q (quartet), dd (doublet of doublet), Tβ F p-and -pTβ F represent the 5'-end and the 3'-end nucleoside residues, respectively, of 3 and 11. 1 H NMR and 13 C NMR chemical shifts recorded in D 2 O were calibrated using residual solvent signal (δ H 4.80 ppm) and from dioxane (δ C 67.8 ppm), respectively. For CDCl 3 δ H and δ C are reported from internal TMS (δ 0.00). 31 P NMR and 19 F NMR chemical shifts were reported from an external capillary standard of 85 % phosphoric acid (δ P 0.00 ppm) and of TFA (δ F -77.00 ppm), respectively. Superscript a and b labels denote interchangeable assignments within a compound. High Resolution Mass Spectra (HRMS) were determined on a Liquid Chromatography Time-of-flight (LCT) spectrometer. Electron Spray Ionization Mass Spectra (ESI-MS) were measured on a LC/MS mass spectrometer. 5'-O-(dimethoxytrityl)-2'-β-fluorothymidine (5). This compound was prepared as previously reported [START_REF] Kois | [END_REF] in 85% yield from 4. 1 H NMR (300 MHz, CDCl 3 ): δ 9.23 (1H, br s), 7.50-7.15 (10H, m), 6.83 (4H, d, J = 3.3 Hz), 6.26 (1H, dd, J = 3.3, 19.5 Hz), 5.05 (1H, dd, J = 3.0, 52.2 Hz), 4.44 (1H, dd, J = 3.7,19.4 Hz), 4.02 (1H, m), 3.78 (6H, s), 3.42 (2H, m), 1.74 (3H, s). 13 C NMR (75 MHz, CDCl 3 ): δ 164. 0,158.9,150.6,144.8,137.2,136.0,135.9,130.3,128.3,128.2,127.3,113.5,110.6,95.5 (d,J = 192 Hz),86.8,84.0 (d,J = 16.5 Hz),83.4,75.9 (d,J = 25.5 Hz),62.9,55.5,12.6. 19

Scheme 1 .

 1 Scheme 1. Synthesis of 3 a .

Figure 1 .

 1 Figure 1. Fractional amount of 3, 2 and of their respective PPs as a function of irradiation time.

  F (282 MHz, CDCl 3 ): δ -197.82. HRMS (ESI + mode) (M+Na) + : calc. for C 31 H 31 O 7 N 2 FNa 585.2013, found 585.2031. 3'-O-Acetyl-2'-β-fluorothymidine (6). This compound was obtained as a white foam in 75% yield from 5 following a similar acetylation/detritylation procedure TOC

ACKNOWLEDGMENT

We thank O. Rakotonandrasana and M.-T. Martin for recording some NMR data and S. Coantic-Castex for confirming NMR values and refining photochemical data. Financial support from the Swedish Natural Science Research Council (Vetenskapsrådet) and EU FP7 programme (JC), and computational time at the National Supercomputer Center (Linköping, Sweden) awarded by the Swedish National Allocation Committee for the High Performance Computing (OP) are gratefully acknowledged, as are the brazilian agencies CNPq, FAPEMIG and CAPES.

ASSOCIATED CONTENT

* Supporting Information NMR and CD spectra, PSEUROT+J HF data, theoretical analysis of CD data, and quantum chemical calculations. These materials are available free of charge via the internet at http://pubs.acs.org.

AUTHOR INFORMATION

Corresponding Author pascale.clivio@univ-reims.fr

Notes

The authors declare no competing financial interest.

For convenience, the name 2'-β-fluorothymidine and those of its derivatives are preferred instead of 1-(2deoxy-2-fluoro-β-D-arabinofuranosyl)thymine. Similarly, the name 2'-α-fluorothymidine and those of its derivatives are preferred instead of 2'-deoxy-2'fluoro-5-methyluridine.


Otman Maalal

Marc Prat

René Peinador

Didier Lasseux

Evaluation of pore size distribution via fluid-fluid displacement porosimetry: The viscous bias

Keywords:

Fluid-fluid displacement porosimetry (FFDP) is a method for determining the pore size distribution of a porous medium from viscous fluid-fluid immiscible displacement (drainage) results, using a simplistic representation of the pore structure by cylindrical parallel tubes. The method leads to unambiguous estimates when the porous medium microstructure is strongly anisotropic and reasonably akin to this representation. The information obtained from this technique is however less straightforward to interpret for more complex microstructures. On the basis of synthetic data obtained from pore network simulations, it is evidenced, in particular, that this method leads to a shift in the distribution toward unrealistic small pore sizes. The shift, that may be very significant, is due to the viscous pressure drop in the displacing fluid. It is shown that significant improvement is achieved when the characterization of the pore size distribution is considered as an inverse optimization problem. The latter is solved using an optimization method presented in a previous work that was restricted to situations where viscous effects are unimportant.

Introduction

Fluid-fluid displacement porosimetry (FFDP) is a technique for characterizing the pore size distribution of a porous medium [START_REF] Morison | A comparison of liquid-liquid porosimetry equations for evaluation of pore size distribution[END_REF][START_REF] Peinador | Characterization of polymeric UF membranes by liquid-liquid displacement porosimetry[END_REF][START_REF] Calvo | Porosimetric characterization of polysulfone ultrafiltration membranes by image analysis and liquid-liquid displacement technique[END_REF][START_REF] Tanis-Kanbur | Membrane characterization via evapoporometry (EP) and liquid-liquid displacement porosimetry (LLDP) techniques[END_REF]. In particular, the method is commonly used for the characterization of thin porous media such as ultrafiltration and microfiltration membranes (see e.g. [START_REF] Peinador | Characterization of polymeric UF membranes by liquid-liquid displacement porosimetry[END_REF][START_REF] Calvo | Porosimetric characterization of polysulfone ultrafiltration membranes by image analysis and liquid-liquid displacement technique[END_REF][START_REF] Tanis-Kanbur | Membrane characterization via evapoporometry (EP) and liquid-liquid displacement porosimetry (LLDP) techniques[END_REF], and is traditionally referred to as the Liquid-Liquid Displacement Porosimetry (LLDP) technique since the two fluids commonly used are liquids. However, the technique can be extended to Gas-Liquid Displacement Porosimetry [START_REF] Islam | Mixed Poiseuille-Knudsen flow model for gas liquid displacement porometry data treatment[END_REF][START_REF] Peinador | Comparison of capillary flow porometry (CFP) and liquid extrusion porometry (LEP) techniques for the characterization of porous and face mask membranes[END_REF] and can be simply referred to as the FFDP technique. In this method, the membrane is saturated by a wetting liquid and another immiscible liquid or gas is forced through the pores. This is accomplished through a series of pressure steps resulting in the immiscible displacement of the wetting liquid by the injected fluid (a drainage process). The injected fluid pressure at the membrane inlet is measured, as well as the injected fluid flow rate through the membrane. The dependence of these two quantities upon the pore size distribution lies in the local capillary pressure invasion threshold which must be overcome for a displacement to occur in a pore, and also in the dependency of the local flow conductances on the pore size. However, as pointed out earlier [START_REF] Morison | A comparison of liquid-liquid porosimetry equations for evaluation of pore size distribution[END_REF], practically all the determinations of the pore size distribution with the FFDP technique have been made considering a simplified representation of the membrane as a bundle of parallel cylindrical tubes. As analyzed in some reported works [START_REF] Mourhatch | Determination of the true pore size distribution by flow permporometry experiments: an invasion percolation model[END_REF][START_REF] Maalal | Determination of the throat size distribution of a porous medium as an inverse optimization problem combining pore network modeling and genetic and hill climbing algorithms[END_REF], this is quite a questionable simplification since the microstructure of most membranes significantly differs from this simplistic representation. Owing to this simplification, the question arises as to whether this technique really allows determining the pore size distribution in thin porous layers which structure departs from a simple system of cylindrical pores in parallel. The main objective of the present article is hence to clarify this question. To this end, the FFDP is simulated numerically for given microstructures. This means that the pore size distribution is perfectly known for each microstructure under consideration, without assuming cylindrical tubes in parallel. Then, the pore size distribution (PSD) is determined from the FFDP numerical simulation, i.e., from the simulated two-phase immiscible displacement in the microstructure of interest, as it would be obtained from a FFDP real experiment. The PSD estimated with this technique and the simplistic pore structure representation can finally be compared to the real PSD in order to evaluate the possible bias and errors associated with the FFDP. According to these simulations, significant discrepancies are found between the PSD extracted from the FFDP and the real PSD. In order to improve the PSD estimation from the FFDP raw data, the PSD extraction model, based on the parallel cylindrical tubes simplification, is abandoned and a more sophisticated approach is developed. It consists in considering the PSD determination as an inverse optimization problem, which is solved using the method presented in [START_REF] Maalal | Determination of the throat size distribution of a porous medium as an inverse optimization problem combining pore network modeling and genetic and hill climbing algorithms[END_REF]. This approach is tested for various PSD, namely, uniform, Gaussian, Log-normal and bimodal. Both steps, i.e., the evaluation and assessment of the commonly used FFDP approach and the inverse optimization approach, are based on numerical simulations of the immiscible fluid-fluid displacement in the pore space of the thin porous medium under consideration. The latter are performed within the classical framework of pore network modeling [START_REF] Blunt | Flow in porous media-pore-network models and multiphase flow[END_REF][START_REF] Blunt | Detailed physics, predictive capabilities and macroscopic consequences for pore-network models of multiphase flow[END_REF][START_REF] Xiong | Review of pore network modelling of porous media: experimental characterizations, network constructions and applications to reactive transport[END_REF][START_REF] Bondino | Relative permeabilities from simulation in 3D rock models and equivalent pore networks: critical review and way forward[END_REF]. Typically, pore network models rest on a yet simplified representation of the pore space as a system of nodes connected by narrower channels. The nodes are local larger volumes in the pore space, which, within the framework of pore network modeling, are referred to as pore bodies. The narrower channels between pores are referred to as pore throats or simply throats. In what follows, simple two-dimensional square pore networks are considered.

As illustrated in Fig. 1, the pore bodies are located at the nodes of a square regular grid and the throats correspond to the segments between nodes. The morphology of pore networks extracted from real porous media microstructures is obviously more complex than a regular grid, (see e.g. [START_REF] Agaesse | Validation of pore network simulations of ex-situ water distributions in a gas diffusion layer of proton exchange membrane fuel cells with X-ray tomographic images[END_REF][START_REF] Raeini | Generalized network modeling: network extraction as a coarse-scale discretization of the void space of porous media[END_REF]. However, it is expected that the crucial features for assessing the FFDP are the pore space interconnectivity (the fact that a pore body is connected to several neighboring pore bodies) and the disordered nature of the pore space (the fact that the local geometrical dimensions vary from one pore body or throat to another). These two characteristics are present in the simple pore networks considered in this study so that the main conclusions reached from this approach are expected to also apply to unstructured networks representative of the pore space of real porous media.

Within the pore network approach, the throats and pore bodies shall be considered as the two microsctructure elements to which a size distribution can be affected, namely the pore body size distribution (PBSD) and the throat size distribution (TSD), respectively. The constrictions in the pore space, which correspond to the throats, play an important role since the menisci tend to get pinned in the constrictions during a drainage process. This results from capillary effects. Moreover, the local viscous pressure drops are highly dependent on the flow through the throats (the narrower the passage, the greater is the pressure drop between two pore bodies for an identical flow rate). For these reasons, the focus is laid here on the determination of the TSD.

The present study can be seen as an extension of the works presented in [START_REF] Mourhatch | Determination of the true pore size distribution by flow permporometry experiments: an invasion percolation model[END_REF] and [START_REF] Maalal | Determination of the throat size distribution of a porous medium as an inverse optimization problem combining pore network modeling and genetic and hill climbing algorithms[END_REF]. In the latter reference, the immiscible displacement was simulated using the standard invasion percolation (IP) algorithm as described for instance in Wilkinson and Willemesen (1983). This corresponds to a situation where the displacement is controlled by capillary effects only. However, in a FFDP experiment, the fluid pressure is not spatially uniform in the displacing phase but decreases on average from the inlet to the outlet. For this reason, the displacement is expected to differ from that predicted by the standard IP algorithm and the flow situation corresponding to the FFDP requires a specific study. Consequently, the viscous pressure drop in the displacing fluid is taken into account, contrary to the study presented in [START_REF] Maalal | Determination of the throat size distribution of a porous medium as an inverse optimization problem combining pore network modeling and genetic and hill climbing algorithms[END_REF]. It is thus important to notice that the situation considered in the present work is much closer to operating conditions of currently available FFDP devices than the situation considered in [START_REF] Maalal | Determination of the throat size distribution of a porous medium as an inverse optimization problem combining pore network modeling and genetic and hill climbing algorithms[END_REF]. As shall be seen, the existence of the viscous pressure drop leads to an important bias in the determination of the pore size distribution if not considered in data interpretation. The identification and the illustration of this viscous bias is the major objective of the present work together with a method for alleviating its impact.

The present article is organized as follows. In Section 2, the method commonly used for determining the TSD from the FFDP raw data is recalled. The pore network model is presented in Section 3. The algorithms used to simulate the displacement and compute synthetic FFDP data is presented in Section 4. Section 5 is dedicated to the assessment of the TSD determination from the commonly used FFDP method using a parallel cylindrical tubes model. In Section 6, a method to determine the TSD from the FFDP curve using PNM simulations combined with an optimization method is briefly presented. Section 7 reports on results obtained with the optimization method when both the FFDP and retention curves are used as target data. A short discussion is presented in Section 8. Conclusions are drawn in Section 9.

Fluid-fluid displacement porosimetry

As indicated above and sketched in Fig. 1, the FFDP is based on the measurement of the pressure, P i , in the displacing fluid at the inlet and the corresponding displacing fluid flow rate Q(P i ). The pressure, P o , at the outlet is fixed and considered to be the same in the displacing and displaced fluids (this point is further discussed in Section 4). The commonly used method [START_REF] Morison | A comparison of liquid-liquid porosimetry equations for evaluation of pore size distribution[END_REF] for determining the TSD from the flow-rate data Q(P i ) in the FFDP explicitly relies on the model of a bundle of parallel cylindrical tubes. On this basis, the flow-rate is computed by applying Poiseuille's law in each tube containing the non-wetting fluid or possibly a mixed Poiseuille-Knudsen flow model in the case of GLDP [START_REF] Islam | Mixed Poiseuille-Knudsen flow model for gas liquid displacement porometry data treatment[END_REF]. In what follows, Knudsen effects are not considered. Moreover, steady, laminar and incompressible flow is assumed and both fluids are supposed Newtonian. This allows writing

Q i = ΔPπ 8μL ∫ +∞ rmin r 4 f Q (r)dr (1)
Here, L is the tube length, ΔP = P i -P o , μ is the displacing fluid dynamic viscosity; the radius of a tube is denoted by r and f Q is the TSD to be determined. For a given inlet pressure P i , the radius r min is defined with the Young-Laplace equation as

r min = 2γ cosθ ΔP ( 2
)
where γ is the interfacial tension and θ the contact angle. In the present O. Maalal et al. International Journal of Multiphase Flow 149 (2022) 103983

work, the displaced fluid is considered as being perfectly wetting and hence, θ = 0. Eq. ( 2) expresses that all the tubes having a radius greater or equal to r min are invaded by the displacing fluid. The greater the inlet pressure, the lower is the threshold radius r min . Combining Eqs. ( 1) and

(2) leads to

Q i = πγ cosθ 4μLr min ∫ +∞ rmin ρ 4 f Q (ρ)dρ (3)
Considering Eq. (3) for r min = r and deriving the relationship with respect to r leads to

dQ i dr = - Q i r - πγ cosθ 4 μL r 3 f Q (r) (4)
and thus to

f Q (r) = - 4 μL πγ cosθ r -3 ( dQ i dr + Q i r ) (5)
Eq. ( 5) can be expressed in discrete form as [START_REF] Morison | A comparison of liquid-liquid porosimetry equations for evaluation of pore size distribution[END_REF] f Qk (r AB ) = -

4 μL πγ cosθ r -3 AB ( Q iB -Q i A r B -r A + Q iAB r AB ) (6)
where Q iA and Q iB are two close values of Q i corresponding to two close values of the inlet pressure P iA and P iB ,

Q iAB = Qi B +Qi A 2
and r AB = (r A + r B ) /2, r A and r B being determined from Eq. ( 2).

The value of f Qk can then normalized as

d Qk (r AB ) = f Qk ∑ n i=1 f Qi (7)
where n is the number of classes considered to compute the throat size distribution. In what follows, water is considered as the wetting fluid and air as the displacing fluid. Thus γ = 0.072 N/m and μ = μ air = 10 -5

Pa.s.

The above equations are similar to Eqs. ( 9)-( 14) in [START_REF] Maalal | Determination of the throat size distribution of a porous medium as an inverse optimization problem combining pore network modeling and genetic and hill climbing algorithms[END_REF]. However, there is a significant difference. In this reference, only quasi-static displacements were considered. As a result, the viscous pressure drop ΔP in Eq. (1) was assumed small compared to the capillary pressure, i.e., the pressure difference between the two fluids. As a consequence, ΔP was independent of r min . The relaxation of the quasi-static hypothesis induces a modification of the radius exponent in Eqs. ( 5) and ( 6), which is now equal to -3 whereas it is equal to -4 in the corresponding equations in [START_REF] Maalal | Determination of the throat size distribution of a porous medium as an inverse optimization problem combining pore network modeling and genetic and hill climbing algorithms[END_REF]. This is because the viscous pressure drop and the capillary pressure are directly correlated as expressed by Eq. ( 2) in the viscous case whereas the pressure drop is independent of the capillary pressure in the QS case considered in [START_REF] Maalal | Determination of the throat size distribution of a porous medium as an inverse optimization problem combining pore network modeling and genetic and hill climbing algorithms[END_REF].

Interconnected model porous medium

For FFDP, immiscible displacement is simulated in the model porous medium depicted in Fig. 1 using the pore network modeling (PNM) approach [START_REF] Blunt | Detailed physics, predictive capabilities and macroscopic consequences for pore-network models of multiphase flow[END_REF]. Obviously, this structure does not assume an underlying cylindrical parallel tubes structure. As illustrated in this figure, pore bodies in the pore network model are spheres of radius r p located at the nodes of a regular square grid. The distance between two adjacent nodes is the lattice spacing a, with a = 350 μm in all the simulations presented below. The size of the network is the number of nodes (pores) in each direction (for example, Fig. 1 shows a 30 × 30 square network). Throats are cylindrical channels of radius r t . The throat size is randomly distributed according to a given probability density function (p.d.f.). Similarly, the pore body sizes are distributed according to a given p.d.f. with the constraint that the size of a pore body is equal to, or greater, than the size of the largest throat to which it is connected. In the remainder of this work, the same type of p.d.f. is considered for both the pore bodies and the throats of a given network. In what follows, four types of p.d.f. are considered, namely uniform, Gaussian, Log-normal and bimodal. Also, it can be noted that the network coordination number, i.e., the number of throats connected to a pore body, is 4 in the 2D square network.

Fluidfluid immiscible displacement simulation

In the FFDP, a wetting fluid is displaced by a non-wetting one, a process called drainage. It is analyzed in details by [START_REF] Lenormand | Numerical models and experiments on immiscible displacements in porous media[END_REF], where it is shown that the resulting flow pattern depends on the competition between capillary and viscous forces. The method for computing the flow is similar to that described in this reference. Since the inlet pressure variation between two pressure steps is small, a simplification is to consider that the pressure variation in the displaced fluid between the two steady states corresponding to the two considered pressure steps can be neglected. The pressure is thus considered as uniform and equal to P o in the displaced fluid. In fact, there is necessarily a flow in the displaced fluid and thus a viscous pressure drop during the transient period between two equilibria (equilibrium refers to the steady-state situation reached after each inlet pressure increment). However, only a small fraction of menisci moves during each transient period since the inlet pressure increments are small. As a result, the transient flow rate in the displaced fluid is small. The approximation that is performed consists in neglecting this flow and the associated displaced fluid viscous pressure drop in the identification of the moving menisci (see below). This means that the viscous flow is only computed in the sub-region of the network occupied by the displacing fluid. The displacing fluid flow computation consists in expressing the mass conservation equation at each network node (pore body) j that is given by

∑ Q i-j,k = 0 (8)
where Q i-j,k is the non-wetting fluid volume flow rate between pore body j and adjacent pore body k. If the throat between pore bodies j and k is occupied by the wetting fluid, then Q i-j,k = 0. If the throat is occupied by the non-wetting fluid, then

Q i-j,k = πr 4 t j,k 8μ nw l j,k ΔP j,k (9)
where r t j,k and l j,k are the radius and length of the throat connecting the two pore bodies, respectively, whereas ΔP j,k is the pressure drop in the displacing fluid between the two neighboring pore bodies. The boundary conditions at the inlet and outlet of the network are P = P i and P = P o = const, respectively. A zero-flux condition is imposed on the lateral faces of the network. This leads to a linear system for the pressure field at the nodes of the network in the non-wetting fluid. This system is numerically solved using the conjugate gradient method. Once the pressure field is obtained, the displacing fluid flow rate Q i is computed using Poiseuille's law applied to all throats saturated by the displacing fluid and connected to the inlet. It can thus be noticed that the pressure is considered to be the same in the wetting fluid and non-wetting fluid at the outlet. This is an approximation. Actually, the pressure in the wetting fluid is expected to be slightly lower than in the non-wetting fluid at the outlet where the non-wetting fluid saturation remains small. The approximation consists in neglecting the pressure difference between the two fluids at the outlet compared to the viscous pressure drop, i.e. the pressure difference between the two fluids at the inlet. The flow computation method described above applies for a given fluid distribution in the network. It must be coupled with an algorithm allowing the determination of the fluid distribution. To this end, a meniscus at the entrance of a throat is defined as stable when the pressure difference between each side of the meniscus is less than the capillary pressure threshold of the throat under consideration (such a throat is referred to as an interfacial throat), which means O. Maalal et al. International Journal of Multiphase Flow 149 (2022) 103983

P i-j,k -P o < P cth (10)
Here, P cth = 2γ cosθ rt and P i-j,k is the pressure in the displacing fluidsaturated pore body adjacent to the throat under consideration, of radius r t . Starting from the fluid distribution obtained at the previous step, the pressure field in the displacing fluid for the new pressure step is computed as indicated above. Then, all unstable menisci are identified. This corresponds to all menisci such that P i-j,k -P o ≥ P cth . The interfacial throat for which (P i-j,k -P o ) -P cth is maximum is invaded with the displacing fluid together with the adjacent pore body if the latter was not already invaded in a previous step. If this pore body was occupied by the displaced fluid, new menisci are then positioned at the entrance of the throats occupied by the displaced fluid and connected to the newly invaded pore body. Thus, the list of interfacial throats is updated so as to remove the invaded throat and add the throats connected to the newly invaded pore body and occupied by the displaced fluid. This procedure is repeated until there is no unstable menisci remaining in the network, i.e., until there exists no interfacial throat for which P i-j,k -P o ≥ P cth .

The resulting fluid distribution corresponds to the equilibrium distribution at the specified value of P i . The above procedure is enforced after the breakthrough (BT). The breakthrough corresponds to the minimum inlet pressure P iBT for which the displacing fluid reaches the outlet. For inlet pressures lower than P iBT , the displacing fluid does not reach the outlet and therefore there is no flow when steady state is reached. The procedure to compute the fluid distribution before BT is similar to the one described above. The important difference is that the pressure is spatially uniform and equal to P i in the displacing fluid since there is no flow before BT when the capillary equilibrium is reached. The invasion rules described above correspond to a drainage process without trapping. An element, pore body or throat, occupied by the displacing fluid is thus assumed to be always connected to the outlet face, possibly through liquid films [START_REF] Zhou | Hydrocarbon drainage along corners of noncircular capillaries[END_REF] present in the pore and throat surface roughness and/or in a subnetwork of small crevices at the pore walls. Implementation of trapping rules is possible (see for instance [START_REF] Joekar-Niasar | Insights into the relationships among capillary pressure, saturation, interfacial area and relative permeability using pore-network modelling[END_REF], but, for the sake of simplicity, trapping is not considered in what follows.

In this regard, it can be recalled that trapping is exacerbated in 2D networks compared to 3D networks [START_REF] Wilkinson | Invasion percolation: a new form of percolation theory[END_REF]. Thus, neglecting trapping can be also seen as an approximate way of considering a situation closer to the 3D case. The crucial point we keep in our approach is that menisci must travel in the pore space and that each meniscus displacement acts as a pore space local geometry probe. This point is further discussed in Section 8.

Examples of fluid distribution in a 30 × 30 2D network obtained with the procedure described above are displayed in Fig. 2 As can be seen from Fig. 2, the displacement occurs in comparatively a larger number of throats and pore bodies in the lower region of the network, that is the region adjacent to the inlet, than in the upper region of the network. This is fully consistent with the existence of the nonlinear viscous pressure drop in the displacing fluid illustrated in Fig. 3. Consider two interfacial throats of the same radius. Such a throat is invaded by the displacing fluid when P i-j,k -P o ≥ P cth . The closer the throat to the inlet, the greater is the local pressure P i-j,k (see Fig. 3). As a result, in the sequence of pressure steps, a throat closer to the inlet is invaded before a throat of comparable radius located further away in the network. This also means that throats of different radii are invaded when the inlet pressure is incremented. For a given inlet pressure, since P i-j,k decreases on average with the increasing distance from the inlet, the critical invasion threshold radius given by the equation P i-j,k -P o = P cth = 2γcosθ rt increases with the increasing distance from the inlet. In other words, results illustrated in Figs. 2 and3 clearly indicate that there is not a unique value of the threshold pore radius given by the Young-Laplace Eq. ( 2) for a given value of the inlet pressure. This is in contrast with the common FFDP assumption, for which a single pore radius is associated to each inlet pressure value. This also differs from the IP approach followed by [START_REF] Maalal | Determination of the throat size distribution of a porous medium as an inverse optimization problem combining pore network modeling and genetic and hill climbing algorithms[END_REF]. In fact, the situation considered in the present article can be analyzed as a situation of invasion percolation in a stabilizing gradient (IPSG) [START_REF] Xu | Invasion percolation with viscous forces[END_REF]. The gradient in the pore occupation probability leads to the formation of a stable invasion front corresponding, in Fig. 2, to the gradual shrinkage of the two-phase zone, which localizes here on the outlet side of the network.

The FFDP curve obtained from the algorithm above described is shown in Fig. 4. Note that this curve corresponds to one realization of the network and does not result from an average over several realizations (contrary to the pressure curves in Fig. 3). This curve typically shows the displacing fluid flow rate as a function of the pressure difference ΔP = P i -P o . Throughout the article, the displacing fluid flow rate, Q i , is normalized. More precisely, the quantity Qi Pi-Po is considered and normalized by Qif P if -Po , Q if and P if being the final displacing fluid flow rate and inlet pressure, respectively. Thus, the quantity represented in Fig. 4 is

χ = Q i ( P if -P o ) Q if (P i -P o ) (12)
It must be noted that χ should not be mistakenly identified as the displacing fluid relative permeability (k r ) employed in the classical theory of two-phase flow in porous media [START_REF] Dullien | Porous Media: Fluid Transport and Pore Structure[END_REF]. Indeed, χ would correspond to k r if the fluid distributions were obtained via quasi-static displacements controlled only by the capillary effects.

Owing to the viscous pressure drop in the displacing fluid, this is obviously not the case in the situation under study. To this extend, χ can be seen as a "dynamic" relative permeability.

TSD extraction from commonly used FFDP procedure

The next step is to apply the TSD extraction procedure commonly used for FFDP, i.e., Eqs. (1) to (7), starting from the computed FFDP curve (Fig. 4). The objective is to analyze how the TSD obtained with this procedure (Eq. ( 6)), referred to as the FFDP TSD, compares with the real TSD. A 2D square 50 × 5, a 2D square 50 × 25 and a 2D square 30 × 30 networks were considered for the comparison. Note that the second value in the network size nomenclature denotes its height in the main flow direction. For each network, the χ(ΔP) curve was computed and the FFDP TSD extracted using Eq. ( 6). For a given network size and distribution, the results were averaged over 100 network realizations so as to eliminate dispersion of the results. For a given network size, the FFDP TSD was extracted from each χ(ΔP) curve yielding 100 FFDP TSD which were finally averaged to obtain the FFDP TSD for the network under concern. The comparison was performed for two types of p.d.f., namely, uniform and Gaussian. The p.d.f. parameters are the throats and pore bodies minimum and maximum radii, r pmin , r pmax and r tmin , r tmax for a uniform distribution, the mean pore and throat radii, r pmoy , r tmoy and standard deviation, σ, in the case of a Gaussian p.d.f.. The p.d.f. parameters considered in this section are reported in Table 1.

The comparison between the real TSD and the FFDP TSD for the various cases is presented in Fig. 5 and Table2. As can be seen, the FFDP TSD is significantly different from the real TSD for all the cases Fig. 3. Displacing fluid slice average pressure profiles in the 30 × 30 network corresponding to the fluid distributions depicted in Fig. 2. Note that the pressure is averaged over 30 realizations whereas only one realization is of course shown in Fig. 2. The inlet and outlet respectively correspond to a normalized height equal to 0 and 1.

Table 1

PBSD and TSD parameters considered in Section 5. O. Maalal et al. International Journal of Multiphase Flow 149 (2022) 103983 considered here. The FFDP procedure introduces a significant positive skewness. In addition to the noticeable asymmetry, it can be seen that the FFDP procedure leads to identify throat sizes smaller than in the real TSD. This is qualitatively consistent with the experimental results reported in [START_REF] Peinador | Comparison of capillary flow porometry (CFP) and liquid extrusion porometry (LEP) techniques for the characterization of porous and face mask membranes[END_REF], showing a shift toward smaller pore sizes between the pore size distribution determined from the FFDP and the one deduced from the capillary pressure curve. The shift can be explained as follows. As illustrated in Fig. 2, the region near the inlet is preferentially invaded due to the greater displacing fluid pressure. When the inlet pressure is sufficiently large, menisci are only present at some distance from the inlet where the pressure in the displacing fluid is less than at the inlet. The minimum size identified from Eq. ( 1) is r min = 2γ cosθ Pi-P0 whereas the actual minimum size must be determined not at the inlet where there is no menisci anymore but where menisci are present, thus r min-real ≈ 2γ cosθ P(z)-P0 where z is the distance from the inlet where the less advanced menisci, i.e. the menisci the closest to inlet, can be found. Since P(z) -P 0 < P i -P 0 , r min . < r min-real . Hence, the throat sizes extracted from Eq. ( 1) are eventually much smaller than the sizes of the throats in which the menisci are actually moving since they move in a region where the displaced fluid pressure is significantly lower than the inlet pressure. The higher the inlet pressure, the greater is the impact since, as discussed before, the two-phase zone, i.e., the zone containing the menisci, shrinks and localizes on the outlet side of the network while increasing the inlet pressure (Fig. 2). The shift is therefore due to the viscous pressure drop and the use of Eq. ( 1) wrongly considering the inlet pressure and not the lower pressures in the region where the menisci actually move.

As reported in Table 2, the errors associated with the parallel cylindrical tube simplification is very significant. In summary, the nature of the distribution, uniform or Gaussian, is not retrieved and the identified sizes are too small. The comparison of the FFDP TSD with the real TSD for the two networks, in Fig. 5, shows that the shift in the distribution toward unrealistic small throat sizes depends on the network height, i.e., the distance between the inlet and the outlet measured in lattice spacing unit corresponding to the number of elements (pores or throats) over the network height. The greater the height, the more pronounced is the shift. This is consistent with the above consideration on the wrong use of Eq. ( 1) since the greater the height the greater the inlet pressure must be to displace the menisci near the outlet. This is an indication that the commonly used FFDP TSD extraction method should not be used, in particular for the characterisation of thick porous materials.

The results presented in Fig. 5 and Table 2 are consistent with the fluid distributions and the pressure profiles depicted in Figs. 2 and3. As discussed previously, the existence of the pressure drop across the network is not consistent with the identification of a single throat size from the Young-Laplace equation, as assumed in the commonly used FFDP procedure. In order to fully invade the network, the pressure difference between the two fluids in the outlet region of the network must be comparable to (slightly greater than) the throat capillary pressure thresholds. Since the displacing fluid pressure is significantly higher in the network inlet region, using the inlet pressure P i together with the Young-Laplace equation leads necessarily to underestimate the throat sizes. Since the viscous pressure drop increases with the network height, greater inlet pressures are necessary to fully invade the network when its height is increased. This is consistent with the increase in the shift toward smaller throat sizes with an increasing network size reported in Fig. 5.

Alternative method for the FFDP TSD extraction

Method of solution

Owing to the serious shortcomings inherent to the commonly used FFDP TSD extraction procedure, the question arises as to whether a better procedure can be devised. As shown in the previous section, the standard and commonly used FFDP TSD extraction procedure does not allow determining the real TSD from the Q i (ΔP) curve. The main reason explaining this unsatisfactory estimation lies in the fact that throats of various sizes are invaded when the pressure is incremented due to the displacing fluid pressure variation across the porous medium whereas only a single throat size is assigned via the Young-Laplace Eq. ( 2) to all the newly invaded throats with the commonly used FFDP procedure, regardless of their actual sizes. This notably introduces a shift toward (unrealistic) small throat sizes.

In order to improve the estimation, the TSD identification problem is addressed as an inverse problem. The direct problem is the determination of the Q i (ΔP) curve for a given microstructure. The inverse problem consists in identifying the TSD from the Q i (ΔP) curve and is treated with an optimization procedure in what follows. The general idea is to gradually modify the microstructure (more precisely the sizes assigned to the throats for a given network structure) so as to retrieve the Q i (ΔP) curve accurately. To this end, the method presented in [START_REF] Maalal | Determination of the throat size distribution of a porous medium as an inverse optimization problem combining pore network modeling and genetic and hill climbing algorithms[END_REF] is used. In this reference, the objective was to identify the TSD from the displacing fluid relative permeability and the retention curve whereas here, the objective is to identify the TSD from the FFDP curve. As in the previously mentioned work, it is assumed that the pore network structure is known. This means that the optimization procedure is performed for a given type of network, i.e., a square network in the present study. Again, the unknown is the TSD. The pore body sizes are not considered in the optimization procedure on the ground that the FFDP flow is essentially controlled by the throats. Thus, the pore size is taken as constant (equal to 80 µm) in what follows. A particular realization of the network, with specified p.d.f. parameters, is used as the reference network whose TSD is hence known. This reference network is the equivalent of the porous medium to be characterized in practice, that is, the one for which the experimental measurement of Q i (ΔP) would have been performed in practice and whose TSD is to be determined.

The optimization procedure is actually fully identical to the one described in [START_REF] Maalal | Determination of the throat size distribution of a porous medium as an inverse optimization problem combining pore network modeling and genetic and hill climbing algorithms[END_REF], considering χ as defined by Eq. ( 12) in place of Q n in the latter reference. For this reason, the optimization algorithm is not described again here. The reader is referred to the above-mentioned article for the details. In essence, the procedure consists in gradually modifying the throat sizes in the network so as the FFDP curve of the optimized network matches the FFDP curve of the reference network using a hybrid genetic and hill climbing algorithm.

Table 2

Comparison between the real TSD and the TSD obtained from the commonly used FFDP procedure (Eq. ( 6)), i.e., the FFDP TSD. Mean values and standard deviations were computed over the 100 realizations.

Results

A reference network, with a Gaussian TSD having a mean radius of 40 μm and a standard deviation of 10 μm, is first considered. The corresponding reference TSD is shown in the inset of Fig. 6a. Note that the initial network population in the optimization process is generated with uniform distributions, thus with a type of p.d.f. different from the Gaussian distribution of the reference network.

As shown in Fig. 6b, the optimization procedure leads to an almost perfect agreement between the χ(ΔP) curves for the reference and optimized networks. Compared to the commonly used FFDP procedure (Fig. 5), the identified TSD is much closer to the reference TSD. There is no shift toward unrealistic small throat sizes and the distribution is much closer to a Gaussian distribution. However, one can note a slight positive skewness. The values for the first four moments of the distribution reported in Table 3 are also in much better agreement with the corresponding values of the real TSD (as can be seen from the comparison with the values reported in Table 2 obtained for the commonly used FFDP procedure).

As reported in Fig. 7 and Table 3, somewhat similar results are obtained when the reference p.d.f. is uniform. In both cases, it can be seen that a few throat sizes greater than the largest throat size in the reference distribution are obtained through the optimization procedure. This is due to the specification of the maximum throat size in the optimization procedure. As explained in [START_REF] Maalal | Determination of the throat size distribution of a porous medium as an inverse optimization problem combining pore network modeling and genetic and hill climbing algorithms[END_REF], the maximum throat radius is approximately specified by multiplying by a factor 1.3 the radius corresponding to the first pressure step at which the displaced fluid saturation is no longer equal to one. The factor 1.3 is introduced because the biggest throat at the inlet can be expected to be smaller than the biggest throat in the whole network. Nevertheless, the fraction of throat sizes in this too large range is quite low. This is an indication that the optimization procedure should lead to still better results with a more refined estimate of the maximum throat radius (via smaller pressure steps for example).

TSD extraction combining FFDP and retention curves

Method of solution

It was shown in [START_REF] Maalal | Determination of the throat size distribution of a porous medium as an inverse optimization problem combining pore network modeling and genetic and hill climbing algorithms[END_REF] that the inversion optimization procedure, considering the displacing fluid relative permeability and the retention curve, led to better results than when only the displacing fluid relative permeability was considered. Thus, one can wonder whether considering the retention curve in addition to the FFDP curve could improve the TSD identification in the case under study here. The retention curve corresponds to a series of capillary quasi-static equilibria within the porous medium and is a classical characterization of the material, also for thin porous media of special interest in the present study (see e.g. [START_REF] Fairweather | The effects of wetproofing on the capillary properties of proton exchange membrane fuel cell gas diffusion layers[END_REF]. However, this implies in practice to perform an additional experiment since the FFDP devices are typically not adapted to measure the retention curve. Examples of specific measurement devices adapted to thin porous media are described for instance in [START_REF] Fairweather | The effects of wetproofing on the capillary properties of proton exchange membrane fuel cell gas diffusion layers[END_REF] and [START_REF] Gostick | Direct measurement of the capillary pressure characteristics of water-air-gas diffusion layer systems for PEM fuel cell[END_REF]. To obtain the drainage retention curve, the procedure is somewhat similar to the FFDP procedure. The porous medium is initially saturated by a wetting fluid. Then a non-wetting fluid is injected through a series of pressure steps. An important difference compared to the FFDP devices is the presence of a porous membrane with smaller pores than in the medium of interest at the outlet, preventing the non-wetting fluid to exit the porous medium. Thus, contrary to the FFDP, there is no displacing fluid flow in the porous medium when an equilibrium is reached. An important difference lies also in the fluid-fluid distribution. In the case of the retention curve, the fluid-fluid distributions are quasi-static, only controlled by capillary effects. By contrast, as illustrated in Fig. 2

Table 3

Comparison between the reference TSD and the optimzed TSD. Case where the reference TSD is the Gaussian TSD depicted in the inset in Fig. 6a. were measured and are available. For the present study, this means that the curve P c (S) can be computed for the reference network as well as for all the networks generated during the optimization procedure. The algorithm for computing the retention curve in a network is described in [START_REF] Maalal | Determination of the throat size distribution of a porous medium as an inverse optimization problem combining pore network modeling and genetic and hill climbing algorithms[END_REF] and therefore is not further again reported here. The optimization procedure is completely analogous to the one described in this reference when both the retention curve and the relative permeability curve is used. The sole modification is the use of χ (Eq. ( 12)) in place of Q n in [START_REF] Maalal | Determination of the throat size distribution of a porous medium as an inverse optimization problem combining pore network modeling and genetic and hill climbing algorithms[END_REF].

Results

In Fig. 8, the results on the retention curve, the FFDP curve and the TSD resulting from the optimization procedure are reported together with the corresponding quantities for the reference network.

As can be seen from the comparison between Figs. 6 and 8, using both the retention and FFDP curves in the optimization process leads to significantly better results. This is also confirmed by the results reported in Table 4. It can be noted in Fig. 8a and b that the retention and FFDP curves resulting from the optimization procedure reproduce very well the corresponding curves of the reference case. As shown earlier [START_REF] Maalal | Determination of the throat size distribution of a porous medium as an inverse optimization problem combining pore network modeling and genetic and hill climbing algorithms[END_REF], it is of course possible to develop the inverse optimization procedure using only the retention curve. For completeness, the results obtained with this strategy are also indicated in Table 4. As can be seen, results when only the retention curve is used are slightly better than when only the FFDP curve is employed. However, the best results are obtained when both the retention and the FFDP curves are used.

In Fig. 9 and Table 5, results obtained for the case of the uniform distribution depicted in the inset of Fig. 9c are reported. A comparison of Fig. 9 with Fig. 7 shows that here, again, significantly better results are obtained when the FFDP curve is used in conjunction with the retention curve in the optimization process. However, results for the uniform distribution can be judged as slightly less satisfactory than for the Gaussian distribution. In particular, it can be noticed that the procedure leads to identify a few throat sizes greater than the greatest throat size in the reference distribution. Again, this can be attributed to the procedure used to estimate the maximum radius of the distribution with the factor 1.3, as mentioned at the end of Section 6.2.

Nevertheless, to further illustrate the capabilities of the inverse optimization procedure as above described, two additional tests were performed considering a reference Log-normal distribution (as depicted in Fig. 10a) and the expected more challenging bimodal distribution depicted in Fig. 11a. For these cases, only the full optimization procedure, i.e., the one combining both the P c (S) and χ(ΔP) curves was used. As can be seen from Figs. 10 and 11, respectively corresponding to the Log-normal and bimodal distributions, and from Table 6, results are quite satisfactory. Here also, the match between the reference capillary pressure and FFDP curves and the curves resulting from the optimization procedure is excellent although they are not reported here for the sake of brevity. Finally, from the four reference TSD tested in this study, the less satisfactory results were obtained for the uniform TSD, mainly because of the uncertainty associated with the specification of the upper bound radius of the distribution.

Discussion

The TSD identification based on the optimization procedure clearly leads to much better results than the conventional FFDP approach. However, several simplifications were introduced. The displacement was computed neglecting the viscous pressure drop in the displaced fluid during the transient period occurring right after the inlet pressure increase until a new steady-state (equilibrium) is reached. As explained before, this approximation was made because the FFDP procedure is based on small inlet pressure increments, corresponding to a few meniscus displacements between two successive pressure steps. Nevertheless, it would be interesting to test the quality of the approximation from comparisons with more advanced models, such as dynamic pore networks for instance, (see e.g., [START_REF] Joekar-Niasar | Non-equilibrium effects in capillarity and interfacial area in two-phase flow: dynamic pore-network modelling[END_REF][START_REF] Joekar-Niasar | Analysis of fundamentals of two-phase flow in porous media using dynamic pore-network models: a review[END_REF]. Trapping was neglected because it is exacerbated in 2D networks compared to 3D networks. In this respect, it would be interesting to consider 3D networks with trapping. It may be expected that some pore sizes are unidentified because of trapping since the trapped regions cannot be visited by the menisci probing the pore space.

Also a network of identical structure, i.e., a square network, was considered for the reference network and the networks generated in the optimization procedure. Obviously, this is a serious limitation. However, the idea was to first test the procedure for this particular case where the network structure is known. The next step would be to test the procedure considering a reference network structure different from the one used in the optimization procedure and explore whether it is still possible to reasonably identify the TSD. In this respect, the connectivity,

Fig. 1 .

 1 Fig. 1. Square pore network and sketch of immiscible displacement in FFDP (displacing fluid in blue, displaced fluid in light grey). The red dashed area corresponds to an exemplary network slice used for spatially averaging the pressure.

  Fig. 2. Example of fluid-fluid distributions, computed in a 30 × 30 2D network with a uniform TSD (r tmin = 30 μm, r tmax = 50 μm) and a uniform PBSD (r pmin = 80 μm, r pmax = 120 μm) for an increasing pressure difference ΔP = P i -P o . The displacing fluid in blue enters the network at the bottom. S is the displaced fluid saturation, i. e., the network volume fraction occupied by the displaced (wetting) fluid. (a) At BT. (b) to (f) After BT. (a) ΔP=3653 Pa, S=0.6 (b) ΔP =3681 Pa, S=0.5 (c) ΔP =3722 Pa, S=0.4. (d) ΔP =3743 Pa, S=0.3 (e) ΔP =3768 Pa, S=0.2 (f) ΔP =3781 Pa, S=0.1

Fig. 4 .

 4 Fig. 4. FFDP curve: displacing normalized fluid flow rate, χ, as a function of ΔP = P i -P o .

Fig. 5 .

 5 Fig. 5. Comparison of the TSD obtained from the commonly used FFDP procedure, Eq. (6), and the real TSD (shown in the insets) for various cases. (a) 50 × 5 network, Gaussian p.d.f., r tmoy = 40 μm, σ = 10 μm. (b) 50 × 5 network, uniform p.d.f., r tmin = 30 μm, r tmax = 50 μm. (c) 50 × 25 network, Gaussian p.d.f., r tmoy = 40 μm, σ = 10 μm. (d) 50 × 25 network, uniform p.d.f., r tmin = 30 μm, r tmax = 50 μm. e 30 × 30 network, Gaussian p.d.f., r tmoy = 40 μm, σ = 10 μm. (f) 30 × 30 network, uniform p.d.f., r tmin = 30 μm, r tmax = 50 μm.

  Fig. 6. Results of the optimization procedure for a case where the reference p.d.f. is Gaussian. (a) TSD of the optimized network (the inset shows the reference TSD). (b) FFDP curve for the reference network and the optimized network.

Fig. 7 .

 7 Fig. 7. Results of optimization procedure for a case where the reference p.d.f. is uniform. (a) FFDP curve for the reference network and the optimized network. (b) TSD of the optimized network (the inset shows the reference TSD).

Acknowledgments

This research was supported by the Project "PEMFC -SUDOE" -SOE1/P1/E0293 which is co-financed by the European Regional Development Fund in the framework of the Interreg Sudoe programme.

International Journal of Multiphase Flow 149 (2022) 103983 i.e., the number of adjacent pores connected to a given pore, appears as an important parameter [START_REF] Meyers | Determination of the pore connectivity and pore size distribution and pore spatial distribution of porous chromatographic particles from nitrogen sorption measurements and pore network modelling theory[END_REF][START_REF] Mourhatch | Determination of the true pore size distribution by flow permporometry experiments: an invasion percolation model[END_REF]. This suggests to include the connectivity as a parameter to be identified in the optimization procedure. Some of the above elements were also discussed in [START_REF] Maalal | Determination of the throat size distribution of a porous medium as an inverse optimization problem combining pore network modeling and genetic and hill climbing algorithms[END_REF] where complementary aspects can be found. In this reference, it was concluded that the combined use of the retention and displacing fluid relative permeability curves leads to better results than the use of the relative permeability curve alone. A similar conclusion is reached here since the combined use of the retention and FFDP curves leads to the best results for all the distributions tested.

The combined use of the retention curve and FFDP curves is thus recommended. For the TSD identification of thin porous media with throats in the micronic range, such as the diffusion medium used in fuel cells for example [START_REF] Fairweather | The effects of wetproofing on the capillary properties of proton exchange membrane fuel cell gas diffusion layers[END_REF][START_REF] Gostick | Direct measurement of the capillary pressure characteristics of water-air-gas diffusion layer systems for PEM fuel cell[END_REF][START_REF] Park | A review of gas diffusion layer in PEM fuel cells: materials and designs[END_REF], obtaining both curves should not be a problem. For thin nanoporous media, the combined approach could be more difficult to use because measurement of the retention curve can be a challenging

Table 4

Comparison between the reference TSD and the TSD obtained from the optimization procedure. Case where the reference TSD is the Gaussian TSD depicted in the inset of Fig. 8c. task. In practice, the retention curve measurement implies to place a porous membrane at the outlet of the porous medium of interest. The pores in the outlet membrane must be smaller than the pores in the porous medium to be characterized in order to prevent the non-wetting fluid to exit the sample. Using a membrane with the required pore size when the pores in the porous medium of interest are in the nanometer range may be a difficulty in practice. This problem is not encountered with the FFDP technique. Also, the FFDP technique can be performed using a pair of immiscible liquids with a relatively low interfacial tension so as to generate the displacement with reasonably low-pressure differences even in the presence of nano-pores. Thus, in that case, the FFDP technique combined with the optimization procedure is an interesting option to characterize the TSD.

Conclusion

Pore network simulations were developed in order to evaluate the fluid-fluid displacement porosimetry (FFDP) technique. It was suggested that two types of microstructure elements, namely pore bodies and pore throats, must be distinguished in the pore space. It was argued that the

Table 5

Comparison between the reference TSD and the TSD obtained from the optimization procedure. Case where the reference TSD is the uniform TSD depicted in the inset in Fig. 9c. distribution commonly obtained with the FFDP technique was the throat size distribution (TSD) rather than the pore body size distribution (PBSD) in accordance with the physical process at play (drainage). Direct Pore network simulations indicate that the commonly used FFDP procedure for identifying the TSD from the FFDP curve introduces a bias toward unrealistic small sizes as a consequence of the viscous pressure drop in the displacing fluid. The thicker the porous medium, the greater is the bias. In practice, this means that the results obtained with the commonly used FFDP procedure should not be regarded as representative of the real TSD when the pore space structure is significantly different from a system of straight non-interconnected quasi-cylindrical pores.

A significantly more sophisticated procedure, combining pore network simulations of the immiscible displacement and an optimization technique based on a hybrid genetic and hill climbing algorithm recently reported in the literature was then tested to improve the TSD identification from the FFDP raw data. This procedure leads to much better results than the usual FFDP interpretation and allows eliminating the bias toward unrealistic small throats. Furthermore, it was shown that the TSD identification could be still improved by considering the retention curve in addition to the FFDP curve in the optimization procedure.

The efficiency of the method proposed here is to the cost of a significantly less simple procedure than the commonly used FFDP interpretation technique. The latter is analytical whereas the

Table 6

Comparison between the reference TSD and the TSD obtained from the optimization procedure making use of both the capillary pressure and FFDP curves as target data. Cases where the reference TSD is the Log-normal TSD depicted in Fig. 10a and the bimodal TSD depicted in Fig. 11a. optimization procedure implies repeated numerical simulations of the immiscible displacements in the pore space. In addition, several simplifications were introduced to test the improved optimization procedure. In particular, identical network structure and size were considered in the reference network and the various networks generated during the optimization procedure. The variability in the pore size was neglected, considering only the throat size distribution. Owing to these simplifications, solving the problem consisting in identifying the TSD (and the PBSD) from macroscopic data such as the FFDP and capillary pressure curves clearly requires further work.

CRediT authorship contribution statement

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.