Yosra Ben Cheikh
email: yosra.ben-cheikh@univ-lehavre.fr

François Massol

Nathalie Giusti-Petrucciani

Marie-Agnès Travers

Impact of epizootics on mussel farms: Insights into microbiota composition of Mytilus species

Keywords: Pathobiome, mussel mortalities, opportunistic bacteria, core microbiome

Outbreaks of marine mussel mortality on French farms could have different aetiologies. One of them implies Vibrio splendidus strains. Beyond the involvement of this pathogen, there is considerable evidence that diseases often result from interactions between several microbes and the host. In this study, we explored the bacterial communities associated with mussel species and the surrounding water collected from a mussel farm affected by mortalities. The microbiota of Mytilus edulis, Mytilus galloprovincialis and their hybrids displayed an abnormal abundance of Proteobacteria, in particular the genera Vibrio, Cobetia and Arcobacter. Despite the dysbiosis, the Mediterranean mussel showed a different microbiota profile with a higher richness and presence of the phylum Bacteroidetes. Bipartite network analyses at the level of bacteria families confirmed this finding and showed that the microbiomes of M. edulis and the hybrids tended to cluster together. In addition, injection of mussels with the virulent V. splendidus induced less mortality rate in M. galloprovincialis compared to the other Mytilus sp.

suggesting a better resistance of the Mediterranean mussel to infection.

Our findings point to a probable aetiology of pathobiome-mediated disease in mussels. To fully understand this phenomenon, more knowledge is needed on the roles of pathobiotic systems and their development during disease establishment.

Introduction

Over the last two decades, European mussel production has been declining. Among the putative causes of this trend, the spread of diseases, algal blooms, lack of spat, predation and low profitability have been mentioned [START_REF] Avdelas | The decline of mussel aquaculture in the European Union: causes, economic impacts and opportunities[END_REF]. Indeed, since 2014, mass mortality events have been recorded in blue mussel stocks in France [START_REF] Bechemin | mars 2014)[END_REF] and since 2016 in the Netherlands [START_REF] Capelle | Observations on recent mass mortality events of marine mussels in the Oosterschelde, the Netherlands[END_REF]. While no cause could explain this phenomena in the Netherlands [START_REF] Capelle | Observations on recent mass mortality events of marine mussels in the Oosterschelde, the Netherlands[END_REF], the involvement of pathogens in mortalities has been demonstrated in French mussel farms [START_REF] Ben Cheikh | Infection dynamics of a V. splendidus strain pathogenic to Mytilus edulis: In vivo and in vitro interactions with hemocytes[END_REF], 2016;[START_REF] Benabdelmouna | Mortality investigation of Mytilus edulis and Mytilus galloprovincialis in France: An experimental survey under laboratory conditions[END_REF][START_REF] Oden | Multilocus sequence analysis of Vibrio splendidus related-strains isolated from blue mussel Mytilus sp. during mortality events[END_REF].

As filter-feeding organisms, marine mussels are constantly exposed to a variety of microorganisms including pathogenic bacteria that can threaten their lives. In addition to external microbes, like all macroorganisms, mussels interact with the microorganisms they harbour, i.e. their microbiota. This dialogue between the host and its microbiota can help maintain their respective integrity, as shown e.g. in leeches [START_REF] Kikuchi | Spatial and Temporal Population Dynamics of a Naturally Occurring Two-Species Microbial Community inside the Digestive Tract of the Medicinal Leech[END_REF][START_REF] Tasiemski | Reciprocal immune benefit based on complementary production of antibiotics by the leech Hirudo verbana and its gut symbiont Aeromonas veronii[END_REF]. In fact, advances in analytical techniques over the last few decades have favoured the expansion of studies on microbial communities, particularly in the human model, highlighting their importance for their host [START_REF] Amato | Incorporating the gut microbiota into models of human and non-human primate ecology and evolution[END_REF][START_REF] Bäckhed | Host-bacterial mutualism in the human intestine[END_REF][START_REF] Martin | Panorganismal Gut Microbiome-Host Metabolic Crosstalk[END_REF][START_REF] Mushegian | Rethinking "mutualism" in diverse host-symbiont communities[END_REF]. Considered the "last organ" in humans [START_REF] Baquero | The microbiome as a human organ[END_REF], a large body of evidence has revealed the beneficial contribution of the microbiota to basic biological processes such as digestive, metabolic and immune functions [START_REF] Callens | Food availability affects the strength of mutualistic host-microbiota interactions in Daphnia magna[END_REF][START_REF] Sirisinha | The potential impact of gut microbiota on your health:Current status and future challenges[END_REF][START_REF] Wang | The Human Microbiota in Health and Disease[END_REF] and its involvement in the ecoevolutionary dynamics of their host [START_REF] Amato | Incorporating the gut microbiota into models of human and non-human primate ecology and evolution[END_REF][START_REF] Kohl | Microbial communities exhibit host species distinguishability and phylosymbiosis along the length of the gastrointestinal tract[END_REF][START_REF] Macke | Life history and eco-evolutionary dynamics in light of the gut microbiota[END_REF][START_REF] Moeller | Cospeciation of gut microbiota with hominids[END_REF][START_REF] Murall | Chapter Five -Invasions of Host-Associated Microbiome Networks[END_REF]. However, like other ecosystems, the composition of the microbiota can be modulated by several factors such as antibiotics and pathogens leading to an altered community structure or dysbiosis [START_REF] Clark | Distinct Shifts in Microbiota Composition during Drosophila Aging Impair Intestinal Function and Drive Mortality[END_REF][START_REF] Costello | The application of ecological theory toward an understanding of the human microbiome[END_REF][START_REF] Lamont | Polymicrobial synergy and dysbiosis in inflammatory disease[END_REF][START_REF] Martin | Panorganismal Gut Microbiome-Host Metabolic Crosstalk[END_REF]. In case of dysbiosis, the mutual harmony between bacteria is disrupted and the healthy balanced microbiota turns into a disease-inducing pathobiota [START_REF] Pitlik | How holobionts get sick-toward a unifying scheme of disease[END_REF][START_REF] Tungland | Chapter 11 -Intestinal Dysbiosis in Obesity, Metabolic Syndrome and Related Metabolic Diseases: Therapeutic Strategies Utilizing Dietary Modification, Pro-and Prebiotics, and Fecal Microbial Transplant (FMT) Therapy[END_REF][START_REF] Wang | The Human Microbiota in Health and Disease[END_REF].

The large number of findings in the human model has sparked research interest in other organisms including bivalves. Numerous studies have investigated microbiota structure in oysters, mussels and to a lesser extent clams, under the influence of factors such as host genotype, environment and pathogens (for review see [START_REF] Paillard | Recent advances in bivalve-microbiota interactions for disease prevention in aquaculture[END_REF]. The structure of the microbiota can be divided into at least two types of communities, namely a stable community of resident bacteria within a host population ("the core microbiome") and a variable microbial community that depends on the surrounding environment [START_REF] Santibáñez | First characterization of the gut microbiome associated with Mytilus chilensis collected at a mussel farm and from a natural environment in Chile[END_REF]. Data on the role of the core microbiome are limited, although some authors suggest a contribution to essential host biological functions [START_REF] Lemanceau | Let the core microbiota be functional[END_REF]. Furthermore, the balance between the core and variable microbiome is still unclear, especially as the dynamics of host-associated communities are directly affected by abiotic and biotic factors. For example, elevated water temperature appears to be a major determinant and strongly influences the composition of the microbiota in mussels and oysters by also favouring pathogens such as Vibrio sp. [START_REF] Li | Temperature elevation and Vibrio cyclitrophicus infection reduce the diversity of haemolymph microbiome of the mussel Mytilus coruscus[END_REF][START_REF] Li | Elevated seawater temperatures decrease microbial diversity in the gut of Mytilus coruscus[END_REF][START_REF] Lokmer | Hemolymph microbiome of Pacific oysters in response to temperature, temperature stress and infection[END_REF]. In addition, water pollution by nanoparticles, nanoplastics or chemicals can affect the diversity of the microbiota as has been demonstrated in Mytilus galloprovincialis [START_REF] Auguste | Impact of nanoplastics on hemolymph immune parameters and microbiota composition in Mytilus galloprovincialis[END_REF][START_REF] Auguste | Exposure to TiO2 nanoparticles induces shifts in the microbiota composition of Mytilus galloprovincialis hemolymph[END_REF][START_REF] Balbi | Insight into the microbial communities associated with first larval stages of Mytilus galloprovincialis: Possible interference by estrogenic compounds[END_REF]. Besides abiotic factors, microbiota imbalance has been linked to diseases [START_REF] King | Oyster disease in a changing environment: Decrypting the link between pathogen, microbiome and environment[END_REF]. Indeed, infection of Crassostrea gigas by OsHV-1 virus induced dysbiosis with subsequent bacteraemia by opportunistic bacteria [START_REF] Lorgeril | Immune-suppression by OsHV-1 viral infection causes fatal bacteraemia in Pacific oysters[END_REF]. Similarly, field studies on oysters affected by mortalities have shown signs of microbiota disruption [START_REF] Lasa | Dynamics of the Pacific oyster pathobiota during mortality episodes in Europe assessed by 16S rRNA gene profiling and a new target enrichment nextgeneration sequencing strategy[END_REF].

For mussels, only a few studies have explored microbiota dynamics in interaction with pathogens. [START_REF] Li | Temperature elevation and Vibrio cyclitrophicus infection reduce the diversity of haemolymph microbiome of the mussel Mytilus coruscus[END_REF] reported the reduction in microbial diversity of M. coruscus haemolymph after exposure to Vibrio cyclitrophicus. Likewise, we recently demonstrated that experimental infection with V. splendidus reduces the diversity of the M. edulis microbiota and favours pathobiontic bacteria [START_REF] Ben Cheikh | Vibrio splendidus infection induces dysbiosis in the blue mussel and favors pathobiontic bacteria[END_REF]. However, the response of the mussel microbiota during disease outbreaks remains unknown.

The aim of this study is to provide new insights into the microbiota of mussels experiencing mortality outbreaks on French farms. For this purpose, the bacterial communities associated with M. edulis, M. galloprovincialis and their hybrids, grown in the same environment, and the surrounding water were analysed. The differences between the microbial communities found in the different species were assessed using ecological network methods to cluster microbiota and link the obtained classification to the type of sampled mussels. The susceptibility of the different types of mussels to infection was then explored by experimental infection with V. splendidus strains.

Material and methods

Mussel collection

45 farmed adults of Mytilus sp. (M. edulis, M. galloprovincialis and hybrid M. edulis x galloprovincialis, 15 individuals from each taxon, size between 4 and 5 cm) were collected from a mussel farm affected by massive mortalities (Bay of Lannion, France 48°44 44.8 N, 3°35 24.9 W) in September 2018 (mean temperature = 17°C). The animals were immediately transported to the laboratory and randomly dispatched in different temperature-controlled (16°C) aerated water tanks filled with seawater sampled from the corresponding site (salinity = 35 PSU). Mussels were taken at random from two different tanks for microbiota analysis or for experimental infection. Seawater from the same site was also sampled for microbiota analyses.

Microbiota analysis

After cleaning, mussels were removed from their shells, flash frozen in liquid nitrogen and stored at -80°C (15 individuals/species). Serial filtrations of 5L of seawater were performed using 8 μm and 5 µm porosity (WhatmanTM NucleoporeTM Polycarbonate membrane filter, 47 mm diameter, UK). The resulting filtered seawater was centrifuged for 15 min at 11000 g.

The membranes and bacterial pellets were stored at -80°C.

DNA extraction, amplicon library construction and Illumina MiSeq sequencing were performed by Biofidal (Vaulx-en-Velin, France). DNA was extracted from frozen mussels ground in liquid nitrogen and from bacteria pellet using Mag Bind® Universal Pathogen DNA Kit (Omega BIO TEK, Georgia, USA) and from the different filters using Quick DNA/RNA MagBead kit (Zymo Research, California, USA) according to the manufacturer's protocol. DNA concentration and quality were checked with the QuantiFluor® dsDNA System (Promega, France). The 16S rRNA gene of bacterial communities was amplified using the 5X HOT BIOAmp ® BlendMaster Mix (Biofidal) and the 341F: CCTACGGGNGGCWGCAG (17 pb) and 805R: GACTACHVGGGTATCTAATCC (21 pb) primers (designed by Biofidal) targeting the variable V3-V4 loops [START_REF] Klindworth | Evaluation of general 16S ribosomal RNA gene PCR primers for classical and next-generation sequencing-based diversity studies[END_REF]. After purification of the PCR products by SPRIselect reagent kit (Beckman Coulter), a second amplification was performed using index adapter oligos (P5/P7). The purified PCR products were normalized to obtain an equimolar library and subsequently pooled. Negative controls and ZymoBIOMICS microbial community standards were included during the library preparation (Supplemental Figure S1).Amplicons were sequenced on the MiSeq platform with V3 reagents producing 2 X 300 bp paired-end reads.

After demultiplexing, primer sequence adapters were trimmed from the reads using cutadapt software version 1.12 [START_REF] Martin | Cutadapt removes adapter sequences from high-throughput sequencing reads[END_REF]. Paired-end sequence reads were collapsed into so-called pseudoreads using sequence overlap with USEARCH version 9.2 [START_REF] Edgar | Search and clustering orders of magnitude faster than BLAST[END_REF]. After chimera and singletons removing, filtered pseudoreads were aligned to the RDP database version 11.5 [START_REF] Cole | Ribosomal Database Project: data and tools for high throughput rRNA analysis[END_REF] with Snap-aligner version 1.0 beta.23 [START_REF] Zaharia | Faster and more accurate sequence alignment with SNAP[END_REF] to identify the taxonomy. Finally, USEARCH was used to create the OTU-table and to calculate alpha diversity. Principal coordinate analyses (PCoA) were performed via xlstat using the Bray-Curtis distance matrix. All FastQ files were deposited in SRA-NCBI under the project accession number PRJNA923705.

Host-microbiota network

To analyse the structure of the host-microbiota association dataset, we used methods borrowed from the analysis of bipartite ecological networks. First, we aggregated data on microbial taxa at the family level and transformed the microbial family abundance table into presence/absence data (i.e., all abundances > 0 were set to 1). We then modelled the obtained bipartite network of hosts by microbial families using a latent block model (LBM) [START_REF] Govaert | Block clustering with Bernoulli mixture models: Comparison of different approaches[END_REF][START_REF] Leger | Clustering methods differ in their ability to detect patterns in ecological networks[END_REF]. The clustering obtained with LBMs is indeed more useful than the one obtained with modularity-search algorithms [START_REF] Newman | Finding community structure in networks using the eigenvectors of matrices[END_REF][START_REF] Traag | From Louvain to Leiden: guaranteeing wellconnected communities[END_REF] as blocks can identify densely connected as well as sparsely connected regions of the network. We tested the correlation between the obtained clustering of hosts and their classification among species using a method described in [START_REF] Massol | A methodological framework to analyse determinants of host-microbiota networks, with an application to the relationships between Daphnia magna's gut microbiota and bacterioplankton[END_REF]: we computed the normalised mutual information (NMI) index between the two classifications of the hosts (by host type and by LBM block); we tested the significance of NMI using randomizations of the bipartite network through the configuration model [START_REF] Strona | A fast and unbiased procedure to randomize ecological binary matrices with fixed row and column totals[END_REF].

Experimental infection

Two strains of V. splendidus, the virulent 10/068 1T1 and the non-virulent 12/056 M24 T1 were used for experimental infection by injection as described in previous studies [START_REF] Ben Cheikh | Infection dynamics of a V. splendidus strain pathogenic to Mytilus edulis: In vivo and in vitro interactions with hemocytes[END_REF], 2016). Briefly, bacteria were grown overnight in Luria Bertani NaCl 20 g.L -1 at 22°C with constant agitation (80 rpm) and resuspended in filtered sterile seawater (FSSW, OD600nm = 1, 2.10 8 CFU mL -1 ) after double washing (FSWW) /centrifugation (1200 g for 10 min).

Anesthetized mussels (2-3 h in magnesium chloride solution at 50 g.L -1 , 1/4: v/v seawater/freshwater) were intramuscularly injected with 100 µL of bacterial suspension or FSSW for negative controls. The animals were transferred to tanks (3 replicate tanks, 10 mussels per tank) filled with 2 L of UV-treated and filtered seawater supplemented with 50 mL of phytoplankton (Isochrysis galbana) and maintained under static conditions at 16°C with aeration. Mortality was monitored over 8 days with removing of newly dead mussels from the tanks.

Statistical analyses

SigmaPlot 12 (Systat Software Inc., Chicago, IL) and R version 4.0.5 [START_REF] Core | A Language and Environment for Statistical Computing[END_REF] were used for statistical analysis. A comparison of OTUs between mussel types and/or seawater was performed. Values were tested for normality (Shapiro-Wilk test) then a one-way ANOVA was performed. Pairwise comparisons between groups were assessed by post-hoc analyses (with Holm-Sidak correction). Statistical significance was accepted for p < 0.05.

Results

Microbial structure in mussels and seawater

The 48 samples (15 M. edulis, 15 M. galloprovincialis, 15 hybrid and 3 seawater samples) sequenced using the 16S rRNA gene generated 4 842 875 raw reads ranging from 70 302 to 126 377 reads per sample with an average of 100 893 reads (Supplemental Table S1). After Chimera filtering, unique pseudoreads averaged 38 597 reads.

The clustering of reads in OTUs revealed a higher richness in the Mediterranean mussel M. galloprovincialis (significant difference with M. edulis, ANOVA p < 0.05) and in seawater with a mean absolute abundance of 80 486 ± 3 214 OTUs and 80 367 ± 808 OTUs respectively. The blue mussel M. edulis and the hybrids had a mean absolute abundance of 68 068 ± 2 322 and 70 637 ± 2 783 OTUs respectively (Figure 1 A). Alpha diversity metrics were higher in M. galloprovincialis compared to the other two mussel species and seawater (non-significant data, supplemental Figure S2A). Beta diversity analysis revealed dissimilarities between mussels and seawater (Supplemental Figure S2B). The compositional structure at the phylum level indicated a similarity between M. edulis and the hybrids (Figure 1 B) with a predominance of Proteobacteria (relative abundance r.a. of 92.3% to 99.4% in M. edulis and 89.9% to 99.9% in hybrids) and the presence of Bacteroidetes (r.a. from 0.3% to 5.5% in M. edulis and from 0 to 5.4% in hybrids). Other rare phyla with an abundance > 1 were detected in some hybrids, e.g., Firmicutes (replicates 2, 6 and 7) and Tenericutes (replicates 1 and 13). Mediterranean mussels harboured a lower proportion of Proteobacteria (r.a. from 47.3% to 99.6%) and a higher proportion of Bacteroidetes (r.a. from 0 to 23.5%). OTUs belonging to Tenericutes, Verrucomicrobia, Fusobacteria and Cyanobacteria were also identified in some mussels with r.a. > 1. The composition of the surrounding seawater was similar to that of the mussels with a notable r.a. of Bacteroidetes (mean 21.1% ± 6.4%) and the predominance of Proteobacteria (mean 77.4% ± 6.9%).

A comparative analysis of the most abundant families (Top 30) and genera revealed a more specific layout for each group although a similar bacterial community was detected in most samples (Figure 2, Supplemental Table S2). Overall, Vibrionaceae, Halomonadaceae, Pseudoalteromonadaceae and Campylobacteraceae were the dominant families in the blue mussel and hybrids (mean r.a. 5-35%) with notably the genera Vibrio, Cobetia, Pseudoalteromonas and Arcobacter but also unclassified Epsilonproteobacteria. Similarly, the Mediterranean mussel harboured Vibrionaceae (Vibrio), Pseudoalteromonadaceae (Pseudoalteromonas) and unclassified Epsilonproteobacteria but also a higher proportion of Campylobacteraceae (Arcobacter) and Flavobacteriaceae (especially Tenacibaculum and Flavobacterium) (mean r.a. 18.2% and 6.8% respectively). In contrast, in addition to unclassified Gammaproteobacteria, the prevailing taxa (r.a. up to 5%) in seawater samples were the families Rhodobacteraceae, Oceanospirillaceae, Pseudoalteromonadaceae,

Crocinitomicaceae and Flavobacteriaceae and the genera Sulfitobacter, Oceanospirillum, Pseudoalteromonas. The majority of Vibrio identified in mussels have been affiliated to the V.

Splendidus group (Supplemental Table S3).

Association between microbial family structure and mussel taxa

The LBM inferred 5 clusters of hosts and 4 clusters of bacterial families (Figure 3A).

Assignment of host individuals to clusters resulted in probabilities of membership at least as high as 88.3%; by contrast, a few bacterial families were only assigned to a particular cluster with a probability of 50.8% (Supplemental Figure S3). interactions can be found in Supplemental Table S4.

The classification of the hosts (48 samples) according to their LBM cluster was significantly congruent with their classification based on sample type (mussel taxon or seawater sample, i.e., 4 groups) according to the NMI (NMI = 0.381, p value < 0.0001 when tested against 10,000 randomizations of the bipartite network). Indeed, a graphical representation of the two classifications evince that M. galloprovincialis projects almost completely onto LBM cluster 2 while LBM clusters 4 and 5 project mostly onto M. edulis and hybrids (Figure 3B). Seawater samples coincide almost perfectly with LBM cluster 3.

Susceptibility of mussels to pathogens

To assess the susceptibility of mussel populations to bacterial infection, M. edulis, M.

galloprovincialis and hybrids were injected with two strains of V. splendidus of contrasting virulence. Exposure of mussels to avirulent strain 12/056 M24 T1 induced limited mortalities in all mussel groups, ranging from 3.3% ± 3.3% in hybrids to 13.3% ± 3.3% in M. edulis at day 8 (Figure 4). In contrast, injection of the virulent strain 10/068 1T1 induced higher mortalities with different rates in the mussel species. The first mortalities appeared 24h post-injection and increased progressively until stabilization on day 8. Hybrids and M. edulis had different patterns at the beginning of the infection (16.6% ± 12% for M. edulis and 46.6% ± 6.6% for hybrids on day 1) but reached 76.6% of mortality at the end of the experiment while M. galloprovincialis had lower mortalities (23.3% ± 3.3% on day 1 and 46.6% ± 12% on day 8). No or negligible mortality was observed in controls.

Discussion

Since 2014, mass mortalities have been reported in French mussel farms. While M. edulis, M.

galloprovincialis and their hybrids represent the main farmed species, the epizootics seem to mostly affect the blue mussel [START_REF] Bechemin | mars 2014)[END_REF]. In general, investigations of mussel disease have so far been limited to the detection and isolation of known pathogens from moribund animals such as V. splendidus [START_REF] Oden | Multilocus sequence analysis of Vibrio splendidus related-strains isolated from blue mussel Mytilus sp. during mortality events[END_REF], combined with experimental infection in the laboratory [START_REF] Ben Cheikh | Infection dynamics of a V. splendidus strain pathogenic to Mytilus edulis: In vivo and in vitro interactions with hemocytes[END_REF], 2016). In this study, we provide an overview of the bacterial communities associated with three mussel species and the surrounding water sampled from a farm affected by mortality epidemics.

The microbiota of the three Mytilus sp. is dysbiotic

According to previous studies, the microbiota of healthy marine bivalves, and mussels in particular is characterized by the presence of different phyla with a predominance of Proteobacteria but also other bacteria belonging to Cyanobacteria, Firmicutes, Planctomycetes, Bacteroidetes, Actinobacteria, Verrucomicrobia, Tenericutes and Chloroflexi [START_REF] Ben Cheikh | Vibrio splendidus infection induces dysbiosis in the blue mussel and favors pathobiontic bacteria[END_REF][START_REF] Musella | Tissue-scale microbiota of the Mediterranean mussel (Mytilus galloprovincialis) and its relationship with the environment[END_REF]. Here, the analysis of whole animal microbiota showed a severe dysbiosis marked by an abnormal dominance of Proteobacteria (Fig. 1, mean r.a. ranging from 86.5% to 96.1%). Microbiota imbalance has been linked to mortality events and infectious diseases in bivalves. Field analysis of the Pacific oyster C. gigas showed that animals experiencing mortality outbreaks at different life stages displayed signs of dysbiosis resulting in a change in the diversity and composition of the microbiota [START_REF] Lasa | Dynamics of the Pacific oyster pathobiota during mortality episodes in Europe assessed by 16S rRNA gene profiling and a new target enrichment nextgeneration sequencing strategy[END_REF]. Similarly, a study on Unionida freshwater mussel during mass mortality events revealed bacterial invasion and shifts in the bacterial microbiome [START_REF] Richard | Mussel mass mortality and the microbiome: evidence for shifts in the bacterial microbiome of adeclining freshwater bivalve[END_REF]. For marine mussels, there are no data available to our knowledge on animal microbiota during mortality events. However, we have recently reported a disruption of the blue mussel microbiota during experimental infection with V. splendidus and demonstrated the implication of some opportunistic pathobiontic bacteria in the infection [START_REF] Ben Cheikh | Vibrio splendidus infection induces dysbiosis in the blue mussel and favors pathobiontic bacteria[END_REF]. Thus, even if the causal origin of mussel mortalities in France still not clear, we can establish a link with pathogenic bacteria, especially as they are detected in Mytilus tissues. Our analysis highlighted the genera Vibrio, Cobetia and Arcobacter as the predominant bacteria in mussel species (Fig. 2). The abundance of Vibrio sp. in the microbiota of mussels is not surprising as these genera are naturally occurring in marine environments and certain strains are often associated with outbreaks of mortality marine bivalves [START_REF] Travers | Bacterial diseases in marine bivalves[END_REF]. In our study, we mainly detected strains belonging to the V. Splendidus clade which includes pathogens of bivalves such as V.

splendidus, V. tasmaniensis and V. crassostreae [START_REF] Ben Cheikh | First evidence for a Vibrio strain pathogenic to Mytilus edulis altering hemocyte immune capacities[END_REF][START_REF] Bruto | Vibrio crassostreae, a benign oyster colonizer turned into a pathogen after plasmid acquisition[END_REF][START_REF] Duperthuy | Use of OmpU porins for attachment and invasion of Crassostrea gigas immune cells by the oyster pathogen Vibrio splendidus[END_REF] but also less studied strains showing virulence for some marine organisms (e.g. V. toranzoniae and V. kanaloae) and environmental species without pathogenic potential (e.g. V. gigantis and V. atlanticus) [START_REF] Lasa | Comparative genomic analysis of two Vibrio toranzoniae strains with different virulence capacity reveals clues on its pathogenicity for fish[END_REF][START_REF] Romalde | New Vibrio species associated to molluscan microbiota: a review[END_REF]. In accordance with our findings, different studies have reported the presence of V. splendidusrelated strains in different French mussel farms affected by mortalities [START_REF] Bechemin | mars 2014)[END_REF]Charles et al., 2020a;[START_REF] Oden | Multilocus sequence analysis of Vibrio splendidus related-strains isolated from blue mussel Mytilus sp. during mortality events[END_REF]. However, these studies are restrictive as they are based on the isolation of culturable bacteria of interest and are therefore not representative of the animal's microflora. In addition to V. splendidus strains, Francisella halioticida has been recently described as potentially involved in blue mussel mortality [START_REF] Cano | First Detection of Francisella halioticida Infecting a Wild Population of Blue Mussels Mytilus edulis in the United Kingdom[END_REF]Charles et al., 2020b). The implication F. halioticida in mussel disease is not confirmed in our study since this strain was detected at very low proportions in some individuals of blue mussels and hybrids (data not shown). In contrast, metagenomic analyses showed the prevalence of bacteria of the genus Cobetia and Arcobacter. While knowledge of Cobetia strains is limited, several reports described the dominance of Arcobacter sp. in oysters and mussels in a context of disease suggesting them as opportunistic bacteria [START_REF] Ben Cheikh | Vibrio splendidus infection induces dysbiosis in the blue mussel and favors pathobiontic bacteria[END_REF][START_REF] Li | Temperature elevation and Vibrio cyclitrophicus infection reduce the diversity of haemolymph microbiome of the mussel Mytilus coruscus[END_REF][START_REF] Lokmer | Hemolymph microbiome of Pacific oysters in response to temperature, temperature stress and infection[END_REF][START_REF] Lorgeril | Immune-suppression by OsHV-1 viral infection causes fatal bacteraemia in Pacific oysters[END_REF]. Nevertheless, the potential role of these bacteria needs further investigation, as some species (e.g. A. cryaerophilus, A. ellisii, A.

nitrofigilis) have virulence genes and the ability to colonise and enter human cells [START_REF] Ferreira | Insights in the pathogenesis and resistance of Arcobacter: A review[END_REF][START_REF] Levican | Adherence to and invasion of human intestinal cells by Arcobacter Species and their virulence genotypes[END_REF].

Microbiota of M. galloprovincialis differs from M. edulis and hybrids

The analysis of microbiomes using a latent block model identified four clusters of bacteria families and five clusters of host individuals (Fig. 3A). Since this clustering proved to be congruent with the one obtained from host types, with a significantly high NMI, this suggests that microbiota are different between mussel species. Indeed, host clusters 4 and 5 belong almost exclusively to M. edulis and hybrids, while M. galloprovincialis individuals are almost all found in host cluster 2. The differences in bacteria family richness between clusters thus confirm that M. galloprovincialis is associated to a richer microbiome than both M. edulis and hybrids. The microbial cluster 3 is more associated with host cluster 2 than with host clusters 4 and 5 (probabilities of interactions: 36.6% vs. 18.1% and 8.7%; Supp. Table S4). Therefore, it probably contains some of the bacterial families peculiar to M. galloprovincialis. Despite these differences, the latent block model also uncovered a potential list of bacterial families that could constitute the 'core microbiome' of all sampled mussels, i.e., microbial cluster 1. This group of 56 bacterial families comprises quite common families such as Aeromonadaceae, Enterobacteriaceae, Pseudomonadaceae, Rickettsiaceae and Vibrionaceae (Supp. Table S5). In particular, this cluster contains the three families comprising the three genera described in Fig.

(Halomonadaceae, Campylobacteraceae and Vibrionaceae).

The shared core microbiota among the three mussel species confirms that the surrounding environment plays a major role in shaping the composition of the microbiome. In agreement with our findings, previous reports have shown that the microbiota of co-cultivated oysters and mussels displayed a high degree of similarity [START_REF] Pierce | Gut Microbiomes of the Eastern Oyster (Crassostrea virginica) and the Blue Mussel (Mytilus edulis): Temporal Variation and the Influence of Marine Aggregate-Associated Microbial Communities[END_REF][START_REF] Vezzulli | Comparative 16SrDNA gene-based microbiota profiles of the Pacific Oyster (Crassostrea gigas) and the Mediterranean Mussel (Mytilus galloprovincialis) from a Shellfish Farm (Ligurian Sea, Italy)[END_REF].

However, host genetics also appears to be a driving factor in the mussel microbiota given the dissimilarity between the microbial communities associated with M. galloprovincialis and the other two species. The specific microbiome of the Mediterranean mussel may confer fitness advantages to the host including resistance and survival from disease as demonstrated in other organisms [START_REF] Kueneman | Probiotic treatment restores protection against lethal fungal infection lost during amphibian captivity[END_REF][START_REF] Madison | A marine probiotic treatment against the bacterial pathogen Vibrio coralliilyticus to improve the performance of Pacific (Crassostrea gigas) and Kumamoto (C. sikamea) oyster larvae[END_REF][START_REF] Zackular | The gut microbiome modulates colon tumorigenesis[END_REF]. Compared to M. edulis and hybrids, M. galloprovincialis showed a lower abundance of Cobetia and Arcobacter species. In addition, despite the high prevalence of V. Splendidus strains, the Mediterranean mussel survived experimental infection with the virulent V. splendidus 10/068 1T1 better (Fig. 4). Resistance of M. galloprovincialis to infectious diseases has been previously reported. [START_REF] Fuentes | Growth, mortality, pathological conditions and protein expression of Mytilus edulis and M. galloprovincialis crosses cultured in the Rı á de Arousa (NW of Spain)[END_REF] described higher warm season mortality of M. galloprovincialis

x M. edulis hybrids than M. galloprovincialis individuals. Similarly, Benabdelmouna et al.

(2018) suggested a higher susceptibility in M. edulis relative to M. galloprovincialis at the spat stage, which is consistent with field mortalities. Thus, the microbiome of M. galloprovincialis may be related to the efficiency of defence mechanisms and resistance to infectious diseases.

Characterization of host-microbiota dynamics at the spatial and temporal scales will be crucial to understanding the roles of the microbiome in maintaining host homeostasis.

Conclusions

In this study, we provide for the first time an overview of the microbiota of M. edulis, M.

galloprovincialis and their hybrids in the context of mortalities in French mussel farm. We reported dysbiosis of the microbiota in the three mussel species and a difference between bacterial communities in seawater and mussel tissues but also between mussel species. Among the pathobiontic bacteria detected, members of Vibrio, Cobetia and Arcobacter genera were the most abundant suggesting their involvement in mussel disease. In addition, M. galloprovincialis showed a different microbiota profile and was also more resistant to infection by V. splendidus.

Thus, we can speculate on the link between the microbiota of Mediterranean mussels and their resistance to infectious disease. Further research is needed to understand the nature of the interactions of mussel microbiota, including the pathobiome, with their host, which may ultimately help to manage and monitor host health to prevent mortalities in mussel farms.

CRediT authorship contribution statement

Yosra Ben Cheikh: Conceptualization, Methodology, Validation, Formal analysis, Investigation, Writing -Original Draft, Writing -Review & Editing, Visualization, supervision.

  The first cluster of hosts consisted in 2 samples, one of M. galloprovincialis and one hybrid, which harboured high numbers of bacterial families (223 families on average). The second cluster of hosts consisted in 20 samples (1 hybrid, 7 M. edulis and 12 M. galloprovincialis) with an average of 141.9 bacterial families per sample. The third cluster included only 4 samples (one replicate of M. galloprovincialis together with the three seawater samples) with an average of 125.5 families per sample. The fourth cluster of hosts consisted in 16 samples (7 M. edulis, 1 M. galloprovincialis, 8 hybrids) with an average of 109 families per sample. The fifth cluster of hosts consisted in 6 samples (1 M. edulis, 5 hybrids) with an average of 69 families per sample. The four clusters of bacterial families harboured 56, 57, 89 and 176 families, respectively, and were found on average in45.46, 30.44, 13.44 and 2.63 host samples respectively. Probabilities of host-bacterial family

François Massol :

 Massol Methodology, Formal analysis, Writing -Original Draft, Writing -Review & Editing, Visualization. Nathalie Giusti-Petrucciani: Investigation. Marie-Agnès Travers: Conceptualization, Writing -Review & Editing.

Figure 1 .

 1 Figure 1. Total OTUs and the relative abundance of bacterial phyla in mussels and seawater. (A) total OTUs presented by box-and-whisker

Figure 2 .

 2 Figure 2. The most abundant bacterial family and genera in Mytilus sp. and seawater. (A-D) heatmaps revealing the top 30 bacterial families

Figure 3 .

 3 Figure 3. Clustering of hosts and bacterial families. (A) Representation of the association between hosts (rows) and bacterial families (columns)

Figure 4 .

 4 Figure 4. Resistance of Mytilus sp. to V. splendidus infection. Mussels were intramuscularly injected with

Acknowledgements

The authors are indebted to Biofidal (Vaulx en Velin, France) for valuable technical assistance.

This study is set within the framework of the « Laboratoire d'Excellence (LabEx) » TULIP (ANR-10-LABX-41).

Data availability statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Declaration of competing interest

The authors declare no conflicts of interest. the highest among all the probabilities of cluster memberships for each host and family). The white circle indicates the value of the median of the distribution. S1. Number of pseudoreads in mussels and seawater before and after chimera filtering. S2. Mean relative abundance of the top 30 bacterial families identified in mussels and seawater, n = 15 for mussels and n = 3 for seawater. Data are percentage ± SEM. S3. Mean absolute abundance of Vibrio, Cobetia and Arcobacter species identified in mussels, n = 15 for mussels. Data are mean ± SEM. Only mean abundance > 100 OTUs in at least one mussel species is considered. S4. Matrix summarizing the estimated probabilities of association between hosts and bacterial families according to the LBM clusters they belong to. For instance, associations between a random host from host cluster 3 and a random family from microbial cluster 4 occur with a probability of 7.88%. S5. Inventory of the bacterial families found in each of the four latent block model clusters for bacteria.

Supplemental

Supplemental Table

Supplemental Table

Supplemental Table

Supplemental Table

Supplemental Table


Haizhou Liu

Lingtao Mao

Xuyang Chang

François Hild

Multiscale damage analyses of red sandstone in uniaxial compression based on advanced digital volume correlation

Keywords: Digital volume correlation, sandstone, X-ray computed tomography, damage initiation and growth, voxel-scale exponent of power-law interpolation 𝐷: damage variable {𝐟}: nodal force stiffness matrix [𝐊 𝐷 𝐞 ]: damaged elementary stiffness matrix ℓ: mesh size ℓ 𝑟𝑒𝑔 : regularization length 𝒖: sought displacement field 𝑣 𝑖 : nodal displacement {𝒗}: displacement column vector dimensionless correlation residual 𝜌 𝑟 : correlation residual 𝑤: mechanical weight 𝜌 𝑚𝑒𝑐 : mechanical residual Acoustic emission COD: crack opening displacement DVC: digital volume correlation FE: Finite element FWHM X-ray computed microtomography

Highlights

• Macroscale, mesoscale, and microscale damage characterization and quantification in uniaxial compression of red sandstone

• Multimesh DVC considering mechanical regularization, damage, and multimesh refinement to measure full-field deformation

• Very fine (sub-voxel) crack opening displacement fields via voxel-scale multimesh DVC to reveal damage initiation and growth

Introduction

Sandstone is a complex geotechnical engineering material featuring naturally occurring damage. It contains cracks and pores, which, in turn, may store gas or water [START_REF] Liu | Investigation of Sandstone Mesostructure Damage Caused by Freeze-Thaw Cycles via CT Image Enhancement Technology[END_REF]. Mechanical failure of sandstone is often accompanied by strain localization, which intrinsically is a multiscale spatiotemporal accumulation process of damage (i.e., crack initiation, propagation, coalescence leading to macroscopic brittle failure [START_REF] Paterson | Experimental rock deformation -The brittle field[END_REF]). In addition, growing fracture networks may alter the mechanical properties of sandstone when approaching failure [START_REF] Faulkner | Slip on "weak" faults by the rotation of regional stress in the fracture damage zone[END_REF]. Conventional laboratory tests such as uniaxial [START_REF] Bagde | Fatigue properties of intact sandstone samples subjected to dynamic uniaxial cyclical loading[END_REF] and triaxial compression [START_REF] Yang | Experimental investigation on the strength, deformability, failure behavior and acoustic emission locations of red sandstone under triaxial compression[END_REF] combined with Acoustic emission (AE) [START_REF] Yang | Spatial acoustic emission evolution of red sandstone during multi-stage triaxial deformation[END_REF][START_REF] Jia | Acoustic Emission Characteristics and Damage Evolution of Coal at Different Depths Under Triaxial Compression[END_REF] monitoring have been developed to investigate the strength, deformation behavior and failure characteristics of sandstone at the macroscopic scale. However, the correspondence between microstructure and damage growth cannot be observed in detail. Therefore, it is also desirable to investigate the mechanical response on smaller scales (e.g., micro-or mesoscales).

Numerous measurement techniques have been developed for the evaluation of rock deformation, such as the casting method [START_REF] Raven | Water flow in a natural rock fracture as a function of stress and sample size[END_REF], and surface topography measurements [START_REF] Pyrak-Nolte | Single fractures under normal stress: the relation between fracture specific stiffness and fluid flow[END_REF]. The growing use of digital cameras and correlation algorithms offered new ways to measure rock deformation, in which 2D/3D digital image correlation (DIC) techniques are developed to quantify non-destructively the deformation and localization of fracture on the 2D/3D surfaces of rock samples [START_REF] Zhang | Experimental investigation of deformation and failure mechanisms in rock under indentation by digital image correlation[END_REF][START_REF] Fakhimi | Insights on rock fracture from digital imaging and numerical modeling[END_REF][START_REF] Xing | Full-field measurement and fracture characterisations of rocks under dynamic loads using high-speed three-dimensional digital image correlation[END_REF]. Electronic speckle pattern interferometry (ESPI) has also been proved useful to estimate rock deformation [START_REF] Haggerty | Observing deformation and fracture of rock with speckle patterns[END_REF]. With the development of experimental techniques, X-ray computed microtomography (XCT) has emerged as one promising tool for observing rock microstructures. The early applications were mainly geared toward extracting relevant statistics on different components (e.g., holes, cracks and high-density minerals) before and after testing [START_REF] Kawakata | The observations of faulting in Westerly granite under triaxial compression by X-ray CT scan[END_REF], including volume fraction, spatial distribution and microstructural features to establish the relationships between structural changes and physical parameters [START_REF] Ramandi | Porosity and permeability characterization of coal: A micro-computed tomography study[END_REF][START_REF] Moscariello | Capillary collapse of loose pyroclastic unsaturated sands characterized at grain scale[END_REF]. At this stage, the image resolution was usually small (e.g., tens of micrometers) while the sample size was large. Rock deformation characteristics were described simply and directly, which was based on image processing (e.g., image segmentation and statistics).

With the improvement of spatial resolution, micro-focus XCT (i.e., μXCT) gives deeper insight into internal microstructures of rocks down to micro-or even nanoscales (Zhao et al., 2020). Damage development in rocks was further characterized by combining in-situ tests. [START_REF] Renard | Microscale characterization of rupture nucleation unravels precursors to faulting in rocks[END_REF] reported quantitative and high-resolution internal damage growth in sandstone undergoing brittle failure in triaxial compression. [START_REF] Huang | Synchrotron X-ray imaging in 4D: Multiscale failure and compaction localization in triaxially compressed porous limestone[END_REF] demonstrated that synchrotron μXCT could be used to directly image in 4 dimensions (3D in space plus time) to analyze the spatiotemporal growth of local porosity and damage at multiple scales in limestone deformed in triaxial compression. In addition, numerical modeling (e.g., pore network modeling method) based on 3D μXCT images is also an attractive way to predict the mechanical properties of sandstone [START_REF] Vajdova | Permeability evolution during localized deformation in Bentheim sandstone[END_REF][START_REF] Han | Numerical computation of elastic properties for porous rocks based on CT-scanned images using direct mapping method[END_REF].

For internal deformation field measurement of rocks, digital volume correlation (DVC) has proven to be effective when combined with in-situ/ex-situ tests [START_REF] Buljac | Digital Volume Correlation: Review of Progress and Challenges[END_REF], especially for the evaluation of mesoscopic and macroscopic strain fields at submillimeter resolution. DVC aims to measure displacement fields by registering reference and deformed volumes [START_REF] Bay | Digital volume correlation: Three-dimensional strain mapping using x-ray tomography[END_REF]. According to the size of the analyzed volume, DVC is divided into local approaches [START_REF] Smith | Digital volume correlation including rotational degrees of freedom during minimization[END_REF] (i.e., subsets independently registered) and global approaches [START_REF] Roux | Three-dimensional image correlation from X-ray computed tomography of solid foam[END_REF]Hild et al., 2011) (i.e., registration over the whole region of interest (ROI) based on finite element discretization). Regularized DVC was further developed to suppress spurious high frequency components contained in measured displacement fields by adding penalty terms (such as Tikhonov regularization [START_REF] Bell | Solutions of Ill-Posed Problems[END_REF], curvature penalization [START_REF] Van Dijk | A global digital volume correlation algorithm based on higher-order finite elements: Implementation and evaluation[END_REF], and mechanical regularization [START_REF] Leclerc | Voxel-Scale Digital Volume Correlation[END_REF][START_REF] Mendoza | Complete mechanical regularization applied to digital image and volume correlation[END_REF]). These two approaches have been applied to study strain concentrations preceding macroscopic failure of sandstone. [START_REF] Tudisco | Full-field Measurements of Strain Localisation in Sandstone by Neutron Tomography and 3D-Volumetric Digital Image Correlation[END_REF][START_REF] Tudisco | An extension of digital volume correlation for multimodality image registration[END_REF] [START_REF] Mao | An Improved Digital Volumetric Speckle Photography Technique with X-ray Microtomography and its Applications to Investigating Strain Localization in Red Sandstone[END_REF][START_REF] Mao | 3-D strain estimation in sandstone using improved digital volumetric speckle photography algorithm[END_REF] investigated the onset and development of strain localization in red sandstone using so-called digital volumetric speckle photography (i.e., local DVC) in different experiments such as indentation and uniaxial compression. [START_REF] Liu | An Adaptive and Reliable Guided Digital Volume Correlation Algorithm for Sandstone Based on 3D Scale-Invariant Feature Transform[END_REF] analyzed the localization of normal and shear strains of red sandstone subjected to Brazilian test preceding macroscopic failure using local DVC with adaptive subset sizes, thereby revealing progressive damage accumulation. [START_REF] Renard | Dynamic In Situ Three-Dimensional Imaging and Digital Volume Correlation Analysis to Quantify Strain Localization and Fracture Coalescence in Sandstone[END_REF] presented incremental strain concentrations and fracture coalescence in sandstone subjected to triaxial compression by using image segmentation and local DVC to reveal damage growth at microscopic and macroscopic scales. [START_REF] Hu | Pore and permeability properties of reservoir sandstone under a uniaxial compression CT test[END_REF] performed a uniaxial compression test with a reservoir sandstone and analyzed the gradual change of axial strains by using global DVC. Despite these successful DVC applications on sandstone to characterize damage growth, damage quantifications remained limited especially for microcrack propagation at very small scales (i.e., sub-voxel resolutions). The main challenge comes from the measurement of very rapid fluctuations of displacement fields combined with small-scale spatial resolutions since the limited available information will give rise to high measurement uncertainties [START_REF] Leclerc | Digital volume correlation: What are the limits to the spatial resolution?[END_REF]. In addition, no damage description was considered in DVC analyses to capture crack opening displacements.

In this work, a multiscale damage analysis of red sandstone has been carried out using in-situ testing and advanced (i.e., multimesh) DVC. The goal is to illustrate the potential of this new DVC approach in contributing to a better understanding of damage development in sandstone from sub-voxel to super-voxel levels. The addressed questions include: when and where do microcracks initiate? How do these microcracks propagate and coalesce with increased loading? What is the proportion and displacement amplitude of microcrack openings preceding failure? The outline of the paper is as follows. In Section 2, uniaxial compression, μXCT system, and red sandstone characteristics are described. In Section 3, the principle of FE-based global DVC is first introduced, then the implementation of multimesh DVC with mechanical regularization, damage variable, and mesh refinement strategy are described. Section 4 first reports the kinematic and damage fields of the sample in its post-failure state, then a back-to-front damage detection strategy is introduced to evaluate damage growth in pre-failure states at different scales (i.e., macro-, meso-, and microscales).

In situ uniaxial compression test

An in-situ uniaxial compression test was performed in a lab tomograph to analyze the damage mechanisms in sandstone. The cylinder (25 mm in diameter and 45 mm in height) made of red sandstone used herein was extracted from Jiang Xi Province, China. The sample mainly contained quartz (61 vol%), feldspar (31 vol%) and clay minerals (8 vol%) determined by X-ray diffraction. The main physical and mechanical properties were as follows: density of 2.07 g/cm 3 , porosity of 15.1 vol%, Young's modulus (E) equal to 3.3 GPa and Poisson's ratio 𝜈 = 0.18. The loading setup for in situ uniaxial compression is shown in Figure 1(a), in which the outer ring was made of polymethyl methacrylate (PMMA). The load and displacement of the upper platen were recorded by the load cell (20 ± 0.002 kN) and grating displacement meter (20 ± 0.001 mm), respectively. The maximum axial load was 15 kN and the stroke speed 0.18 mm/min.

[Insert Figure 1 here]

The hardware parameters of XCT are gathered in Table 1. Nineteen scans were performed, namely, two scans in the reference configuration for measurement uncertainty analyses, sixteen scans in the deformed configurations preceding failure, and the last scan for post-failure investigations, which was considered to define the crack path as discussed in Section 4. The axial stress/strain curve is plotted in Figure 1(b), where the corresponding axial stresses preceding failure were 0, 4. 1, 8.2, 9.2, 10.2, 11.2, 12.2, 13.2, 14.2, 15.3, 16.3, 17.3, 18.3, 19.3, 20.4, 21.4, and 24.0 MPa, respectively. When reaching the peak stress (27.0 MPa), the sample was completely broken. In each scan, the load was maintained for approximately 25 min, a series of 720 radiographs was acquired to reconstruct 3D images of the sample via Feldkamp cone-beam algorithm [START_REF] Yang | Locating of 2π-projection view and projection denoising under fast continuous rotation scanning mode of micro-CT[END_REF].

The load drops correspond to small force relaxations during scan acquisitions.

[Insert Table 1 here] Figure 2(a) shows 3D renderings (by Tomviz (Schwartz et al., 2012) software) of the volume of interest (VOI) of the sample in the reference configuration. The definition is 700 × 700 × 750 vx with a resolution of 46 μm / vx. No significant artifacts and initial defects are observed in the VOI, and the contrast was uniformly distributed except for several bright spots with high gray levels mostly associated with minerals. Figure 2(b) shows the corresponding gray level histogram with the skewness. The dynamic range of the histogram was approximately 4000 GL (3000-7000 GL), which was limited compared to the 16-bit digitization. The skewness is positive, thus indicating that the histogram was not an ideal Gaussian distribution and the overall trend expanded to the right.

[Insert Figure 2 here]

To compute the equivalent speckle size, the normalized autocorrelation [START_REF] Minguet | On the validation of a priori estimates of standard displacement uncertainties in T3-stereocorrelation[END_REF] was computed by fast Fourier transform for the VOI in the reference configuration. Figure 3(a) shows the normalized autocorrelation of the VOI (zoom about the correlation peak). The speckle sizes (d) were then determined as the mean full width at half maximum (FWHM) along the three perpendicular directions (Figure 3(b)). Along the xand y-directions, they have a similar level (i.e., ≈2.7 vx), but in the z-direction, it is slightly lower (i.e., 2.4 vx), which may be related to rock compaction under gravity. Another way to evaluate the speckle quality is through their contrast, namely, the root-mean-square (RMS) gray level gradients [START_REF] Leclerc | Digital volume correlation: What are the limits to the spatial resolution?[END_REF]. The values in the three directions are 135, 134 and 138 GL / vx, respectively. According to [START_REF] Leclerc | Digital volume correlation: What are the limits to the spatial resolution?[END_REF], a larger RMS image gradient gives rise to smaller measurement uncertainties. Therefore, the measurement uncertainty along z is expected to be slightly lower than that in the other two directions (i.e., 𝜎 𝑧 < 𝜎 𝑥 ≅ 𝜎 𝑦 ).

[Insert Figure 3 here]

For a better understanding about damage features, 3D renderings of the VOI and middle slices along three perpendicular directions in the post-failure configuration are shown in Figure 4. Severe large-scale cracks are observed in the post-failure scan. In particular, a small piece fell off on the right edge of the sample, which is challenging for DVC measurements. It has to be emphasized that no visible cracks appear on the 3D images prior to macroscopic (brittle) failure, which not only points out the rapid damage process but also the challenges of damage detection and quantification.

[Insert Figure 4 here]

Multimesh DVC analyses

1 Principle of regularized DVC

Finite element DVC [START_REF] Roux | Three-dimensional image correlation from X-ray computed tomography of solid foam[END_REF][START_REF] Gomes-Perini | A multigrid PGD-based algorithm for volumetric displacement fields measurements[END_REF] was applied to measure displacement fields from volumetric images acquired during the in situ experiment. Considering the gray level conservation between the volume in the reference state 𝑓 and deformed configuration 𝑔, the correlation residual 𝜌 𝑟 reads

𝜌 𝑟 (𝒙, {𝒗}) = 𝑓(𝒙) -𝑔[𝒙 + 𝒖(𝒙, {𝒗})] (1)
where the sought displacement field 𝒖 is discretized using 8-noded rectangular hexahedra (i.e., C8 elements

𝒖(𝒙) = ∑ 𝚿 𝑖 𝐶8 (𝒙)𝑣 𝑖 𝑖
, where 𝚿 𝑖 𝐶8 are the shape functions corresponding to the nodal displacements 𝑣 𝑖 ). The unknowns then become the nodal displacements 𝑣 𝑖 that are gathered in the column vector {𝒗} . A Gauss-Newton [START_REF] Hild | Toward 4D mechanical correlation[END_REF] iteration scheme was used to search for the minimum of the L2-norm of the global correlation residual (i.e., Φ 𝑐 2 ({𝒗})) until reaching a convergence criterion written in terms of the norm of corrections to {𝒗}. In regularized DVC, linear elasticity is considered in an FE sense by introducing a penalty term based on the equilibrium gap [START_REF] Claire | A finite element formulation to identify damage fields: The equilibrium gap method[END_REF])

[𝐊]{𝒗} = {𝐟} (2)
where [𝐊] is the rectangular stiffness matrix restricted to inner nodes and traction-free surfaces, and {𝐟} the nodal force vector that must vanish in the absence of body forces. The L2-norm of the force residuals then reads

Φ 𝑚 2 ({𝒗}) = {𝒗} 𝐓 [𝐊] 𝐓 [𝐊]{𝒗} (3)
Taking into account these two cost functions leads to the minimization of their weighted sum

{𝒗 𝒎 } = arg min ∑ [Φ 𝑐 2 ({𝒗}) + 𝑤Φ 𝑚 2 ({𝒗})] VOI (4)
The cost functions are made dimensionless by considering a trial displacement field associated with a plane wave of vector 𝜅 [START_REF] Leclerc | Voxel-Scale Digital Volume Correlation[END_REF]. The regularization weight then is proportional to the so-called regularization length ℓ 𝑟𝑒𝑔 raised to the power 4 [START_REF] Leclerc | Digital volume correlation: What are the limits to the spatial resolution?[END_REF] 𝑤 ∝ (2𝜋|𝜅|ℓ 𝑟𝑒𝑔 ) 4

(5)

Within such regularization scheme, there are two characteristic length scales. The first one is associated with the element size ℓ of the mesh, and the second one is ℓ reg . When ℓ reg < ℓ, the displacement fluctuations (and their corresponding uncertainties) are controlled by the finite element discretization.

Conversely, when ℓ reg ≥ ℓ, the displacement fluctuations are required to be mechanically admissible over a domain whose size is proportional to ℓ reg . Thanks to regularized DVC, it was shown that voxel-scale registrations could be performed [START_REF] Leclerc | Voxel-Scale Digital Volume Correlation[END_REF]2012).

2 Damage model

The previous mechanical regularization scheme leads to good approximations when linear elasticity is relevant. However, it no longer applies for elements traversed by cracks [START_REF] Hild | Damage measurements via DIC[END_REF]. The corresponding mechanical model is then corrected by introducing a damage variable D ranging from 0 to 1 (i.e.,

[𝐊 𝐷 𝐞 ]=[𝐊 𝐞 ](1 -𝐷)
, where [𝐊 𝐞 ] and [𝐊 𝐷 𝐞 ] are the elementary stiffness matrices before and after correction), which reduces the regularization weight by a factor √1 -𝐷. Ideally, for brittle damage of sandstone, D should be equal to 1 to fully disable the regularization within individual elements. However, due to limited contrast within very small elements, a low regularization weight was still needed to allow for DVC convergence.

Therefore, D tends to be close to 1 as much as possible. In the present case, D was greater than 0.95, which not only ensured DVC convergence, but also allowed high displacement gradients to be captured. Damage zones thus need to be identified. Gray level residuals (GLRs) [START_REF] Hild | Damage measurements via DIC[END_REF][START_REF] Vrgoč | Damage characterization in fiber reinforced polymer via Digital Volume Correlation[END_REF] were shown to be an effective indicator, especially for voxel level crack openings.

3 Refinement and damage relocalization

Considering computational efficiency and measurement uncertainty, a coarse mesh (see cyan mesh in Figure 5) was first used to obtain a rough damage identification. To capture finer damage zones (i.e., morphology of cracks), refinement of the coarse mesh is desirable [START_REF] Sciuti | Pyramidal adaptive meshing for Digital Image Correlation dealing with cracks[END_REF], see blue mesh in Figure 5.

Then, the damage detection quantity (i.e., GLR) was used again in the fine mesh to obtain better resolved crack paths (see red mesh in Figure 5).

[Insert Figure 5 here]

4 Multimesh regularized DVC with damage

Once the above analyses were completed, multimesh and regularized DVC with damage was run. A Master-Slave elimination method [START_REF] Liu | An Adaptive and Reliable Guided Digital Volume Correlation Algorithm for Sandstone Based on 3D Scale-Invariant Feature Transform[END_REF] was used to deal with hanging nodes, which are defined as nodes with different numbers of connectives located on a surface between two elements (black triangle in Figure 5).

The flowchart of the advanced DVC analysis is shown in Figure 6, which includes three steps:

Step 1-global analyses: regularized DVC (coarse mesh) is performed and a damage variable (D) associated with GLR-based damage detection is considered.

Step 2-local analyses: mesh refinement and damage relocalization are carried out for each damaged element.

Step 3-global analyses: multimesh DVC combined with damage is finally run to analyze deformation fields in the presence of cracks.

[Insert Figure 6 [Insert Figure 7 here] Similar trends are observed for the fitted power laws (𝜎 ∝ ℓ 𝑟𝑒𝑔 𝑏 ), namely, the larger the regularization length, the lower the standard displacement and strain uncertainties [START_REF] Leclerc | Digital volume correlation: What are the limits to the spatial resolution?[END_REF]. This general trend was expected because larger regularization lengths correspond to higher mechanical regularization weights, thereby filtering out more high frequency components in deformation fields. The power law exponents of 𝜎 𝑢 and 𝜎 𝜖 are approximately -1.7 and -2.2, respectively, which are close to theoretical predictions (i.e., 𝑏 = -3 2 ⁄ for 𝜎 𝑢 and 𝑏 = -5 2 ⁄ for 𝜎 𝜖 ) [START_REF] Leclerc | Voxel-Scale Digital Volume Correlation[END_REF]. It is worth emphasizing that too large regularization lengths may give rise to excessive boundary effects [START_REF] Leclerc | Digital volume correlation: What are the limits to the spatial resolution?[END_REF]. In addition, the uncertainty level along the z-direction is smaller than the other two due to better contrast, which is consistent with analyses in terms of RMS gray level gradients. The difference for the three principal strains shows that the maximum (or minimum)

principal strains are always more fragile, thus have higher uncertainties.

Registration quality

Two residuals are utilized to evaluate the quality of image registration and mechanical equilibrium, namely, (5)

𝜌 𝑚𝑒𝑐 = √{𝒗} T [𝐊] T [𝐊] {𝒗} (6)
Figure 8 shows the changes of these two residuals for different pre-failure states (except the post-failure state for which the 'relaxation' process stopped at ℓ 𝑟𝑒𝑔 = 60 vx due to severe damage). Similar trends are observed from these results, namely, (1) an increase of ℓ 𝑟𝑒𝑔 reduced 𝜌 𝑚𝑒𝑐 because more weight is put on mechanical regularization, thus the kinematic description is more constrained leading to larger 𝜌 𝑑𝑣𝑐 ;

(2) both residuals increase with the axial stress for the same regularization length ℓ 𝑟𝑒𝑔 especially when ℓ 𝑟𝑒𝑔 > 60 vx. Because the linear elastic assumption is gradually less satisfied by the true complex deformation process, the registration quality decreases accordingly;

(3) a distinct inflection point is observed for each applied stress, whose regularization length (i.e., ℓ 𝑟𝑒𝑔 = 60 vx) is considered a suitable choice to balance the two residuals. For lower regularization lengths, the DVC residuals are approximately independent of ℓ 𝑟𝑒𝑔 . Conversely, for large regularization lengths, the two residuals have opposite trends (i.e., significant decrease in 𝜌 𝑚𝑒𝑐 induces a sharp increase of 𝜌 𝑑𝑣𝑐 ); (4) the curve for an axial stress of 4.1 MPa is essentially vertical, which suggests that the sample underwent almost pure elastic deformations without too much deviation from mechanical admissibility.

The residual 𝜌 𝑑𝑣𝑐 changes remain very limited. It is worth noting that the deformation of smaller elements is more likely to be linear and elastic (if undamaged). Therefore, a larger mechanical regularization is expected to act effectively as a low-pass filter; ℓ 𝑟𝑒𝑔 = 60 vx is still used in analyses with refined meshes.

[Insert Figure 8 here]

Coarse mesh

A first coarse mesh (with a mean length ℓ = 21 vx) is considered for regularized DVC according to the sample geometry, as shown in Figure 9(a,b). To help convergence, a 'relaxation' process [START_REF] Tomičevć | Mechanics-aided digital image correlation[END_REF] associated with the regularization length (i.e., ℓ 𝑟𝑒𝑔 = 540, 180, 60, 20, 10, 0 vx) is followed in the DVC calculations. The measured x-direction displacement field of the last scan before failure with a regularization length of 60 vx is shown in Figure 9(c). The corresponding displacement uncertainty is less than 0.01 vx (or 0.28 μm) according to the uncertainty analysis of Section 4.1.1 (Figure7(a)). A displacement amplitude of about ±50 μm (or 1.1 vx) occurs in the symmetrical edge region, which illustrates that the deformation before failure is very limited, and damage detection is challenging in the current deformed state.

[Insert Figure 9 here]

Damage identification in post-failure state

Gray level residuals in coarse mesh

Given the fact that no visible cracks are seen in the pre-failure scans, the crack path is determined with the post-failure scan in which cracks fully developed. First, the highest gray level residuals label cracks and highlight where the DVC registration failed. Second, the registration in ROI #1 is deemed trustworthy because it has the lowest gray level residuals except for regions near cracks (due to regularization). On the contrary, higher gray level residuals are observed in ROIs #2 and #3, where untrustworthy registrations occur. It is worth noting that the chip (i.e., right of the blue dotted line) was removed in the subsequent DVC calculations because no grayscale information was available.

[Insert Figure 10 here]

Segmented DVC in coarse mesh

A segmented DVC strategy is considered to deal with the registrations of ROIs #2 and #3. ) is then evaluated based on the mesh size and the maximum principal strain (i.e., COD ≈ 𝜖 1 ℓ). The quantified mean crack opening is approximately 900 μm (or 19.6 vx). However, only a rough match is observed between the crack opening displacement field and gray level residual field due to limitations of the coarse mesh.

[Insert Figure 11 here]

Gray level residual with fine mesh

The detected damaged elements with the coarse mesh are further refined to 5 vx elements, and local regularized DVC is performed to obtain a finer damage localization and initial nodal displacements. Figure 12(a)

shows the histogram of gray level residuals for the fine mesh, and the thresholding method (i.e., threshold = 1300 GL) is utilized to detect new damaged elements (Figure 12(b)). The threshold used here is higher than that used in the coarse mesh (i.e., threshold = 850 GL) due to significant crack openings in the fine mesh.

[Insert Figure 12 here]

Multimesh DVC Analysis

Construction of multimesh

The new finite element table of connectivity is generated according to the mesh information (i.e., coarse mesh, fine mesh, and damage mesh) for multimesh DVC calculations (ℓ 1 = 21 vx, ℓ 2 = 5 vx and ℓ 𝑟𝑒𝑔 = 60 vx, D = 0.96). It is worth noting that the damage variable D cannot be set to 1 due to poor contrast in the fine elements. To reduce the computation cost, a contiguous sub-volume containing the cracked surfaces is considered, in which three types of nodes are observed, namely, master (red) nodes in the coarse mesh, slave (blue) nodes, and damaged (green) nodes in the fine mesh (Figure 13). Since damaged nodes are splitable, they should be completely surrounded by slave nodes.

[Insert Figure 13 here] [Insert Figure 14 here]

Crack opening displacement fields

Multiscale damage analyses before failure

Three different characteristic scales are distinguished to describe the details of damage growth in red sandstone under uniaxial compression. At the macroscale, the ROI approaches the entire specimen with a volume of tens of cm 3 , and mean strains are used to evaluate the overall damage process. At the mesoscale, the ROI is based on the FE-mesh with a volume scale of tens of mm 3 , and element-wise strains are estimated from crack opening displacement fields are deduced. At the microscale, the spatial resolution is generally close to the voxel-scale (i.e., of the order of 50 μm). Therefore, very detailed deformations of the considered body are available.

Macroscale study

Macroscale damage is estimated with the mean strains (i.e., maximum principal strain 𝜖 1 , minimum principal strain 𝜖 3 , volumetric strain 𝜖 v ) over the whole ROI of size 32.2 × 32.2 × 34.5 mm 3 (or 700 × 700 × 751 vx), in which the volumetric strain is defined as the sum of the three principal strains 1 . Figure 15 displays the three mean strains measured by regularized DVC (ℓ = 21 vx and ℓ 𝑟𝑒𝑔 = 60 vx) as functions of the applied stress 𝜎 z . The volumetric strain reveals a similar trend with conventional results (Závacký, et al., 2021) and four deformation stages are distinguished. In stage 1 (0 < 𝜎 z < 16.3 MPa), 𝜖 v decreases as 𝜎 z increases, which demonstrates the sample was gradually compressed and 𝜖 v is mainly controlled by the minimum (i.e., negative) principal strain 𝜖 3 . The two principal strain curves approximately exhibit a linear trend, which indicates that no extensive damage occurred in the sample. For stage 2 (16.3 MPa < 𝜎 z < 20.4 MPa), 𝜖 v gradually increases but still remains negative and the two principal strain curves continue to grow but with a nonlinear trend. Such trend indicates the presence of microcracks inside the sample, and 𝜖 1 contributes more to 𝜖 v . In stage 3 (20.4 MPa < 𝜎 z < 24.0 MPa), 𝜖 v becomes positive, which illustrates the growth of microcracks that leads to volume expansion (i.e., dilatancy) of the sample. Nonlinear trends are observed for the two 1 The following convention is used. Positive strains refer to expansion, and negative strains to contraction. principal strain responses. Stage 4 corresponds to the post-failure (unloaded) state of the sample for which the maximum principal strains and volumetric strain experience large increments due to the development of large macrocracks.

[Insert Figure 15 here]

Mesoscale analysis

Mesoscale damage is related to element-wise (i.e., 0.23 mm in the fine mesh) crack opening displacement fields measured by multimesh DVC with the same settings as those used in the post-failure state (Figure 13). There are two reasons for using the post-failure state to assess the damage state. First, it is very difficult to detect microcracks in pre-failure stages especially for small spatial resolutions due to the brittleness of rocks.

Therefore, back-to-front damage detection in time sequence can accurately characterizes the initiation and propagation of microcracks. Second, crack openings are expected to be monotonic with increasing loading.

Therefore, the complete crack path can be provided by the post-failure state in which the macrocracks are easily detected as they remained open. Back-to-front damage detection combined with measurement uncertainties enables damage growth to be quantified at earlier stages using small spatial resolutions.

The COD threshold is defined as a factor k times the standard COD uncertainty (i.e., 𝜎 𝐶𝑂𝐷 = 𝑘𝜎 𝐶𝑂𝐷 𝑠𝑡𝑑 , with 𝜎 𝐶𝑂𝐷 𝑠𝑡𝑑 = 0.03 vx or 1.4 μm) assessed with the repeated scan of the reference configuration. The entire crack propagation process is revealed by analyzing crack opening displacement fields after thresholding with different COD thresholds (i.e., k = 3, 4, 5, 6, see Figure 16). It has to be emphasized that the total CODs in the cracked regions should be estimated by the sum of CODs of adjacent damaged elements. The larger the COD threshold, the fewer microcracks appear for the same deformed state. In particular, for the last scan before failure, two microcracked surfaces are observed when factor k is less than 5. However, when factor k is greater than 5, only one microcracked surface appears. Given the fact that the beginning of stage 2 (Figure 15) corresponds to the initiation of microcracks, it should not contain too many microcracks. Therefore, k = 5 is selected to assess the COD detection threshold (i.e., 𝜎 𝐶𝑂𝐷 = 0.15 vx or 7 μm).

[Insert Figure 16 here]

The mean and standard deviation of CODs (i.e., μ 𝐶𝑂𝐷 and S 𝐶𝑂𝐷 ), and the damaged volume (V 𝐷 ) and its volume fraction (V 𝐶𝑂𝐷 ) defined as the percentage of damaged volume are gathered in Table 2. When loaded to 17.3 MPa, scattered microcracks initiate at the bottom of the sample with a mean COD equal to 20.8 μm and a volume faction of 0.22% (which is the beginning of stage 2 shown in Figure 15). These microcracks gradually coalesce and propagate vertically with a mean COD of 24.1 μm and a volume fraction of 6.3% when reaching an axial stress of 20.4 MPa (i.e., end of stage 2, see Figure 15). Beyond this stress level, the microcracks propagate rapidly with a mean COD ≈ 28.2 μm and a volume fraction ≈ 36.5%. Only one front microcrack is observed but without the rear one, which presumably initiated during final failure.

[Insert Table 2 here]

Microscale analysis

While mesoscale DVC calculations were carried out successfully to evaluate damage growth of the main cracks, the assessment of the fine cracks is not sufficiently detailed. Therefore, it was necessary to run voxel-scale multimesh DVC to analyze microscale damage that is related to voxel-wise (i.e., 1 voxel ≡ 46 μm) previous ones due to the COD uncertainty (0.17 vx or 7.8 μm) and mesh size limitation. In addition, even though the morphology of cracks at the mesoscale is rather fine, it can be further refined via microscale (or voxel-scale) analyses.

GLRs
[Insert Figure 18 here]

The detailed implementation steps of voxel-scale multimesh DVC are as follows:

Step Step 2: A new FE-mesh of size 5 vx is produced on the ROI. Regularized DVC (ℓ = 5 vx, ℓ 𝑟𝑒𝑔 = 60 vx, 𝐮 0 ← 𝐮 1 ) is run with the initial guess 𝐮 1 to obtain nodal displacements {𝐮 2 } and the corresponding GLR field (i.e., GLR 2 ) as shown in Figure 19(b).

Step 3: Damage detection is carried out according to GLR 2 (red mesh in Figure 19(c)). Regularized DVC is rerun with damage (ℓ = 5 vx, ℓ 𝑟𝑒𝑔 = 60 vx, D = 0.99, 𝐮 0 ← 𝐮 2 ) to obtain 𝐮 3 and GLR 3 .

Step 4: A new FE-multimesh is generated on the ROI, namely, 5 vx-mesh for undamaged zone and 1 vx-mesh (Figure 19(d)) for damaged elements. Multimesh nodal displacements 𝐮 4 are interpolated based on 𝐮 3 .

Step 5: Generate reduced FE-multimesh by taking out the 1 vx-elements belonging to any damaged zone.

Run voxel-scale multimesh DVC (ℓ 1 = 5 vx, ℓ 2 = 1 vx, ℓ 𝑟𝑒𝑔 = 60 vx, 𝐮 0 ← 𝐮 4 ) to calculate 𝐮 5 and GLR 5 . It is observed that this DVC calculation leads to trustworthy results according to low gray level residuals (Figure 19(e)).

Step 6: Grow one layer of 1 vx-elements along both sides of the cropped zone. Rerun voxel-scale multimesh DVC (ℓ 1 = 5 vx, ℓ 2 = 1 vx, ℓ 𝑟𝑒𝑔 = 60 vx, 𝐮 0 ← 𝐮 5 ) to calculate 𝐮 6 and GLR 6 (Figure 19(f)).

Low gray level residuals are observed except for several limited regions, which indicates that almost no damage occurred in the considered meshes and the DVC calculations are trustworthy.

Step 7: Continue growing one layer of 1 vx-elements as above and repeat Step 6 but with the initial guess of 𝐮 6 to calculate 𝐮 7 and GLR 7 (Figure 19(g)). This new DVC calculation converged after only a few iterations. Higher gray level residuals occur in the added layer, which indicates that it does contain damage. In particular, the cropped region is very fine and close to 1 vx width. Therefore, Step 7 can be stopped.

Step 8: Full FE-multimesh displacement fields 𝐮 8 are interpolated based on 𝐮 7 , and thus GLR 8 is calculated in which finer damage localization is observed (Figure 19(h)).

Step 9: Damage detection is carried out again according to GLR 8 (red mesh in Figure 19(i)). Multimesh DVC is run with damage (ℓ 1 = 5 vx, ℓ 2 = 1 vx, ℓ 𝑟𝑒𝑔 = 60 vx, D = 0.99, 𝐮 0 ← 𝐮 8 ) to assess the final nodal displacements {𝐮 𝑓 } and the corresponding residuals GLR 𝑓 .

[Insert Figure 19 here]

The same DVC settings as described in Figure 19(i) are applied to the pre-failure scans to assess damage growth. The COD threshold is equal to 4.6 μm (or 0.1 vx), which is defined as 5 times the standard COD uncertainty. Table 3 lists the results of the corresponding damage quantification. In addition, thresholded COD fields are plotted around the COD curve vs. axial stress (Figure 20). When loaded to 18.3 MPa, a few scattered microcracks initiate with a mean COD of 5.1 μm (or 0.11 vx), which is very close to the detection threshold, and only a damage volume fraction of 2% is reached. Two microcracks (μ 𝐶𝑂𝐷 = 6.0 μm or 0.13 vx, V 𝐶𝑂𝐷 = 14%) are observed when reaching an axial stress of 20.4 MPa. These two microcracks propagate and gradually coalesce into a unique mesocrack when loaded to 24.0 MPa. At this point, the damage volume fraction is 48% with a mean COD of 7.5 μm (or 0.16 vx). For the post-failure state, the voxel-level cracks with opening of 40.9 μm (or 0.98 vx) appear. The voxel-scale COD fields reveal damage initiation and propagation very clearly.

[Insert Table 3 here] [Insert Figure 20 here]

Conclusion and perspective

In this work, an in situ uniaxial compression test applied to red sandstone was carried out based on μXCT imaging. Mechanically regularized multimesh DVC considering damage was used to evaluate damage growth at three different scales. In particular, voxel-scale multimesh DVC was implemented to study the early stages of damage inception.

From the macroscopic point of view, the overall deformation of the sample was assessed via mean strain curves (i.e., 𝜖 1 , 𝜖 3 , 𝜖 v ) measured over the whole ROI. It was shown that damage had set in for an applied stress equal to 61% (i.e., 16.3 MPa) of the ultimate compressive strength. Microcracks gradually developed leading to dilatancy when the applied stress was greater than 76% (i.e.,20.4 MPa) of the ultimate strength. For the post-failure state (i.e., beyond the ultimate stress), macroscale damage occurred suddenly in the form of two longitudinal cracks.

The estimation of element-wised crack opening displacements and crack morphology was carried out and gave more insight into damage growth at the mesoscale, especially for sub-voxel crack opening displacements (CODs). Microcracks initiated from the bottom of the sample, then gradually coalesced and propagated vertically. Only one cracked surface located in the front of the sample was observed for the last scan before failure with a mean COD ≈ 28.2 μm and a volume fraction ≈ 36.5%.

The analysis of the gray level residuals and COD fields at the microscale by voxel-scale multimesh DVC enabled damage growth to be characterized and quantified in a very fine way. It provided further understanding of the damage chronology. This study thus demonstrates the potential of multimesh DVC for the quantification of damage at very small scales (i.e., down to the voxel scale).

The significance of this study lies in two points. First, only visible cracks (i.e., CODs greater than one voxel) are generally extracted from CT images. Therefore, sub-voxel openings were hard to capture before the peak loading because of the brittleness of sandstone. However, sub-voxel crack initiation and propagation were clearly measured at earlier stages by the proposed DVC procedure, which is very rare for this type of research.

Second, the kinematic and damage fields of sandstone were quantified for a very fine mesh (i.e., at the voxel-scale), which is impossible in conventional DVC calculations with a spatial resolution of tens of voxels at the least.

The limitation of this study stems from two aspects. First, uniaxial compression is not a strictly continuous process because the applied force is maintained for each scan. Therefore, fast damage growth of sandstone cannot be captured. Second, the computational cost and speed of multimesh DVC are too large and slow for the complete ROI once finer discretization (e.g., at the voxel scale) was considered for cracked regions. Given the above limitations, future work will try to design and perform in situ tests with continuous loading and on-the-fly acquisitions. The corresponding measurement technique will need to be further improved (e.g. projection-based DVC, [START_REF] Kosin | A projection-based approach to extend digital volume correlation for 4D spacetime measurements[END_REF] to measure bulk deformations and damage in rocks.

Figure Captions

  determined 3D strain fields of sandstone during triaxial experiments using local DVC based on neutron tomography and XCT. Their results clearly demonstrated strain localization in the cracked zones. Stappen et al. (2022) used the same algorithm to visualize grain-scale deformation mechanisms of sandstone in triaxial compression with 4D μXCT images.

  The standard displacement and strain uncertainties were estimated based on two scans in the reference configuration (i.e., scans ref and #0, see Figure1(b)), which is assessed as the standard deviation 𝜎 of the nodal displacements (i.e., 𝜎 u 𝑥 , 𝜎 u 𝑦 , 𝜎 u 𝑧 ) and principal strains (i.e., 𝜎 𝜖 1 , 𝜎 𝜖 2 , 𝜎 𝜖 3 ) measured by regularized DVC with a mesh size ℓ of 21 voxels (Figure7(a,b)).

  the dimensionless correlation residual 𝜌 𝑑𝑣𝑐 and mechanical residual 𝜌 𝑚𝑒𝑐

  Figure 10(a) shows the voxel-wise gray level residual field measured by regularized DVC (ℓ = 21 vx and ℓ 𝑟𝑒𝑔 = 60 vx). Figure 10(b) displays the corresponding histogram with two peaks. It also shows the gray level residuals in the xOy middle slice (z = 16 mm), which is divided into three smaller regions of interest (i.e., ROI #1, ROI #2 and ROI #3) separated by cracks. A ternary operation is performed on the same slice based on thresholding (i.e., GLR <350 → 0, 350≤GLR <850 → 1, GLR ≥850 → 2, Figure 10(c)). The trinary slice clearly reveals the registration quality of regularized DVC.

  Figure11(b,c). Compared to Figure10(b,c), the gray level residuals in these three regions are reduced except for the cracked regions. In particular, cracks are further recovered and more damage details appear (e.g., two cracks inside the white ellipse in Figure11(b) while only one is present in Figure10(b)). The element-wise crack opening displacement (COD) field (see Figure11(d)) is then evaluated based on the mesh size and the

Figure 14

 14 Figure 14 displays different fields characterizing damage (i.e., displacement field in the x-direction, gray level residual field and COD field). Very high displacement gradients (Figure 14(a)) are observed in the cracked regions. Local displacement jumps are greater than 800 μm (or 17.4 vx). The localized zones in the gray level residuals (Figure 14(b)) correspond to the damaged regions. The crack shape is better and finely captured compared to Figure 11(a), which illustrates the benefit of multimesh DVC. The element-wise crack opening displacements (Figure 14(c)) quantify the damaged zones. The total opening in cracked regions is estimated

  and voxel-wise CODs. According to the post-failure GLR field (Figure17(a)) assessed in Section 4.2.2, a small ROI of size 120 × 175 × 10 vx (or 5.5 × 8.1 × 0.5 mm 3 ) containing one fine crack is chosen (Figure17(b,c)) considering the computational cost.[Insert Figure17here] Thresholded mesoscale COD fields in this small ROI measured by regularized DVC (ℓ = 5 vx, ℓ 𝑟𝑒𝑔 = 60 vx and D = 0.99) are shown in Figure 18. Localized zones are observed only in the last three scans but not in

1 :

 1 Reference ROI and corrected deformed ROI in the post-failure state are shown in Figure 19. The displacement field 𝐮 1 for the deformed volume correction is measured by regularized DVC (ℓ = 21 vx, ℓ 𝑟𝑒𝑔 = 60 vx, D = 0.98).

Fig. 1 .Fig. 5 Fig. 6 Fig. 8 Fig. 9

 15689 Fig. 1. In situ setup for uniaxial compression (a) and axial stress/strain curve (b) depicting the 19 scans acquisitions of the in situ test Fig 2. (a) 3D renderings of the VOI in the reference configuration. (b) Corresponding gray level histogram with skewness Fig 3. (a) Normalized autocorrelation of the VOI in the reference configuration. (b) Characteristic speckle sizes along three perpendicular directions determined by its FWHM Fig 4. 3D renderings of the post-failure configuration. (a) VOI. and (b) middle slices (i.e., xOy plane for z = 16 mm, yOz plane for x = 15.5 mm, xOz plane for y = 15.5 mm) Fig. 5 Example of mesh refinement and damage relocalization Fig. 6 Flowchart of multimesh DVC analyses in the presence of cracks Fig 7. Standard displacement (a) and strain (b) uncertainties as functions of regularization length. The solid lines depict power-law interpolation with exponent b Fig. 8 Mechanical residual 𝜌 𝑚𝑒𝑐 as a function of dimensionless correlation residual 𝜌 𝑑𝑣𝑐 for different applied loads and regularization lengths Fig. 9 Coarse mesh shown in xOy (a) and xOz (b) planes for the reference configuration. (c) Measured displacement field (ℓ 𝑟𝑒𝑔 = 60 vx) in the x-direction for the last pre-failure state Fig. 10 Gray level residual fields in the post-failure state. (a) 3D rendering (values less than 850 GL are transparent). (b) Gray level residuals in the z-direction middle slice (see red dotted circle in subfigure(a)) and its histogram with thresholding values. (c) Trinary image after thresholding Fig. 11 Segmented DVC results. (a) Segmentation and DVC analyses in two individual regions (i.e., ROIs #2 and #3). (b) Gray level residuals in the z-direction middle slice measured by regularized DVC (D = 0.98). (c) Trinary image after thresholding. (4) Crack opening displacement field in the same slice.

Fig. 12

 12 Fig. 12 Damage detection in the fine mesh. (a) Histograms of elementary GLRs and selected threshold. (b) 3D rendering of damaged elements Fig. 13 Schematic view of nodes in multimesh DVC. (a) 3D rendering. (b) Middle xOy slice (z = 15 mm). The red nodes depict masters, blue nodes are slaves, and damaged nodes are green Fig. 14 Damage related fields measured by multimesh DVC. The top row shows 3D renderings and the bottom row the corresponding middle xOy slice. (a) Displacement along x-direction (in μm). (b) Gray level residuals. (c) Crack opening displacement (in μm) Fig. 15 Mean measured strains as functions of axial stress and the four deformation stages Fig. 16 Thresholded mesoscale COD fields for different axial stresses and with four different COD thresholds

  Fig. 1

Fig. 12

 12 Fig. 12

  Fig. 13

Fig. 15

 15 Fig. 15

Fig. 16

 16 Fig. 16

Fig. 17

 17 Fig. 17

Fig. 18

 18 Fig. 18

Fig. 19

 19 Fig. 19

Fig. 20

 20 Fig. 20

Table 2 .

 2

	axial stress (MPa)	17.3 18.3 19.3 20.4 21.4 24.0
	V 𝐷 (mm 3 )	61	107 315 1219 3129 7339
	V 𝐶𝑂𝐷 (%)	0.22 0.50 2.0	6.3	14.5 36.5
	μ 𝐶𝑂𝐷 (μm)	20.8 21.2 23.8 24.1 25.1 28.2
	𝑆 𝐶𝑂𝐷 (μm)	2.8	2.8	3.2	4.1	4.2	6.1

Funding

This work was financially supported by the National Key Research and Development Program of China (2022YFC2904102), Science Fund for Creative Research Groups of the National Natural Science Foundation of China (52121003), Fundamental Research Funds for the Central Universities of China (2022YJSMT04), and the China Scholarship Council.

Availability of data and material

The raw/processed data required to reproduce these findings will be available upon request.

Code availability

The code will be available upon request.

Declarations

Conflicts of interest

The authors declare no conflict of interest.

Author Contributions