Sylvie Vandoolaeghe

Alessandra Blaizot

Danie Boudiguet

Valérie Bougault

Eduardo Dei Cas

Benoît Foligne

Anne Goffard
email: anne.goffard@univ-lille2.fr

Hélène Lefranc

Bénédicte Oxombre

Thomas Trentesaux

Bernard Vandenbunder

Isabelle Wolowczuk

Laurence Delhaes
email: laurence.delhaes@pasteur-lille.frˆdeceased

A charter for biomedical research ethics in a progressive, caring society

Keywords: Charter, Ethics, Research, Life science, Integrity

Background: Given that advances in research continuously raise new ethical issues, a multidisciplinary working group of investigators involved in biomedical research has gathered to discuss and compare ethical viewpoints in their daily practice. Methods: The working group has drafted a Charter for Ethics in Biomedical Research that encompasses all the steps in the research process, i.e. from the initial idea to analysis and publication of the results.

Results: Based on key principles for ethically responsible research, the Charter may serve as a tool for performing research, discussing research issues and training researchers.

Introduction

Over the past few decades, outstanding progress in the field of biotechnology has increased awareness of ethical issues. Researchers, national bodies and international organizations have therefore been prompted to develop guidelines on research ethics and integrity [START_REF] Fluss | The evolution of research ethics: the current international configuration[END_REF].

A well-known chapter in the history of research with human subjects opened after World War II, with the establishment of the Nuremberg Code [2]. This was followed by the Declaration of Helsinki, which constitutes the basis for all subsequent national and international ethical guidelines [START_REF] Medical | World medical association declaration of helsinki: ethical principles for medical research involving human subjects[END_REF]. However, the ethical concepts in these important international guidelines do not cover aspects of daily research practice. In the field of biotechnology, a range of questions concern aspects of a researcher's day-to-day work: training, practice, societal and professional roles, objectivity, credibility, evalu-ation and constraints. In 1983, France was the first country to create a National Consultative Committee on Ethics in the Health and Life Sciences (Comité Consultatif National d'Ethique pour les sciences de la vie et de la santé, CCNE) to "provide advice on ethical issues and social issues raised by the progress of knowledge in the fields of biology, medicine and health". Soon after the foundation of the CCNE, France introduced legislation on the creation of ethics forums and "Think Tanks" in French hospitals. Hence, Lille University Hospital and the Lille Faculty of Medicine launched the Espace Ethique Hospitalier et Universitaire de Lille (EEHU; http://www.eehu-lil le.fr/) in 2010. Since 2012, the "Ethics and Research" working group has met within the framework of the EEHU (under the impetus of Dr Eduardo Dei Cas [START_REF] Boury | Biomedical research: the debate on the reduction and emergence concepts[END_REF][START_REF] Boury | Current bioethical issues in parasitology[END_REF][START_REF] Aliouat-Denis | In memoriam: Eduardo Dei-Cas (1945-2014)[END_REF]) to consider the role of life and health science researchers in a university hospital in the 21st century (from fundamental research to clinical applications). After a number of multidisciplinary debates on issues generated by the research, progress in research and the representations of research that ensue, the group came to feel that the fruit of its discussions could be usefully summed up as a guide to ethical issues for researchers.

Materials and methods

Constitution of an "Ethics and Research" working group

In 2011, individual and institutional stakeholders in the life and health sciences based in and around the city of Lille (northern France) were invited to form a working group. Notably, people with expertise in medicine, pharmacy, dentistry, nursing, sports science, biology, chemistry, physics, literature, philosophy and languages were invited to participate. Ultimately, a group of 14 people from a wide range of areas of biomedical research met to discuss research ethics (Table 1). The group meets every two months to discuss issues and share its members' experiences. It combines critical thinking and ethical analysis within a framework of discourse ethics [START_REF] Habermas | Justification and Application: Remarks on Discourse Ethics[END_REF]. The group's mission goes beyond normative ethics, although it does question the sources, grounds, and legitimacy of ethical matters.

Methods

The

Results

The charter for biomedical research ethics

The Charter for Biomedical Research Ethics is based on debate, open-mindedness and dialogue [START_REF] Boury | Biomedical research: the debate on the reduction and emergence concepts[END_REF][START_REF] Boury | Current bioethical issues in parasitology[END_REF][START_REF] Aliouat-Denis | In memoriam: Eduardo Dei-Cas (1945-2014)[END_REF][8][9][START_REF][END_REF][11][START_REF] Elliott | Science, policy, and the transparency of values[END_REF]. Rather than imposing additional constraints on researchers, the Charter's objective is to put forward common points The Charter describes how these principles can be translated into actions during the different steps in the research process, from the initial working hypothesis to publication and dissemination of the results and consideration of opportunities for further research. The Charter addresses the following five aspects of research:

Defining the context and objectives of research

As an essential prerequisite to a research project, the researcher must determine the "state of the art". This enables her/him to (i) present and acknowledge the field's results and current limitations (integrity and honesty), (ii) cite the corresponding authors, and (iii) think of ways of advancing the state of the art (critical thinking and open-mindedness).

The researcher should precisely define her/his working hypothesis by delimiting the field of investigation (rigorousness). She/he should take account of what is known from previous research, while seeking to add value. Accordingly, she/he will justify her/his efforts to make further progress (creativity and altruism).

Selecting materials and setting up methods

Methodology (integrity, rigorousness and verifiability)

The researcher should select the methodology that is best suited to answering the research question. Methods should not be chosen or adapted in order to obtain predefined results. The researcher must check that the chosen materials and methods are appropriate and can be used optimally under the research project's various constraints. She/he must seek to reduce experimental bias as much as possible, while considering the limitations of her/his approach.

From the beginning of the work onwards, all raw data should be stored so as to ensure rigorousness and transparency in accordance with current legislation. This requirement covers the research hypothesis, communication with various organizations or journals, bibliographic references, the methods used, the data collected, and the structuring, analysis, publication and communication of results. Compliance with these criteria helps to structure the research and enables research processes to be verified and replicated.

Commitments, interests and working relationships

The researcher must not overrate her/his work or mislead funding bodies, research institutions, colleagues, the media or the general public (integrity, honesty, and accuracy).

The researcher must state the identity of her/his sponsors and funding sources and must describe the nature of relationship(s) she/he has with them (e.g. acting as a consultant or an expert advisor). She/he must declare any personal, professional or financial interests in order to maintain the necessary level of trust in her/his own research and (as an extended responsibility) in research in general. In the event of a potential conflict of interest, the researcher can only agree to participate if her/his impartiality is beyond doubt.

When a researcher is asked to engage in a research activity, she/he must always honestly estimate her/his ability to successfully complete that activity. She/he should commit to maintaining and upgrading her/his skill level and knowledge and to achieving proficiency in the techniques required to advance her/his research (competence).

The researcher must carry out her/his research in accordance with the current legislation and deontological and institutional guidelines, rules and regulations (legality).

She/he must respect the colleagues and students which whom she/he works by treating them fairly and without subjecting them to discrimination or any kind of abuse (respect, courtesy and loyalty in working relationships within the group, between groups, towards her/ his colleagues and hierarchy). She/he should seek to manage tensions within a group as best she/he can and must acknowledge the constraints related to employment contracts.

Participants (caution, carefulness, respect, confidentiality and reliability)

The researcher and her/his institution are ethically obliged to take account of the risk-benefit ratio associated with research work. In clinical research, the "subjects" are human beings and so their dignity, freedom and privacy must be protected. The researcher must minimize the risk of harm to others, and should act with caution and foresight. She/he should take specific precautions to protect vulnerable persons. When the research is performed in a country other than her/his own, she/he should respect the local culture and environment. She/he should admit to and correct (to the greatest extent possible) any mistakes, omissions or damage caused. Moreover, the researcher must be cautious and take precautions if her/his experiments are likely to have an environmental impact.

For biomedical research involving human subjects, the sponsor must obtain approval from the relevant independent ethics committee/investigational review board.

If animals are used, the researcher must treat them carefully and respectfully. She/he must limit the animals' suffering, in line with the "3Rs" (Reduction, Replacement, and Refinement). She/he must not implement an unnecessary or poorly designed animal experiment or any animal experiment that has not been approved by the relevant institutional animal care and use committee.

Analysing results

When analysing research results, the researcher must not disclose confidential personal data (respect, carefulness and confidentiality). The researcher should not be let her/his preconceived ideas influence the analysis or interpretation of her/his results. She/he should differentiate between her/his own opinions and the scientific experiment (impartiality). She/he must also distinguish between hypotheses and speculations on one hand and the results and conclusions of her/his research on the other (accuracy). The researcher should avoid bias when analysing and interpreting her/his own data. According, peer reviews or reviews of a colleague's work should be objective (integrity and objectivity).

The researcher must seek to avoid negligence (rigorousness). She/he should maintain a critical view of her/his own work and that of her/his peers -particularly during collaborative work (critical thinking). In interdisciplinary projects, she/he should place the requisite trust in colleagues and pay attention to the different views that they may bring to her/his area of expertise (open-mindedness).

A researcher's analysis of the results and her/his subsequent establishment of the conclusions must be driven by a requirement for truth and knowledge. The research must not be transferred too quickly to the media or business stakeholders outside the scientific community. These results must not be falsified by constraints, obligations or pressure. In contrast, results should not be set aside for too long in an effort to "protect" them.

When the researcher evaluates the research of her/his peers, she/he should carry out the mission thoroughly, fairly and without disclosing confidential information. She/he shall review papers and grant applications in a timely manner and from a neutral viewpoint, and must be guided solely by scientific considerations. Reviewing is an important part of her/his overall responsibility to the discipline of research. When supervising students, she/he shall offer appropriate advice and promote competence and responsibility (responsible tutoring).

Communicating and publishing research results

The researcher must comply with patent rights, copyright and other intellectual property rights. She/he acknowledges other people's work and contributions. She/he must not present other people's work as her/his own. She/he must respect her/his colleagues' ideas and never plagiarize them [8] (honesty and accuracy).

The researcher accepts responsibility for her/his contributions as a group member, co-author or author. When publishing, the author should honestly distinguish between (i) co-authors who meet the criteria for authorship and (ii) colleagues who should solely be acknowledged (i.e. minor contributors).

The principal investigator should appropriately manage her/his group members' respective contributions by precisely defining their respective roles (rigorousness).

The researcher should communicate her/his results in a timely and open manner, in order to establish ownership of intellectual property and establish priority. Raw data must be stored, so that they can be validated and/or replicated (transparency, rigorousness, precision and verifiability). The researcher takes account of all her/his results -including unwanted or unexpected onesand should not filter them arbitrarily (accuracy). She/he should not disclose confidential information on patients, participants or other people involved in research before, during or after publication of the results (confidentiality). However, on request, the results should be made available to the participants.

Defining follow-up actions and perspectives for further research

By performing high-quality research with integrity, the researcher creates and maintains confidence in her/his work. In this respect, she/he must be alert to any lack of integrity within her/his research group (responsibility).

In order to enable and promote responsible behaviour by the researcher, her/his research environment should:

promote competence via access to continuing education and training. set clear, realistic working rules that are compatible with the research project's time-frame and constraints. implement appropriate actions to correct or retract publications, if required. define fair criteria for career advancement and thus enable both individual and collective development.

Discussion

Scientific research is rooted in the quest for truths and knowledge that serve humanity; it fits into a vision of a progressive, caring society. Researchers must have freedom of speech and intellectual independence with regard to their science, from hypotheses through to implementation and conclusions. In this respect, funding agreements are drawn up to precisely define the rights and obligations of each party (including those of the research institution) and to manage intellectual property, communication and exploitation of the results. However, the research process must always leave some room for serendipity.

Given that (i) today's societies have many points of reference, (ii) scientific and biomedical practices are constantly evolving and (iii) the use of scientific and biomedical practices may create dilemmas and raise obstacles, we need to find a way of thinking about research ethics that enables us to question the legitimacy of the means employed and the research's possible consequences. This approach is necessarily placed in a specific context -here and now, with circumstances, constraints, a past and a future. The present Charter for Biomedical Research Ethics is not a tool for accrediting a given research project as being "ethical" or not; in fact, it incites each researcher to think about the objectives she/he defines, the means she/he uses to achieve them, and the risks and limitations related to her/his research. Accordingly, the researcher takes responsibility not only for her/his own research but also for research in general.

Research in the life and health sciences generates representations of the living world [START_REF] Elliott | Science, policy, and the transparency of values[END_REF] via (a) frameworks for conceptualizing research hypotheses, (b) biomedical practices (in which a patient is considered as an object, e.g. in biotechnology), (c) policies on precaution and prevention, (d) education and training [START_REF] Schopper | Innovations in research ethics governance in humanitarian settings[END_REF][START_REF] Carrese | The Essential Role of Medical Ethics Education in Achieving Professionalism: The Romanell Report[END_REF], (e) a rationale for political decisions on biological research, (f ) perceptions and images of the living world conveyed by the media and in the cultural and artistic domains.

Biomedical research deals with emergence of the complexity of the living world in its universality and singularity, which requires open-mindedness, critical thinking, multiple analytical strategies, creativity and unceasing conscientiousness. In view of the complexity and social impacts of research and its mutually stimulating disciplines (biomedical science, the human and social sciences, etc.), it is essential to question research's resources, limitations and meaning. What are we aiming to achieve? What can be expected? What should be avoided? How can science and society move forward? As its title suggests, the present Charter for Ethics in Biomedical Research is primarily aimed at researchers in the fields of life and health sciences worldwide. However, it could be useful for other stakeholders in the paramedical professions and/or other fields of science. Indeed, when disciplines meet or are confronted with one another, a new vision of the living worldperhaps better adjusted to the complexity of lifemay emerge. Along these lines, today's research programmes are growing at the interfaces between nanoscience, computer science and biology. Some researchers recommend a global, integrated approach for strengthening the links between human health, animal health and environmental care (the "One Health" concept, for example [START_REF] Kahn | In memoriam: Great 21st century physician "One Health" leader dies -Ronald M. Davis, MD, Past President[END_REF]). Another example is the interdisciplinary workgroups (including representatives of the human and social sciences) organized by the North American "Global Health" consortium [START_REF]Health Engagement: Sharpening a Key Tool for the Department of Defense[END_REF]. The European Union's research directives for the "Horizon 2020" program [START_REF][END_REF] also encourage various types of multidisciplinary and/ or international cooperation; these correspond to an innovative, globally coherent vision of issues that require international guidelines for research integrity [START_REF] Resnik | International standards for research integrity: an idea whose time has come?[END_REF].

In conclusion, we believe that present Charter for Ethics in Biomedical Research (based on common principles) will evolve as a function of the context and scientific progress. Ethical issues must be continually questioned because they reflect the infinite complexity of the living world. We consider that the living world is focused on organisms and their futureorganisms that interact with each other and with their environment, compatible with an ongoing world, simultaneously complex, historical and open. Nevertheless, certain fundamental principles will always be applied, regardless of contextual and conceptual changes.

for

  Infection and Immunity of Lille (CIIL), Institut Pasteur de Lille (IPL), F-59019 Lille, France. 10 UMR CNRS 8402, F-59021 Lille, France. 11 CIIL, IPL, Lactic Acid Bacteria and Mucosal Immunity Research Group, Center for Infection and Immunity of Lille (CIIL), Institut Pasteur de Lille (IPL), F-59019 Lille, France. 12 University Lille, CNRS, Inserm, CHU Lille, Institut Pasteur de Lille, U1019 -UMR 8204 -CIIL -Centre d'Infection et d'Immunité de Lille, F-59000 Lille, France. 13 Délégation à la Recherche Clinique et à l'Innovation, Maison Régionale de la Recherche Clinique Hospitalière et Universitaire, Lille University Hospital, F-59037 Lille, France. 14 LIRIC -UMR 995, F-59000 Lille, France. 15 Interdisciplinary Research Institute, USR 3078 CNRS, University of Lille 1, F-59658 Villeneuve d'Ascq, France. 16 Lung Infection and Innate Immunity (LI3), BP 245, F-59019 Lille, France. Received: 28 April 2015 Accepted: 5 November 2015

  The present Charter for Biomedical Research Ethics in a Progressive, Caring Society is intended to be simple, open to challenge and revision, and accessible to all researchers and stakeholders involved in health and life science research. The Charter provides a number of practical tools for researchers:

	-A tool for improving research and ethical practices:
	Enhancing the quality and legitimacy of research,
	emphasizing the researcher's individual responsibility
	and promoting personal consideration and collective
	debate on research ethics.
	-A tool for debate and dialogue within research
	groups: Facilitating interactions within the research
	group, between research groups, and between
	research groups and institutions.
	-A training tool for students and young researchers:
	Facilitating learning and counselling, encouraging
	the proper supervision of young researchers,
	fostering their awareness of ethics and introducing
	them to research ethics.

Table 1

 1 The members of the multidisciplinary "Ethics and Research" working group principles for action for those working in the fields of life and health science.The Charter is based on the following principles:

	First name and family name	Qualifications/position/field of research	Affiliation
	Alessandra Blaizot	DDS, MSc/Ethics	CHRU a of Lille, University Hospital, University
			of Lille
	Danie Boudiguet	Director of Nursing, Nurse anaesthetist (I.A.D.E b	CHRU of Lille, University Hospital
		school) RNA/Teaching and training
	Valérie Bougault	PhD/Physiology & Sport sciences	University of Lille
	Eduardo Dei Cas	MD, PhD/Mycology & Parasitology	CHRU of Lille, University Hospital, University
			of Lille
	Laurence Delhaes	MD, PhD/Mycology & Parasitology	CHRU of Lille, University Hospital, University
			of Lille
	Benoît Foligné	PhD/Microbiology & Immunology	Institut Pasteur de Lille, University of Lille
	Anne Goffard	MD, PhD/Virology	CHRU of Lille, University of Lille
	Hélène Lefranc	BSc, Hospital Engineer at DRCI c, , Clinical Research	CHRU of Lille University Hospital
		Project Manager/Research Management
	Bénédicte Oxombre-Vanteghem	PhD/Neuro-immunology	University of Lille
	Thomas Trentesaux	DDS, PhD/Ethics	CHRU of Lille, University Hospital, University of Lille

Bernard Vandenbunder Vandenbunder PhD/Physics -Cell Biology -Biophotonics Centre National de la Recherche Scientifique Sylvie Vandoolaeghe MSc, Assistant Coordinator of the Lille EEHU $$ /Philosophy CHRU of Lille, University Hospital Isabelle Wolowczuk PhD/Immunology and metabolism Centre National de la Recherche Scientifique a Centre Hospitalier Régional et Universitaire. b Infirmier Anesthésiste Diplômé d'Etat. c Délégation à la Recherche Clinique et à l'Innovation. $$ RNA Registered Nurse Anaesthetist, Centre National de la Recherche Scientifique. EEHU Espace Ethique Hospitalier et Universitaire of reference and

Acknowledgments

The authors thank Catherine Brossaud (Senior Health Manager, IFCS, CHRU Lille) for her unceasing interest in the working group and its discussions; Nico Decock (Nurse Anaesthetist Trainer at the School of Nurse Anaesthetists, Lille-2 University) for his attentive and critical reading of the manuscript; and Dr David Fraser (Biotech Communication SARL, Damery, France) for editorial assistance.

In memoriam

The parasitologist Dr Eduardo Dei-Cas (1945-2014) was one of the leading experts on the opportunistic fungus Pneumocystis. In 2000, he initiated a debate on ethics in biology, medicine and research at Lille Universities, France. He was particularly interested in the debate concerning emergence and reductionism. Dr Dei-Cas notably participated in seminars as a member of the Centre for Medical Ethics (CME, Catholic Institute of Lille, Lille, France), organized a symposium on ethical and epistemological questions at the EMOP X (Paris, 2008), and shared his expertise during a course on "Dealing with the ethics of the living world" (Philosophy Research Unit, University of Lille 3 -Professor F. Worms) (2007-2012). In 2012, Dr Eduardo Dei-Cas launched the present "Ethics and Research" working group in the Ethics Department (Espace Ethique Hospitalier et Universitaire de Lille, EEHU) at Lille University and University Hospital.

Competing interest

The authors declare that they have no conflicts of interest.

Authors' contributions

All authors (SV, AB, DB, VB, EDC, BF, AG, HL, BO, TT, BV, IW and LD) collected and reviewed the literature, contributed to discussions on ethics and reviewed the manuscript for critical content. SV, AG, BF, BV, IW and LD wrote the manuscript. All authors read and approved the final manuscript.

Author details 1 Departement of Ethics, Espace Ethique Hospitalier et Universitaire de Lille (EEHU), Lille University Hospital, Faculty of Medicine, F-59000 Lille, France.

2 University of Lille, F-59000 Lille, France. 3 School of Dentistry, Lille University, F-59000 Lille, France. 4 University of Paris Descartes Laboratory of Medical Ethics and Forensic Medicine, EA 4569, F-75008 Paris, France. 5 School of Anaesthetist Nursing, Lille University Hospital, F-59000 Lille, France. 6 EA 4488 "Activité Physique, Muscle, Santé", Faculty of Sports Sciences, F-59790 Ronchin, France. 7 Department of Microbiology, Lille University Hospital, Faculty of Medicine, F-59000 Lille, France. 8 Biology and Diversity of Emerging Eukaryotic Pathogens Research Group (BDEEP), F-59019 Lille, France. 9 Center


F Xing

R Masson

S Lopez
email: s.lopez@brgm.fr

Parallel Vertex Approximate Gradient discretization of hybrid dimensional Darcy flow and transport in discrete fracture networks

This paper proposes a parallel numerical algorithm to simulate the flow and the transport in a discrete fracture network taking into account the mass exchanges with the surrounding matrix. The discretization of the Darcy fluxes is based on the Vertex Approximate Gradient finite volume scheme adapted to polyhedral meshes and to heterogeneous anisotropic media, and the transport equation is discretized by a first order upwind scheme combined with an Euler explicit integration in time. The parallelization is based on the SPMD (Single Program, Multiple Data) paradigm and relies on a distribution of the mesh on the processes with one layer of ghost cells in order to allow for a local assembly of the discrete systems. The linear system for the Darcy flow is solved using different linear solvers and preconditioners implemented in the PETSc and Trilinos libraries. The convergence of the scheme is validated on two original analytical solutions with one and four intersecting fractures. Then, the parallel efficiency of the algorithm is assessed on up to 512 processes with different types of meshes, different matrix fracture permeability ratios, and different levels of complexity of the fracture network.

Introduction 1.Hybrid dimensional flow and transport models

This article deals with the simulation of the Darcy flow and transport in fractured porous media for which the fractures are modeled as interfaces of codimension one. In this framework, the d -1 dimensional flow and transport in the fractures is coupled with the d dimensional flow and transport in the matrix leading to the so called hybrid dimensional Darcy flow and transport model.

For the Darcy flow model, we focus on the particular case where the pressure is continuous at the interfaces between the fractures and the matrix domain. This type of Darcy flow model introduced in [START_REF] Alboin | Modeling fractures as interfaces for flow and transport in porous media[END_REF], [START_REF] Amir | Décomposition de domaine et préconditionnement pour un modèle 3D en milieu poreux fracturé Proceeding[END_REF] corresponds physically to pervious fractures for which the ratio of the normal permeability of the fracture to the width of the fracture is large compared with the ratio of the permeability of the matrix to the size of the domain. Note that it does not cover the case of fractures acting as barriers for which the pressure is discontinuous at the matrix fracture interfaces (see [START_REF] Flauraud | Domain Decomposition for an asymptotic geological fault modeling[END_REF], [START_REF] Martin | Modeling fractures and barriers as interfaces for flow in porous media[END_REF], [START_REF] Angot | Asymptotic and numerical modelling of flows in fractured porous media[END_REF] for discontinuous pressure models). It is also assumed in our model that the pressure is continuous at the fracture intersections. It corresponds to the assumption that the ratio between the permeability at the fracture intersections and the width of the fracture is large compared to the ratio between the tangential permeability of each fracture and its length. We refer to [START_REF] Formaggia | A reduced model for Darcy's problem in networks of fractures[END_REF] and [START_REF] Schwenck | Dimensionally reduced flow models in fractured porous media: crossings and boundaries[END_REF] for more general reduced models taking into account discontinuous pressures at fracture intersections in dimension d = 2.

The hybrid dimensional transport model is derived in [START_REF] Alboin | Modeling fractures as interfaces for flow and transport in porous media[END_REF] in the case of a convection diffusion flux for the matrix and fracture concentration. In this work, a purely advective model is considered. It requires the specification of the transmission conditions at the matrix fracture interfaces and at fracture intersections which, to our knowledge, have not been done so far at the continuous level.

The discretization of the hybrid dimensional Darcy flow model with continuous pressures has been the object of several works. In [START_REF] Karimi-Fard | An efficient discrete-fracture model applicable for general-purpose reservoir simulators[END_REF] a cell-centred Finite Volume scheme using a Two Point Flux Approximation (TPFA) is proposed assuming the orthogonality of the mesh and isotropic permeability fields. Cell-centred Finite Volume schemes can be extended to general meshes and anisotropic permeability fields using MultiPoint Flux Approximations (MPFA) following the ideas of [START_REF] Tunc | A model for conductive faults with non matching grids[END_REF], [START_REF] Sandve | An efficient multi-point flux approximation method for discrete fracture-matrix simulations[END_REF], [START_REF] Ahmed | Control-volume distributed multi-point flux approximation coupled with a lower-dimensional fracture model[END_REF] and [START_REF] Ahmed | Threedimensional control-volume distributed multi-point flux approximation coupled with a lower-dimensional surface fracture model[END_REF]. In [START_REF] Alboin | Modeling fractures as interfaces for flow and transport in porous media[END_REF] and [START_REF] Hoteit | An efficient numerical model for incompressible two-phase flow in fracture media[END_REF] a Mixed Finite Element (MFE) method is proposed, and Control Volume Finite Element Methods (CVFE) using nodal unknowns have been introduced for such models in [START_REF] Bogdanov | Two-phase flow through fractured porous media[END_REF], [START_REF] Reichenberger | A mixed-dimensional finite volume method for multiphase flow in fractured porous media[END_REF], [START_REF] Monteagudu | Control-volume model for simulation of water injection in fractured media: incorporating matrix heterogeneity and reservoir wettability effects[END_REF], [START_REF] Matthai | A: Finite element-node-centred finite-volume two-phase-flow experiments with fractured rock represented by hybrid-element[END_REF], [START_REF] Geiger | Massively parallel sector scale discrete fracture and matrix simulation[END_REF]. The Hybrid Finite Volume and Mimetic finite difference schemes, belonging to the family of Hybrid Mimetic Mixed Methods [START_REF] Droniou | Gradient schemes: A generic framework for the discretisation of linear, nonlinear and nonlocal elliptic and parabolic equations[END_REF], have been extended to hybrid dimensional models in [START_REF] Faille | Model reduction and discretization using hybrid finite volumes of flow in porous media containing faults[END_REF], [START_REF] Antonietti | Mimetic finite difference approximation of flows in fractured porous media[END_REF] as well as in [START_REF] Brenner | Gradient discretization of hybrid dimensional Darcy flows in fractured porous media[END_REF], [START_REF] Brenner | Gradient discretization of hybrid dimensional Darcy flows in fractured porous media with discontinuous pressure at matrix fracture interfaces[END_REF] in the more general Gradient Discretization framework [START_REF] Droniou | The Gradient discretization method: A framework for the discretization of linear and nonlinear elliptic and parabolic problems[END_REF]. Non-matching discretizations of the fracture and matrix meshes are studied in [START_REF] D'angelo | A mixed finite element method for Darcy flow in fractured porous media with non-matching grids[END_REF], [START_REF] Fumagalli | A reduced model for flow and transport in fractured porous media with non-matching grids[END_REF], [START_REF] Berrone | An optimization approach for large scale simulations of discrete fracture network flows[END_REF] and [START_REF] Schwenck | Dimensionally reduced flow models in fractured porous media: crossings and boundaries[END_REF].

Regarding the hybrid dimensional advective transport model, an explicit first order upwind scheme combined with the MPFA Darcy fluxes is used in [START_REF] Ahmed | Control-volume distributed multi-point flux approximation coupled with a lower-dimensional fracture model[END_REF], [START_REF] Ahmed | Threedimensional control-volume distributed multi-point flux approximation coupled with a lower-dimensional surface fracture model[END_REF], and [START_REF] Sandve | An efficient multi-point flux approximation method for discrete fracture-matrix simulations[END_REF]. At fracture intersections, the authors neglect the accumulation term and the concentration unknown is eliminated using the flux conservation equation in order to avoid severe restrictions on the time step caused by the small volumes. A CVFE method is used in [START_REF] Reichenberger | A mixed-dimensional finite volume method for multiphase flow in fractured porous media[END_REF] with a first order upwind approximation and a fully implicit time integration of the two phase flow model to avoid small time steps. Higher order methods have also been developed in the CVFE method of [START_REF] Matthai | A: Finite element-node-centred finite-volume two-phase-flow experiments with fractured rock represented by hybrid-element[END_REF] using a MUSCL type second order scheme for the saturation equation and also in [START_REF] Hoteit | An efficient numerical model for incompressible two-phase flow in fracture media[END_REF] where a Discontinuous Galerkin method is used for the transport saturation equation with an Euler implicit time integration in the fracture network and an explicit time integration in the matrix domain. In [START_REF] Haegland | Streamline methods with application to flow and transport in fractured media[END_REF], a streamline method is developed in 2D based on the hybrid dimensional Darcy flow velocity field. The solution is very accurate for purely advective transport but this method requires that the fractures be expanded and seems difficult to extend to the case of a complex 3D network in practice.

Content and objectives of this work

In this work, we focus on the Vertex Approximate Gradient (VAG) scheme introduced in [START_REF] Eymard | Small-stencil 3D schemes for diffusive flows in porous media[END_REF] for diffusion problems and extended in [START_REF] Brenner | Vertex approximate gradient scheme for hybrid dimensional two-phase Darcy flows in fractured porous media[END_REF], [START_REF] Brenner | Gradient discretization of hybrid dimensional Darcy flows in fractured porous media[END_REF], [START_REF] Brenner | Gradient discretization of hybrid dimensional Darcy flows in fractured porous media with discontinuous pressure at matrix fracture interfaces[END_REF] to hybrid dimensional Darcy flow models. The VAG scheme uses nodal and fracture-face unknowns in addition to the cell unknowns which can be eliminated without any fill-in. Thanks to its essentially nodal feature, it leads to a sparse discretization on tetrahedral or mainly tetrahedral meshes. The VAG scheme has the major advantage, compared with the CVFE methods of [START_REF] Bogdanov | Two-phase flow through fractured porous media[END_REF], [START_REF] Reichenberger | A mixed-dimensional finite volume method for multiphase flow in fractured porous media[END_REF], [START_REF] Monteagudu | Control-volume model for simulation of water injection in fractured media: incorporating matrix heterogeneity and reservoir wettability effects[END_REF] or [START_REF] Matthai | A: Finite element-node-centred finite-volume two-phase-flow experiments with fractured rock represented by hybrid-element[END_REF], that it avoids the mixing of the control volumes at the fracture matrix interfaces, which is a key feature for its coupling with the transport model. As shown in [START_REF] Brenner | Vertex approximate gradient scheme for hybrid dimensional two-phase Darcy flows in fractured porous media[END_REF] for two phase flow problems, the VAG scheme allows for a coarser mesh size at the matrix fracture interface for a given accuracy. For the discretization of the transport hybrid dimensional model, we will use in this work a simple first order upwind scheme with explicit time integration. The extension to second order MUSCL type discretization will be considered in a future work. Our main objective in this paper is to develop a parallel algorithm for the VAG discretization of hybrid dimensional Darcy flow and transport models, and to assess the parallel scalability of the algorithm.

Starting from the hybrid dimensional Darcy flow model of [START_REF] Brenner | Vertex approximate gradient scheme for hybrid dimensional two-phase Darcy flows in fractured porous media[END_REF] and [START_REF] Brenner | Gradient discretization of hybrid dimensional Darcy flows in fractured porous media[END_REF], we first derive the hybrid dimensional transport model for a general fracture network taking into account fracture intersections and the coupling with the matrix domain. Then, the VAG discretization of the Darcy flow model is recalled and the VAG Darcy fluxes are used to discretize the transport model with an upwind first order discretization in space and an Euler explicit time integration. A key feature of this discretization is the definition of the control volumes which is adapted to the heterogeneities of the porous medium. This can be achieved thanks to the fact that, on the one hand, the VAG scheme keeps the cell unknowns and, on the other hand, the VAG Darcy fluxes are constructed independently of the definition of the control volumes. In particular, the control volumes are constructed in such a way that, at matrix fracture interfaces, the volume is taken only in the fracture. Otherwise, the fracture will be enlarged artificially and the front velocity will not be accurately approximated in the fractures as it it the case for usual CVFE methods. Note also that we do not eliminate the concentration unknowns at fracture intersections as was done in [START_REF] Sandve | An efficient multi-point flux approximation method for discrete fracture-matrix simulations[END_REF], [START_REF] Ahmed | Control-volume distributed multi-point flux approximation coupled with a lower-dimensional fracture model[END_REF] and [START_REF] Ahmed | Threedimensional control-volume distributed multi-point flux approximation coupled with a lower-dimensional surface fracture model[END_REF] for cell centred discretizations. In the case of a nodal discretization like the VAG scheme, this elimination is not possible since these unknowns are connected to the matrix and it is not needed since the size of the control volumes at fracture intersections is roughly the same as the size of any control volume located at the matrix fracture interface.

Our parallelization of the hybrid dimensional flow and transport numerical model is based on the SPMD (Single Program, Multiple Data) paradigm. It relies on a distribution of the mesh on the processes with one layer of ghost cells in order to allow for a local assembly of the discrete systems. The linear system for the Darcy flow is solved using different linear solvers and preconditioners implemented in the PETSc and Trilinos libraries.

In order to validate the convergence of the scheme, two analytical solutions are constructed for the hybrid dimensional flow and transport model. We consider the case of a single nonimmersed fracture as well as the case of four intersecting fractures. The analytical solutions for the transport model are obtained by integration of the matrix and fracture equations along the characteristics of the velocity field taking into account source terms coming from the matrix fracture transmission conditions. Then, we study the parallel scalability of the Darcy flow and transport solvers on up to 512 processes. Our numerical investigation includes different levels of complexity of the fracture network with a number of fractures ranging from a few to a few hundreds. It covers different types of meshes namely hexahedral, tetrahedral and prismatic meshes as well as a large range of permeability ratios between the fracture network and the matrix domain. In addition, the influence of the choices of the linear solver and of the preconditioner is also studied for the solution of the Darcy flow equation.

The paper is organized as follows. Section 2 recalls the geometrical and functional setting introduced in [START_REF] Brenner | Gradient discretization of hybrid dimensional Darcy flows in fractured porous media[END_REF] for a general 2D fracture network immersed in a surrounding 3D matrix. Then, the hybrid dimensional Darcy flow and transport models are introduced. In Section 3, the VAG discretization is recalled for the Darcy flow model and extended to the transport model. The parallel implementation of the scheme is detailed in section 4. Section 5 is devoted to the description of the test cases including the analytical solutions and to the numerical investigation of the parallel scalability of the algorithm. Let Ω denote a bounded domain of R d , d = 2, 3 assumed to be polyhedral for d = 3 and polygonal for d = 2. To fix ideas the dimension will be fixed to d = 3 when it needs to be specified, for instance in the naming of the geometrical objects or for the space discretization in the next section. The adaptations to the case d = 2 are straightforward.

We consider the asymptotic model introduced in [START_REF] Alboin | Modeling fractures as interfaces for flow and transport in porous media[END_REF] where fractures are represented as interfaces of codimension 1. Let I be a finite set and let Γ = i∈I Γ i and its interior Γ = Γ \ ∂Γ denote the network of fractures Γ i ⊂ Ω, i ∈ I, such that each Γ i is a planar, polygonal, simply connected, open domain included in an oriented plane P i of R d . It is assumed that the angles of Γ i are strictly smaller than 2π and that Γ i ∩ Γ j = ∅ for all i = j. For all i ∈ I, let us set

Σ i = ∂Γ i , Σ i,j = Σ i ∩ Σ j , j ∈ I \ {i}, Σ i,0 = Σ i ∩ ∂Ω, Σ i,N = Σ i \ ( j∈I\{i} Σ i,j ∪ Σ i,0
), and Σ = (i,j)∈I×I,i =j Σ i,j . It is assumed that Σ i,0 = Γ i ∩ ∂Ω. We will denote by dτ (x) the d -1 dimensional Lebesgue measure on Γ. On the fracture network Γ, we define the function space

L 2 (Γ) = {v = (v i ) i∈I , v i ∈ L 2 (Γ i ), i ∈ I}, endowed with the norm v 2 L 2 (Γ) = i∈I v i 2 L 2 (Γ i ) . Its subspace H 1 (Γ) is defined as the space of functions v = (v i ) i∈I such that v i ∈ H 1 (Γ i ), i ∈ I
with continuous traces at the fracture intersections. The space H 1 (Γ) is endowed with the norm

v 2 H 1 (Γ) = i∈I v i 2 H 1 (Γ i )
and its subspace with vanishing traces on Σ 0 = i∈I Σ i,0 is denoted by H 1 Σ 0 (Γ). Let us also consider the trace operator γ i from H 1 (Ω) to L 2 (Γ i ) as well as the trace operator γ from H 1 (Ω) to L 2 (Γ) such that (γv) i = γ i (v) for all i ∈ I.

On Ω, the gradient operator from H 1 (Ω) to L 2 (Ω) d is denoted by ∇. On the fracture network Γ, the tangential gradient ∇ τ acting from H 1 (Γ) to L 2 (Γ) d-1 is defined by

∇ τ v = (∇ τ i v i ) i∈I ,
where, for each i ∈ I, the tangential gradient ∇ τ i is defined from H 1 (Γ i ) to L 2 (Γ i ) d-1 by fixing a reference Cartesian coordinate system of the plane P i containing Γ i . We also denote by div τ i the divergence operator from H div (Γ i ) to L 2 (Γ i ).

The function spaces arising in the variational formulation of the hybrid dimensional Darcy flow model are

V = {v ∈ H 1 (Ω) such that γv ∈ H 1 (Γ)},
and its subspace

V 0 = {v ∈ H 1 0 (Ω) such that γv ∈ H 1 Σ 0 (Γ)}.
The space V 0 is endowed with the following Hilbertian norm

v V 0 = ∇v 2 L 2 (Ω) d + ∇ τ γv 2 L 2 (Γ) d-1 1/2 .
Let Ω α , α ∈ A denote the connected components of Ω \ Γ, with A being the set of connected components of Ω \ Γ. Let us define the space H div (Ω \ Γ) = {q m = (q m,α ) α∈A | q m,α ∈ H div (Ω α )}. Using the orientation of P i we can define the two sides ± of the fracture Γ i , for all i ∈ I. For all q m ∈ H div (Ω \ Γ), let γ ± n,i q m denote the normal trace of q m on the side ± of Γ i with the normal oriented outward from the side ±. Let us define the Hilbert function space

H(Ω, Γ) = { q m = (q m,α ) α∈A , q f = (q f,i ) i∈I | q m ∈ H div (Ω \ Γ), q f,i ∈ L 2 (Γ i ) d-1 , div τ i (q f,i ) -γ + n,i q m -γ - n,i q m ∈ L 2 (Γ i ), i ∈ I}, and its closed Hilbert subspace W (Ω, Γ) = {(q m , q f ) ∈ H(Ω, Γ) | α∈A Ωα (q m,α • ∇v + div(q m,α )v)dx + i∈I Γ i (q f,i • ∇ τ i γ i v + (div τ i (q f,i ) -γ + n,i q m -γ - n,i q m )γ i v)dτ (x) = 0 ∀ v ∈ V 0 }. (1)
The last definition corresponds to imposing in a weak sense the conditions i∈I γ n,Σ i q f,i = 0 on Σ \ Σ 0 and γ n,Σ i q f,i = 0 on Σ i,N , i ∈ I, where γ n,Σ i is the normal trace operator on Σ i (tangent to Γ i ) with the normal oriented outward from Γ i , and using the extension of γ n,Σ i q f,i by zero on Σ \ Σ i .

Hybrid dimensional Darcy Flow Model

In the matrix domain Ω \ Γ (resp. in the fracture network Γ), let us denote by

Λ m ∈ L ∞ (Ω) d×d (resp. Λ f ∈ L ∞ (Γ) (d-1)×(d-1)
) the permeability tensor such that there exist λ m ≥ λ m > 0 (resp.

λ f ≥ λ f > 0) with λ m |ξ| 2 ≤ (Λ m (x)ξ, ξ) ≤ λ m |ξ| 2 for all ξ ∈ R d , x ∈ Ω, (resp. λ f |ξ| 2 ≤ (Λ f (x)ξ, ξ) ≤ λ f |ξ| 2 for all ξ ∈ R d-1 , x ∈ Γ).
We also denote by µ the fluid viscosity which is assumed constant and by d f ∈ L ∞ (Γ) the width of the fractures assumed to be such that there exist

d f ≥ d f > 0 with d f ≤ d f (x) ≤ d f for all x ∈ Γ.
Given ū ∈ V , the strong formulation of the hybrid dimensional Darcy flow model amounts to: find u ∈ V and (q m , q f ) ∈ W (Ω, Γ) such that u -ū ∈ V 0 and

         div(q m,α ) = 0 on Ω α , α ∈ A, q m,α = -Λm µ ∇u on Ω α , α ∈ A, div τ i (q f,i ) -γ + n,i q m -γ - n,i q m = 0 on Γ i , i ∈ I, q f,i = -d f Λ f µ ∇ τ i γ i u on Γ i , i ∈ I. (2)
The weak formulation of (2) amounts to: find u ∈ V such that u -ū ∈ V 0 and the following variational equation is satisfied for all v ∈ V 0 :

Ω Λ m (x) µ ∇u(x) • ∇v(x)dx + Γ d f (x) Λ f (x) µ ∇ τ γu(x) • ∇ τ γv(x)dτ (x) = 0. (3)
The existence and uniqueness of the solution to (3) derives from the Lax Milgram theorem and a Poincaré inequality stated in [START_REF] Brenner | Gradient discretization of hybrid dimensional Darcy flows in fractured porous media[END_REF].

Hybrid dimensional transport model

Let γ n be the normal trace operator on ∂Ω with the normal oriented outward from Ω. Let us define

∂Ω -= {x ∈ ∂Ω | γ n q m (x) < 0}, Σ - i,0 = {x ∈ Σ i,0 | γ n,Σ i q f,i (x) < 0}, i ∈ I, as well as the following subset of Σ \ Σ 0 : Σ -= {x ∈ Σ \ Σ 0 | i∈I |γ n,Σ i q f,i (x)| = 0}.
We consider a linear, purely advective model with velocity q m in the matrix domain and q f in the fracture network. The matrix concentration is denoted by c m (c m,α in each connected component Ω α , α ∈ A) and the fracture concentration, representing the average concentration in the fracture width, is denoted by c f (c f,i in each fracture Γ i , i ∈ I). The 2D equation in the fracture network is as usual obtained by integration of the 3D advection equation in the width of the fractures. For a purely advective equation, the transmission condition at the matrix fracture interfaces states that the input normal flux in the matrix is obtained using the upwind fracture concentration c f . At the fracture intersection Σ -, an additional unknown c f,Σ -must be introduced and the transmission conditions state that the normal fluxes sum to zero and that the input normal fluxes are obtained using the upwind concentration c f,Σ -.

Let be given the input boundary conditions cm

∈ L ∞ (∂Ω -), cf,i ∈ L ∞ (Σ - i,0 ), i ∈ I, and the initial conditions c 0 m ∈ L ∞ (Ω\Γ), c 0 f ∈ L ∞ (Γ).
Let us denote by φ m (x) the porosity in the matrix and by φ f (x) the porosity in the fracture network. The transport hybrid dimensional model

amounts to: find c m ∈ L ∞ (Ω\Γ)×(0, T ) , c f ∈ L ∞ Γ×(0, T ) , and c f,Σ -∈ L ∞ Σ -×(0, T ) , such that:                                φ m ∂ t c m,α + div(c m,α q m,α ) = 0 on Ω α × (0, T ), α ∈ A φ f d f ∂ t c f,i + div τ i (c f,i q f,i ) = γ + n,i c m q m + γ - n,i c m q m on Γ i × (0, T ), i ∈ I, (γ ± n,i c m q m ) -= c f (γ ± n,i q m ) - on Γ i × (0, T ), i ∈ I, (γ n,Σ i c f,i q f,i ) -= c f,Σ -(γ n,Σ i q f,i ) - on (Σ i \ Σ i,0 ) × (0, T ), i ∈ I, j∈I γ n,Σ j c f,j q f,j = 0 on (Σ \ Σ 0 ) × (0, T ), (γ n c m q m ) -= cm (γ n q m ) - on ∂Ω × (0, T ), (γ n,Σ i c f,i q f,i ) -= cf,i (γ n,Σ i q f,i ) - on Σ i,0 × (0, T ), i ∈ I, c m = c 0 m on (Ω \ Γ) × {t = 0}, c f = c 0 f on Γ × {t = 0}, (4)
where the notations a + = max(a, 0) and a -= min(a, 0) are used for all a ∈ R.

3 Vertex Approximate Gradient Discretization (VAG)

VAG discretization of the Darcy flow model

In the spirit of [START_REF] Eymard | Small-stencil 3D schemes for diffusive flows in porous media[END_REF], we consider generalized polyhedral meshes of Ω in the sense that the edges of the mesh are linear but its faces are not necessarily planar. Roughly speaking, each face is assumed to be defined by the union of the triangles joining each edge of the face to a so called face centre. This definition has the advantage to include in particular hexahedral cells with non planar faces.

Let M be the set of cells which are disjoint open polyhedral subsets of Ω such that K∈M K = Ω. For each K ∈ M, it is assumed that there exists x K ∈ K \ ∂K, the socalled "centre" of the cell K, such that K is star-shaped with respect to x K . We then denote by F K the set of interfaces of non zero d -1 dimensional measure among the interior faces K ∩ L, L ∈ M \ {K}, and the boundary interface K ∩ ∂Ω, which possibly splits in several boundary faces. Let us denote by

F = K∈M F K
the set of all faces of the mesh. Note that the faces are not assumed to be planar, hence the term "generalized polyhedral mesh". For σ ∈ F, let E σ be the set of interfaces of non zero, d -2 dimensional measure among the interfaces σ ∩ σ , σ ∈ F \ {σ}. Then, we denote by

E = σ∈F E σ
the set of all edges of the mesh. Let V σ = (e,e )∈E 2 σ ,e =e e ∩ e be the set of nodes of σ. For each K ∈ M we define V K = σ∈F K V σ , and we also denote by

V = K∈M V K
the set of all nodes of the mesh. It is then assumed that for each face σ ∈ F, there exists a socalled "centre" of the face x σ ∈ σ \ e∈Eσ e such that x σ = s∈Vσ β σ,s x s , with s∈Vσ β σ,s = 1, and β σ,s ≥ 0 for all s ∈ V σ ; moreover the face σ is assumed to be defined by the union of the triangles T σ,e defined by the face centre x σ and each edge e ∈ E σ .

The mesh is also supposed to be conforming w.r.t. the fracture network Γ in the sense that for each i ∈ I there exists a subset

F Γ i of F such that Γ i = σ∈F Γ i
σ. We will denote by F Γ the set of fracture faces i∈I F Γ i . The following notations will be used for convenience:

M s = {K ∈ M | s ∈ V K }, M σ = {K ∈ M | σ ∈ F K }, F Γ,s = {σ ∈ F Γ | s ∈ V σ }, and F Γ,K = F K ∩ F Γ .
This geometrical discretization of Ω and Γ is denoted in the following by D.

The VAG discretization was introduced in [START_REF] Eymard | Small-stencil 3D schemes for diffusive flows in porous media[END_REF] for diffusive problems on heterogeneous, anisotropic media. Its extension to the hybrid dimensional Darcy model is based on the following vector space of unknowns:

X D = {v K , v s , v σ ∈ R, K ∈ M, s ∈ V, σ ∈ F Γ },
and its subspace with homogeneous Dirichlet boundary conditions on ∂Ω:

X 0 D = {v ∈ X D | v s = 0 for s ∈ V ext }.
where V ext = V ∩ ∂Ω denotes the set of boundary nodes, and V int = V \ ∂Ω denotes the set of interior nodes.

A finite element discretization of V is built using a tetrahedral sub-mesh of M and a second order interpolation at the face centres x σ , σ ∈ F \ F Γ defined by the operator I σ :

X D → R such that I σ (v) = s∈Vσ β σ,s v s .
The tetrahedral sub-mesh is defined by T = {T K,σ,e , e ∈ E σ , σ ∈ F K , K ∈ M} where T K,σ,e is the tetrahedron joining the cell centre x K to the triangle T σ,e (see Figure 2 for examples of such tetrahedra). For a given v ∈ X D , we define the function π T v ∈ V as the continuous piecewise affine function on each tetrahedron of T such that π

T v(x K ) = v K , π T v(s) = v s , π T v(x σ ) = v σ , and π T v(x σ ) = I σ (v) for all K ∈ M, s ∈ V, σ ∈ F Γ , and σ ∈ F \ F Γ .
The nodal basis of this finite element discretization will be denoted by η

K , η s , η σ , for K ∈ M, s ∈ V, σ ∈ F Γ .
The VAG discretization of the hybrid dimensional Darcy flow model ( 2) is based on its weak formulation [START_REF] Alboin | Modeling fractures as interfaces for flow and transport in porous media[END_REF]. Given ūs , s ∈ V ext , it amounts to: find u D ∈ X D with u s = ūs for all s ∈ V ext and such that for all v D ∈ X 0 D one has

Ω Λ m (x) µ ∇π T u D (x) • ∇π T v D (x)dx + Γ d f (x) Λ f (x) µ ∇ τ γπ T u D (x) • ∇ τ γπ T v D (x)dτ (x) = 0. (5)
Following [START_REF] Brenner | Gradient discretization of hybrid dimensional Darcy flows in fractured porous media[END_REF], this Galerkin Finite Element formulation (5) can be reformulated in terms of discrete conservation laws using the following definition of the VAG fluxes. For all v D ∈ X D , the VAG matrix fluxes connect the cell K ∈ M to its nodes and fracture faces ν ∈ V K ∪ F Γ,K :

F K,ν (v D ) = - K Λ m (x) µ ∇π T v D (x) • ∇η ν (x)dx = ν ∈V K ∪F Γ,K a ν K,ν (v K -v ν ) (6) with a ν K,ν = K Λm(x) µ ∇η ν (x) • ∇η ν (x)dx.
The VAG fracture fluxes connect the face σ ∈ F Γ to its nodes s ∈ V σ : Then, the Galerkin Finite Element formulation ( 5) is equivalent to: find u D ∈ X D satisfying the following set of discrete conservation equations and Dirichlet boundary conditions:

F σ,s (v D ) = - σ d f (x) Λ f (x) µ ∇ τ γπ T v D (x) • ∇ τ γη s (x)dτ (x) = s ∈Vσ a s σ,s (v σ -v s ) (7) with a s σ,s = σ d f (x) Λ f (x) µ ∇ τ γη s (x) • ∇ τ γη s (x)dτ (x).
                   s∈V K F K,s (u D ) + σ∈F Γ,K F K,σ (u D ) = 0, K ∈ M s∈Vσ F σ,s (u D ) + K∈Mσ -F K,σ (u D ) = 0, σ ∈ F Γ K∈Ms -F K,s (u D ) + σ∈F Γ,s -F σ,s (u D ) = 0, s ∈ V int , u s = ūs , s ∈ V ext .

First order upwind discretization of the transport model 3.2.1 Definition of control volumes

The VAG discretization of the hybrid dimensional transport model combines the VAG matrix and fracture fluxes ( 6), [START_REF] Balay | PETSc Users Manual[END_REF] with the following definition of the control volumes based on partitions of the cells and of the fracture faces. These partitions are respectively denoted, for all K ∈ M, by

K = ω K ∪ s∈V K \(Vext∪V Γ ) ω K,s
and, for all σ ∈ F Γ , by

σ = ω σ ∪ s∈Vσ\Vext ω σ,s .
Then, the control volumes are defined by ω K for all cells K ∈ M, by ω σ for all fracture faces σ ∈ F Γ , and by ω s = K∈Ms ω K,s , for all nodes s ∈ V int \ V Γ , and by

ω s = σ∈F Γ,s ω σ,s , for all nodes s ∈ V Γ \ V ext .
Note that this definition avoid the mixing of the fracture and matrix rocktypes at the control volumes s ∈ V Γ \ V ext and σ ∈ F Γ . This is exhibited in Figure 4 in comparison with an alternative choice mixing the matrix and fracture rocktypes which artificially enlarges the fractures. We refer to [START_REF] Brenner | Gradient discretization of hybrid dimensional Darcy flows in fractured porous media[END_REF] for numerical comparisons on a two phase flow model of these two types of choices of the control volumes. The same idea is applied for all nodes located at different rocktype interfaces. In practice, for such a node s ∈ V int \ V Γ (resp. s ∈ V Γ \ V ext ), the set ω K,s (resp. ω σ,s ) should be non empty only for the cell(s) K (resp. fracture face(s) σ) with the largest permeability among those around the node s (see [START_REF] Eymard | Vertex centred discretization of compositional multiphase Darcy flows on general meshes[END_REF] for details).

In practice, the above partitions of the cells and fracture faces does not need to be built. It is sufficient to define the matrix volume fractions

α K,s = ω K,s dx K dx , s ∈ V K \ (V ext ∪ V Γ ), K ∈ M,
constrained to satisfy α K,s ≥ 0, and s∈V K \(Vext∪V Γ ) α K,s ≤ 1, as well as the fracture volume fractions

α σ,s = ωσ,s d f (x)dτ (x) σ d f (x)dτ (x) , s ∈ V σ \ V ext , σ ∈ F Γ ,
such that α σ,s ≥ 0, and s∈Vσ\Vext α σ,s ≤ 1. Then, the porous volumes of the control volumes are set to

φ K = (1 - s∈V K \(V Γ ∪Vext) α K,s ) K φ m (x)dx, K ∈ M φ σ = (1 - s∈Vσ\Vext α σ,s )d f,σ σ φ f (x)dτ (x), σ ∈ F Γ , φ s = σ∈F Γ,s α σ,s d f,σ σ φ f (x)dτ (x), s ∈ V Γ \ V ext , φ s = K∈Ms α K,s K φ m (x)dx, s ∈ V \ (V ext ∪ V Γ ), with d f,σ = σ d f (x)dτ (x) σ dτ (x)
.

Time integration

For N ∈ N * , let us consider the time discretization

t 0 = 0 < t 1 < • • • < t n-1 < t n • • • < t N = T of the time interval [0, T ].
We denote the time steps by ∆t n = t n+1 -t n for all n = 0,

• • • , N -1.
Given cs , s ∈ V ext with arbitrary values on the set of ouput boundary nodes andc

V + ext = {s ∈ V ext | F K,s (u D ) ≥ 0 ∀K ∈ M s and F σ,s (u D ) ≥ 0 ∀σ ∈ F Γ,s },
                       φ K c n+1 K -c n K ∆t n + s∈V K H K,s (c n D ) + σ∈F Γ,K H K,σ (c n D ) = 0, K ∈ M, φ σ c n+1 σ -c n σ ∆t n + s∈Vσ H σ,s (c n D ) - K∈Mσ H K,σ (c n D ) = 0, σ ∈ F Γ , φ s c n+1 s -c n s ∆t n - K∈Ms H K,s (c n D ) - σ∈F Γ,s H σ,s (c n D ) = 0, s ∈ V int , c n+1 s = cs , s ∈ V ext ,
with the following explicit upwind two point fluxes

H K,ν (c n D ) = c n K F K,ν (u D ) + + c n ν F K,ν (u D ) - H σ,s (c n D ) = c n σ F σ,s (u D ) + + c n s F σ,s (u D ) -. (8)
The solution of this explicit upwind scheme classically satisfies the following maximum principle

m ≤ c n+1 µ ≤ M for all µ ∈ V ∪ F Γ ∪ M \ V + ext , with M = max µ∈V∪F Γ ∪M\V + ext c 0 µ and m = min µ∈V∪F Γ ∪M\V + ext c 0 µ ,
provided that the following Courant-Friedrichs-Lewy (CFL) condition

∆t n ≤ min(∆t M , ∆t F Γ , ∆t V ), (9)
is satisfied with

               ∆t M = min K∈M φ K s∈V K F K,s (u D ) + + σ∈F Γ,K F K,σ (u D ) + , ∆t F Γ = min σ∈F Γ φ σ s∈Vσ F σ,s (u D ) + + K∈Mσ (-F K,σ (u D )) + , ∆t V = min s∈V int φ s K∈Ms (-F K,s (u D )) + + σ∈F Γ,s (-F σ,s (u D )) + .

Parallel implementation in ComPASS

The hybrid dimensional Darcy flow and transport discrete model is implemented in the framework of the code ComPASS (Computing Parallel Architecture to Speed up Simulations) [START_REF] Dalissier | ComPASS: a tool for distributed parallel finite volume discretizations on general unstructured polyhedral meshes[END_REF], which focuses on parallel high performance simulation (distributed memory, Message Parsing Interface -MPI) adapted to general unstructured polyhedral meshes (see [START_REF] Eymard | High performance computing linear algorithms for two-phase flow in porous media[END_REF]).

Mesh non overlapping and overlapping decompositions

Let us denote by N p the number of MPI processes. The set of cells M is partitioned into N p subsets M p , p = 1, ..., N p using the library METIS [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF]. The partitioning of the set of nodes V and of the set of fracture faces F Γ is defined as follows: assuming we have defined a global index of the cells K ∈ M let us denote by K(s), s ∈ V (resp. K(σ), σ ∈ F Γ ) the cell with the smallest global index among those of M s (resp. M σ ). Then we set

V p = {s ∈ V | K(s) ∈ M p }, and
F p Γ = {σ ∈ F Γ | K(σ) ∈ M p }. The overlapping decomposition of M into the sets M p , p = 1, ..., N p ,
is chosen in such a way that any compact finite volume scheme such as the VAG scheme can be assembled locally on each process. Hence, as exhibited in Figure 5, M p is defined as the set of cells sharing a node with a cell of M p . The overlapping decompositions of the set of nodes and of the set of fracture faces follow from this definition:

V p = K∈M p V K , p = 1, • • • , N p , and
F p Γ = K∈M p F K ∩ F Γ , p = 1, • • • , N p .
Figure 5: Example of mesh decomposition.

The partitioning of the mesh is performed by the master process (process 1), and then, each local mesh is distributed to its process. Therefore, each MPI process contains the local mesh (M p , V p , F p Γ ), p = 1, 2, ..., N p which is splitted into two parts:

own mesh: (M p , V p , F p Γ ), ghost mesh: (M p \M p , V p \V p , F p Γ \F p Γ ).
We now turn to the parallel implementation of the discrete hybrid dimensional Darcy flow model (2) and transport model (4).

Parallelization of the discrete hybrid dimensional Darcy flow

On each process p = 1, ..., N p , the local problem of the discrete hybrid dimensional Darcy flow (2) is defined by the set of unknowns

u µ , µ ∈ V p ∪ F p Γ ∪ M
p and the set of equations

                   s∈V K F K,s (u D ) + σ∈F Γ,K F K,σ (u D ) = 0, K ∈ M p , s∈Vσ F σ,s (u D ) + K∈Mσ -F K,σ (u D ) = 0, σ ∈ F p Γ , K∈Ms -F K,s (u D ) + σ∈F Γ,s -F σ,s (u D ) = 0, s ∈ V int ∩ V p , u s = ūs , s ∈ V ext ∩ V p . (10)
Note that this includes the equations of the nodes s ∈ V p , of the fracture faces σ ∈ F p Γ and of the cells K ∈ M p , both those own cells in K ∈ M p and the ghost cells K ∈ M p \ M p . The set of equations can be rewritten as the following rectangular linear system

  A p vv A p vf A p vc A p f v A p f f A p f c A p cv A p cf A p cc     U p v U p f U p c   =   b p v b p f b p c   (11)
where

U p v ∈ R V p , U p f ∈ R F p Γ and U p c ∈ R M p
denote the vector of process p own and ghost unknowns at nodes, fracture faces and cells respectively. The above matrices have the following sizes

A p vv ∈ R V p ×V p , A p vf ∈ R V p ×F p Γ , A p vc ∈ R V p ×M p , A p f v ∈ R F p Γ ×V p , A p f f ∈ R F p Γ ×F p Γ , A p f c ∈ R F p Γ ×M p , A p cv ∈ R M p ×V p , A p cf ∈ R M p ×F p Γ , A p cc ∈ R M p ×M p . and b p v ∈ R V p , b p f ∈ R F p Γ , b p c ∈ R M p
denote the corresponding right hand side vectors. The matrix A p cc is a non singular diagonal matrix and the cell unknowns can be easily eliminated without fill-in leading to the following Schur complement system

A p U p v U p f = b p (12)
with

A p := A p vv A p vf A p f v A p f f - A p vc A p f c (A p cc ) -1 A p cv A p cf , b p := b p v b p f - A p vc A p f c (A p cc ) -1 b p c , and
U p c = (A p cc ) -1 (b p c -A p cv U p v -A p cf U p f ). (13)
The linear system ( 12) is built locally on each process p and transfered to the parallel linear solver library PETSc [START_REF] Balay | PETSc Users Manual[END_REF] or Trilinos [START_REF] Heroux | Trilinos Users Guide[END_REF]. The parallel matrix and the parallel vector in PETSc or Trilinos are stored in a distributed manner, i.e. each process stores its own rows. We construct the following parallel linear system

AU = b, (14)
with

A :=      A 1 R 1 A 2 R 2 . . . A Np R Np      process 1 process 2 . . . process N p , U :=        U 1 v U 1 f U 2 v U 2 f . . .        process 1 process 2 . . . , b :=      b 1 b 2 . . . b N p      process 1 process 2 . . . process N p
where R p , p = 1, 2, ..., N p is a restriction matrix satisfying

R p U = U p v U p f . The matrix A p R p ,
where

U :=         U 1 v U 1 f U 2 v U 2 f . . .        
is the vector of own and ghost node and fracture-face unknowns on all processes. The matrix S, containing only 0 and 1 entries, is assembled once and for all at the beginning of the simulation. Finally, thanks to (13), the vector of own and ghost cell unknowns U p c is computed locally on each process p.

In conclusion, the parallel implementation of the discrete hybrid dimensional Darcy flow can be summarized as: Algorithm 1: Parallel implementation of the discrete hybrid dimensional Darcy flow 1: Assemble locally on each process the rectangular linear system (11), 2: Compute locally on each process the Schur complement ( 12) of ( 11), 3: Construct the parallel linear system [START_REF] Chen | Well flow models for various numerical methods[END_REF] in PETSc or Trilinos, 4: Solve the parallel linear system ( 14) to obtain the solution at own nodes and fracture faces, 5: Communicate the solution at ghost nodes and fracture faces from (15), 6: Compute locally on each process the solution at own and ghost cells from (13).

Parallelization of the discrete hybrid dimensional transport model

The parallel implementation of the transport model ( 4) with an explicit upwind discretization of the fluxes consists of the following four steps.

1. Compute the Darcy matrix and fracture fluxes defined by ( 6) and ( 7).

2. Compute the maximum time step ∆t satisfying the CFL condition ( 9) and set ∆t n = ∆t for all n = 0, • • • , N -2, and ∆t

N -1 = T -(N -1)∆t with N = T ∆t . 3. For each time step n = 0, 1, ..., N -2, 3a. Compute c n+1 s , c n+1 σ and c n+1 K , s ∈ V int ∩V p , σ ∈ F p Γ , K ∈ M p solution of the following explicit equations                        φ K c n+1 K -c n K ∆t + s∈V K H K,s (c n D ) + σ∈F Γ,K H K,σ (c n D ) = 0, K ∈ M p , φ σ c n+1 σ -c n σ ∆t + s∈Vσ H σ,s (c n D ) - K∈Mσ H K,σ (c n D ) = 0, σ ∈ F p Γ , φ s c n+1 s -c n s ∆t - K∈Ms H K,s (c n D ) - σ∈F Γ,s H σ,s (c n D ) = 0, s ∈ V int ∩ V p , c s = cs , s ∈ V ext ∩ V p . ( 16
)
3b. Get the node and fracture-face ghost unknowns

c n+1 s , c n+1 σ , s ∈ V int ∩ (V p \V p ), σ ∈ F p Γ \F p
Γ using the PETSc or Trilinos matrix vector product with the matrix S defined in [START_REF] Dalissier | ComPASS: a tool for distributed parallel finite volume discretizations on general unstructured polyhedral meshes[END_REF], as was done for the Darcy flow solution U :

        C 1 v C 1 f C 2 v C 2 f . . .         = S        C 1 v C 1 f C 2 v C 2 f . . .       
Thanks to our mesh decomposition, step 1 and step 3a are performed locally on each process. For step 2, the maximum time step ∆t p is computed locally on each process p, then the time step ∆t is obtained using the MPI reduce operation.

In the following numerical experiments the parameters are set to tan(θ) = 1 2 , Λ f = 20 and d f = 0.01. The mesh is a topologically Cartesian n x × n x grid. Figure 6 shows an example of the mesh with n x = 20 as well as the analytical solution in the matrix obtained at time t f = 0.5 chosen as the final time of the simulation. The time step is defined by the maximum time step allowed by the CFL condition [START_REF] Bogdanov | Two-phase flow through fractured porous media[END_REF]. Figure 8 exhibits the convergence of the relative L 1 errors between the analytical solution and the numerical solution at time t f both in the matrix domain and in the fracture as a function of the grid size n x = 100, 200, 400, 800, 1600. Figure 7 shows the analytical solution and the numerical solutions obtained at time t f along the fracture. In both cases, we can observe the expected convergence of the numerical solution to the analytical solution with a higher order of convergence in the fracture due to the fact that at time t f the analytical solution in the fracture is continuous as exhibited in Figure 7. (

) 18
where

L 1 = ∂ t + k 1 ∂ x + β 1 and L 2 = ∂ t -k 2 ∂ y + β 2 .
This system can also be integrated analytically along the characteristics of the matrix and fracture velocity fields, but it leads to complex computations. It is much easier to obtain the stationary solution of this system which is defined in the fractures by

c f,1 (x) = e - β 1 k 1 x , c f,4 (y) = e - β 2 k 2 (1-y) , c 0 = e - β 1 k 1 x 0 + re - β 2 k 2 (1-y 0 ) r + 1 , c f,2 (x) = e - β 1 k 1 x c 0 e β 1 k 1 x 0 + β 1 k 1 r 1 e (r 1 x- 3β 2 4k 2
) -e

(r 1 x 0 - 3β 2 4k 2 ) , c f,3 (y) =    e β 2 k 2 y c 0 e - β 2 k 2 y 0 + β 2 k 2 r 2 e (- r 2
4 - β 1 4k 1 tan(θ 1 ) ) -e (-r 2 y 0 - β 1 4k 1 tan(θ 1 ) ) , if y < 1 4 , e β 2 k 2 y c 0 e - β 2 k 2 y 0 + β 2 k 2 r 2 e (-r 2 y- β 1 4k 1 tan(θ 1 ) ) -e (-r 2 y 0 - β 1 4k 1 tan(θ 1 ) ) if y > 1 4 , with r 1 = β 1 k 1 + β 2 k 2 tan(θ 1 ) and r 2 = β 2 k 2 + β 1 k 1 tan(θ 1 )
, and in the matrix by

c m,1 (x, y) = 0, c m,2 (x, y) = c f,4 (y), c m,3 (x, y) = 0 if y < 1 4 , c f,1 y-1 4 tan(θ 1 ) if y > 1 4 , c m,4 (x, y) = c f,3 (y) if y < y 0 , c f,2 y-1 4 tan(θ 1 )
if y > y 0 .

In the following numerical experiments the parameters are set to tan(θ 1 ) = 5 8 , tan(θ 2 ) = 1 4 , Λ f,1 = 200, Λ f,3 = 400, and d f,1 = d f,3 = 0.01. The mesh is, as for the previous test case, a topologically Cartesian n x × n x grid exhibited in Figure 9 for n x = 20. Figure 9 also shows the stationary analytical solution in the matrix. The time step is again defined by the maximum time step allowed by the CFL condition ( 9) and the simulation time is chosen large enough to obtain the numerical stationary solution.

Figure 10 exhibits the convergence of the relative L 1 errors between the stationary analytical and the numerical solutions both in the matrix domain and in the fracture as a function of the grid size n x = 100, 200, 400, 800. We can again observe the expected convergence of the numerical solution to the analytical solution with a higher order of convergence in the fracture network due to the fact that the solution is continuous on each individual fracture as exhibited in Figure 10. This property is always true when looking at the solution at the matrix time scale and could be exploited in a future work by using an implicit time integration in the fracture coupled to an explicit time integration in the matrix domain with a higher order discretization in space in the spirit of what has been done in [START_REF] Hoteit | An efficient numerical model for incompressible two-phase flow in fracture media[END_REF].

Fracture network with hexahedral meshes

The objective of this subsection and of the next subsection is to investigate the parallel scalability of the algorithms described in Section 4. In this subsection we consider a topologically Cartesian mesh of size n x × n x × n x of the cubic domain (0, 1) 3 exhibited in Figure 11 for n x = 32. The mesh is exponentially refined at the interface between the matrix and the fracture network exhibited in Figure 12. The permeabilities are isotropic and set to Λ f = 20 in the fracture network and to Λ m = 1 in the matrix. The porosities are set to φ m = φ f = 1 and the fluid viscosity is set to µ = 1. The initial concentration is set to 0 both in the matrix domain and in the fracture network and a value of 1 of the concentration is injected on the bottom boundaries of the matrix and of the fracture network. The pressure is fixed to u = 1 and γu = 1 on the bottom boundary and to u = 0 and γu = 0 on the top boundary. The remaining lateral boundaries are considered impervious. Figure 12 exhibits the tracer concentrations obtained with the mesh n x = 128 at times t = 0, t = 0.2, t = 0.4 and at the final simulation time t f = 0.5. Table 1 presents the numbers of linear solver iterations for the stationary pressure solution for a number of MPI processes ranging from N p = 2 to N p = 512 and with the mesh size n x = 128 corresponding to roughly 2.1 × 10 6 cells, 2.1 × 10 6 nodes and 5.2 × 10 4 fractures faces. Both the GMRES and BiCGStab linear solvers from the PETSc library are tested combined with either the Boomer AMG preconditioner from the Hypre library [START_REF]Hypre -Parallel high performance preconditioners[END_REF], the Aggregation AMG preconditioner from the Trilinos library [START_REF] Heroux | Trilinos Users Guide[END_REF] or the block Jacobi ILU(0) preconditioner from the Euclid library. No restart is used for the GMRES linear solver. Table 2 shows the corresponding computation times both for the setup phase of the preconditioner and for the solve phase of the linear solver.

According to these tables, the GMRES and the BiCGStab linear solvers combined with the Boomer AMG preconditioner are good choices for a number of processes N p 128, while the BiCGStab linear solver combined with the block Jacobi ILU(0) preconditioner is more efficient for this test case for N p = 256 and N p = 512. This was expected since the Boomer AMG preconditioner requires a sufficiently large number of unknowns per core to maintain a good parallel scalability due to the high level of communications in particular in the setup phase of the preconditioner. For this linear system, the number of unknowns per MPI process is roughly 4100 for N p = 512 which is too small for this type of preconditioner while the block Jacobi preconditioner still maintains a good parallel scalability for such a number of unknowns per MPI process. On the other hand, as shown in Table 3, Boomer AMG exhibits an optimal scalability while ILU(0) is not scalable in terms of iteration count with respect to the ratio

Λ f
Λm between the fracture and matrix permeabilities. The same remark also holds in terms of scalability with respect to the mesh size which means that the ILU(0) preconditioner is only advantageous for small size and mildly heterogeneous problems. Tables 1 and2 also clearly show that the BiCGStab linear solver outperforms the GMRES linear solver for cases requiring a large number of iterations due to the fact that the cost of the orthogonalization procedure increases with the dimension of the Krylov subspace. The Aggregation AMG preconditioner yields a larger number of iterations compared with the Boomer AMG preconditioner but has a much lower setup time resulting for this test case in a total lower CPU time. However, this implementation of the Aggregation AMG preconditioner seems to lack robustness with respect to the matrix fracture permeability ratio as exhibited in Table 3. Next, Figure 13 plots the total (Darcy flow and transport models) computation time and the computation time for the transport model only as a function of the number of processes. In these runs the GMRES linear solver is used combined with the Boomer AMG preconditioner for N p 128 and with the ILU(0) preconditioner for N p = 256, 512. For the range 2 -512 of the number of processes, it appears that the computation time of the Darcy flow linear system solution remains small compared with the transport model computation time. This can be checked by comparison of the computation times in Table 2 and in Figure 13. This explains the good scalability obtained for both the total and transport computation times thanks to the explicit nature of the time integration scheme.

Fracture network with tetrahedral meshes

This test case considers tetrahedral meshes of the cubic domain (0, 1) 3 conforming to the fracture network. An example of tetrahedral mesh showing both the matrix domain and the fracture network is exhibited in Figure 14. All the physical parameters, initial and boundary conditions are the same as for the previous test case. The mesh used in this subsection contains about 6.2 × 10 6 cells, 9.7 × 10 5 nodes and 7.1 × 10 4 fracture faces. Figure 15 exhibits the tracer concentrations obtained with this tetrahedral mesh at times t = 0, t = 0.2, t = 0.4 and at the final simulation time t f = 0.5.

As for the previous test case, Tables 4, 5 and 6 investigate the performance of the Darcy flow system linear solution for both the GMRES and BiCGStab linear solvers and for the same three preconditioners as in the previous test case. The conclusions are basically the same as for the hexahedral mesh test case. The Boomer AMG preconditioner exhibits an optimal robustness with respect to the matrix fracture permeability ratio

Λ f
Λm . On the other hand it requires a rather high number of unknowns per MPI process to maintain a good parallel scalability due to the high level of communications in particular in the setup phase. The ILU(0) preconditioner can be an interesting alternative but only for small size and midly heterogeneous problems. The aggregation AMG preconditioner from the Trilinos library used in our test cases seems to lack robustness and we did not manage to make it work better through tuning its parameters.

Figure 16 plots the total (Darcy flow and transport models) computation time and the computation time of the transport model only as a function of the number of processes. In these runs the GMRES linear solver is used combined with the Boomer AMG preconditioner for N p 128 and with the ILU(0) preconditioner for N p = 256, 512. Compared with the previous subsection, an even better parallel scalability of the transport model computation time is observed in the right Figure 16. This can be explained by the ratio of roughly 6 between the number of cells and the number of nodes typical of a tetrahedral mesh. For a topologically Cartesian mesh, this ratio is roughly 1. Since the cell concentrations are computed locally in each process, this explains the better scalability observed for this tetrahedral mesh compared with the previous hexahedral mesh. On the left Figure 16, it is observed that the linear system solution computation time is no longer small compared with the transport computation time for N p = 256 and 512. Hence, it significantly reduces the parallel efficiency of the simulation for a large number of processes, say N p = 256, 512 in this test case.

Application to a complex fracture network

In this subsection, our algorithm is applied to a complex fracture network kindly provided by M. Karimi-Fard and A. Lapène from Stanford University and TOTAL. Figure 17 exhibits the mesh of the domain Ω = (0, 5888.75) × (0, 3157.5) × (0, 250) (m) which contains about 1.2 × 10 7 prismatic cells, 6.5 × 10 6 nodes and 5.13 × 10 5 fracture faces. This 3D mesh is defined by the tensor product of a triangular 2D mesh with a uniform vertical 1D mesh with 24 intervals. The fracture network exhibited in Figure 18 contains 581 connected components. It is a set of 21376 × 24 faces of the 3D mesh defined by the tensor product of a subset of 21376 edges of the triangular 2D mesh with the 1D vertical mesh. The 2D triangular mesh contains 517540 cells and is refined in the neighbourhood of the fracture network down to an average size of 3.5 m. Figure 18 also shows the location of the injection well and of the two production wells. Each well is vertical of radius r w = 0.1 m and its centre in the horizontal plane is located at the middle of a fracture edge in the 2D triangular mesh. In the vertical direction, only the 12 fracture faces at the center of the 1D mesh are perforated. The permeabilities are isotropic and set to Λ f = 10 -11 (m 2 ) in the fracture network and to Λ m = 10 -15 (m 2 ) in the matrix domain. The porosities are set to φ m = φ f = 0.1, the fracture width to d f = 1 m and the fluid viscosity to µ = 10 -3 Pa.s -1 .

The initial concentration is set to 0 both in the matrix domain and in the fracture network. A total volume of 5.0×10 6 m 3 is injected in one year at the injector well with a tracer concentration of 1. The pressures of each perforated fracture face σ of the producer wells are fixed to p w = 0 and the flow rates are given by the Peaceman model

q σ = W I σ (p σ -p w ),
where p σ is the pressure in the fracture face and W I σ the well index of the fracture face. This well index is computed following Peaceman methodology [START_REF] Peaceman | Interpretation of well-block pressures in numerical reservoir simulation[END_REF], [START_REF] Peaceman | Interpretation of well-block pressures in numerical reservoir simulation with nonsquare grid blocks and anisotropic permeability[END_REF], [START_REF] Chen | Well flow models for various numerical methods[END_REF] by expanding the fracture face as a box of size dx × d f × dz. The analytical pressure solution obtained for a vertical well with the well pressure p w , the well radius r w and the well flow rate q w per unit length is imposed at the 8 corners of the box. Then, the flow rate q w dz is imposed at the box center and the pressure p c at the box center can be computed analytically using the VAG scheme. We deduce the well index W I = q w dz p c -p w leading in this simple case to the analytical formula the tracer concentration after one year of injection and at final time both in the matrix domain and in the fracture network. Figure 23 shows the total computation times with different number of MPI processes N p = 16, 32, 64, 128, 256, 512. It is observed that the total computation time exhibits a rather good scalability. In addition, the linear solver (GMRES+Boomer AMG) for the pressure converges in no more than 25 iterations whatever the number of MPI processes. Also the comparison of the total and transport computation times in Figure 23 shows that the time for the pressure solution remains small compared with the transport computation time up to N p = 512.

Conclusion

This paper introduced a parallel VAG scheme for the simulation of a hybrid dimensional Darcy flow and transport model in a discrete fracture network taking into account the mass exchanges with the matrix. The convergence of the scheme was validated on two original analytical solutions for a flow and transport model that includes fractures. The parallel efficiency of the algorithm was studied for different complexities of fracture networks, and a large range of matrix fracture permeability ratios and different type of meshes. The numerical results exhibit a very good parallel strong scalability as expected from the explicit nature of the time integration of the transport model with a better result on tetrahedral meshes thanks to the communication free computation of the cell unknowns. The Darcy flow solution is remarkably robust using the Boomer AMG preconditioner on all types of fracture networks, meshes and for all permeability ratios that have been tested. On the other hand, it requires as usual a rather high number of unknowns per process to maintain a good parallel scalability. Future work includes the extension of the parallel algorithm to hybrid dimensional multiphase flow models and the use of a more accurate second order MUSCL scheme for the transport model.

Figure 1 :

 1 Figure 1: Example of a 2D domain with 3 intersecting fractures Γ 1 , Γ 2 , Γ 3 and 2 connected components Ω 1 , Ω 2 .

Figure 2 :

 2 Figure 2: For a cell K with one fracture face σ in bold: cell unknown v K located at x K , fractureface unknown v σ located at x σ , node unknowns v s , v s , face centre x σ of face σ , triangle T σ,e (convex hull of e and x σ ), triangle T σ ,e (convex hull of e and x σ ) and tetrahedron T K,σ ,e (convex hull of x K , x σ and e ).

Figure 3 :

 3 Figure 3: Matrix fluxes (in blue) and fracture fluxes (in red) inside a cell K with a fracture face σ (in bold). The matrix fluxes F K,ν connect the cell K to its nodes and fracture faces ν ∈ V K ∪ F Γ,K . The fracture fluxes F σ,s connect the face σ to the nodes s ∈ V σ of σ.

Figure 4 :

 4 Figure 4: Example of choices of the control volumes at cells, fracture face, and nodes, in the case of two cells K and L splitted by one fracture face σ (the width of the fracture has been enlarged in this figure). The left figure exhibits the good choice with no mixing of fracture and matrix rocktypes while the right figure exhibits the bad choice enlarging artificially the fracture.

  0 D ∈ X D such that c 0 s = cs for all s ∈ V ext , the transport discrete model amounts to find c n+1 D ∈ X D for all n = 0, • • • , N -1 satisfying the following discrete conservation laws and Dirichlet input conditions

  b p are stored in process p. The parallel linear system[START_REF] Chen | Well flow models for various numerical methods[END_REF] is solved using the GMRES or BiCGStab algorithm preconditioned by different type of preconditioners as discussed in the numerical section. The solution of the linear system provides on each process p the solution vector U and fracture-face unknowns. Then, the ghost node unknowns u µ , µ ∈ (V p \V p ) and the ghost fracture-face unknowns u µ , µ ∈ (F p Γ \F p Γ ) are recovered by a synchronization step with MPI communications. This synchronization is efficiently implemented using a PETSc or Trilinos matrix vector product U = SU

Figure 6 :

 6 Figure 6: Left: example of mesh with n x = 20 where the red line is the fracture. Right: analytical solution of (17) at time t f = 0.5.

Figure 7 :

 7 Figure 7: Left: analytical solution and numerical solutions along the fracture at time t f with n x = 100, 200, 400, 800, 1600. Right: zoom view of left figure around x = 0.5.

Figure 9 :

 9 Figure 9: Left: example of mesh with n x = 20 where the red lines account for the four fractures. Right: stationary analytical solution of (18).

Figure 10 :

 10 Figure 10: Left: stationary analytical solution in the four fractures as a function of the x coordinate. Right: relative L 1 errors and references in the matrix domain and in the fracture network between the stationary analytical and numerical solutions as a function of the grid size n x = 100, 200, 400, 800.

Figure 11 :

 11 Figure 11: Hexahedral mesh in the matrix domain (left) and in the fracture network (right) with n x = 32.

Figure 12 :

 12 Figure 12: Concentration in the matrix domain and in the fracture network obtained at different times with the mesh n x = 128. A threshold concentration of 0.2 is used in the matrix domain.

Figure 13 :

 13 Figure 13: Total computation time (left) and computation time for the transport model (right) vs. the number of MPI processes for the mesh size n x = 128.

Figure 14 :

 14 Figure 14: Example of tetrahedral mesh of the matrix domain (left) conforming to the fracture network (right).

5 Figure 15 :

 515 Figure 15: Concentration in the matrix domain and in the fracture network obtained at different times for the tetrahedral mesh. A threshold of 0.2 is used in the matrix domain.

Figure 16 :

 16 Figure 16: Total computation time (left) and computation time for transport model (right) vs. number of MPI processes with tetrahedral mesh.

  ), D = 0.5 dx 2 + d 2 f .The production lasts 8 years.

Figure 17 :

 17 Figure 17: Prismatic mesh of the domain Ω defined by the tensor product of a vertical 1D uniform mesh with a 2D triangular mesh.

Figure 19

 19 Figure19plots the mean tracer concentration in each well as a function of time as well as the total volume of tracer as a function of time in the matrix, in the fracture network and their sum. Figure20exhibit the pressure solution in the matrix domain and Figures21 and 22shows the tracer concentration after one year of injection and at final time both in the matrix domain and in the fracture network. Figure23shows the total computation times with different number of MPI processes N p = 16, 32, 64, 128, 256, 512. It is observed that the total computation time exhibits a rather good scalability. In addition, the linear solver (GMRES+Boomer AMG) for the pressure converges in no more than 25 iterations whatever the number of MPI processes. Also the comparison of the total and transport computation times in Figure23shows that the time for the pressure solution remains small compared with the transport computation time up to N p = 512.

Figure 20

 20 Figure19plots the mean tracer concentration in each well as a function of time as well as the total volume of tracer as a function of time in the matrix, in the fracture network and their sum. Figure20exhibit the pressure solution in the matrix domain and Figures21 and 22shows the tracer concentration after one year of injection and at final time both in the matrix domain and in the fracture network. Figure23shows the total computation times with different number of MPI processes N p = 16, 32, 64, 128, 256, 512. It is observed that the total computation time exhibits a rather good scalability. In addition, the linear solver (GMRES+Boomer AMG) for the pressure converges in no more than 25 iterations whatever the number of MPI processes. Also the comparison of the total and transport computation times in Figure23shows that the time for the pressure solution remains small compared with the transport computation time up to N p = 512.

  Figure19plots the mean tracer concentration in each well as a function of time as well as the total volume of tracer as a function of time in the matrix, in the fracture network and their sum. Figure20exhibit the pressure solution in the matrix domain and Figures21 and 22shows the tracer concentration after one year of injection and at final time both in the matrix domain and in the fracture network. Figure23shows the total computation times with different number of MPI processes N p = 16, 32, 64, 128, 256, 512. It is observed that the total computation time exhibits a rather good scalability. In addition, the linear solver (GMRES+Boomer AMG) for the pressure converges in no more than 25 iterations whatever the number of MPI processes. Also the comparison of the total and transport computation times in Figure23shows that the time for the pressure solution remains small compared with the transport computation time up to N p = 512.

Figure 18 :

 18 Figure 18: Fracture network showing the location of the single injection well and of the two production wells.

Figure 19 :

 19 Figure 19: Mean tracer concentration in both production wells as a function of time (left) defined as the ratio between the well tracer flow rate and the well fluid flow rate (equal in our case to the well fracture-face tracer concentration). Volume of tracer as a function of time in the matrix domain, in the fracture network and their sum (right).

Figure 20 :

 20 Figure 20: Pressure on the matrix domain.

Figure 21 :

 21 Figure 21: Tracer concentration after one year of injection in the matrix domain (left) and in the fracture network (right).

Figure 22 :

 22 Figure 22: Tracer concentration at final time in the matrix domain (left) and in the fracture network (right).

Figure 23 :

 23 Figure 23: Total computation time in hours (left) and computation time for transport model (right) vs. number of MPI processes with prismatic mesh.

Table 1 :

 1 Number of linear solver iterations vs. the number of MPI processes obtained with different linear solvers and preconditioners for the mesh size n x = 128.

	N p	2	4	8	16	32	64 128 256 512
	GMRES + Boomer AMG	15	15	15	15	15	16	15	15	15
	GMRES + Aggregation AMG 59	78	65	39	65	54	73	53	62
	GMRES + ILU(0)	751 707 655 644 648 634 633 624 613
	BiCGStab + Boomer AMG	9	9	9	9	9	10	9	9	10
	BiCGStab + ILU(0)	508 476 484 503 473 513 491 487 484

Table 2 :

 2 Linear solver setup phase and solve phase computation times vs. the number of MPI processes obtained with different linear solvers and preconditioners for the mesh size n x = 128.

	N p		2	4	8	16	32	64	128 256 512
	GMRES	Setup 34.1	20.1	16.3	11.7 11.3 11.9 12.1 19.2 29.3
	Boomer AMG	Solve 26.3	15.6	14.8	7.2	5.2	3.8	2.5	5.2	9.6
	GMRES	Setup	4.7	1.9	1.6	1.5	2.3	2.9	4.4	6.6 11.3
	Aggregation AMG	Solve 45.1	20.9	17.0	9.7	5.2	2.5	2.3	1.5	3.2
	GMRES	Setup 16.9	21.3	16.3	23.2 14.6 11.0 9.7	6.0	4.8
	ILU(0)	Solve 672.3 590.9 281.6 163.9 71.4 30.7 16.7 8.3	4.0
	BiCGStab	Setup 38.0	23.3	15.0	10.3	9.1	9.4 12.8 14.8 23.8
	Boomer AMG	Solve 37.1	21.3	11.5	7.4	4.1	2.9	2.5	4.4 10.0
	BiCGStab	Setup 18.9	19.9	16.5	22.1 14.3 12.4 9.4	5.8	3.9
	ILU(0)	Solve 179.4 111.7 86.0	59.9 27.9 15.4 8.0	4.2	2.2

Table 3 :

 3 Number of linear solver iterations vs. the matrix fracture permeability ratioΛm for n x = 128 and N p = 2, 128.

	Λ f

Table 4 :

 4 Number of linear solver iterations vs. the number of MPI processes obtained with different linear solvers and preconditioners for the tetrahedral mesh.

	N p	2	4	8	16	32	64 128 256 512
	GMRES + Boomer AMG	11	12	12	12	12	12	12	12	12
	GMRES + Aggregation AMG	38	78	40	39	52	-	35	-	52
	GMRES + ILU(0)	1003 725 717 682 667 656 644 629 612
	BiCGStab + Boomer AMG	8	7	8	8	8	8	8	8	8
	BiCGStab + ILU(0)	565 513 527 544 535 483 489 483 473
	-: The relative residual norm stagnates after a few iterations.
	-: Some future investigations are necessary.

Table 5 :

 5 Linear solver setup phase and solve phase computation times vs. the number of MPI processes obtained with different linear solvers and preconditioners for the tetrahedral mesh.

	N p		2	4	8	16	32	64	128 256 512
	GMRES	Setup time 12.4	7.8	4.9	5.0	4.3	6.2 7.2 13.5 22.4
	Boomer AMG	Solve time	8.0	5.5	2.9	1.7	1.1	0.9 1.4 3.1	6.9
	GMRES	Setup time	3.7	1.9	1.2	1.8	2.1	1.6 2.9 3.3	4.7
	Aggregation AMG	Solve time 19.7	20.9	5.1	2.7	2.0	-	1.5	-	3.0
	GMRES	Setup time	5.7	7.4	7.2	5.6	4.7	5.2 3.4 2.8	1.8
	ILU(0)	Solve time 560.6 254.4 150.0 66.5 30.1 15.2 7.7 4.1	2.8
	BiCGStab	Setup time 21.5	14.5	9.9	6.5	5.3	5.9 8.2 12.4 19.7
	Boomer AMG	Solve time 24.0	10.2	6.1	3.5	1.8	1.5 2.1 4.3	9.5
	BiCGStab	Setup time	5.8	6.4	6.4	5.4	4.7	5.0 3.4 2.6	1.8
	ILU(0)	Solve time 110.4 63.0	39.0 19.2 11.6 5.4 2.8 1.4	1.2

-: The residual norm stagnates after a few iterations.

Table 6 :

 6 Number of linear solver iterations vs. the matrix fracture permeability ratio

	Λ f Λm for

-:

The solver doesn't converge in 1200 iterations.

Acknowledgments

This work is supported by a joint project between INRIA and BRGM Carnot institutes (ANR, INRIA, BRGM). This work was also granted access to the HPC and visualization resources of "Centre de Calcul Interactif" hosted by University Nice Sophia Antipolis. This work was partially funded by ADEME BRGM research project Orbou (convention 1305C0131). We would like also to thank Mohammad Karimi-Fard and Alexandre Lapène for kindly providing us the mesh of Section 5.5.

Numerical experiments

All the numerical tests have been implemented on the Cicada cluster of the University Nice Sophia Antipolis consisting of 72 nodes (16 cores/node, Intel Sandy Bridge E5-2670, 64GB/node). We always fix 1 core per process and 16 processes per node. The communications are handled by OpenMPI 1.8.2 (GCC 4.9) and PETSc 3.5.3.

The first two test cases are designed in order to validate the Darcy fluxes and the convergence of the transport model discretization on two analytical solutions including one fracture for the first test case and four intersecting fractures for the second test case. In the remaining test cases, the parallel scalability of our Darcy flow and transport solvers is assessed with different types of fracture networks and meshes and different matrix fracture permeability ratios. In particular, the last test case applies our algorithm to a complex fracture network with hundreds of fractures.

Numerical convergence for an analytical solution with one fracture

Let us set Ω = (0, 1) 2 , and denote by (x, y) the Cartesian coordinates of x. We then define

x tan(θ)}, and Ω 2 = Ω \ Ω 1 . We consider a single fracture defined by Γ = (x 1 , x 2 ) = ∂Ω 1 ∩ ∂Ω 2 with tangential permeability Λ f > 0, and width d f > 0. The matrix permeability is isotropic and set to Λ m = 1, the matrix and fracture porosities are set to φ m = φ f = 1, and the fluid viscosity is set to µ = 1. The pressure solution is fixed to u(x, y) = 1 -x. In this case, the transport model (4) reduces to the following system of equations which specifies our choice of the boundary and initial conditions:

where

It is assumed that k > 1. This system can be integrated along the characteristics of the matrix and fracture velocity fields leading to the following analytical solution:

100 200 400 800 1600

Numerical convergence for an analytical solution with four intersecting fractures

Let Ω = (0, 1)

4 -tan(θ 2 ), 1), and the intersection of x 1 x 2 and x 3 x 4 equal to

It is assumed that θ 1 , θ 2 ∈ (0, arctan( 3 4 )). We consider the four fractures Γ 1 = (x 1 , x 0 ), Γ 2 = (x 0 , x 2 ), Γ 3 = (x 3 , x 0 ), Γ 4 = (x 0 , x 4 ), with tangential permeabilities Λ f,1 = Λ f,2 > 0, and Λ f,3 = Λ f,4 > 0, and with widths d f,1 = d f,2 > 0, and

The fractures partition the domain Ω in the following four subdomains

. It is assumed that k 1 > 1 and k 2 tan(θ 2 ) > 1. The matrix and fracture porosities are set to φ m = φ f = 1 and the fluid viscosity is set to µ = 1.

The pressure solution is set to u(x, y) = 1 -x. In that case, the transport model (4) reduces to the following system of equations which specifies our choice of the boundary and initial conditions: find c m,α (x, y, t), α = 1, • • • , 4, c f,1 (x, t), c f,2 (x, t), c f,3 (y, t), c f,4 (y, t), and c 0 (t) such that 1 4 + x tan(θ 1 ), t) = 0 on (x 0 , 1) × (0, T ), L 2 c f,3 (y, t) -β 2 c m,3 ( 3 4 -y tan(θ 2 ), y, t) = 0 on (0, y 0 ) × (0, T ), L 2 c f,4 (y, t) -β 2 c m,1 ( 3 4 -y tan(θ 2 ), y, t) = 0 on (y 0 , 1) × (0, T ), c f,2 (x 0 , t) = c f,3 (y 0 , t) = c 0 (t) on (0, T ), (r + 1)c 0 (t) -c f,1 (x 0 , t) -rc f,4 (y 0 , t) = 0 on (0, T ), c f,1 (0, t) = c f,4 (1, t) = 1 on (0, T ), c f,1 (x, 0) = 0 on (0, x 0 ), c f,2 (x, 0) = 0 on (x 0 , 1), c f,3 (y, 0) = 0 on (0, y 0 ), c f,4 (y, 0) = 0 on (y 0 , 1),