Dat Tri

Ngo

André Fourno

Benoit Noetinger

Modeling of transport processes through large-scale discrete fracture networks using conforming meshes and open-source software

Keywords: DFN, Conforming mesh, Transport processes, Vertex-centered finite volume method

Most industrial and field studies of transport processes in Discrete Fracture Networks (DFNs) involve strong simplifying assumptions, especially at the meshing stage. High-accuracy simulations are therefore required for validating these simplified models and their domain of validity. The present paper proposes an efficient workflow based on open-source software to obtain transport simulations. High-quality computational meshes for DFNs are first generated using the conforming meshing approach FraC. Then, a tracer transport model implemented in the open-source code DuMux is used for simulating tracer transport driven by the advection-dispersion equation. We adopt the box method, a vertex-centered finite volume scheme for spatial discretization, which ensures concentration continuity and mass conservation at intersections between fractures. Numerical results on simple networks for validation purposes and on complex realistic DFNs are presented. An a-posteriori convergence study of the discretization method shows an order of convergence O(h) for tracer concentration with h the mesh size.

Introduction

During recent decades, modeling of transport processes in fractured porous media has received increasing attention from the geoscientist community because of its wide range of applications in geothermal energy [START_REF] Aquilina | A tracer test at the Soultz-Sous-Forets Hot Dry Rock geothermal site[END_REF], petroleum exploration and production [START_REF] Bourbiaux | Fractured reservoir simulation: A challenging and rewarding issue[END_REF], CO 2 geological storage and nuclear waste disposal [START_REF] Grenier | Assessment of Retention Processes for Transport in a Fractured System at Äspö (Sweden) Granitic Site: from Short-Time Experiments to Long-Time Predictive Models[END_REF][START_REF] Fourno | Development and qualification of a smeared fracture modelling approach for transfers in fractured media[END_REF]. Fracture networks might consist of pre-existing natural or pressure-induced fractures in the case of hydraulic fracturing [START_REF] Delorme | Unconventional production forecast needs integration of field hydraulic stimulation data through fracture model calibration and optimized numerical scheme[END_REF].

Solute transport in fractured media should be studied using either continuous [START_REF] Svensson | A continuum representation of fracture networks. Part I: Method and basic test cases[END_REF][START_REF] Karimi-Fard | Generation of coarse-scale continuum flow models from detailed fracture characterizations[END_REF][START_REF] Fourno | A continuum voxel approach to model flow in 3d fault networks: A new way to obtain up-scaled hydraulic conductivity tensors of grid cells[END_REF] or discrete models [START_REF] Erhel | Flow simulation in three-dimensional discrete fracture networks[END_REF][START_REF] Hyman | Conforming delaunay triangulation of stochastically generated three dimensional discrete fracture networks: a feature rejection algorithm for meshing strategy[END_REF][START_REF] Berrone | Non-stationary transport phenomena in networks of fractures: effective simulations and stochastic analysis[END_REF]. When fractures are densely distributed and the length of fractures is of the same order of magnitude, continuous models can be used safely, because there is a representative elementary volume (REV) [START_REF] Long | Porous media equivalents for networks of discontinuous fractures[END_REF][START_REF] Neuman | A proposed conceptual framework and methodology for investigating flow and transport in swedish crystalline rocks[END_REF]. At scales larger than the REV, the contribution of fractures may be captured in equivalent parameters that are obtained via upscaling processes [START_REF] Karimi-Fard | Generation of coarse-scale continuum flow models from detailed fracture characterizations[END_REF][START_REF] Matthai | Upscaling two-phase flow in naturally fractured reservoirs[END_REF][START_REF] Jourdain | Upscaling permeability for fractured concrete: meso-macro numerical approach coupled to strong discontinuities[END_REF]. On the other hand, for sparsely distributed fracture networks, possibly no REV exists; the solute is carried along connected preferential pathways of high permeability. In that situation, in order to capture the effects related to the fracture sparseness, the domain should be modeled by discrete approaches. The basic idea is to keep an explicit representation of fracture networks. In this work, we focus on solute transport in discrete fracture networks (DFNs) assuming that flow only occurs through fracture planes, the matrix surrounding fractures is disregarded.

Fractures are generally considered as lower-dimensional objects, i.e. 1D or 2D objects embedded in a 3D domain. Pipe-network approaches, which are based on the assumption of flow channeling within DFNs, model fractures as 1D-in-3D objects. With this strong topological simplification, pipe-network models allow to significantly reduce the computational cost for large-scale reservoirs. The accuracy of such models has been investigated for steady-state flow and transport [START_REF] Cacas | Modeling fracture flow with a stochastic discrete fracture network: calibration and validation: 1. The flow model[END_REF][START_REF] Dershowitz | Derivation of equivalent pipe network analogues for three-dimensional discrete fracture networks by the boundary element method[END_REF][START_REF] Gylling | The channel network model-a tool for transport simulations in fractured media[END_REF]. [START_REF] Noetinger | A quasi steady state method for solving transient Darcy flow in complex 3D fractured networks[END_REF] recently proposed a quasi steady-state method for solving transient Darcy flow in complex 3D fractured networks only. This work is then extended for accounting for transfers arising from the matrix [START_REF] Noetinger | A quasi steady state method for solving transient Darcy flow in complex 3D fractured networks accounting for matrix to fracture flow[END_REF]. At the lowest order of approximation, this method needs only one unknown that is assigned at each identified intersection location between fractures. The pipe network model is then recovered. Nevertheless, further studies should focus on validation of the pipe-models for transient flow and transport.

Another simplifiying assumption is to consider fractures as planar objects. High-accuracy simulations of flow and transport within DFNs can be obtained on conforming or non-conforming meshes. By using nonconforming meshes, advanced numerical approaches, such as the PDE-constrained optimization approach [START_REF] Berrone | A pde-constrained optimization formulation for discrete fracture network flows[END_REF][START_REF] Benedetto | The virtual element method for discrete fracture network simulations[END_REF][START_REF] Benedetto | A globally conforming method for solving flow in discrete fracture networks using the virtual element method[END_REF][START_REF] Berrone | Non-stationary transport phenomena in networks of fractures: effective simulations and stochastic analysis[END_REF] and the Mortar method [START_REF] Erhel | Flow simulation in three-dimensional discrete fracture networks[END_REF][START_REF] Pichot | A mixed hybrid mortar method for solving flow in discrete fracture networks[END_REF][START_REF] Pichot | A generalized mixed hybrid mortar method for solving flow in stochastic discrete fracture networks[END_REF], allow to alleviate the meshing-related difficulties. Nevertheless, applying such a method may lead to a high additional computational effort and requires adaptations of classical numerical schemes to non-conforming meshes. Therefore, although creating high-quality conforming meshes for such networks is a non-trivial task due to geometrical complexities, it is still being dealt with by many recent studies [START_REF] Hyman | Conforming delaunay triangulation of stochastically generated three dimensional discrete fracture networks: a feature rejection algorithm for meshing strategy[END_REF][START_REF] Hyman | dfnWorks: A discrete fracture network framework for modeling subsurface flow and transport[END_REF][START_REF] Huang | A numerical method for simulating fluid flow through 3-d fracture networks[END_REF][START_REF] Fourno | An original and useful approach to mesh a discrete fracture network using a delaunay triangulation: Application on flow and transport upscaling from characterization scale to reservoir scale[END_REF]. [START_REF] Hyman | Conforming delaunay triangulation of stochastically generated three dimensional discrete fracture networks: a feature rejection algorithm for meshing strategy[END_REF] developed the feature rejection algorithm for meshing (FRAM) to generate threedimensional conforming grids for DFNs. This approach is based on constraining generations of DFNs to exclude arbitrarily small features from the network, only features greater than a user-prescribed minimum length scale h will be generated. Nonetheless, by removing small fractures accuracy of simulations can be affected due to the loss of fracture connectivity within the resulting mesh. [START_REF] Huang | A numerical method for simulating fluid flow through 3-d fracture networks[END_REF][START_REF] Fourno | An original and useful approach to mesh a discrete fracture network using a delaunay triangulation: Application on flow and transport upscaling from characterization scale to reservoir scale[END_REF] recently proposed a simple yet efficient meshing approach for DFNs, the so-called "Fracture Cut Method for Meshing" (FraC). The main idea of the approach is to decompose each fracture into a set of connected closed contours, between which common segments will be discretized in a conforming manner.

Even tiny fractures are taken into account to conserve the DFN connectivity, and a strategy for moving or adding intersecting points is applied to ensure an acceptable quality of the final mesh.

Furthermore, the resulting conforming meshes obtained from these meshing approaches have often a complex "non-manifold" topology [START_REF] Sander | The Dune FoamGrid implementation for surface and network grids[END_REF]. Handling such a situation usually requires the introduction of additional data fields and logic with respect to standard manifold grids. Some finite element codes using Mixed and Hybrid Finite Element Methods (MHFEM), for instance Cast3M (2017) [START_REF] Fourno | An original and useful approach to mesh a discrete fracture network using a delaunay triangulation: Application on flow and transport upscaling from characterization scale to reservoir scale[END_REF], allow to avoid this extra implementation effort. On the other hand, standard finite volume codes often only allow for manifold topologies that are predominant in realistic applications and do not address non-manifold issues. Currently, there are few finite volume simulators for subsurface flow and transport, e.g. PFLOTRAN [START_REF] Lichtner | PFLOTRAN User Manual[END_REF][START_REF] Hammond | Evaluating the performance of parallel subsurface simulators: An illustrative example with PFLOTRAN[END_REF] and DuMux [START_REF] Flemisch | Dumu x : DUNE for multi-{phase, component, scale, physics,...} flow and transport in porous media[END_REF], able to deal with non-manifold topologies. Resulting meshes from the FRAM method were used to obtain single-phase flow solutions using PFLOTRAN, followed by Lagrangian-based particle transport modeling [START_REF] Hyman | Conforming delaunay triangulation of stochastically generated three dimensional discrete fracture networks: a feature rejection algorithm for meshing strategy[END_REF][START_REF] Makedonska | Particle tracking approach for transport in three-dimensional discrete fracture networks[END_REF][START_REF] Karra | Effect of advective flow in fractures and matrix diffusion on natural gas production[END_REF][START_REF] Hyman | dfnWorks: A discrete fracture network framework for modeling subsurface flow and transport[END_REF]. The objective of this paper is to describe an efficient workflow that provides high-accuracy solute transport simulations using conforming meshes. The proposed workflow uses the FraC approach for meshing and DuMux for simulation frameworks in order to profit from its Euclidean-based transport modeling capability.

It is noteworthy that other advanced conformal mesh generators such as dfnGen [START_REF] Hyman | Conforming delaunay triangulation of stochastically generated three dimensional discrete fracture networks: a feature rejection algorithm for meshing strategy[END_REF] could be employed for the network generation and meshing, the FraC approach was chosen simply because of its user-friendliness. The workflow is able to handle with non-trivial DFN configurations, which remain a major challenge for classical methods. Furthermore, the numerical framework can easily be extended to consider multiphase, multi-component flow, transient single-phase flow or reactive transport as well as other hydrological modelling issues using DFN models.

The contribution is organized as follows. Section 2 provides the mathematical formulation and numerical methods for transport processes using the DuMux [START_REF] Flemisch | Dumu x : DUNE for multi-{phase, component, scale, physics,...} flow and transport in porous media[END_REF] code. Section 3 is focused on the methodology of the FraC approach used to obtain conforming meshes for DFNs. Numerical results of solute transport on simple and complex DFNs are shown in Section 4, both for model validation and demonstration purposes. For validation purpose, solutions of solute transport obtained with DuMux on a simple singlefracture DFN are compared to analytical solutions. Another validation on a more complex 7-fracture DFN is also proposed using numerical results from [START_REF] Ahmed | Three-dimensional control-volume distributed multi-point flux approximation coupled with a lower-dimensional surface fracture model[END_REF] as reference. The second part of the section is dedicated to case studies of transport processes on a synthetic benchmark 33-fracture network and on a complex and realistic DFN in order to assess the reliability of the approach. Finally, we conclude the paper in Section 5 by giving some directions for further studies.

Mathematical formulation and numerical methods

The motion of a non-reactive solute in the subsurface can be described by the advection-diffusion equation, [START_REF] Bear | Flow and contaminant transport in fractured rock[END_REF]as follows:

ω ∂C ∂t = ∇.(D.∇C -Cv) + q, (1)
where ω (-), D (L 2 /T) are the porosity and the dispersion coefficient tensor for the solute in water, respectively. C (M/L 3 ) stands here for the solute concentration, v (L/T) for the velocity field and q (M/(L 3 .T))

for the external sink and/or source term. The velocity field v is given by the Darcy's law:

v = - K µ (∇p + ρg), (2)
in which p (M/(L.T 2 )) denotes the water pressure, g the gravity acceleration vector, K (L 2 ) the permeability tensor, µ (M/(L.T)) the dynamic viscosity of water and ρ (M/L 3 ) is its density. The water flow is assumed to be in the steady-state, therefore Equation (1) can be recast as:

ω ∂C ∂t = D∇ 2 .C -v.∇C + q.
(3)

The one-dimensional form of Equation (3) reads as follows

ω ∂C ∂t = D L ∂ 2 C ∂x 2 -v. ∂C ∂x + q. ( 4
)
Although several physical models available in DuMux allow to study a wide range of subsurface problems, spanning from single-phase single-component isothermal to three-phase three-component non-isothermal flow, non-reactive solute transport is however not implemented in the code. Fortunately, there is the 1p2c v is constant, therefore, the continuity equation of the 1p2c model becomes the advection-diffusion equation.

The non-manifold issue of the mesh topology is solved in DuMux (2017) via the DUNE grid manager DUNE::FoamGrid [START_REF] Sander | The Dune FoamGrid implementation for surface and network grids[END_REF]. A vertex-centered finite volume scheme, also known as the box method (BM) (see, for example, [START_REF] Huber | Node-centered finite volume discretizations for the numerical simulation of multiphase flow in heterogeneous porous media[END_REF]), is chosen for solving flow and transport in porous media. Concretely, the computational domain T h is discretized by a finite element grid, the so-called "primary mesh", where h is the mesh size. Since fractures are considered as surface objects, the unit cell could be either triangles or quadrilaterals. A "secondary mesh" which describes the finite volumes is constructed based on the finite element grid. A finite volume B i , also known as a box, is obtained by connecting the gravity center of element i with associated edge midpoints (see Figure 1). Thus, each finite volume is associated with a node of the primary mesh. Furthermore, the model variables are defined for control-volumes but not finite elements, they are thus assigned to the grid nodes. Moreover, local physical properties of porous media, such as permeability and porosity, are specified element-wise, the corresponding values assigned to primary mesh secondary mesh node j node i

B i
Figure 1: Spatial discretization of the box method including two meshes: "primary mesh" in dashed lines and "secondary mesh" in continuous lines. The gray polygon represents the box Bi corresponding to the node i. et al. (2015); [START_REF] Makedonska | Particle tracking approach for transport in three-dimensional discrete fracture networks[END_REF]. We refer to the online version of the paper for the colors mentioned in this caption.

each control volume are determined by a weighted average of the properties of the sub-control volumes.

The BM inherits advantages from both the finite element and finite volume method: the finite element property makes it possible to deal with unstructured grids and the finite volume method ensures mass conservation. In addition, for flow and transport simulations on non-manifold conforming meshes, the BM unconditionally guarantees pressure continuity and mass conservation at the intersections between planar surfaces. Assume that we are considering two intersecting fractures (the red and blue ones) as in Figure 2. Since the mesh is conforming along the intersection line, the control volumes associated to the node i within each fracture share common edges (e i ij and e i ij ) of the sub-control volumes, on which the fluxes are computed. Then, mass conservation at the fracture intersection can naturally be ensured, and the box method is therefore chosen as discretization method for the following flow and transport simulations.

Moreover, although the temporal discretization may be done using an explicit or an implicit Euler scheme in DuMux, only implicit models are retained in our applications. The time-step size is controlled in an adaptive manner based on the number of iterations needed for the Newton's method to converge for the last time integration [START_REF] Flemisch | Dumu x : DUNE for multi-{phase, component, scale, physics,...} flow and transport in porous media[END_REF].

Mesh generation

Description of the meshing method

The FraC method involves two primary stages, i.e. pre-meshing and meshing stages (see Figure 3. The pre-meshing stage includes two steps that are all implemented using the Computational Geometry Algorithms Library (CGAL): (i) data reading/ importing (S 1 ); (ii) fracture intersections and cutting (S 2 ). The triangulation (S 3 ) and merging closed contours (S 4 ) processes gather together in the meshing stage, which is done by using Los Alamos Grid Toolbox (2013) because CGAL is unfortunately not able to store nonmanifold data. The main features in each stage can be summarized as follows.

(i) S 1 : Data reading/importing

The input for FraC comes generally from fractured site characterizations that provide distributions of fracture forms, orientations, dimensions, and spatial locations. Fractures are then representatively approximated by concave polygons. Each polygon F i will then be discretized with a user-defined characteristic length h i , that is subsequently used as the target edge length for F i in the meshing stage. It means that the edge length of the mesh triangles is of the same order of magnitude as h i . Hence, resolution of the final mesh can be easily controlled by adjusting h i ; high values of h i result in coarse meshes and reducing h i yields mesh refinement. The characteristic length h i is fairly similar to the minimum length scale in the FRAM approach [START_REF] Hyman | Conforming delaunay triangulation of stochastically generated three dimensional discrete fracture networks: a feature rejection algorithm for meshing strategy[END_REF][START_REF] Hyman | dfnWorks: A discrete fracture network framework for modeling subsurface flow and transport[END_REF].

(ii) S 2 : Fracture intersection and cutting

The next step involves finding intersections between the polygons associated with the fractures. An octreelike approach [START_REF] Khvoenkova | An optimal method to model transient flows in 3D discrete fracture network[END_REF][START_REF] Hyman | Conforming delaunay triangulation of stochastically generated three dimensional discrete fracture networks: a feature rejection algorithm for meshing strategy[END_REF] is used to optimize the localization of fractures in the space. More precisely, an axis-aligned minimum bounding box BB i is created around F i where BB i is defined by the minimal and maximal value of the corresponding coordinate of the polygon vertices. Then a preliminary test is performed to check intersection between the bounding boxes BB i and BB j of two fractures F i and F j respectively. The number of intersection tests between bounding boxes performed for each fracture can be reduced using a three-dimensional octree structure based on fracture locations in the corridor. This approach works well if the size of the fractures is almost of the same order of magnitude.

If BB i and BB j are not touching in space, obviously F i and F j do not intersect. Otherwise, if BB i and BB j intersect, the intersection test between F i and F j is then performed. The intersection line L ij between the F i and F j planes is first determined. Let ∂F α be the boundary of F α with α ∈ {i, j}, then finding intersections between the fractures becomes finding intersections between ∂F i , ∂F i with L ij . We denote n (α) ij the number of intersection points between ∂F α with L ij . So far, two cases can occur:

-

n (i) ij = 0 or n (j) ij = 0: The fractures do not intersect. -n (i) ij , n (j) ij > 0: In most cases, four intersection points p (1) i , p (2) i , p (1) j , p (2) j (n (i) ij = n (j) ij = 2) may be found, where p (1,2) i lie on F i and p (1,2) j lie on F j . Degenerated cases of tangent fractures where n (i) ij or n (j)
ij or both of them equal to one may also exist. The intersection points are then sorted to find out endpoints. Let I ij be the segment connecting two endpoints, I i be the portion of I ij which lies only on F i but not on F j , I j be the portion which lies only on F j but not on F i , and I ij be the portion which lies on both F i and F j (see Figure 4). Finally, the fractures intersect only if

I ij = Ø.
The fracture cutting step expresses the key idea of the FraC approach, it is solely applied for intersecting fractures. Each of two intersecting fractures F i and F j will then be cut along the extended intersection I ij to obtain interconnected closed contours, subsequently referred to as "sub-fractures". To guarantee the homogeneity of the final mesh, the sub-fractures and more precisely, the segment I ij should also be discretized according to the mesh characteristic length of the original fracture. For that, special focus is put on I ij which is also sub-divided in a homogeneous or adaptive manner. In the first case we define the target edge length h ij = min(h i , h j ) common for I ij and its sub-segment including I i , I j , I ij . Here, h i and h j are the target edge length of F i and F j respectively. In the second case, only I ij will be discretized with h ij , the target edge length for I i and I j remains h i and h j respectively. At the endpoints of I i , I j and I ij , two strategies are flexibly applied: depending on neighborhood area, discretization points of ∂F α are either moved (Type I) or created (Type II) to model the intersection points (Figure 5a). These strategies allow avoiding too short edges resulting in grid cells with extremely high aspect ratio in the final mesh.

I i (F i ) I ij I ij I j (F j )
Once all edges of sub-fractures are discretized, new intersection tests between them and other fractures/ sub-fractures are carried out. By a series of cutting steps, the original fractures are decomposed into a large number of connected closed contours (Figure 5b) that will be used as input for the triangulation step S 3 .

The main benefit of series of cutting steps is that multiple-intersection issues vanish. For example, a problem of triple intersections between three fractures F i , F j , F k , is transformed to standard problems of finding intersections between F k and the sub-fractures of F i and F j . The drawback of this strategy is that it increases the number of intersection tests and, consequently, the computational time. The consequence of this feature of the method is however beyond the scope of the present paper and will be addressed in another work.

(iii) S 3 : Triangulation

In this step, each sub-fracture will be meshed according to its own target edge length h i (Figure 5c).

LaGriT uses a conforming Delaunay triangulation algorithm to ensure that the line of intersection between any two sub-fractures is preserved. More details can be found in [START_REF] Murphy | A point-placement strategy for conforming delaunay tetrahedralization[END_REF] or [START_REF] Hyman | Conforming delaunay triangulation of stochastically generated three dimensional discrete fracture networks: a feature rejection algorithm for meshing strategy[END_REF]. meshed closed contours from step S 3 will be merged into a single triangulation data structure, in which there will probably be duplicate vertices along the intersection lines (Figure 5d). These duplicates are detected and removed, resulting in a final fully connected mesh, that can be saved to *.inp files which can be read by the AVS graphics program, both in binary or ASCII format. A mesh conversion of *.inp files to *.dgf file (DUNE Grid Format, DGF) is then required to make output mesh files compatible with the DuMux code.

I ij I i (F i ) (F j )
A detailed description of the DGF and its utilization can be found on the homepage of DUNE.

Meshing examples

In this section, two applications of the FraC approach for meshing complex DFNs will be shown. The first example is a benchmark DFN, which has also been used by many previous studies [START_REF] Khvoenkova | An optimal method to model transient flows in 3D discrete fracture network[END_REF][START_REF] Fourno | An original and useful approach to mesh a discrete fracture network using a delaunay triangulation: Application on flow and transport upscaling from characterization scale to reservoir scale[END_REF], and the second one is a semi-synthetic DFN where the network's properties are closely based on the geological data from the fictional but realistic Bloemendaal reservoir [START_REF] Verscheure | Joint inversion of fracture model properties for co2 storage monitoring or oil recovery history matching[END_REF].

3.2.1. Benchmark 33-fracture DFN (DFN 1 )
The benchmark fracture network contains 33 disc-shaped fractures inside a 3 m × 3 m × 3 m cube.

A common characteristic length h is set for all fractures. Several meshing realizations are performed with decreasing h to obtain a set of meshes with different degrees of refinement. Concretely h ranges from 9.6 × 10 -1 m to 3.75 × 10 -3 m resulting in meshes that contains from 2 × 10 3 to 8.0 × 10 6 grid cells. In DFN 1a corresponding to h = 0.12 m (left) and a fine grid DFN 1b corresponding to h = 0.03 m (right). This

shows the ability of our code to generate not only extremely fine meshes but also very coarse meshes of the DFN while preserving the fracture network topology. The inset in this figure exhibits the mesh conformity at a multiple-intersection location.

Semi-synthetic DFN of the Bloemendaal reservoir (DFN 2 )

The semi-synthetic DFN consists of more than 50,000 rectangular-shaped fractures that lie on a domain of 12 km × 15 km × 1.4 km. The flowing zone depth of the reservoir is about 200 m (Figures 7a and7b).

We use the algorithm of the FracaFlow code (Beicip-FranLab, 2017) to generate the fracture network that includes two principal fracture sets: one set primarily oriented along the x-axis (red fractures, Figure 7a)

with an average fracture spacing of 10.0 m and another primarily oriented along the y-axis (blue fractures, Figure 7a) with an average fracture spacing of 7.5 m. For the above colors, we refer the reader to the online version of the paper.

In practice, simulations using the DFN concept rarely model the whole fracture network of the reservoir.

For example, in the oil field, the most common tests based on DFN conceptual models, such as well-tests, flow-meters or interference tests, consider the influence zone of a few kilometers around the wells [START_REF] Bourbiaux | An integrated workflow to account for multi-scale fractures in reservoir simulation models: implementation and benefits[END_REF]. Therefore, in this case we focus solely on a part of the semi-synthetic DFN: only fractures inside a bounding box of L 0 × L 0 × 1.4 km (red box, Figure 7b) are taken into account, with L 0 initially equal to 1 km; the DFN contains around 5 × 10 2 fractures. An example of resulting mesh is illustrated in Figure 7c.

We then expand the DFN by increasing L 0 to ten kilometers. The number of fractures and consequently the intersections between them inside the DFN are thus increased. Figure 7d grids.

Numerical results

In this section, numerical results are presented to evaluate the ability of DuMux to simulate solute transport in DFNs with increasing levels of topological complexity. Three validation tests are carried-out: For simplicity, the DFNs in all examples are assumed to be homogeneous with a constant fracture aperture of b = 10 mm for Example 3 and b = 1 mm for the others. An isotropic permeability for all fractures K = 10 3 mD is applied for Example 2 to be consistent with that of [START_REF] Ahmed | Three-dimensional control-volume distributed multi-point flux approximation coupled with a lower-dimensional surface fracture model[END_REF] while, for other cases, the fractures are assigned a scalar permeability tensor K = 10 4 mD. The fluid viscosity is constant and equals µ = 10 -3 Pa.s.

Example 1: Single inclined fracture

The single fracture is an inclined planar rectangle of length 20 m (along the y-axis) and of width √ 5 m.

The angle θ between the fracture normal vector and the z-axis satisfies tanθ = 2.

A concentration C(t = 0) = 0 is applied throughout the entire domain as initial condition. Flow occurs along the y-axis by adopting a pressure gradient ∆P = 10 8 Pa along the y-direction. Tracer is continuously injected into the domain with an imposed concentration C 0 (relative concentration C = C /C0 = 1), the prescribed boundary conditions are summarized as follows:

           C = C 0 on y = 0 m, C = 0 on y = 20 m, Φ = 0 otherwise. ( 5
)
The transport analytical solution of this problem is given by [START_REF] Ogata | A solution of the differential equation of longitudinal dispersion in porous media: fluid movement in earth materials[END_REF]:

C(x, t) = C 0 2 erf c x -vt 2 √ D L t + exp vx D L erf c x + vt 2 √ D L t (6)
where C 0 is the imposed concentration on the inlet and erf c is the complementary error function. When the ratio vx/D L (i.e. the Péclet number P e) is large enough, the second term in equation ( 6) is negligibly (P e2 = 10 2 ). Mesh size of the coarsest mesh hmax = 0.4m. Discontinous line on (a) represents the approximate solution from Equation ( 7) while the one on (b) illustrates the exact solution from Equation ( 6).

small and often disregarded [START_REF] Ogata | A solution of the differential equation of longitudinal dispersion in porous media: fluid movement in earth materials[END_REF][START_REF] Bodin | Simulation and analysis of solute transport in 2d fracture/pipe networks: The solfrac program[END_REF], one obtains the approximate solution as follows:

C(x, t) = C 0 2 erf c x -vt 2 √ D L t . (7)
According to [START_REF] Ogata | A solution of the differential equation of longitudinal dispersion in porous media: fluid movement in earth materials[END_REF], only a maximum error of less than 3% between the exact and approximate solutions in Equations ( 6) and ( 7) is shown for P e > 500.

Different meshes are used to model this problem with h decreasing from h max = 0.4m, the mesh size of the coarsest mesh, to h = h max /8. Figure 8 shows a mesh example of the fracture for h = 0.4m. Two cases according to two dispersion coefficient values are studied. We consider in the first case a convection-dominant process with a low dispersion coefficient D 1 = 1 × 10 -9 m 2 .s -1 and the mean fluid velocity v = 5 × 10 -4 m.s -1 resulting in a Péclet number P e 1 = L × v/D 1 = 10 7 . In the second case, v remains unchanged but the dispersion coefficient is increased up to D 2 = 1 × 10 -4 m 2 .s -1 , the corresponding Péclet number is thus equal to P e 2 = 10 2 .

Figure 9 depicts concentration profiles obtained from numerical simulations on different meshes with variable h for both cases at t = 2 × 10 4 (s). Analytical solutions from Equation ( 7) for P e 1 = 10 7 500 and from Equation ( 6) for P e 2 = 10 2 < 500 are also illustrated. In the first case, there are offset between the analytical and numerical curves mainly at the front location (around y = 10 m). This discrepancy is possibly caused by both the error of the approximation in Equation ( 7) and numerical diffusion effects. It can be noted that, when the grid is refined, the numerical curves tend to the analytical solution. The numerical diffusion effect vanishes for the second case when the physical diffusion effect is stronger. Numerical solutions, particularly for the solution on the finest mesh, i.e. h = h max /8 = 0.05 m, are in perfect agreement with the analytical one.

The exact solution from the second case D 2 = 1×10 -4 m 2 .s -1 is now used for convergence study. Indeed, to estimate the BM convergence order, a L 2 norm of tracer concentration error e C is defined as:

e C = Σ Bi∈T h |B i | × (C ex,Bi -C Bi ) 2 Σ Bi∈T h |B i | × C 2 ex,Bi 1 2 . ( 8
)
Here, C Bi and C ex,Bi stand for the discrete and analytical concentration solution of the box B i respectively while |B i | denotes the box area.

For k ≥ 2, the convergence order col2(k) of the L 2 norm of the error on the solution between grid level k and k -1 can be evaluated as:

col2(k) = log(e C (k)/e C (k -1)) log(h(k)/h(k -1)) . ( 9
)
Figure 10 displays the convergence rate for the concentration of the box method, showing an order of convergence of O(h).

Example 2: 7-fracture DFN

In this benchmark exercise, solutions of DuMux are compared to those from [START_REF] Ahmed | Three-dimensional control-volume distributed multi-point flux approximation coupled with a lower-dimensional surface fracture model[END_REF], who use the cell-centred control-volume distributed multi-point flux approximation (CVD-MPFA) approach as discretization method. The studied DFN lies inside a 200 × 150 × 15 m 3 system. For validation purposes, simulation on 2D model from [START_REF] Ahmed | Three-dimensional control-volume distributed multi-point flux approximation coupled with a lower-dimensional surface fracture model[END_REF] is chosen as a reference. Our mesh (Figure 11a) has the same resolution as that used by [START_REF] Ahmed | Three-dimensional control-volume distributed multi-point flux approximation coupled with a lower-dimensional surface fracture model[END_REF] (figure 11b). The conforming Delaunay mesh generated by the FraC approach contains 4554 grid cells while the Ahmed's mesh contains 4662 triangular cells. Figure 11 shows the locations of the injector I and that of the two producers P 1 and P 2 . Fluid is injected through the injector at the rate of 2.739 × 10 -3 m 3 /day. The imposed pressure at both producers is set to 10 bar. Both gravity and diffusion effects are neglected. Figure 12a illustrates the result of DuMux for the tracer concentration t = 5 years, which matches perfectly with the reference result shown in Figure 12b. An excellent agreement between concentration plots at the producers P 1 w.r.t time in Figure 12c also validates our solutions.

Example 3: 33-fracture benchmark synthetic DFN

The third validation is undertaken on the DFN 1 based on previous studies by [START_REF] Fourno | A continuum voxel approach to model flow in 3d fault networks: A new way to obtain up-scaled hydraulic conductivity tensors of grid cells[END_REF][START_REF] Fourno | An original and useful approach to mesh a discrete fracture network using a delaunay triangulation: Application on flow and transport upscaling from characterization scale to reservoir scale[END_REF].

Firstly, equivalent permeabilities are computed from the numerical solutions of steady-state single-phase flow and compared with the published results. Secondly, solute transport is modeled and analyzed from a convergence point of view.

To determine the diagonal entries of the effective permeability tensor K ii , i ∈ {x, y, z}, a pressure gradient ∆P is applied on two corresponding opposite boundaries of the DFN and other boundaries are assumed to be impermeable. The equivalent permeability of a DFN can be calculated by applying the inverse Darcy's method when the fluid flow is at steady state [START_REF] Zimmerman | Effective transmissivity of two-dimensional fracture networks[END_REF].

K ii = µ × Q ii A × ∆P , (10)
where Q ii is the corresponding macroscopic flux and A is the cross-sectional area.

Note that the off-diagonal terms of the equivalent permeability tensor, i.e. K ij (i, j ∈ {x, y, z}, i = j), can also be computed with a different setup of the boundary conditions. However, the calculation of off-diagonal equivalent permeability is out-of-scope of this work, we only investigate the principal terms K xx , K yy and K zz ; the solutions will be used for validation purposes. The DFN is assumed to be homogeneous and isotropic with the intrinsic permeability is set to 10 4 mD.

Figure 13 compares the numerical results of the equivalent permeabilities obtained from simulations with different codes including DuMux, FracaFlow [START_REF] Khvoenkova | An optimal method to model transient flows in 3D discrete fracture network[END_REF], MP-Frac (De Dreuzy et al., 2013) and Cast3M [START_REF] Fourno | A continuum voxel approach to model flow in 3d fault networks: A new way to obtain up-scaled hydraulic conductivity tensors of grid cells[END_REF]. One can note that DuMux's results for the effective permeabilities along all three principal directions are similar to those obtained by MP-Frac and Cast3M. The results obtained by FracaFlow are slightly different from the others, but remain of the same order of magnitude.

These discrepancies are caused by the loss of the DFN connectivity in the FracaFlow meshes. Indeed, [START_REF] Khvoenkova | An optimal method to model transient flows in 3D discrete fracture network[END_REF] modeled disc-shaped fractures by planar octagons, therefore some DFN intersections cannot be captured resulting in a consistent underestimation of the equivalent permeabilities along the y and z directions. On the contrary, it appears that the fracture connectivity in the DFN is accurately captured by FraC into resulting meshes, that are eligible for flow and transport simulations using standard simulators.

In the experiment of transport through the DFN 1 , solute is continuously injected from an inflow boundary of the DFN. Flow occurs through the DFN along the x-axis by applying a pressure gradient ∆P = 10 5 Pa along this axis. The dispersion coefficient is isotropic, D = 10 -9 m 2 .s -1 . Simulations are run on different meshes with the grid size ranging from 1×10 5 elements (coarsest mesh) to 8×10 6 elements (finest mesh). The characteristic length h varies from h = h max = 0.03m for the coarsest mesh to h = h max /8 = 3.75 × 10 -3 m for the finest mesh, Figure 14a illustrates the pressure field and iso-value lines of the pressure while Figure 14b displays the norm of flow velocity reconstructed from numerical results using a weighting of the sub control volume velocities. In Figure 14b, blue colors represent low values (down to 1.877 ×10 -17 m.s -1 ); warmer colors indicate higher values (up to 3.527 ×10 -3 m.s -1 ). Color scale is set maximum at 1 ×10 -4 m.s -1 to distinguish high-and low-flow regions. We can easily notice from Figure 14 narrowly-spaced pressure iso-value lines for high-flow regions (for instance, in fractures close to the inlet) and widely-spaced ones for low-flow areas. Figure 15 shows the tracer concentration distribution of the numerical simulation obtained through DuMux. It reports the concentration evolution in the DFN at different timescales: after a short time of injection, t 1 = 1 hour (Figure 15a), an intermediate time t 2 = 5 hours (Figure 15b) and a relatively long time t 3 = 25 hours (Figure 15c) and t 4 = 30 hours (Figure 15d) when the solute has reached the outflow boundary. We notice that for long times, i.e. after 25 hours, the solute appears nearly in the whole DFN; in addition, the solute displacement occurs slowly in the DFN, which is expressed through the similar solute distributions on Figures 15c and15d. Clearly, this slow motion is the result of both low diffusion effect and low velocity: the distribution of low concentration areas fits correctly with low-flow regions in Figure 14b.

Figure 16 exhibits the concentration evolution at the observation point A of coordinate (2.35m, 1.35m, 0m) obtained from simulations run on 4 grids, the mesh refinement ratio is a factor of 2. Obviously, the arrival time for coarse meshes is earlier than for fine meshes due to numerical diffusion. The observed curves obtained from the finest meshes, i.e. the third and fourth grid levels, overlap each other (discrepancy less than 1%). For instance, the relative concentration at this point at t = 15h for the two coarsest meshes is about 0.85 and 0.79 respectively, and for both two finest mesh is roughly equal to 0.73. Through this figure, the convergence of the solution with respect to mesh refinement is also judged.

Example 4: realistic DFN of the Bloemendaal reservoir

We study in the last example solute transport within a complex realistic DFN of the Bloemendaal reservoir with L 0 = 3 km. The network contains 3088 fractures and the mesh consists of about 1.1 million grid cells.

As mentioned above, the fracture network is well connected and seems to be homogenizable.

Single-phase flow simulations to compute effective permeabilities of the DFN block as described in subsection 4.3 are also carried out, giving K xx = 4.58 × 10 -18 m 2 and K yy = 3.0 × 10 -18 m 2 . It is worth noticing that, although there are more y-axis-oriented fractures than x-axis-oriented ones within the DFN, the effective permeability of DFN along the x-axis is surprisingly about 1.5 times greater than the y-axis one. This can be explained by the fact that the fracture network is more connected and therefore more conductive in the x-direction than in the y-direction.

Fluid is injected into the domain through an injection well located in the middle of the DFN, the injection rate is 5 × 10 -2 m 3 /day. A Dirichlet condition P = 10 bar is applied for all lateral boundaries while a no-flow Neumann condition is applied for top and bottom boundaries of the domain. As in Example 3, we consider the case where the dispersion coefficient is isotropic, D = 10 -9 m 2 .s -1 . Figure 17 shows solute concentration distribution at t = 10 years (Figure 17a), t = 100 years (Figure 17b) and t = 250 years (Figure 17c).

We note from these figures that there are low-concentration regions which exist locally within the solute plume. However, disregarding local phenomena, the DFN behaves as a continuous homogeneous medium at macroscopic scale: the concentration distribution has an elliptical form.

In practice, the overall development of a solute plume within a 2D or 3D heterogeneous medium may be precisely measured through a moment analysis. This approach is an appropriate tool for examining spatial effects of the domain on developing solute plumes, for identifying loss of mass during flow process and for assessing large-scale dispersive behavior of the flowing domain itself [START_REF] Freyberg | A natural gradient experiment on solute transport in a sand aquifer: 2. spatial moments and the advection and dispersion of nonreactive tracers[END_REF][START_REF] Tompson | Numerical simulation of solute transport in three-dimensional, randomly heterogeneous porous media[END_REF].

Let us remind that the ijk-th moment of the concentration distribution in space, M ijk is defined as [START_REF] Aris | On the dispersion of a solute in a fluid flowing through a tube[END_REF])

M ijk (t) = ∞ -∞ ∞ -∞ ∞ -∞ ωC(x, y, z, t)x i y j z k dxdydz. ( 11
)
where ω [-] and C [-] denote the porosity and the relative tracer concentration, respectively. The zerothorder moment M 000 (t) is the total mass of the solute in the domain. The first-order moment about the origin, normalized by M 000 , describes the propagation of the center of mass O c (x c (t), y c (t), z c (t)) of the solute plume:

x c (t) = M 100 /M 000 ; y c (t) = M 010 /M 000 ; z c (t) = M 001 /M 000 . (12)
The centered second-order tensor defines a spatial covariance tensor reflecting the spreading of mass about O c :

σ =      σ xx σ xy σ xz σ yx σ yy σ yz σ zx σ zy σ zz      (13)
where

σ xx = M 200 M 000 -x 2 c , σ yy = M 200 M 000 -y 2 c , σ zz = M 200 M 000 -z 2 c , σ xy = σ yx = M 110 M 000 -x c y c , σ xz = σ zx = M 101 M 000 -x c z c , σ yz = σ zy = M 011 M 000 -y c z c .
Diagonalization of σ provides the orientation of principal axes (x , y , z ) as well as the eigenvalues of the spatial covariance tensor σ x x , σ y y , σ z z .

In addition, the derivative over time of the covariance is proportional to the dispersion tensor [START_REF] Aris | On the dispersion of a solute in a fluid flowing through a tube[END_REF][START_REF] Freyberg | A natural gradient experiment on solute transport in a sand aquifer: 2. spatial moments and the advection and dispersion of nonreactive tracers[END_REF], one could expect that the square root of the ratio between the eigenvalues of the centered second-order moments, σ x x /σ y y , reflects the anisotropy of the elliptic solute plume for the homogenized media.

Figure 18a reports the variation versus time of the centered second-order moments σ xx , σ yy , σ zz and that of the eigenvalues of the σ tensor, i.e. σ x x , σ y y , σ z z . We note that along the x-and y-axis, the centered second-order moments and the corresponding eigenvalue match correctly, only a difference between σ zz and σ z z is observed. However, since both σ zz and σ z z are roughly two orders of magnitude less than those along the x-and y-axis, changing of the z -axis with respect to the z-axis can also be negligible. This proves that the principal axes for the developing direction of the solute plume x , y , z and the original orientations

x, y, z nearly coincide.

Figure 18b illustrates the time-evolution of the σ x x /σ y y ratio. One can notice that σ x x /σ y y is fairly correlated with that of the effective permeabilities K xx /K yy ; a discrepancy exists between them but it is still acceptable and generally remains stable over time. At early times when the amount of solute in the domain remains low, σ x x /σ y y is a little greater than the average value because local fracture connectivity could have significant effects on the ensemble solute distribution. After 50 years, local effects vanish and the curve becomes relatively stable until 100 years. However, σ x x /σ y y increases slightly again after 150 years because the solute has already reached lateral boundaries of the DFN.

Summary and discussion

The purpose of this work is to obtain high-accuracy simulations of solute transport within large, threedimensional Discrete Fracture Networks (DFNs). A workflow is described, including both the mesh generation for DFNs using the conforming mesh method FraC and the simulation framework using the DuMux code.

The meshing approach FraC, initially developed by [START_REF] Fourno | An original and useful approach to mesh a discrete fracture network using a delaunay triangulation: Application on flow and transport upscaling from characterization scale to reservoir scale[END_REF], relies on the idea of decomposing fractures along extended intersection segments between fractures to obtain a set of interconnected closedW contours, that will then be discretized with the same point density on shared segments between them to ensure the conformity of the final mesh. Conforming meshes were then used for transport simulations with DuMux using a vertex-centered discretization method, or box method. A number of tests were conducted for both validation and demonstration purposes. Validation tests gave satisfactory results: (i) numerical solutions of solute transport on a DFN of single inclined fracture are in agreement with analytical solutions by [START_REF] Ogata | A solution of the differential equation of longitudinal dispersion in porous media: fluid movement in earth materials[END_REF], an a-posteriori convergence analysis showed a convergence rate of O(h) for concentration, where h is the mesh size; (ii) results on 7-fracture DFN obtained by DuMux using the box method match perfectly with those obtained by the cell-centred control-volume distributed multi-point flux approximation (CVD-MPFA) method [START_REF] Ahmed | Three-dimensional control-volume distributed multi-point flux approximation coupled with a lower-dimensional surface fracture model[END_REF]; and (iii) simulations on a benchmark 33-fracture DFN showed convergence of solution with increasing mesh refinement. Finally, we considered development of solute plume injected into a DFN produced from the Bloemendaal reservoir data that consists of thousands of fractures. This test case reveals that our method is able to deal with large-scale DFN containing a moderate fracture number.

High-accuracy simulations obtained through this work appear to be very useful for validation of earlier theoretical studies on simplified models, for instance pipe-networks models [START_REF] Noetinger | A quasi steady state method for solving transient Darcy flow in complex 3D fractured networks[END_REF][START_REF] Noetinger | A quasi steady state method for solving transient Darcy flow in complex 3D fractured networks accounting for matrix to fracture flow[END_REF], as well as for benchmarking purposes to test the reliability and performance of commercial or academic simulators. These results also motivate many directions for further works, both for the mesh generation and at the numerical stage. Objectively speaking, numerical simulations on large grids using DuMux could lead to prohibitive CPU times. In order to overcome this issue, several measures could be examined, either adopting the adaptive mesh refinement (AMR) approach [START_REF] Sander | The Dune FoamGrid implementation for surface and network grids[END_REF] or distributing the workload on different resources by using parallel computing.

It has also been shown that one can only run massively parallel simulations in DuMux if the DUNE grid manager allows it. Nevertheless, to the best of the authors' knowledge, DUNE::FoamGrid is not yet parallelizable. On the contrary, the AMR approach is supported by DUNE::FoamGrid, i.e. simulations can be run on time-evolving meshes with local mesh refinement/coarsening inside DuMux. On the other hand, the AMR approach can be included into FraC frameworks itself without major difficulties. In principle, only vicinity areas of intersections should be refined and a coarse meshing is applied for the remaining. This work can easily be carried out within LaGriT as in the dfnWorks framework [START_REF] Hyman | dfnWorks: A discrete fracture network framework for modeling subsurface flow and transport[END_REF].

Furthermore, although all DFNs in this paper are assumed to be homogeneous for simplicity, the framework is, however, able to handle heterogeneity. More precisely, the step S 4 of the FraC approach provides flexibility for both fracture-to-fracture and in-fracture heterogeneity. Heterogeneous parameters such as fracture aperture, transmissivity, permeability and porosity are assigned to each mesh element. In the case of fracture-to-fracture heterogeneity, grid cells in a same fracture are assigned a common material index matID, the heterogeneous parameters are then distributed according to matID. To model inside-fracture heterogeneity, spatially correlated random fields for fracture aperture are first generated using the R package gstat [START_REF] Pebesma | Multivariable geostatistics in S: the gstat package[END_REF]. The correlation length of random fields λ is proportional to fracture size. Sequential unconditional Gaussian simulations are carried out to generate random fields for fracture aperture constrained by a simple kriging variogram. The permeability and transmissivity fields are then derived from the aperture distribution using the well-known transmissivity-aperture cubic law [START_REF] Adler | Fractured porous media[END_REF]. Heterogeneous data are then directly attached into *.dgf grid files and will be read in the DuMux via the DUNE grid interface DUNE-grid [START_REF] Bastian | A generic grid interface for parallel and adaptive scientific computing. Part II: Implementation and tests in DUNE[END_REF]. Preliminary results show that simulations on heterogeneous DFN are in general more time-consuming than on homogeneous ones, hence once again, this highlights the importance of parallel computing.

In addition, for the purpose of various applications, more sophisticated problems including, for instance, transient single-phase flow, multi-phase multi-component flow, sorption or reactive transport could also be examined.

Finally, another direction to extend our works is to account for the matrix to fracture flow. That may pose challenges for both the mesh generation and simulation procedure. For instance, meshing dense-distributed DFNs with matrix could yield poor quality tetrahedrons that would dramatically degrade numerical performances. Furthermore, it requires advanced discretization methods able to deal with the exchange at matrix-fracture interfaces. For that, Vertex Approximate Gradient (VAG) Scheme [START_REF] Brenner | Vertex approximate gradient scheme for hybrid dimensional two-phase darcy flows in fractured porous media[END_REF] or cell-centred control-volume distributed multi-point flux approximation (CVD-MPFA) method [START_REF] Ahmed | Three-dimensional control-volume distributed multi-point flux approximation coupled with a lower-dimensional surface fracture model[END_REF] could be excellent candidates.

  model simulating single-phase flow of a compressible fluid consisting of two components where the fluid properties such as the density and the viscosity depend upon the pressure. The primary variables of the 1p2c model are the pressure p and the mole or mass fraction of dissolved components x κ . Basing on this model, transport processes can be simulated in DuMux with a slight modification of the code, in which the pressure-dependency of the fluid properties is removed. The steady-state fluid flow are considered, i.e.

Figure 2 :

 2 Figure 2: Illustration of the flux distribution for intersecting control volumes. The orange and yellow polygons belong to the red fracture, and the green and gray polygons belong to blue fracture. In this example, inflow occurs from the orange polygon and outflow occurs through the remaining polygons. This figure is reproduced from Hyman

Figure 3 :

 3 Figure3: General flow chart from 3D mesh generation of DFNs basing on geological input datasets to numerical simulations of flow and transport processes. The framework consists of two major parts: (I) Mesh generation using the FraC approach and (II) Numerical simulations in DuMux. For the first part, the implementation of the FraC approach is carried out by using CGAL and LaGriT for pre-meshing and meshing stage, respectively. Detailed steps of each stage are also listed.

Figure 4 :

 4 Figure 4: An example of two intersecting fractures Fi (blue) and Fj (gray). The extended intersection Iij (bold line) between them involves Iij (bold line), Ii and Ij (dashed lines).

  (iv) S 4 : Merging meshed sub-fractures and post-treatment This step bridges the gap between the meshing part and the simulation part of our workflow. All resulting

Figure 5 :

 5 Figure 5: Illustration of the meshing steps: step S1 -fracture intersection (a) and fracture cutting (b); step S3polygon triangulation (c); and step S4 -merging meshed polygon (d).

Figure 6 ,Figure 6 :

 66 Figure6, we show, as an example, the DFN mesh generated by FraC with two values of h: a coarse grid

Figure 7 :

 7 Figure 7: Semi-synthetic DFN: entire fracture network in lateral view (a) and in top-down view (b). The red box in (b) indicates a part of the DFN that is being meshed. An example of resulting mesh for the case L0 = 1 km where the DFN contains around 5 × 10 2 fractures is shown in (c), and (d) reports the variation of the fracture number in the bounding box (in red) w.r.t. the box size.

Figure 8 :

 8 Figure 8: A snapshot of the mesh for a single inclined planar fracture embedded in 3D domain with h = 0.4m. The green rectangular cuboid indicates the bounding box around the fracture.

Figure 9 :

 9 Figure 9: Concentration profile for the test case of continuous injection of tracer within a single fracture at t = 2×10 4(s): (a) v1 = 5 × 10 -4 m.s -1 , D1 = 1 × 10 -9 m 2 .s -1 (P e1 = 10 7 ), and (b) v2 = 5 × 10 -4 m.s -1 , D2 = 1 × 10 -4 m 2 .s -1

Figure 10 :

 10 Figure 10: Log2-plot of the discrete concentration L 2 errors for the case D2 = 1 × 10 -4 m 2 .s -1 showing an order of convergence O(h) for the concentration.

  Figure 11: 7-fracture DFN: mesh generated using FraC approach (a) and that used in Ahmed et al. (2015) (b).

Figure 12 :

 12 Figure 12: Concentration plots for the 7-fracture DFN.

Figure 13 :

 13 Figure 13: Mesh validation test of the benchmark 33-fracture DFN: equivalent permeabilities computed from the simulation results of single-phase flow at steady-state.

Figure 14 :

 14 Figure 14: 33-fracture DFN: pressure field and iso-value lines of pressure (a) and the norm of the velocity field (b) from the simulation on a mesh of h = hmax = 0.03 m.

Figure 15 :

 15 Figure 15: Snapshots of solute concentration distribution over time on the mesh DFN 1b of mesh size h = hmax = 0.03 m in Figure 6. Point A illustrates an observation point located at (2.35m, 1.35m, 0m).

Figure 16 :

 16 Figure 16: Concentration evolution at the observation point A.

Figure 17 :Figure 18 :

 1718 Figure 17: Solute concentration distribution over time on a mesh of L0 = 3 km.

Acknowledgements

The authors would like to thank the Computational Earth Science Group of LANL for providing us the source code of LaGriT and dfnWorks.

Counts of words

6560+78+569 (19/18/374/13) Total text+headers+captions

Counts of characters


Monica Rolando

Pedro Escoll

Carmen Buchrieser

Legionella pneumophila restrains autophagy by modulating the host's sphingolipid metabolism

Keywords: bacterial nutrition, Legionella pneumophila, lipid metabolism, xenophagy sphingosine-1 phospahte lyase

Sphingolipids are bioactive molecules playing a key role as membrane components, but they are also central regulators of many intracellular processes including macroautophagy/autophagy. In particular, sphingosine-1-phosphate (S1P) is a critical mediator that controls the balance between sphingolipid-induced autophagy and cell death. S1P levels are adjusted via S1P synthesis, dephosphorylation or degradation, catalyzed by SGPL1 (sphingosine-1-phosphate lyase 1). Intracellular pathogens are able to modulate many different host cell pathways to allow their replication. We have found that infection of eukaryotic cells with the human pathogen Legionella pneumophila triggers a change in the host cell sphingolipid metabolism and specifically affects the levels of sphingosine. Indeed,

L. pneumophila secretes a protein highly homologous to eukaryotic SGPL1 (named LpSPL).

We solved the crystal structure of LpSPL, showed that it encodes lyase activity, targets the host's sphingolipid metabolism and plays a role in starvation-induced autophagy during

L. pneumophila infection to promote intracellular survival.

Xenophagy is a cellular mechanism allowing an infected cell to rapidly degrade the invading bacteria by relying on the core autophagy machinery. The important role of autophagy in restricting the replication of pathogens has been demonstrated for several intracellular bacteria: the invaders are surrounded by autophagosome-like structures that allow fusion with lysosomes and consequently the degradation of the pathogen. Furthermore, it has been shown that the induction of autophagy suppresses the survival of intracellular bacteria. Thus the subversion of autophagy is a key strategy employed by pathogens to block the host response and to promote their survival.

Legionella pneumophila is one of those intracellular pathogens that interferes with the host autophagy machinery. L. pneumophila is a Gram-negative bacterium that is naturally found in aquatic environments where it replicates in protists, but that can also cause a severe pneumonia in humans called Legionnaires' disease. The high conservation of many signaling pathways in human macrophages and protists allows L. pneumophila to invade and to replicate in human cells. Once the host cell is infected, L. pneumophila is able to delay its delivery to lysosomes and to build up a specialized vacuole where it replicates until nutrients are depleted. The ability of L. pneumophila to subvert host defenses to set up its intracellular cycle relies on its uniqueness to encode over 300 secreted effector proteins that interfere with diverse cellular pathways. Our analyses of L. pneumophila genomes showed for the first time that many of these effectors share high similarity with eukaryotic proteins, and are thus proteins never or only rarely found in prokaryotic genomes, a finding that led to the hypothesis that these proteins had been acquired through horizontal gene transfer from its hosts.

One example is the lpp2128/spl gene that encodes a protein highly similar to eukaryotic SGPL1, named LpSPL for Legionella peumophila SGPL1. Phylogenetic analyses of the LpSPL sequence show that it branches within the eukaryotic clade of SGPL1, closest to the one of Acanthamoeba castellanii, further supporting the idea that the LpSPL-encoding gene had been acquired from a protist host. In order to characterize its function and role in host-pathogen interaction, we resolved its crystal structure. This showed that LpSPL exhibits a dimeric multidomain architecture like eukaryotic SGPL1 and revealed a high level of structural conservation in its active site composition with its eukaryotic counterparts. Eukaryotic SGPL1 catalyzes the irreversible cleavage of S1P, a product of the sphingolipid degradation pathway. By using a fluorogenic homolog of S1P we showed that LpSPL exhibits lyase activity in both transfected and infected cells. Furthermore, mutations of the residues reported to be involved in eukaryotic SGPL1 activity also abrogated the lyase activity of LpSPL.

Sphingolipids are a family of metabolites essential for membrane biogenesis, but they are also important signaling molecules modulating numerous physiological processes. The central component of the sphingolipid degradation pathway is ceramide that, hydrolyzed from sphingomyelin or synthetized de novo, undergoes successive conversions to various other sphingolipid intermediates like ceramide-1-phosphate, sphingosine and S1P. S1P levels are tightly fine-tuned through synthesis and degradation because once S1P is present in the intracellular milieu it functions as a cellular mediator (Figure 1A).

To understand the role LpSPL plays in the host cell, we measured host cell sphingolipid levels during infection by using mass spectrometry techniques. We observed that infection of macrophages with L. pneumophila leads to an overall reduction of bioactive sphingolipids such as sphingomyelin and ceramide, showing for the first time that L. pneumophila infection affects host cell sphingolipid metabolism. Interestingly, cells infected with a L. pneumophila strain that does not express LpSPL reveal an accumulation of sphingosine compared to cells infected with wild-type L. pneumophila. Thus, the bacterium specifically secretes LpSPL to prevent an increase of sphingosine levels in the infected host cell. This observation correlates with the fact that ectopically expressed LpSPL localizes at the endoplasmic reticulum, like eukaryotic SGPL1.

Given the role of sphingolipids as regulators of autophagy, we further explored whether secretion of LpSPL into the host cell affects the autophagic response to infection. Indeed, L. pneumophila is already known to interfere with the autophagy machinery: 2 secreted effectors have been previously identified to target the autophagy machinery (RavZ and LegA9).

By using multiple assays, including SQSTM1/p62 and LC3 western blots we found that LpSPL is a third effector secreted by L. pneumophila that modulates autophagy, as LpSPL restrains starvation-induced autophagy by acting on autophagosome biogenesis (Figure 1B). Importantly, this effect depends on its enzymatic activity as a catalytically inactive mutant fails to decrease host authophagy. By using high-content image-based analyses we quantified the formation of LC3 puncta in cells infected with wild type or a mutant strain deleted for the spl gene, showing that LpSPL limits the autophagic response during L. pneumophila infection to counteract antibacterial response by the host cell.

Taken together, our work discovered that the modulation of the sphingolipid metabolism of the host cell is a sophisticated strategy used by L. pneumophila to exploit the host autophagy machinery and to evade the host cell response.

Figure 1 .

 1 Figure 1. Legionella pneumophila modulates sphingolipid metabolism and autophagy during infection. (A) Simplified cartoon of sphingolipid metabolism: i) sphingomyelin from plasmamembranes is turned into ceramide by the enzyme sphingomyelinase (SMase). ii) Ceramide is then metabolized by ceramidase into sphingosine that iii) is phosphorylated by SPHK (sphingosine kinase) into sphingosine-1-phosphate (S1P), iv) a reaction that can be reverted by SGPP (sphingosine-1-phosphate phosphatase). v) S1P is secreted into the extracellular milieu, where it might act in a paracrine or autocrine way, or vi) may by broken down to hexadecenal and ethanolamine-1-phosphate (EA1P) by sphingosine-1-phosphate lyase (SGPL1). L. pneumophila decreases sphingolipid levels during infection (red arrows). It encodes 3 eukaryotic-like proteins

Acknowledgments

CB's lab was funded by the Institut Pasteur, the French Region Ile de France (DIM Malinf), ANR-10-LABX-62-IBEID and Infect-ERA project EUGENPATH (ANR-13-IFEC-0003-02); PE was funded by the Fondation pour la Recherche Médicale grant N° DEQ20120323697.


Luiz Antonio

Ricardo Paschoeto
email: ricardopas@bol.com.br

Santos

ICT and Financial Inclusion in the Brazilian Amazon

Keywords: Financial Inclusion, Development, ICT for Development, Agência Barco, Marajó Island

The challenge of providing the infrastructure of public services in the less developed regions of Brazil has mobilized the Brazilian government in the quest for new and creative approaches that can reduce the major inter-regional disparities in the country. One of the initiatives implemented include access to the financial system, since, by way of example, such access is almost non-existent on Marajó Island in the state of Pará in the Brazilian Amazon. To change this reality, an innovative e-government project is the itinerant bank branch installed in a boat, named Agência Barco, to serve the riverine populations of regions with low population density, transportation difficulties and limitations in access to information and communication technology (ICT). Thus, the main objective of this research is to identify how the financial inclusion indicators have been influenced by the work of Agência Barco on Marajó Island from the ICT standpoint. The results obtained led to the conclusion that Agência Barco has been able to attend the needs of access to financial products and services demanded by the population of Marajó Island, as well as identify opportunities for broadening financial education and inclusion through this e-government venture.

Introduction

Studies of the World Bank [START_REF]World Development Indicators[END_REF] and the Brazilian Institute of Geography and Statistics [START_REF]Informações Sociais, Demográficas e Econômicas[END_REF], among others, point to a continuous improvement of social and economic development indicators in Brazil in the last two decades, with a reduction of inequalities between income classes. However, when analyzing the Brazilian Municipal Human Development Index, regional inequalities are perceived, and it is possible to identify a broad variance in the opportunities available to Brazilians [START_REF]Índice de desenvolvimento humano municipal brasileiro[END_REF]. These inequalities also prevail in access to financial services, which can be seen in the indicators of the Financial Inclusion Reports of 2010 and 2015, which reveal a wide dispersion between the Units of the Federation [START_REF] Bcb | Relatório de Inclusão Financeira N. 2. Brasília[END_REF][START_REF] Bcb | Atendimento bancário no país -Distribuição do quantitativo de municípios por Região e UF[END_REF], as in March 2015 when 240 municipalities were without banking services of any kind (branches, service outlets or ATMs). Furthermore, in March 2015 a total of 1922 Brazilian municipalities did not have bank branches, i.e. 34.5% of Brazilian municipalities [START_REF] Bcb | Atendimento bancário no país -Distribuição do quantitativo de municípios por Região e UF[END_REF]. Therefore, there is still vast asymmetry in Brazil in access to and use of banking services [START_REF] Crocco | Exclusão financeira no Brasil : uma análise regional exploratória[END_REF].

Thus, this article seeks to investigate an ICT-equipped itinerant bank in vessels -an egovernment project developed by Caixa Econômica Federal (CAIXA 1 ) and named Agência Barco2 -which seeks to attend the populations of riverine regions with access to banking services, in cities where there is not a single local bank branch. In specific terms, the universe of this work includes the Agência Barco that serves Marajó Island 3 in the state of Pará -a region with a low human development index (HDI), including the city of Melgaço, which has the lowest HDI in Brazil [START_REF]Índice de desenvolvimento humano municipal brasileiro[END_REF] -aiming to identify how the Agência Barco, by means of ICT, influences the financial inclusion of the Marajó Island region.

2

Theoretical Background

ICT and Financial Inclusion

For the Brazilian Central Bank, financial inclusion is the process of effective access and use by the population of financial services suited to their needs, thereby contributing to their quality of life. However, not all individuals and companies have access to the financial system, either due to the lack of availability of services and products for given sectors of society or the lack of a financial culture among the people. Despite progress in the relationship between citizens and the financial system and the increased presence of financial institutions in almost all Brazilian municipalities due to the success of the Banking Correspondent (BC) model [START_REF] Jayo | Um mapeamento descritivo dos modelos de gestão de redes de correspondentes bancários no Brasil[END_REF][START_REF] Leonardi | Multiplex Appropriation in Complex Systems Implementation: The Case of Brazil's Correspondent Banking System[END_REF], the country still has less than half of the number of bank branches per capita existing in developed countries [START_REF] Bader | Logística da distribuição bancária: tendências, oportunidades e fatores para inclusão financeira[END_REF].

Moreover, by adopting a holistic approach, studies on ICT4D -Information and Communication Technology for Development -seek to analyze how ICT can promote development while respecting the complexity of local, national and international conditions [START_REF] Heeks | Do Information and Comunication Technologies (ICTs) Contribute to Development?[END_REF]. Thus, ICT is the latest enabler of financial inclusion, since its mass diffusion is the most significant technological change in low-income communities in recent years, leading to the emergence of inclusive financial services, especially those related to mobile communication technology [START_REF] Heeks | Do Information and Comunication Technologies (ICTs) Contribute to Development?[END_REF].

The Dynamic Info-Inclusion Model (2iD)

In an attempt to overcome the lack of research in the ICT4D field in Brazil, Joia [START_REF] Joia | Bridging the digital divide: some initiatives in Brazil[END_REF][START_REF] Joia | Inclusão digital no Brasil: um modelo heurístico de natureza dinâmica[END_REF] proposed the dynamic info-inclusion model (2iD), which evaluates digital inclusion encompassing both political, technical, educational and social aspects, as well as the dynamics of a virtuous cycle of participation and empowerment, as shown in Figure 1.

In this context, Joia [START_REF] Joia | Bridging the digital divide: some initiatives in Brazil[END_REF][START_REF] Joia | Inclusão digital no Brasil: um modelo heurístico de natureza dinâmica[END_REF] suggests that the sustainability of digital inclusion should include not only financial and economic factors, but also factors that reflect the government's concern for the continuity of public policies for digital inclusion. Furthermore, according to the author, the "education" component in the model should go beyond the mere training of individuals, incorporating awareness of the opportunities generated by ICTs for socio-economic change. Moreover, the author stresses the need to consider the environment and the context, in order to create specific content that meets the expectations and needs of the location of the individuals in which it develops the infoinclusion project.

On the other hand, in the dynamic process of implementation of actions, there occurs increased awareness of individuals, particularly the awakening of interest in issues relating to the use of ICT. From that moment onwards, these individuals begin to demand content, education, services and access to ICT. The feedback, coupled with the broadening of the empowerment cycle arising from this dynamic, generates the implementation of new initiatives for digital inclusion, as seen in Figure 1.

The 2iD Model Adapted to Financial Inclusion (2iDf)

For the development of this research, the 2iD dynamic info-inclusion model created by Joia [START_REF] Joia | Bridging the digital divide: some initiatives in Brazil[END_REF][START_REF] Joia | Inclusão digital no Brasil: um modelo heurístico de natureza dinâmica[END_REF] was adapted to evaluate financial inclusion, as shown in Figure 2. For this purpose, the theoretical framework that supports the components of the 2iD model adapted to financial inclusion, called 2iDf, is presented in Figure 3.

Thus, in the adapted model of Figure 2, the following aspects of financial inclusion are addressed:  Infrastructure and access -This deals with the creation of individual and collective conditions for the population of the locations to access financial services involving:

Infoinclusion Local Content Education

Infrastructure and Access

Economic, Political, Ethical, & Legal Sustainability

Increase in awareness

Feedback and growth

Increase in demand

Implementation

o The presence of banking facilities and technological resources; o Availability of services and products; o Accessible cost for access to the Agência Barco and the products offered.



Financial Education -This deals with the training of people to use the products offered in all its possibilities, comprising: o Presentation of the characteristics of products and services that enhance the financial resources of the citizen; o Presentation of the characteristics of products and services to support entrepreneurial activities that can generate more income in the locations; o Promotion of the use of electronic transactions, thereby replacing cash.



Services and Products -This involves offering of a portfolio of financial products and services that takes into consideration the reality of the financial inclusion in each location served.  Sustainability -This involves the maintenance and updating of financial products and services offered in the locations, including economic aspects -logistic costs, personnel and ICT -as well as aspects related to public policies and legal and ethical issues. Just as with the 2iD model of Joia [START_REF] Joia | Bridging the digital divide: some initiatives in Brazil[END_REF], the dynamic process of the 2iDf model perceives the government as being responsible for implementing the financial inclusion initiatives. In this dynamic process, through the implementation of financial inclusion actions, there is an increase in the awareness of individuals, particularly the awakening of interest in questions relating to the use of financial products and services. Thereafter, these individuals begin to demand more financial services and products, more educational activities and, consequently, more access to technologies that enable the use of financial

Infoinclusion

Services and Products

Financial

Education

Infrastructure and Access

Economic, Political, Ethical,& Legal Sustainability

Increase in awareness

Feedback and growth

Increase in demand

Implementation

services. The feedback coupled with the broadening of the empowerment cycle resulting from this dynamic leads to the implementation of new financial inclusion initiatives.

Component of The

Methodological Procedures

This article used the case study method [27], with data collection by means of interviews and direct observation. For the processing and analysis of data, the content analysis technique [START_REF] Bardin | Análise de Conteúdo[END_REF] and the application of the dynamic info-inclusion model adapted to financial inclusion (2iDf) were used. Thus, a single case study is investigated with one unit of analysis [27] represented by a service outlet of the Marajó Island Agência Barco of CAIXA, with data obtained from document research, a questionnaire and interviews. A directed sample is used, in which individuals are selected on the basis of certain characteristics regarded as relevant by the researchers and participants [START_REF] Gil | Como Elaborar Projetos de Pesquisa[END_REF], and employing techniques of content analysis for data analysis [START_REF] Silva | Análise de Conteúdo: Exemplo de Aplicação da Técnica para Análise de Dados Qualitativos[END_REF].

Data Collection

To achieve the objectives proposed in this case study a literature review was initially conducted on the subject in question. Documents were also gathered from the financial institution under analysis, IBGE, UNDP and the Central Bank.

Open interviews were conducted [31] with employees of the financial institution to assist in providing an in-depth description of the case in order to identify relevant information not available in documentary sources. However, in order to grasp the perception of the Agência Barco users regarding the services and aspects of financial inclusion, semi-structured interviews were staged [32]. The interviews were conducted during the visit of one of the researchers to the Agência Barco in Marajó. A total of 23 clients and five employees from CAIXA and five servants of organizations involved participated in these interviews, which were recorded and later transcribed verbatim for analysis. Moreover, assystematic non-participant direct observation was also conducted [31].

All the above data were obtained over the course of five days in July 2015 during a trip made by one of the researchers in the Agência Barco in the Marajó Island in Brazilian Amazon.

Data Analysis

Interviews conducted with Agência Barco clients and employees were handled with content analysis techniques, with a priori categorization based on constructs of the 2iDf model, with alphanumeric coding and grouping by frequency of occurrence, i.e. by repeating of contents common to the majority of respondents [START_REF] Bardin | Análise de Conteúdo[END_REF]. The categorization followed the mixed model [START_REF] Laville | A construção do saber: Manual de metodologia da pesquisa em ciências humanas[END_REF], that is, it made it possible to add new categories as registration units were regrouped.

To assist in the analysis of content, lexical analysis (which applies statistical methods to the description of the vocabulary) was applied before the analysis of content. This was done to ensure that the data analysis therefrom was fully implemented, encompassing several possibilities that might arise or emerge [START_REF] Freitas | Análise Léxica e Análise de Conteúdo: Técnicas complementares, sequenciais e recorrentes para exploração de dados qualitativos[END_REF].

Registration units for clipping of excerpts were defined by words and expressions that alluded to the static and dynamic components of the 2iDf model. Lexical analysis led to the initial identification of 864 words and phrases which, in turn, were grouped into 235 initial categories. The recurring process of lexical analysis and analysis of content led to a new categorization phase, totaling 97 categories with 821 occurrences (see Figure 4). From this phase of the review process onwards, the intermediary categories were also coded according to their influence -positive or negative -in relation to the components of the 2iDf model. The data were interpreted by means of comparison with the constructs highlighted in the theoretical framework, which supported the components of the 2iDf model, seeking to identify the impact of the Agência Barco on financial inclusion of locations served, and to identify which of the components of the 2iDf model generated opportunities for actions that could contribute to regional development.

The "Agência Barco" Case

CAIXA is a Brazilian public bank having experience in the operation of social programs of the Federal Government and of banking inclusion to expand attendance to populations still without access to banking services, especially through Banking Correspondents. However, in many areas with difficult access by land, or very far from municipalities with dynamic economic activity, there are difficulties for the business model of the Banking Correspondents to attend the needs of the populations. This situation is especially relevant in the riverine locations of the states of the Amazon region, where river trips between towns and larger cities can take more than one day.

In this scenario in which barriers to the physical presence of banking institutions are well-nigh insuperable, the 'Itinerant Riverine CAIXA Service Units' project, which became known as Agência Barco, was conceived.

The Agência Barco consists of boats designed and built exclusively for operation of a bank agency where bank staff perform their activities like any other bank branch of the institution. In addition, the boats are equipped with infrastructure for support to partnerships with other agencies of the Federal and State Government. These include, among others, the Ministry of Health, Ministry of Culture, Ministry of Labor, Ministry of Education, Secretariat on Policies for Women, as well as Courts of Justice to carry out public and institutional policies in the locations. The Marajó Island Agência Barco was inaugurated in January 2014 to serve ten cities: Soure, Salvaterra, Ponta de Pedras, Muana, São Sebastião da Boa Vista, Curralinho, Bagre, Breves, Melgaço and Portel. It is a vessel with three decks, with a total area of one thousand and seventy-seven square meters, with capacity for seventy-six people seated in the service area of one hundred and forty square meters, plus a further twenty passengers who remain on the boat throughout the trip, namely five employees of CAIXA, five employees of the organizations involved, four security guards and six crew members.

One of the prerequisites for opening the Agência Barco imposed by the regulatory bodies was that no cash would be stored or transacted on the boat for security reasons. Thus, the business strategy has the support of a lottery office in each city served by the boat, such that the amounts in cash can be handled by these service outlets, also operated locally by CAIXA.

The Agência Barco makes a monthly trip, referred to as a cycle, remaining two days on average in each location, working during normal banking hours and offering all the services of a normal branch, except for cash transactions. The displacement between cities is usually at the end of the day in the early evening and night, depending on the sailing conditions and the weather (wind, rain, etc.). The main services offered are: opening accounts, microcredit operations, financing of building materials, security bonds, life insurance, direct consumer credit, registration and resetting passwords on federal government citizen cards, registration and regularization of the social integration program, release of length of service pension funds, release of unemployment benefits, registration and regularization of social security cards, family allowance benefits, among others.

The technological solution for the Agência Barco includes satellite communication to connect to the CAIXA Datacenter located in Brasilia, which is structured to support real time data and voice applications with autopointing and autotracking mechanisms. These functions enable execution and transmission of back-end processes of the Agência Barco, even with the vessel in transit, and reduce the occurrence of communication failures due to movement of the boat on its moorings. The satellite connection is concentrated in the teleport of the telecommunications services operator and forwarded by terrestrial circuits to the CAIXA Datacenter located in Brasilia. In addition to this, the vessel is equipped with direct access to the Internet via cell phone networks.

Results

Of the ten municipalities served by the Agência Barco on Marajó Island, eight are classified as having low or very low human development and only two are classified as having medium human development, which is well below the Brazilian global HDI [START_REF]Índice de desenvolvimento humano municipal brasileiro[END_REF]. Moreover, only four of the ten municipalities served by the Agência Barco have bank branches, and only the municipality of Breves has a CAIXA branch, which is the financial institution responsible for putting the public policies of the federal government into operation.

The Financial Inclusion Index of the State of Pará was analyzed in the second Financial Inclusion Report of the Brazilian Central Bank [START_REF] Bcb | Relatório de Inclusão Financeira N. 2. Brasília[END_REF], with Marajó Island featuring the lowest indicators of the state, which confirms the lack of the availability of banking services in the Agência Barco operating region.

It was revealed that the majority of the customers interviewed (47.8%) had little schooling (incomplete primary education) and a little over half of the sample received social benefits from the federal government -with only one respondent stating that this is the main source of family income.

Of the twenty-three customers interviewed, approximately 50% knew of the existence of the Agência Barco through third parties, i.e. a neighbor or relative saw the boat in the harbor and passed on the information, and six customers interviewed came to the branch by boat, coming from tributaries or creeks in the regions surrounding the cities.

In addition to this, there has been a marked increase in cell phone penetration, with all of the customers confirming that they owned cell phones, although around half of them did not have any form of Internet access.

Of the twenty-three respondents, 40.9% reported having no formal relationship with financial institutions through checking or savings accounts, though only three said they had never had a bank account. In addition, more than half of the respondents reported their preference for full withdrawal of money deposited in a bank account.

With respect to financial education, there was little experience or awareness regarding other products such as, for example, loans and investments among those interviewed.

The content analysis conducted resulted in tables that summarize the opinions of respondents regarding the presence of each of the components of the 2iDf model. They identified the intermediate categories and their frequency of occurrence, as well as the influence of each intermediate category in the final category, namely a positive influence or presence of aspects related to a given component of the model or a negative influence or no aspects related to a given component of the model.

The perception of the existence of the static elements of the 2iDf model proves to be more accentuated than the dynamic components, as evidenced by a higher frequency of comments -attaining more than 93% of all remarks made by respondents. Thus, the low frequency of remarks related to the dynamic aspects of the model suggests that the financial education process failed to generate awareness of the potential opportunities offered by financial inclusion in the daily lives of individuals [START_REF] Joia | Inclusão digital no Brasil: um modelo heurístico de natureza dinâmica[END_REF].

To assess the relative impact of the presence of each static component of the 2iDf model, a scale was created with three levels to represent the perception of the presence of a given component of the model and its positive aspects -Figure 5 and Table 1.

Perception of presence of the component Caption

High> 75%

25% <= Average <= 75%

Low < 25%

Based on the above scale, a graphical representation of the results analyzed is shown in Figure 6, associated with the perceptions of the presence of positive aspects of the static components of the 2iDf model.

Conclusions

This research identified the need for investment in ICT infrastructure to meet the commercial agreements in order to provide a wider range of services, as well as enable partnerships with other state and federal agencies for services related to issuance of documents and processing retirement benefits, for example. In addition, the need to create incentives for dissemination of ICT infrastructure in regions with lower population density was identified, such as for example, the creation of basic cell phone packages with access to transaction services (cell phone banking, credit and debit transactions, cell phone payment).

It is also clear that without minimum financial education of the customers, it will be difficult for Agência Barco to offer more than the most basic and simple financial services to the population served.

Lastly, it was perceived that the riverine population in Marajó Island is still not even aware of the potential benefits accrued from being financially included as the dynamic components of the 2iDf model were not considered important by the local population.

In sum, financial inclusion is still a challenge in Brazil and tackling this situation depends on technological innovations, business models and public policies that can provide faster inclusion of the population currently excluded from financial system in the country.

Financial Education

Infrastructure and Access

Increase in awareness

Feedback and growth

Increase in demand

Figure 1 -

 1 Figure1-Dynamic Info-inclusion Model (2iD). Source: Adapted from Joia[11, pp 308]

Figure 2 .

 2 Figure 2. 2iD Model adapted for financial inclusion (2iDf)

Figure 3 .

 3 Figure 3. Theoretical framework for the 2iD model adapted for financial inclusion (2iDf)

Figure 4 .

 4 Figure 4. Categories of the content analysis

Figure 5 .

 5 Figure 5. Classification of the components of the 2iDf model

Figure 6 .

 6 Figure 6. Final representation of the 2iDf model

Table 1 .

 1 Scale of the presence of the static components of the 2iDf model

CAIXA is a Brazilian federal bank in charge of implementing social programs in the country.

Agência Barco means Riverboat Branch in English.

More information about Marajó Island is available at https://en.wikipedia.org/wiki/Maraj%C3%B3, retrieved on January 19, 2017.