Zexuan Xu
email: zexuanxu@lbl.gov

Nicolas Massei

Ingrid Padilla

Andrew Hartmann

Bill Hu

Characterization, modeling, and remediation of karst in a changing environment

Keywords: Karst, Data, Modeling, Changing Climate

This introductory editorial paper provides a review and prospective outlook of the achievements and challenges in karst research under a changing environment. A brief discussion of the past and future karst research has been focused on: (1) data and new technologies; (2) modeling of karst flow and reactive transport; (3) responses of karst hydrosystems to climate variability and changes across scales.

Introduction

Karst is a term used to describe a special landscape con-taining caves and extensive underground water systems that is developed by dissolution in the carbonate rocks such as limestone and gypsum [START_REF] Ford | Karst hydrogeology and geomorphol-ogy[END_REF]. Karst terrains undergo considerable dissolution to create and enlarge joints, fractures, bedding planes, and other open-ings in which groundwater flows. These terrains show dis-tinctive surface and subsurface features associated with sinkholes, springs, caves, and sinking, losing, and gaining streams, which are characterized by welldeveloped conduit networks having high connectivity (e.g., Xu et al. 2015a, b) and/or highly transmissive zones (e.g., [START_REF] Kuniansky | Simulating groundwater flow in karst aquifers with distributed parameter models-comparison of porous-equivalent media and hybrid flow approaches[END_REF]. These characteristics make some karst groundwater sys-tems highly productive and become important freshwater resources for human consumption and ecological integrity of streams, wetlands, and coastal zones [START_REF] Padilla | Fate, transport, and exposure of emerg-ing and legacy contaminants in karst system:state of knowledge and uncertainty[END_REF]White 2018). Karst regions underlie approximately 10% of the ice-free continental area, and karst aquifers par-tially supply drinking water resources to almost a quarter of the world's population [START_REF] Ford | Karst hydrogeology and geomorphol-ogy[END_REF]. Stress on water resources in karst aquifers has increased significantly in recent decades, associated with the limited quantity of water supply and deterioration of water quality [START_REF] Hartmann | Karst water resources in a changing world: review of hydrologi-cal modeling approaches[END_REF]. Available water resources have become the major constraint for many karst regions, for example, in Southwest China [START_REF] Lu | Karst water resources and geo-ecology in typical regions of China[END_REF], where management of land, people and water resources is a critical environmental chal-lenge in rocky desertification areas [START_REF] Jiang | Rocky desertification in Southwest China: impacts, causes, and restoration[END_REF]). In addition, many geohazards in karst areas, for instance, sinkhole collapse and land subsidence, are strongly related to over-utilization of groundwater, inappropriate limestone mining, and other human factors [START_REF] Kuniansky | The current status of mapping karst areas and availability of public sinkhole-risk resources in karst terrains of the United States[END_REF].

The characteristics of large connectivity and void spaces, which make karst aquifer systems highly productive, also make them highly vulnerable to contamination [START_REF] Göppert | Solute and colloid transport in karst conduits under low-and high-flow conditions[END_REF]. Surface karst features and large per-meability in the subsurface enhance the ability to transport contaminants over long distances, particularly through well-developed conduit networks. The porous matrix of the karst rocks and the significant amount of sediments trapped in karst formations provide high storage capacity for contami-nants that can be slowly released for long periods of time. As a result, karst systems can serve as an important route for contaminant exposure to humans and ecosystems [START_REF] Padilla | Historical contamination of groundwater resources in the north coast karst aquifers of Puerto Rico, dimensión[END_REF]. Groundwater quality impacted from point and nonpoint sources of contamination have been widely stud-ied in karst areas (e.g., [START_REF] Padilla | Fate, transport, and exposure of emerg-ing and legacy contaminants in karst system:state of knowledge and uncertainty[END_REF]Xu et al. 2015b;[START_REF] Yu | Spatiotemporal changes of CVOC concentrations in karst aquifers: analysis of three decades of data from Puerto Rico[END_REF].

Karst systems are highly dynamic and heterogeneous, with complex hydrological, geochemical, and biological processes that affect the transport of water, solutes, and non-aqueous entities over a wide range of spatial and temporal scales. Fundamental knowledge exists on many of these pro-cesses and several predictive and characterization models have been developed and applied to karst systems. Tremen-dous challenges and uncertainties are faced when trying to characterize, predict exposure, implement remedial actions, and manage contaminated systems, particularly in a chang-ing environment, including climate change, human impacts, etc. A need exists to advance knowledge and technology for understanding the hydrological response of karst systems to those changes, and to develop appropriate management and remediation approaches.

Interdisciplinary studies, including collaborations and expertise from a full range of scientific and engineering disciplines, are in high demand for assessment of water resources availability, chemical dynamics, and potential con-taminant exposure in a karst aquifer. For example, advanced monitoring sensors and instruments for physical and chemi-cal hydrology can enhance data availability in field geol-ogy, particularly in deep subsurface karst conduit networks. Interactions between biogeosciences, atmospheric sciences, and marine and oceanographic sciences help fully address earth system processes in karst studies. Statistical, math-ematical, scientific computing, and data sciences knowledge contributes to build the state-of-the-art physical-based and/ or data-driven modeling approaches to the study of karst hydrogeology. Collaborative works between engineers, hydro-scientists, and public health experts are mapping the way to assess direct links between contamination in karst regions and adverse health outcomes [START_REF] Cordero | Team science applied to environmental health research: karst hydrogeology and preterm birth in Puerto Rico[END_REF].

Major achievements and challenges in karst research

Papers in this thematic issue help to summarize the knowl-edge learned in the past, to review the progress of ongoing research, and to predict major challenges in the future. This section is organized into three perspectives: (1) data and new technologies; (2) modeling of karst flow and reactive transport; (3) responses of karst hydrosystems to climate variability and changes across scales.

Data and new technologies

Limited laboratory and field observational data have been a major constraint in karst research. Laboratory experiments have been extensively developed to study the groundwater and solute transport (e.g., [START_REF] Faulkner | Laboratory analog and numerical study of groundwater flow and solute transport in a karst aquifer with conduit and matrix domains[END_REF][START_REF] Anaya | Estimating preferential flow in karstic aquifers using statistical mixed models[END_REF], and reactive transport (e.g., [START_REF] Geyer | Multitracer test approach to characterize reactive transport in karst aqui-fers[END_REF], in the interaction among karst conduits, fractures, and porous medium. Unlike the alluvial sandstone aquifers, karst fea-tures, including sinkholes, karst windows and conduit networks, are difficult to be successfully explored by traditional hydrogeological field methods. For example, the subsurface karst conduit networks are found at nearly 100 m below land surface in the Woodville Karst Plain, North Florida [START_REF] Werner | Determination of groundwater flow patterns from cave exploration in the Woodville Karst Plain, Florida[END_REF], which led to difficulties in field investigation. Tracer tests have been widely applied in the studies of karst aquifers to resolve the active and often inaccessible part of cave and conduit networks, and to evaluate the hydraulic connectivity in karst aquifers [START_REF] Goldscheider | Tracer tests in karst hydrogeology and speleology[END_REF][START_REF] Werner | Determination of groundwater flow patterns from cave exploration in the Woodville Karst Plain, Florida[END_REF]. Geochemistry data were widely collected to study the global carbon cycle, including carbon flux data for biogeochemical cycle [START_REF] Gill | Belowground carbon flux links biogeochemi-cal cycles and resource-use efficiency at the global scale[END_REF], aquatic carbonate concentration for the contribution of carbonate weathering to atmospheric CO 2 [START_REF] Huang | Inorganic carbon flux and its source in the karst catchment of Maocun, Guilin, China[END_REF], and karst-related carbon sink [START_REF] Zeng | Hydrologically-driven vari-ations in the karst-related carbon sink fluxes: insights from high-resolution monitoring of three karst catchments in South-west China[END_REF], from watershed to global scales.

Utilization of advanced instruments and technologies is critical to provide more data and detailed descriptions to gain insights into karst research, particularly in the deep subsurface conduit networks considering the difficulties and complexities of karst exploration. [START_REF] Xu | Long distance seawater intrusion through a karst conduit network in the Woodville Karst Plain[END_REF] reported a flow speed, flow direction, and electrical conductivity moni-toring work over a long-term period, with the help of cave divers at the 100-m deep conduit in the Woodville Karst Plain, Florida. This study discovered potential seawater intrusion through subsurface conduit networks in the coastal karst aquifer, which indicated a trend of rising sea levels, and analyzed its impacts on groundwater contamination. It is also crucial to capitalize on long-term data, where available, and accordingly to develop and maintain long-term observa-tories of karst hydrosystems. On the other hand, investiga-tion of conduit network geospatial distribution and physical geometry is in high demand in karst research by deploy-ing advanced and new monitoring techniques. In recent years, [START_REF] Kenney | Neutrally buoyant sensor apparatus and method for mapping a water pathway[END_REF] developed a neutrally buoyant sensor method for mapping a water pathway in the direc-tion of subsurface conduit flow. In [START_REF] Weidner | Underwater cave mapping using stereo vision[END_REF], a systematic approach for the 3-D mapping of underwater caves utilized the intersection of video light with the cave boundaries to describe the outline of the cave. Further works are required to apply these newly developed techniques in field application, improve the data resolution and quality, and transform measurements to actionable information.

Modeling of karst flow and reactive transport

Modeling karst groundwater flow and contaminant trans-port in karst systems pose particular challenges because of their highly heterogeneous nature, wide spectrum of flow regimes (ranging from laminar to turbulent flow), lack of proper data for accurate characterization, and limited applicable models with full physical processes, among others. Numerical models are often used to describe flow and trans-port trends, as long as they are applied within the limits of their realm. Several modeling approaches exist for simu-lating groundwater flow and contaminant transport in karst aquifers, for example, the lumped and distributed models of varying complexities [START_REF] Ghasemizadeh | Review: groundwater flow and transport modeling of karst aquifers, with particular reference to the North Coast Limestone aquifer sys-tem of Puerto Rico[END_REF]Hart-mann et al. 2014;[START_REF] Scanlon | Can we simulate regional groundwater flow in a karst system using equivalent porous media models? Case study, Barton Springs Edwards aquifer[END_REF]. Dual-permeability hydrogeological system is one of the major challenges in simulating groundwater flow in a karst aquifer, which con-sists of both high-permeability conduits and fractures, and low-permeability surrounding porous media. Different from the porous media in the alluvial aquifer, groundwater flow within the highpermeability conduits and fractures can eas-ily become non-Darcian and even beyond the upper criti-cal Reynolds numbers. The non-Darcian conduit flow and Darcian porous medium flow are computed and coupled in the discrete-continuum hybrid models [START_REF] Liedl | Simulation of the development of karst aquifers using a cou-pled continuum pipe flow model[END_REF][START_REF] Shoemaker | Documentation of a conduit flow process (CFP) for MOD-FLOW-2005[END_REF], which have been widely used to simulate variable-saturation flow conditions [START_REF] Reimann | Effects of dynamically variable saturation and matrix-conduit coupling of flow in karst aquifers[END_REF]) and groundwater contamination transport issues (Xu et al. 2015a) in karst aquifers. The recharge reservoir model is developed and coupled with discrete-continuum models to simulate surface and subsurface flow in karst systems under various rainfall conditions [START_REF] Chang | Modeling the hydrological behavior of a karst spring using a nonlinear reservoir-pipe model[END_REF][START_REF] Chen | A new approach to evalu-ate spatiotemporal dynamics of controlling parameters in distrib-uted environmental models[END_REF]. A densitydependent discrete-continuum model has been developed in [START_REF] Xu | Development of a discrete-continuum VDFST-CFP numerical model for simulating seawater intrusion to a coastal karst aquifer with a conduit system[END_REF] as a new mod-eling method to study seawater intrusion through conduit networks in a karst aquifer, providing a way to simulate the impacts of sea level rise under changing climate. Future work of more accurate parameterizations is highly desired and deterministic in karst modeling applications, particularly the exchange flux between conduit and porous medium, the dispersion coefficient, and friction factor within the con-duit. In addition, high-performance computing and advanced numerical methods are desired to solve the karst hydrology models more efficiently and accurately, with fully coupled physical-based processes in a finer resolution.

At regions with limited data availability, simplified karst modeling approaches are applied that represent the domi-nant dual permeability processes of karst systems by varying combination of buckets that represent the important stor-ages of the karst system: the soil, epikarst, karst matrix, and conduits [START_REF] Hartmann | Karst water resources in a changing world: review of hydrologi-cal modeling approaches[END_REF]. In recent years, research focused on quantifying the uncertainty that results from those simplified karst modeling approaches [START_REF] Chen | A new approach to evalu-ate spatiotemporal dynamics of controlling parameters in distrib-uted environmental models[END_REF][START_REF] Adinehvand | A step-wise semi-dis-tributed simulation approach to characterize a karst aquifer and to support dam construction in a data-scarce environment[END_REF][START_REF] Xu | Numerical modeling and sensitivity analy-sis of seawater intrusion in a dual-permeability coastal karst aqui-fer with conduit networks[END_REF] and quantify-ing the information content of auxiliary data as hydrochemi-cal information (Hartmann et al. 2017a, b) or geophysical measurements [START_REF] Mazzilli | On the inclusion of ground-based gravity measurements to the calibration process of a global rainfalldischarge reservoir model: case of the Durzon karst sys-tem (Larzac, southern France)[END_REF]). Future research is nec-essary to further explore the value of auxiliary data and to link the advantages of these simplified modeling approaches and the above-mentioned physically based approaches.

The subsurface conduit network is formed by a positive feedback mechanism between dissolution by water under-saturated with respect to calcite and the amount of water percolating through a carbonate rock [START_REF] Worthington | Self-organized permeability in car-bonate aquifers[END_REF]. Physical-based reactive transport modeling (RTM) in hydrology is important to investigate the evolution of frac-tures and conduits, and to compute carbon mass balance in multiple temporal and global spatial scales. For example, [START_REF] Rehrl | Conduit evolution in deep-seated settings: conceptual and numerical models based on field observations[END_REF] developed a discrete-continuum RTM to couple aquatic chemical reaction processes with ground-water flow and transport models to examine the interrelation between the hydrogeological environment and conduit development in a deep-seated setting. In a recent study, de Rooij and Graham (2017) presented a hydrochemical model that can be used to generate plausible karstic conduit networks that honor what is known about geology, hydrology, and topography of a karst system. Conduit evolution modeling, however, is still an oversimplification of reactive transport processes that requires further investigations in the limited number of chemistry processes, coarse spatial resolution, and lack of validation. On the other hand, sto-chastic and statistical modeling in karst conduit evolutions have been widely applied to help build the conduit networks in hydrologic numerical modeling [START_REF] Jaquet | Sto-chastic discrete model of karstic networks[END_REF]Pardo-lguzquiza et al. 2012). Prior information, sampling methods, computational cost, and lack of real conduit network valida-tion data are still unresolved issues and open questions for stochastic simulations of karst networks.

Response of karst hydrosystems to climate variability and change across scales

To understand the impacts of global climate change on water resources in karst aquifers, the scalability of hydro-logic models becomes critical, from pore scale to regional scale. Owing to the very heterogeneous nature of the karst, such hydrosystems may display particular hydrological behaviors that need to be thoroughly analyzed and modeled. The interaction between watershed-scale hydrologic models and regional, even global-scale climate models has been investigated in a number of studies, such as the US National Water Model (currently under development, http:// water.noaa.gov/about/nwm) and European Flood Awareness System [START_REF] Kauffeldt | Technical review of large-scale hydrological models for imple-mentation in operational flood forecasting schemes on continental level[END_REF]. None of the above, how-ever, address the complexity of karst aquifers, or have been applied in either karst watersheds or on a large regional scale of karst land mass. Some recent work showed that the reli-ability of presently applied large-scale simulation tools is strongly limited [START_REF] Hartmann | A large-scale simulation model to assess karstic ground-water recharge over Europe and the Mediterranean[END_REF](Hartmann et al. , 2017)). In addition, global-scale biogeochemical simulations are essential to evaluate the roles of karst regions in global carbon cycles, particular carbon sink effects in carbonate rock dissolution and weathering. The upscaling of current hydrology models in karst aquifers from the dualpermeability watershed scale to global scale becomes important to carbon budget estima-tion. The hydrologic and aquatic geochemistry processes in a karst aquifer, including interaction of non-linear and linear flow, and reactive transport in carbonate rock dissolution and precipitation processes, have potential impacts on global carbon cycle and should be included in future earth system modeling studies.

On the other hand, most karst distributed models are still developed based on some traditional hydrology models that only solve the Darcy equation in the aquifer, for example, MODFLOW-CFP [START_REF] Shoemaker | Documentation of a conduit flow process (CFP) for MOD-FLOW-2005[END_REF], with very sim-plified land surface processes. The hydrological interactions between atmosphere and land have been addressed in several integrated hydrology models, including PF.WRF [START_REF] Maxwell | Development of a coupled groundwater-atmosphere model[END_REF]) and WRF-Hydro [START_REF] Gochis | The WRF-Hydro model technical description and user's guide, version 1.0[END_REF]; however, they have rarely been integrated with dual-permeability hydrology models and applied in karst regions. These inte-grated hydrology models are advanced with physical-based land surface models to study the responses of increasing temperature, dynamic vegetation, precipitation variability, and human factors. For example, early snowmelt processes under changing climate control groundwater table fluctua-tion and then affect carbonate dissolution and water chemis-try [START_REF] Winnick | Snowmelt controls on concentration-discharge relationships and the balance of oxidative and acid-base weath-ering fluxes in an alpine catchment, East River, Colorado[END_REF]. Integration of land surface models and dual-permeability hydrology models becomes the key for a successful simulation of hydrological and biogeochem-ical responses in karst aquifers under climate change.

Preview of the papers expected to be accepted

The scientific community of earth science has become even more aware of the significant and grand challenges fac-ing water resources and ecosystems in karst aquifers. An increasing number of karst-related research studies have been initiated and developed worldwide in recent years. For example, the number of abstracts submitted to the karst session in AGU has increased in the past 3 years (22 in 2015, 28 in 2016 and 37 in 2017). Among these, a significant number of studies have been cross-listed with biogeoscience, global environmental change, atmospheric, and oceanography.

In response to the increasing amount of karst research, this thematic issue on "Characterization, Modeling, and Remediation of Karst in a Changing Environment (AGU)" with Environmental Earth Sciences (EES) will be composed of papers on research related to the characterization, moni-toring, and modeling of hydrogeology, morphology, biogeochemistry, geohazards, and remediation of carbonate and karst water systems at multiple temporal and spatial scales.


Luc Brémaud

Xiran Cai
email: xirancaicaixr@shanghaitech.edu.cn

Renald Brenner
email: renald.brenner@sorbonne-universite.fr

Quentin Grimal
email: quentin.grimal@sorbonne-universite.fr

Maximum effect of the heterogeneity of tissue mineralization on the effective cortical bone elastic properties

Keywords: bone mineralization, porosity, elasticity, homogenization, mesoscale

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

About 10% of the adult skeleton is renewed every year by the remodeling process which, in cortical bone, produces secondary osteons. Bone cells deposit a collagen matrix which is progressively mineralized resulting in an increase of the average degree of mineralization of bone (DMB) [START_REF] Boivin | Changes in bone remodeling rate influence the degree of mineralization of bone[END_REF]. As a consequence of this continuous remodeling process, cortical bone mineralized extracellular matrix (or shortly, matrix) is an assembly of portions of tissue of various ages characterized by different mineralization levels. The oldest part of this tissue is found in the interstitial tissue, which is composed of remnants of old osteons which have been partially remodeled. The difference of DMB between intact osteons and interstitial tissue is directly observed in microradiographs or indireclty as differences of elasticity with nanoindentation [START_REF] Zysset | Elastic modulus and hardness of cortical and trabecular bone lamellae measured by nanoindentation in the human femur[END_REF], or acoustic impedance (a proxy for tissue elasticity) with scanning acoustic microscopy [START_REF] Raum | Microelastic imaging of bone[END_REF] as illustrated in Fig. 1. Typically, osteons with the lowest mineral density contain 70% to 75% of the mineral content of the most highly calcified tissue parts [START_REF] Boivin | Changes in bone remodeling rate influence the degree of mineralization of bone[END_REF]. For instance, Lefevre et al. found that DMB measured from quantitative microradiography in elderly adults was 1.260 (± 0.029) g cm -3 in interstitial tissue vs. 1.079 (± 0.037) g cm -3 in osteonal tissue [START_REF] Lefèvre | Compositional and mechanical properties of growing cortical bone tissue: a study of the human fibula[END_REF].

Fig. 1: Scanning acoustic microscopy (SAM) maps of acoustic impedance (MRayl) of human femoral cortical bone in a plane transverse to osteons. The relative differences in impedance reflect to a large extent the variations of mineralization levels in different portions of the tissue, a larger impedance (or elasticity) corresponding to more mineralized, and older, tissue [START_REF] Raum | Derivation of elastic stiffness from site-matched mineral density and acoustic impedance maps[END_REF]. The main image is a region of 5 × 5mm 2 scanned at 50 MHz with a resolution of about 30µm. The inset is a SAM image at 200 MHz with a resolution of about 8µm [START_REF] Granke | Change in porosity is the major determinant of the variation of cortical bone elasticity at the millimeter scale in aged women[END_REF]. The vascular porosity appears in dark blue.

The effective elastic properties of cortical bone at the mesoscale (Grimal et al. 2011a), or millimeter scale, are determined by the elastic properties of the mineralized tissue (matrix) and the properties of the pore network. Micromechanical models have been used for several decades to calculate the effective elastic properties for given representative volume elements (RVE) of bone. These models have been extensively used to investigate the relationships between material composition, microarchitecture and elastic properties. The heterogeneity of the elastic properties of the matrix, i.e., the differences between osteonal and interstitial tissues, has been considered in some studies [START_REF] Crolet | Compact bone: numerical simulation of mechanical characteristics[END_REF][START_REF] Dong | Prediction of cortical bone elastic constants by a two-level micromechanical model using a generalized self-consistent method[END_REF][START_REF] Grimal | Derivation of the mesoscopic elasticity tensor of cortical bone from quantitative impedance images at the micron scale[END_REF][START_REF] Hamed | Multiscale modeling of elastic properties of cortical bone[END_REF]). However, many other authors have considered an homogeneous matrix [START_REF] Hellmich | can the diverse elastic properties of trabecular and cortical bone be attributed to only a few tissue-independant phase properties and their interactions[END_REF][START_REF] Sevostianov | Impact of the porous microstructure on the overall elastic properties of the osteonal cortical bone[END_REF][START_REF] Parnell | The influence of mesoscale porosity on cortical bone anisotropy. investigations via asymptotic homogenization[END_REF][START_REF] Baumann | The relative influence of apatite crystal orientations and intracortical porosity on the elastic anisotropy of human cortical bone[END_REF][START_REF] Gagliardi | Tissue mineral density measured at the sub-millimetre scale can provide reliable statistics of elastic properties of bone matrix[END_REF]Cai et al. 2019a). The errors on predicted mesoscale properties when the differences in elasticity between osteonal and interstitial tissues are disregarded have not been quantified. Noteworthy, these errors may strongly depend on the underlying microstructure (porosity and pore organization).

The objective of the present paper is to quantify the range of errors on predicted mesoscale effective elastic properties due to the assumption of an homogeneous mineralized matrix. We compare the effective elasticity of a RVE assuming (1) different elastic properties in osteonal and interstitial tissue and (2) average elastic properties. In order to account for the variability of bone microstructure, we use a collection of high resolution three-dimensionnal images of the pore network to build RVEs of cortical bone. For each RVE, we assume a constant osteonal wall thickness and we artificially vary this thickness between 35 and 140 µm to create RVEs with different amounts of osteonal tissue. With this parametric study, we estimate a maximum effect of the heterogeneity This work provides data to quantitatively discuss the relevancy of important assumptions of cortical bone models. The results are helpful to guide the interpretation of effective elastic properties at the millimeter scale. These effective properties are involved in finite element models for the prediction of bone strength (van Rietbergen and Ito 2015; [START_REF] Engelke | FEA to Measure Bone Strength: A Review[END_REF] or bone remodeling (Martínez-Reina and Pivonka 2019), and for the in vivo assessment of porosity by inverse approaches [START_REF] Minonzio | Ultrasound-Based Estimates of Cortical Bone Thickness and Porosity Are Associated With Nontraumatic Fractures in Postmenopausal Women: A Pilot Study[END_REF]

Model

Specimens and bone microstructure modelling

We have used a collection of synchrotron radiation micro-computed tomography (SR-µCT) 3D images of volumes of interest (VOI) from a previous study (Cai et al. 2019a) (Fig. 2). Briefly, cortical bone specimens were harvested from the mid-diaphysis of the left femur of 25 human cadavers. Among the donors, 13 were females and 12 were males 77 ± 11.5,mean±SD). The nominal size of the specimens was 3 × 4 × 5 mm 3 in radial (x 1 ), circumferential (x 2 ) and axial (x 3 ) directions of the bone, respectively.

The femurs were provided by the Départment Universitaire d'Anatomie Rockefeller (Lyon, France) through the French program on voluntary corpse donation to science. The tissue donors or their legal guardians provided informed written consent to give their tissue for investigations, in accord with legal clauses stated in the French Code of Public Health.

Images of the specimens were obtained with isotropic voxel size of 6.5 µm performed at the European Synchrotron Radiation Facility (ESRF, Grenoble, France), Figure 3. The 3D volume of each specimen was slightly rotated with the image processing software Fiji [START_REF] Schindelin | Fiji: An open-source platform for biological-image analysis[END_REF]) using bilinear interpolation so that Fig. 2: A representative reconstructed cross-sectional slice from the SR-µCT data (pixel size 6.5 µm) before segmentation.

the geometric coordinates coincide with the material coordinates defined by the faces of the specimen. As perpendicularity and parallelism errors were about 1 • [START_REF] Cai | Quantification of stiffness measurement errors in resonant ultrasound spectroscopy of human cortical bone[END_REF], the orientation correction was small. Thereafter, axis x 3 was approximately along the direction of osteons and axes x 1 and x 2 were perpendicular to osteons. The images were then binarized by simple thresholding to obtain two phases : pores and mineralized matrix (Cai et al. 2019a). In each specimen, a VOI of approximately 2.8 × 3.9 × 4.8 mm 3 corresponding to one representative volume element (RVE) for the homogenization was selected manually. The VOI size was much larger than the minimum RVE size for cortical bone (Grimal et al. 2011a). After a convergence study (Cai et al. 2019a), the voxel size was then increased to 35 µm to reduce the computational cost, which lead to about 1.2 million voxels per specimen.

We then built RVEs made of three phases: the pores and two phases within the mineralized matrix to model its heterogeneity, that is, an assembly of osteonal (O) and interstitial (I) tissues. Images were processed to obtain segmentation masks for the osteonal and interstitial phases based on the observation that the Haversian canals run roughly along an osteon's center line [START_REF] Maggiano | Three-dimensional reconstruction of Haversian systems in human cortical bone using synchrotron radiation-based micro-CT: Morphology and quantification of branching and transverse connections across age[END_REF]). The osteonal tissue was defined to be within a ring of a given thickness e around pores. Each slice of the 3D image stack, which was roughly perpendicular to the osteon axis, was processed independently. The ring was obtained by erosion with a disk as structuring element (Figures 3 and4).

Notation for the elasticity tensor

We use the classical two-indices Voigt notation for the stiffness tensor components C ij . For an orthotropic material with principal directions aligned with the frame axes (x 1 , x 2 , x 3 ), the stiffness tensor C has nine independent moduli : C ii (i = 1, 2, 3) denote the longitudinal stiffnesses, C ii (i = 4, 5, 6) denote the shear moduli and only three non-diagonal terms are different from zero: C 12 , C 13 and C 23 . For a transversely Fig. 3: Three-dimensional representative volume element (RVE) with vascular pores in black. Left: two-phase model where the matrix (white) is considered to be a mixture of interstitial and osteonal tissue. Right: three-phase model where the matrix is made of two phases (white: osteonal tissue; dark gray: interstitial tissue). The direction of the bone diaphysis axis (and approximate osteon axis) is given by x 3 .

Mineralized matrix elasticity model

Following [START_REF] Grimal | Derivation of the mesoscopic elasticity tensor of cortical bone from quantitative impedance images at the micron scale[END_REF]Grimal et al. ( , 2011b)), the mineralized matrix elastic properties are modeled in a transversely isotropic framework as functions of the mineral volume fraction f ha . We use the "mineral foam matrix with collagen inclusions" micromechanical Fig. 5: Elastic coefficients of the mineralized matrix as a function of mineral volume fraction f ha [START_REF] Grimal | Derivation of the mesoscopic elasticity tensor of cortical bone from quantitative impedance images at the micron scale[END_REF](Grimal et al. , 2011b)). model introduced by Hellmich et al. ( 2004) which is based on two idealizations: (1) at a length scale of 100 nm, hydroxyapatite crystals and ultrastructural water with non-collagenous organic material constitutes a mineral foam;

(2) at a length scale of 5 to 10 microns, collagen fibers are embedded into the mineral foam. Volume fractions of the collagen fibers and water-filled nanoscale pores are respectively denoted f col and f w . The model assumes fixed elastic properties of the constituents (mineral, collagen, and water) and uses self-consistent and Mori-Tanaka homogenization schemes to derive the effective properties. Since f ha + f col + f w = 1 the model has in fact only two independent parameters. Furthermore, using an empirical relationship between f col /f w and f ha [START_REF] Raum | Derivation of elastic stiffness from site-matched mineral density and acoustic impedance maps[END_REF][START_REF] Broz | Material and compositional properties of selectively demineralized cortical bone[END_REF], the mineralized matrix stiffness tensor C m is obtained as a function of a single parameter which we choose to be f ha . The evolution of matrix elastic coefficients as a function of f ha is plotted in Figure 5 for illustration. Different values of f ha are assumed for the osteonal and interstital tissues as described below.

Choice of RVE parameters

Each RVE is defined by a specific microstructure θ (i.e, a 3D image of the vascular porosity), the transversely isotropic elastic properties of osteonal and interstitial tissues as well as osteon thickness e. As the elastic properties of the matrix phases are functions of f ha , as explained above, the RVE effective elasticity tensor is a function C = f (θ, e, C m;O , C m;I ) with C m;O and C m;I the elastic stiffness tensors of osteonal and interstitial tissues, respectively.

Interstitial tissue is assumed to be fully mineralized, that is, C m;I = C m (f I ha = 0.43) and we consider a mineral volume fraction of osteonal tissue f O ha = 0.40,

C m;O = C m (f O ha = 0.4
). These constitutive assumptions correspond to a longitudinal stifness C m;O 33 which is 19.75% smaller than C m;I 33 (Figure 6). These values were selected based on the elasticity differences observed in nanoindentation measurements of osteonal and interstitial tissues. For instance, [START_REF] Rho | Elastic properties of human cortical and trabecular lamellar bone measured by nanoindentation[END_REF] measured dehydrated tibia specimens embedded in epoxy resin and reported Young's moduli measured along the osteon axis of 22.5 ± 1.3GPa and 25.8 ± 0.7GPa in osteonal and interstitial tissues respectively, that is, a Young's modulus 12.8% smaller in osteonal tissue. [START_REF] Zysset | Elastic modulus and hardness of cortical and trabecular bone lamellae measured by nanoindentation in the human femur[END_REF] reported Young's moduli measured along the osteon axis in diaphyseal wet femoral bone specimens of 19.1 ± 5.4GPa in osteonal and 21.2 ± 5.3GPa in interstitial For the purpose of comparison, we also calculated the stiffness tensor C m;hom of an equivalent homogeneous matrix. The elasticity tensor is obtained from the mineralized matrix elasticity model by considering the average mineral content f ha , that is

C m;hom = C m ( f ha ), f ha = ν I f I ha + ν O f O ha ν I + ν O (1)
where ν I and ν O are the volume fractions of interstitial and osteonal phase, respectively. Due to the quasi-linear variation of the elastic moduli in the considered range of mineral volume fraction (Figure 5), it almost coincides with the Voigt bound

C m ( f ha ) ≈ C m . ( 2
)
Osteon thickness e was chosen according to literature data reporting dimensions of the Haversian canals and of the osteons. [START_REF] Britz | The relation of femoral osteon geometry to age, sex, height and weight[END_REF] conducted a twodimensional histomorphometric analysis on microradiographs of femoral diaphyseal tissue from 88 donors (male and female 17 to 97 years old from the Melbourne Femur Research Collection (https://dental.unimelb.edu.au/research/melbourne-femurresearch-collection). They reported osteon diameter in the range [155 -325] µm with average of 220 ± 28 µm. [START_REF] Cooper | Age-dependent change in the 3D structure of cortical porosity at the human femoral midshaft[END_REF] analyzed the properties of the vascular porous network in three-dimensions from Micro-CT images of 79 donors (males and females 18 to 92 years old, also from the Melbourne Femur Collection). They reported pore diameters in the range [56.0 -456.3] µm with average of 117.2 ± 89.6 µm. [START_REF] Gauthier | 3D analysis of the osteonal and interstitial tissue in human radii cortical bone[END_REF] et al. conducted a SR-µCT study on samples from the radius of eight female donors (70.3 ± 13.7 y.o.). Osteonal canal diameter was (mean, SD) 64.7 ± 23.1µm and osteon diameter was 184.0 ± 13.3µm. Based on these mean values of osteon and pore diameters, the average thickness of the osteonal wall is between 50 and 60 µm with an expected large interval of variation which cannot be calculated from the experimental data of different sources.

We used e = 35, 70, 105 and 140µm. The smallest value corresponds to the typical (average) size of an osteon. The other values were chosen in order to model a maximum effect of the matrix heterogeneity on effective elastic properties at the mesoscale.

Pores were assumed to be filled with a fluid whose elastic properties are similar to those of water, that is a bulk modulus equal to 2.2 GPa and a null shear modulus. The parameters of the RVE are summarized in Table 1.

Overall elasticity inferred from numerical homogenization

To investigate the influence of the bone tissue heterogeneity on the effective elastic properties, we have performed unit-cell computations on the RVE samples with the microstructural assumptions previously described. The size of the constitutive heterogeneities (i.e pores and osteons) being much smaller than the sample volume, the Hill-Mandel macro-homogeneity condition is fulfilled. The overall properties can thus advantageously be determined by assuming periodic-boundary conditions on the microstructural unit-cell Ω. The local problem to be solved thus reads, ∀x ∈ Ω,

curl(curl T ε(x)) = 0, div σ(x) = 0, σ(x) = C(x) : ε(x),
(3) with periodicity conditions on the boundary ∂Ω [START_REF] Suquet | Homogenization Techniques for Composite Media[END_REF]. ε and σ are respectively the linearized strain and Cauchy stress. C is the tensor field of the local elastic moduli which values correspond to the elastic properties of the pores, the osteonal or the interestitial tissue (Table 1) depending on the coordinate. To solve the problem (3), we have used a fast Fourier transform (FFT) based numerical scheme [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF][START_REF] Moulinec | Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials[END_REF]. The principle of this method is to solve iteratively the implicit integral equation for the strain field ε(x)

     ε(x) = E + Ω Γ (0) (x -x ) : τ (x )dx , τ (x) = (C(x) -C (0) ) : ε(x) (4)
with E the macroscopic strain and Γ (0) denoting the strain Green operator corresponding to a reference homogeneous medium with elasticity C (0) . This numerical method, which is widely used in engineering material mechanics, has been recently used for cortical and trabecular bones [START_REF] Colabella | Calculation of cancellous bone elastic properties with the polarization-based fft iterative scheme[END_REF][START_REF] Gagliardi | Estimation of the effective boneelasticity tensor based on µCT imaging by a stochastic model. A multi-method validation[END_REF]Cai et al. 2019a[START_REF] Cai | Anisotropic elastic properties of human cortical bone tissue inferred from inverse homogenization and resonant ultrasound spectroscopy[END_REF]. It can be mentioned that it allows to perform calculations directly on the microstructure digital image. The effective elastic moduli tensor C is classically defined by

σ Ω = C : ε Ω (5)
where . Ω denotes a volume average over the unit-cell. No symmetry assumptions have been made on C and six independent loadings were performed, for each sample, to determine the whole set of overall elastic coefficients (21). The reader is referred to Cai et al. (2019a) for a detailed convergence analysis.

Results and discussion

The mineral content of the interstitial and osteonal tissues being chosen, the effective elasticity tensor C is solely a function of the specific pore network microstructure θ and the considered osteon thickness e. Unit-cell FFT computations have been performed for each parameter pairs (θ, e) by considering either an heterogeneous matrix (pores coated with osteons) or an equivalent homogeneous matrix. Numerical results are reported in Figure 7 for the 25 different pores network microstructures with porosities ranging from about 3 to 20%. The relative errors of elastic moduli (i.e., the relative difference between effective elastic moduli of heterogenous matrix mineralization and homogeneous matrix mineralization) as a function of porosity are shown in Figure 8 for the largest osteon thickness. The average error values are given in Table 2 for each osteon thickness.

The relative error between effective elasticity coefficients using a homogeneous versus a heterogeneous matrix depends on the porosity and the osteonal thickness. The error weakly depends on the elastic modulus and varies from ∼ 1% for e = 35 µm to ∼ 7% for e = 140 µm. For e = 140 µm, it can be observed that the difference is almost constant for a porosity between 5 and 10%. Nevertheless, fluctuations of the error of about 2% are observed between samples of similar porosity but different microstructure, evidencing an influence of the distribution of pores on the error. For lower or higher porosity values, the difference decreases. More bone samples with a high porosity would be necessary for a more quantitative analysis.

The equivalent homogeneous matrix elasticity (Eq. 1) is almost defined by the Voigt upper bound because the 'homogeneous' mineral content is calculated as the arithmetic mean. Accordingly, accounting for the matrix morphological heterogeneity in the homogenization is expected to lead to lower effective moduli, which is observed in our results. Furthermore, for a given microstructure, this softening effect consistently increases with the osteon thickness, except for high porosity values. At a larger scale, by modeling the global stiffness of a mouse femur from high-resolution CT, Blanchard et al. ( 2013) also observed that calculating apparent properties from the average mineral content overestimates the stiffness compared to a model incorporating the heterogeneity of the mineral distribution in the femur.

In our work, the calculated error should be considered as an upper bound because we have designed the study to model a maximum effect of the heterogeneity. We selected values of f ha for the osteonal and interstitial phases leading to differences in elastic properties in the upper range of the differences observed experimentally. We have also considered that all pores are surrounded by softer tissue (young osteon) and that the thickness of the osteons' wall around canals is constant inside each RVE. In reality many pores are surrounded by older tissue, reducing the difference in mechanical properties with interstitial tissue.

We used an idealized model of cortical bone that incorporates the main features of the heterogeneity and that was parametrized in order to observe a range of effects. More realistic and sample-specific models could be designed using calibrated SR-µCT data providing DMB for each voxel. The specific elastic heterogeneity could be reconstructed after converting DMB to elasticity for each voxel. A major drawback of this approach is that the law to convert DMB to a stiffness tensor at this scale is not known and likely dependent on the location in the tissue. Another sample-specific approach would be to perform a segmentation of osteons based on the differences of gray level between osteons and interstitial tissue [START_REF] Gauthier | 3D analysis of the osteonal and interstitial tissue in human radii cortical bone[END_REF]) and then allocate a different stiffness per tissue type. However, even with high quality acquisitions this would have required a prohibitive amount of work for the samples of the present study. Considering more realistic models of tissue heterogeneity would likely lead to errors smaller than the errors we calculated. Interestingly, the error is found to be the largest in the most common porosity interval for non-pathological bone (∼5-10%). For high porosities we have observed a reduction of the error that could be due to the fact that in these cases, the osteonal tissue covers a very large fraction of the matrix volume resulting in a quasi-homogeneous matrix. However, in the case of high porosities, this error might be somewhat underestimated : some of the osteons are in reality more mineralized to an heterogeneity. Modeling cortical bone elastic properties, as any modeling approach, is a matter of idealizations consisting in neglecting certain details of the organization of the tissue. The results of the present study complete a detailed analysis of some idealizations relevant for cortical bone mesoscale elasticity modeling conducted during the last decade. The accuracy of more or less refined models has been evaluated by comparing model predictions to experimental values of anisotropic elastic coefficients obtained from samples of human donors [START_REF] Granke | Change in porosity is the major determinant of the variation of cortical bone elasticity at the millimeter scale in aged women[END_REF][START_REF] Baumann | The relative influence of apatite crystal orientations and intracortical porosity on the elastic anisotropy of human cortical bone[END_REF]Cai et al. 2019b). The samples were typically imaged with high-resolution µCT (as for the present work) to obtain a representation of the microstructure and the matrix mineral content. Then, the question arises as to whether an accurate estimation of elastic properties can be obtained from this information, that is, whether the variability of mesoscale elastic properties can be explained from the mineral content of the matrix and an image of the microstructure.

Porosity (volume fraction of pores) alone is known to explain the largest part of elasticity variations [START_REF] Granke | Change in porosity is the major determinant of the variation of cortical bone elasticity at the millimeter scale in aged women[END_REF]Cai et al. 2019b). Indeed, a popular simple model for cortical bone is a homogeneous transversely isotropic matrix with fixed properties pervaded by elongated pores with a single orientation (e.g., [START_REF] Hellmich | can the diverse elastic properties of trabecular and cortical bone be attributed to only a few tissue-independant phase properties and their interactions[END_REF][START_REF] Parnell | Analytical methods to determine the effective mesoscopic and macroscopic elastic properties of cortical bone[END_REF][START_REF] Baumann | The relative influence of apatite crystal orientations and intracortical porosity on the elastic anisotropy of human cortical bone[END_REF][START_REF] Gagliardi | Tissue mineral density measured at the sub-millimetre scale can provide reliable statistics of elastic properties of bone matrix[END_REF]). This simple model, which disregards the details of the pore structure and the interindividual variations of matrix elasticity, accounts well for the variations of elastic properties with porosity [START_REF] Hellmich | can the diverse elastic properties of trabecular and cortical bone be attributed to only a few tissue-independant phase properties and their interactions[END_REF][START_REF] Granke | To what extent can cortical bone millimeter-scale elasticity be predicted by a two-phase composite model with variable porosity?[END_REF]Cai et al. 2019b), or with the average orientation of mineral crystals [START_REF] Baumann | The relative influence of apatite crystal orientations and intracortical porosity on the elastic anisotropy of human cortical bone[END_REF], with errors compared to experimental data typically below 10% (average error for a collection of samples). Nevertheless this model has several potential biases and its accuracy for predicting the elasticity of a given specimen may not be as good as 10%. As shown by comparing this simple model to an FFT model accounting for the details of the pore network structure (similar to the present work), disregarding the complexity of the pore network and model it as a collection of cylindrical pores leads to an overestimation of effective elasticity of 1% to 15% depending on the elastic coefficient C ij and porosity level (Cai et al. 2019a).

Apart from the porosity, the average matrix mineral content explains a part of the experimental variation of effective mesoscale elasticity: taken together, these two variables were found to explain 76 to 91% of the variability of the different elastic coefficients (Cai et al. 2019b). In the elderly population without specific bone pathology, the coefficient of variation of the average mineral content was found to be 2% (Cai et al. 2019b). Using inverse homogenization based on mesoscale anisotropic stiffness data for the same samples, [START_REF] Cai | Anisotropic elastic properties of human cortical bone tissue inferred from inverse homogenization and resonant ultrasound spectroscopy[END_REF] calculated a coefficient of variation of matrix elastic coefficients between 3 and 7% depending on the elastic constant, which was correlated to the mineral content. (Note that above reported percentages are for coefficients of variations, and that the actual ranges of variation are somewhat larger). Accordingly, an accurate model of cortical bone at the mesoscale should also consider an accurate modeling of the anisotropic elastic properties of the matrix, which would reflect the mineral content of the bone matrix. Then, the question arises as to whether such a model of matrix elasticity should account for the heterogeneous distribution of the mineral at the scale of 10-100 µm, e.g., lower mineral content in osteons compared to interstitial tissue. This is precisely the contribution of the present work where we have shown that disregarding the heterogeneity of the matrix and assuming that its elasticity can be derived from the average matrix mineral content leads to an overestimation of effective elasticity of a maximum of about 7% in extreme cases.

The main limitation of the present study is that we have used a relatively simple micromechanical model to obtain the elasticity tensor with a single parameter f ha associated to the mineral content. Due to the coarse assumptions of the micromechanical model, we argue that this parameter should not be seen as the actual volume fraction of mineral content but rather a lumped parameter driving the matrix stiffness. The specific values of this parameter for osteonal and interstitial tissues in the present study were defined based on the calculated elastic parameters and corresponding experimental data. An alternative could be to use empirical relationships between each matrix elastic coefficient and mineral content [START_REF] Cai | Anisotropic elastic properties of human cortical bone tissue inferred from inverse homogenization and resonant ultrasound spectroscopy[END_REF] or more sophisticated models (Hellmich et al. 2004;[START_REF] Hamed | Multiscale modeling of elastic properties of cortical bone[END_REF]. The rationale for modeling the elasticity tensor with a single parameter is that experimentally, in a population of individuals without documented bone pathology, the collagen properties or the properties of the mineral crystals (such as cristallinity) are not correlated with elastic properties contrary to the mineral content (Cai et al. 2019b).

Another potential modeling issue is the bias introduced by the scheme to calculate effective properties. This, however, can be ruled out as we have shown that for a sufficiently large RVE (typically larger than (1.5) 3 mm 3 ) the upper (kinematic uniform boundary conditions) and lower (static uniform boundary conditions) bounds of apparent properties are very close (Grimal et al. 2011a).

Conclusion

For typical experimental values of osteon thickness and realistic porous microstructures, the error on effective elasticity due to neglecting the mineralized matrix heterogeneity (e.g., not distinguishing osteonal and interstitial tissue) is about 7% at the maximum. This suggests that matrix heterogeneity can likely be disregarded in cortical bone models in most practical cases. In particular, it seems less important to model this heterogeneity compared to (1) the individual-specific matrix elastic properties scaled to the average matrix mineral content and (2) the effect of the complex pore architecture which is not fully captured when the porosity is modeled as a collection of ellipsoids or cylinders. These conclusions should hold for bone tissue without a specific pathology. The effect of matrix heterogeneity may however turn out to be more important in some pathologies [START_REF] Roschger | Bone mineralization density distribution in health and disease[END_REF].

This work provides data to quantitatively discuss the handling of mineralized matrix heterogeneity in cortical bone elasticity models. The results should be helpful to guide the interpretation of effective elastic properties at the millimeter scale. These are involved in finite element models for the prediction of bone response to loads for the assessment of strength (van Rietbergen and Ito 2015; [START_REF] Engelke | FEA to Measure Bone Strength: A Review[END_REF] or remodeling [START_REF] Martelli | Strain energy in the femoral neck during exercise[END_REF]. Experimentally-derived effective elastic coefficients can be used to determine mineralized matrix properties by inverse approaches [START_REF] Sanz-Herrera | Multiscale characterisation of cortical bone tissue[END_REF][START_REF] Cai | Anisotropic elastic properties of human cortical bone tissue inferred from inverse homogenization and resonant ultrasound spectroscopy[END_REF]. Finally, effective elasticity can also be probed in vivo with ultrasound, from which porosity can be estimated by an inverse approach [START_REF] Minonzio | Ultrasound-Based Estimates of Cortical Bone Thickness and Porosity Are Associated With Nontraumatic Fractures in Postmenopausal Women: A Pilot Study[END_REF].

Fig. 4 :

 4 Fig. 4: A transverse cross-section of one cortical bone RVE showing the vascular pores (blue) and the two phases within the mineralized matrix (light orange: osteonal tissue; dark orange: interstitial tissue) for different values of osteon thickness: e = 35µm (left), e = 70µm (center), and e = 140µm (right)

Fig. 6 :

 6 Fig. 6: Difference in elastic properties between the interstitial and osteonal tissues, calculated as 100 × C m;I 33 -C m;O 33 C m;I 33

Fig. 7 :

 7 Fig. 7: Evolution of the effective elastic coefficients C ij with porosity for different osteon thicknesses: (a) e = 35µm, (b) e = 70µm and (c) e = 140µm (solid symbols: two-phase bone matrix composed of interstitial tissue and osteons; open symbols: equivalent homogeneous matrix)

Fig. 8 :

 8 Fig. 8: Relative difference for thickness e = 140µm. Note that the denominator represents the effective elastic moduli considering homogeneous matrix mineralization C = C hom. matrix and ∆ C = | C hom. matrix -C het. matrix |.

Table 1 :

 1 Summary of parameters of the representative volume element. Elasticity tensors of osteons C m;O and interstitial tissue C m;I are calculated with a micromechanical model for a fixed volume fraction of mineral f ha .

	e (µm)	35	70	105	140
	∆ C 11 / C 11 (%)	1.05	2.94	4.75	5.84
	∆ C 33 / C 33 (%)	1.1	3.22	5.42	6.97
	∆ C 23 / C 23 (%)	0.94	2.51	3.87	4.78
	∆ C 44 / C 44 (%)	1.34	3.77	5.92	7.4
	∆ C 66 / C 66 (%)	1.24	3.64	5.16	6.64

Table 2 :

 2 Relative difference averaged over all microstructural samples with porosities ranging from 3 to 20% for different osteon thicknesses e (See the definition of ∆ C/ C in the legend of Fig.8)

Acknowledgements

The authors would like to thank ESRF for the access of beamline at ID 19 and 17 and the help from Cécile Olivier and Françoise Peyrin (CREATIS, CNRS 5220, INSERM U1206, Lyon) for performing SR-µCT experiments.

This work has received financial support from Engineering Department of Sorbonne Université (UFR 919).

Conflict of interest

The authors declare that they have no conflict of interest.


Fabio Coppini

Yiye Jiang

Sonia Tabti

Predictive models on 1D signals in a small-data environment

Keywords: small-data, machine learning, inferring techniques, data augmentation, imputation, Variational Autoencoder, LSTM, Fast Fourier Transform, Discrete Wavelet Transform, time series, class imbalance

This report is concerned by one project done during the Semaine d' Études Mathématiques et Entreprises (SEME) at the Institut de Mathématiques de Bordeaux. The subject, proposed by the company FieldBox.ai, concerns the use of machine learning algorithms and data augmentation techniques, applied to small datasets composed of 1D signals measurements. The target variable is supposed to be continuous, i.e., a regression problem. By first reviewing the literature and existing methods on data augmentation, we propose two procedures to tackle this problem: one allows to create synthetic observations for a specific range of target values and it is based on a perturbation method in Fourier/Wavelet space; the other is based on neural networks and uses a particular version of the Variational Autoencoder known as LSTM-VAE. Our methods are applied to an open dataset, available at the UCI repository, and show encouraging results for a common class of machine learning algorithms.

1 Introduction

Data Science and Machine Learning have gained an exponentially growing attention in the last fifteen years and a lot of effort has been put in deploying algorithms capable of crunching massive datasets. Nowadays, much of the attention is focused on big data and deep neural networks that take days to be trained; however, from a pure industrial viewpoint, most of the common situations that a small to medium-size company needs to face, concern datasets with relatively few observations and classical machine learning algorithms. The aim of this work is to address such small-data framework for a particular kind of data input.

Small-data environment

Suppose D is a dataset consisting of n observations and d features, i.e.,

D = (X 1 , . . . , X d , y) (1.1)
where X i ∈ R n represents the n observations related to the i-th feature, for i = 1, . . . , d, and y ∈ R n stands for n realizations of the target variable. We use the following notation X = (X 1 , . . . , X d ) ∈ R n×d for the n × d matrix of observations and thus (X, y) for the dataset D.

A big-data environment is usually modeled by the fact that when the number of observations n tends to infinity, the number of features d is fixed or potentially grows slower than n. It is more difficult to define a relation for n and d when talking about small-data. Conditioned by industrial applications, we consider a dataset to fit a small-data environment whenever the number of observations per feature, i.e., n/d, is approximately constant. In our setting, we mostly consider the situation n ≈ 500 and d ≈ 10.

1D signals/time-series input

We further restrict our analysis to the case where each single feature is a 1-dimensional (1D) signal, i.e., for i = 1, . . . , d, X i (•) = X i (t) ∈ R with t > 0 and the vector of observations satisfies

X i = (X i (t 1 ), . . . , X i (t n )) , for t 1 < t 2 < • • • < t n , (1.2)
where t i > 0 are different timestamps at which the signal X i is measured.

We consider two different scenarios for the target variable: either it does not explicitly depend on the time variable, i.e., it is a pure 1D signal, or it explicitly depends on it and can be treated as a time-series on its own. The proposed solutions in sections 2 and 4 work under the first and the second assumption respectively.

1D signal

The target variable y(•) is supposed to be a function of the features only. In particular, we suppose that R y(t) = f (X 1 (t), . . . , X d (t)), t > 0.

(1.3)

Remark 1.1. The previous assumption assures that y is constant whenever the signals are so. In particular, this hypothesis apriori excludes the tools available from timeseries analysis such as autoregression models and stationary analysis.

Time-series

The target variable y(•) is supposed to be a function of the features and the time. In particular, we suppose that

R y(t) = f (t, X 1 (t), . . . , X d (t)), t > 0.
(1.4)

Problem formulation

We would like to answer the following question:

does the 1D signal nature of the problem gives more information that can help building more reliable models?

The approach suggested by FieldBox.ai relies on data augmentation techniques, i.e., in generating new synthetic observations in the train set while maintaining a correct target prediction. One of the main goals of this procedure is to reduce overfitting. Within this perspective, the main question can be reformulated in the following points:

1. Is it possible to infer new synthetic observations by exploiting the intrinsic time nature of the dataset?

2. How can one adapt the data augmentation techniques to deal with a continuous target variable?

3. Does augmenting the dataset size with new observations help in improving robustness of the model, e.g., reducing overfitting, while only slightly increasing the bias?

Existing results on data augmentation (see next subsection for the known literature) mainly focus on classification problems: from a conceptual viewpoint, the second question represents a challenging issue. The last question poses an interesting problem already known in imputation statistics. This report proposes a few ideas to tackle previous points, in particular 1. and 2.. Experimental results confirming the choice of this mathematical framework are also provided, we refer to our GitHub repository for the source code, see https: //github.com/fdesmond/seme-ts.

Brief state of the art on data augmentation

There is a fast growing literature on data augmentation techniques [START_REF] Dyk | The Art of Data Augmentation[END_REF][START_REF] Devries | Dataset Augmentation in Feature Space[END_REF][START_REF] Kuchnik | Efficient Augmentation via Data Subsampling[END_REF], yet most of it focuses on classification problems and deep learning algorithms, see, e.g., [START_REF] Shorten | A survey on Image Data Augmentation for Deep Learning[END_REF][START_REF] Wen | Time Series Data Augmentation for Deep Learning: A Survey[END_REF][START_REF] Guennec | Data Augmentation for Time Series Classification using Convolutional Neural Networks[END_REF], with only few results on the regression case and classical machine learning methods, e.g., [START_REF] Dubost | Hydranet: Data Augmentation for Regression Neural Networks[END_REF][START_REF] Krell | Data Augmentation for Brain-Computer Interfaces: Analysis on Event-Related Potentials Data[END_REF][START_REF] Hooker | DARE: Data-Augmented Regression for Extrapolation[END_REF][START_REF] Wang | Data Augmentation for EEG-Based Emotion Recognition with Deep Convolutional Neural Networks[END_REF].

Data augmentation is usually sought for one of the following reasons:

• to reduce class imbalance, especially in classification problems whenever one or more labels are highly underrepresented [START_REF] Chawla | SMOTE: Synthetic Minority Over-sampling Technique[END_REF];

• to avoid overfitting in neural networks when the training dataset is not large enough, e.g., in image recognition [START_REF] Shorten | A survey on Image Data Augmentation for Deep Learning[END_REF];

• to improve robustness of certain classifiers, by interpolation and/or extrapolation of observations in the training set [START_REF] Devries | Dataset Augmentation in Feature Space[END_REF];

• to reduce the curse of dimensionality, e.g., by adding points in sparse regions [START_REF] Hooker | DARE: Data-Augmented Regression for Extrapolation[END_REF].

From an abstract viewpoint, data agumentation techniques fall into the two main categories: domain independent techniques, i.e., methods that apply to roughly any problem; and domain dependent techniques, where the dataset domain (e.g., images, time-series, etc) is taken into account for a tailor-made augmentation.

Examples of domain independent techniques are given by bootstrapping techniques, with our without adding a noise perturbation, see, e.g., [START_REF] Chawla | SMOTE: Synthetic Minority Over-sampling Technique[END_REF][START_REF] Buciluǎ | Model compression[END_REF], but also the use of autoencoders, e.g., [START_REF] Devries | Dataset Augmentation in Feature Space[END_REF][START_REF] Ohno | Auto-encoder-based generative models for data augmentation on regression problems[END_REF]. In such cases, the idea is to create a new observation by slightly perturbing an existing one and to assign it the same label: the perturbation is usually not performed on the input data itself, but on a transformed version of it in some abstract space. This space, usually called feature space, can be again domain independent (e.g., in the case of autoencoders, dimensionality reduction techniques, etc.), or use information from the data domain (e.g., the Fourier space for time-series).

Looking at domain dependent techniques, we've found only two references citing 1D signals [START_REF] Wang | Data Augmentation for EEG-Based Emotion Recognition with Deep Convolutional Neural Networks[END_REF][START_REF] Wen | Time Series Data Augmentation for Deep Learning: A Survey[END_REF]. Most of the results concern data augmentation for time-series classification problems [START_REF] Guennec | Data Augmentation for Time Series Classification using Convolutional Neural Networks[END_REF][START_REF] Forestier | Generating synthetic time series to augment sparse datasets[END_REF][START_REF] Fawaz | Data augmentation using synthetic data for time series classification with deep residual networks[END_REF][START_REF] Smyl | Data Preprocessing and Augmentation for Multiple Short Time Series Forecasting with Recurrent Neural Networks[END_REF][START_REF] Bandara | Improving the Accuracy of Global Forecasting Models using Time Series Data Augmentation[END_REF][START_REF] Oh | Time-Series Data Augmentation based on Interpolation[END_REF] as shown below.

Time-series data augmentation

Following [START_REF] Wen | Time Series Data Augmentation for Deep Learning: A Survey[END_REF], we recognize the following approaches that are used for generating new time-series from existing ones:

• perturb the signal in a suitable feature space (e.g., time domain, frequency domain, time-frequency domain);

• decompose the signal as a time-series with trend, seasonality and residual and bootstrap the residual;

• fit a statistical model on the signal (e.g, ARMA, ARCH, etc.) and perturb its parameter space;

• use learning methods (Autoencoder, Generative Adversarial Network, Reinforcement Learning, etc.).

The new time-series are usually associated to the same label of the parent.

As already anticipated in Section 1.3, a key difficulty of our case is the fact that we are dealing with a regression problem and that there is no label to group the inputs into similar objects. Moreover, much of the cited results assume the time-series to be an observation on itself, while for us, the signals represent the features of the input data, i.e., a column and not a row in the dataset D! Fortunately, we can still manage to reproduce some of the ideas present in the literature, we refer to Section 2 (in particular to subsection 2.3) and Section 4.

Remark 1.2. Interpolation is another possible approach: to add observations by interpolating points in each signal X i and in y. However, this method requires the signals, and the output, to be regular, an hypothesis that we do not want to assume. We refer to [START_REF] Oh | Time-Series Data Augmentation based on Interpolation[END_REF] for some (not so satisfactory) result in this direction.

Proposed approaches

To tackle this problem, we propose two solutions that are presented in sections 2 and 4. The experimental results of the first methodology are discussed right after in Section 3, while we refer to GitHub page4 for the second one.

2 Quantile-based augmentation in feature space

Augmentation procedure

The main principle of the proposed method is the following one: define a time window, i.e., a subset of {t 1 , . . . , t n }, and randomly select some of the features among {X 1 , . . . , X d }; on the time window, replace the chosen features with a suitable perturbed version of them, but keep the other features and the target variable unchanged. Finally, append perturbed and non-perturbed features (together with the target variable) to the original dataset.

To further augment the dataset, on can iterate this procedure with another time window and new features.

Algorithm 1 General procedure Denote by d ∈ {1, . . . , d} the number of features to be perturbed and by σ > 0 the strength of the perturbation one wants to apply. return D ∪ D aug concatenate the observations 5: end procedure Observe that this procedure can be iterated to further augment the dataset by choosing new d and σ. For such implementation, we refer to the GitHub page.

The general procedure dfaug is summarized in Algorithm 1 with the main function aug in Algorithm 2. The functions slice and distort, corresponding to the slicing and perturbation part respectively, are discussed in the following subsections 2.2 and 2.3.

Algorithm 2 Augmentation function

For D = (X 1 , . . . , X d , y) a dataset and d ∈ {1, . . . , d}. The coefficient σ > 0 stands for the distortion magnitude applied on the selected time-series. for k = 1, . . . , d do 4:

X i k ← distort(X i k , σ) perturb the selected feature 5:
end for 6:

return (X 1 , . . . , X d , y) output the synthetic observations 7: end function Observe that if d = d, then all the features are perturbed.

Slicing

Slicing consists in selecting a subset of the rows in the input dataset: in our case, this is equivalent to retaining certain timestamps, e.g., by considering { t1 , . . . , tn } ⊂ {t 1 , . . . , t n } for ñ ≤ n, one obtains the following dataset D sliced = X 1 ( ti ), . . . , X d ( ti ), y( ti ) i=1,...,ñ .

A naive method to slice a dataset is to consider a specific time window, i.e., a subset of consecutive timestamps {t j for j ∈ J ⊂ {1, . . . , n}}. However, this procedure did not yield significant improvements in our setting and we implemented a more sophisticated procedure that we explain hereafter.

In the classification case, one can apply the augmentation procedure separately on each class, in the regression case this is clearly not possible since the target variable apriori takes an infinite number of values and no class is defined. However, if we group observations for which the target variable is not varying too much, we can treat them as belonging to a same abstract class and augment this class by keeping the target values within the range defined by the class itself. There are many different ways to implement this idea, we propose the following procedure.

Select the timestamps t i for which the target variable is varying within a certain range, i.e., for some constant c ∈ R and ε > 0, let J ε c ⊂ {1, . . . , n} be such that

|y(t i ) -c| < ε, for i ∈ J ε c . (2.1)
By choosing different c and , one can get a partition of {1, . . . , n} in terms of intervals as J ε c . Equivalently, one can partition the range of y with real intervals I ε c ⊂ R and obtain J ε c as the indices of the timestamps t i for which y(t i ) ∈ [c -ε, c + ε] =: I ε c . This can be obviously done for a general interval I ⊂ R, see Algorithm 3.

Observe that, if the intervals are suitably chosen, e.g., select I to be a quantile interval of the empirical distribution of y, this procedure allows to create fictitious time-series for which the output is approximately constant.

Algorithm 3 Slicing function

For I a chosen interval range in R. The function dfaug is designed for one specific interval, but if one iterates it on the whole partition, then the new augmented dataset keeps the empirical distribution of the output vector unchanged.

A straightforward way to partition the dataset with respect to the target variable, is to compute the quantile intervals associated to y and to later apply the augmentation procedure to each interval. In the case where the target variable distribution is not uniform, one can apply the augmentation procedure more times to the sliced datasets which correspond to the observations having a lower empirical frequency, e.g., related to rare events, thus providing a way to reduce -some sort of -imbalance in the dataset.

Remark 2.1. Observe that, while the time order in the sliced dataset is maintained, the sampled t j are not necessarily consecutive as in the naive time-window selection. It is possible to break the time structure since we are supposing that y does not explicitly depend on the time variable.

Distortion

To perturb the selected features, we have decided to work in a suitable feature space. In order to exploit the 1D signal nature of the inputs, we have focused on two classical representation space for signals:

1. Fourier space; 2. Wavelet space.

In the frequency domain, we have tested different perturbations both at the level of the phase and at the level of the amplitude. The resulting signal, which has been obtained by inverting back the one perturbed in feature space, shows non-local distortions depending on the most affected modes.

Concerning the Wavelet space, the perturbation has been applied on the first halfplane only, showing the most encouraging results. During hyperparameter tuning, we have found that a mix of these three procedures yields the best results.

We refer to Figure 1 for a example with the three perturbations.

3 Experimental setting

Appliances energy prediction Data Set

A suitable dataset satisfying the hypothesis in subsection 1.2 is given by Appliances energy prediction Data Set. The dataset is used to create regression models of appli-ances energy use in a low energy house. It consists, among other features, of several measurements of temperature, humidity and pressure in the different rooms.

We select 25 out of the 29 features in the original dataset, which are lights, T1, RH 1, T2, RH 2, T3, RH 3, T4, RH 4, T5, RH 5, T6, RH 6, T7, RH 7, T8, RH 8, T9,RH 9, T out, Press mm hg, RH out, Windspeed,Tdewpoint, Appliances. Additionally, we construct two temporal variables: daytime (t d ), dayofweek (t w ), which mark the time point when the observation is recorded in a day and the day it is recorded in a week respectively. To capture the nonlinear trends of time series, we also add the polynomial p(t d , t w , m) regressor in the predictive model to capture the nonlinear trend, which is motivated from the polynomial regression. The order m is set as 12.

The target variable is the energy appliance (energy use in Wh). The entire dataset has 19735 observations, denoted as D A . We refer to the GitHub page related to this project for more details: https://github.com/fdesmond/seme-ts.

Subsampling

We simulate the small-data environment. To this end, we subsample the 5% of D A , this yields the small-data set D B . Train-test split is performed with a proportion of 0.75 : 0.25, so that the small input dataset on which we will train the machine learning algorithms is given by 740 observations. Remark 3.1. The train-test split procedure is performed at the level of D A so that the test set (train set respectively) in D B is actually a 5%-sampling of the test set (train set respectively) in D A .

Hyperparameters

The augmented dataset, which we call D C , is obtained by firstly partitioning the observations into different subsets corresponding to the target interval quantiles; and then by applying the augmentation procedure to different portions of the dataset and with different distortion levels.

We choose to partition the dataset into 9 intervals corresponding to the quantile intervals of y: the augmentation procedure is thus applied to vectors of dimension5 approximately 80.

Model fitting and evaluation

The same machine learning algorithm (e.g., k-NN, Linear Regression, etc.) is fit to the train sets in D A and D B and on D C . For each algorithm, we thus obtain three regression models A, B and C corresponding to the three training datasets.

The evaluation is done by comparing the models B and C on the test set of D A : we choose to use the root mean squared error (RM SE) and the R 2 score as metrics. results can be found in the notebook DL aug.ipynb present in https://github.com/ fdesmond/seme-ts/.

More specifically, we use the long short-term memory variational autoencoder (LSTM-VAE) introduced in [START_REF] Bowman | Generating sentences from a continuous space[END_REF]. We decide not to break the temporal order during train-test-split as well as the small data sampling, see again Section 4.2 for the details. Three folds of facts motivate this choice. Firstly, autoencoders areunsupervised models which is well adapted to our problematic. Secondly, the temporal order of observations can be furthermore exploited; many network architectures have been designed for this sequential data type, with one of the most widely used being LSTM. Lastly, Variational AutoEncoder has sampling scheme, from which we can profit, to generate the restored inputs as augmented data, yet of different noisy levels. We will elaborate this point in the following section.

Data augmentation using LSTM-VAE

The main difference of variational autoencoder compared to the deterministic autoencoder is that, instead of learning a deterministic mapping from each input to its latent vector, it learns the mapping from the input to the posterior distribution of its latent vector. The posterior distributions are always assumed to be multivariate Gussian with diagonal covariance.

To reconstruct an input, VAE samples a latent vector from its posterior distribution, and pass it to the decoder. Thus even for the same input, VAE's ouptut will be different across separate runs. This property can be employed to generate the artificial data associated with the input. Furthermore, for a trained VAE, we can multiply a scale to all the covariance mapped from the inputs, in order to have the smoother/noisier reconstructed signals. These smoother/noisier reconstructions of inputs can also be used as augmented data. We call the scale smooth level, which is the hyperparameter of this approach.

Specific to sequential data, we use LSTM layer in the beginning of the encoder, accordingly at the end of the decoder. The network input is a matrix of dimension timesteps×(d + 1), where d denotes the number of individual time series in equation (1.3), and timesteps is another hyperparameter which is the common length of individual sequences. Note that, we include the response variable as well in the input sequences. We aim to let the network analyze its complex relation with all other regressor sequences, and reconstruct the multivariate sequences which can follow this relation. This property is to be explored more closely in the future. For the network architecture used in the illustrative experiment, we refer to the code on GitHub and the notebook cited at the beginning of the section.

Sampling details and latent space visualization

Since we have LSTM modules in the network, we simply take the first 75% of D A as the training set, leaving the rest as test set. For the construction of small data, we take the last 5% of the training set of D A as the training set of D B , while the first 5% of the big test set as the small test set. We train the network on the small training set to generate the augmented data for ML model fitting.

Before using the trained LSTM-VAE to synthesize artificial data, we would like to know how much information it has learnt from the small training set. In addition to examine the reconstruction performance, we propose to inspect the pattern of latent representation of entire training set of D B . Even though we are not in the classification problem as the common application of VAE, we can still justify the network usefulness by attaching each latent mean with the day of week of its input multivariate sequences6 as classification labels, and verify if the means cluster in the latent space according to their labels.

Conclusions

We did not push the analysis further but we believe that these methods can be improved and give even better results.

The main ideas behind our work can be summed up into the following points:

• (procedure 1) Slice the original dataframe based on target-dependent classes;

• (procedure 1) Perturb the signals at the level of the frequency-domain (i.e., in feature space).

• (procedure 2) Exploit the latent space distribution of VAE to draw samples whose distribution is close to the one of the samples in the original dataset.

1 . 3 1 . 4 1 . 5 6 2 2 . 1 9 3

 1314156219 -data environment . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 1D signals/time-series input . . . . . . . . . . . . . . . . . . . . . . . . 4 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Brief state of the art on data augmentation . . . . . . . . . . . . . . . 5 Proposed approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . Quantile-based augmentation in feature space 7 Augmentation procedure . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 Slicing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3 Distortion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Experimental setting 9 3.1 Appliances energy prediction Data Set . . . . . . . . . . . . . . . . . . 9 3.1.1 Subsampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.1.2 Hyperparameters . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.2 Model fitting and evaluation . . . . . . . . . . . . . . . . . . . . . . . . 10 3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 4 LSTM-VAE based model 11 4.1 Data augmentation using LSTM-VAE . . . . . . . . . . . . . . . . . . 12 4.2 Sampling details and latent space visualization . . . . . . . . . . . . . 13 5 Conclusions 13

1 :D

 1 procedure dfaug(D, d, σ) aug ← aug( D, d, σ) obtain synthetic observations 4:

1 :

 1 function aug(D, d, σ) 2: i 1 , . . . , i d ∈ {1, . . . , d} sampled without replacement 3:

1 :

 1 function slice(D, I) 2:D ← D[y ∈ I]only select the observations for which y(t i ) ∈ I

Figure 1 :

 1 Figure 1: An example of the three different distortions applied to a same signal.

Figure 3 :

 3 Figure 3: Predictions comparison. Comparing a small time-window in D A of the true target (in red) with the two predictions by B (in purple) and C (in blue) respectively.

https://github.com/fdesmond/seme-ts

The dimension here is given by: # number of observations train set / # number of quantile interval.

We use the average day of week if some sequences across different days in a week.

Acknowledgments

The team is indebted towards several persons and structures. In random order, we would like to thank: professor Edoardo Provenzi and the Institut de Mathématiques de Bordeaux 1 for the organization, the kind hospitality and the support the global situation notwithstanding; the company FieldBox.ai 2 and, especially, the lead data scientist Sonia Tabti not only for proposing us the subject, but also for insightful discussions during and after the week; the agency AMIES 3 for making this experience possible and for its effort in merging mathematicians and real world problems.

We also would like to thank all the other participants for the pleasant time we had together and the fruitful exchanges.

This report not only concerns the work done during the SEME in Bordeaux, where the PhD student Julien Bidan also contributed, but also, and mainly, on successive investigations and tries. Further ideas came up during discussions between the first two authors who decided to write them down and who equally contributed in writing this report.

Results

The results are given in Table 3.3. See also Figures 2 and3.

RMSE

LSTM-VAE based model

In this section, we propose another framework for data augmentation, which is based on deep learning technique. The related codes and notebook with the illustrative