Maria Repina

Didier Renard

Frédéric Bouyer
email: frederic.bouyer@cea.fr

Vincent Lagneau

Coupling image analysis and thermo-mechanical simulation results to produce a model of the fracture network in a nuclear glass canister

Keywords: Nuclear waste, fracture network, image analysis, thermo-mechanical simulation, equivalent model, nuclear glasses

Introduction

Vitrification of radioactive waste has been acknowledged as a feasible and reliable application to condition high-level long-lived nuclear fission products in geological disposal facilities on time scales of up to millions of years. This process involves calcination of liquid fission product solutions followed by 2 a melting stage at around 1100 °C in a metallic vessel in which the molten glass and the calcinated fission product solution are directly heated by electromagnetic induction before being poured into a metallic container. As the glass cooling and solidification occur, the stress increase stemming from the mechanical equilibrium between the internal liquid glass and external solidified glass causes the development of a fracture network. The thermal gradient is very important at the beginning of the cooling along the walls of the stainless steel container that results in the creation of a dense fracture network. After, when the gradient fades, a less dense radially growing fracture network creates in the center of the block.

The chemical performance and alteration of glass under leaching has been extensively studied during the last thirty years [START_REF] Pierce | Modeling Interfacial Glass-Water Reactions: Recent Advances and Current Limitations[END_REF][START_REF] Carl | Micro-Continuum Approaches for Modeling Pore-Scale Geochemical Processes[END_REF][START_REF] Gin | An international initiative on long-term behavior of high-level nuclear waste glass[END_REF]. The results of these experiments have provided a better understanding of the mechanisms of glass alteration at the micro level and facilitated elaboration of the glass reactivity with allowance for the alteration layer (GRAAL) model [START_REF] Frugier | SON68 Nuclear glass dissolution kinetics: Current state of knowledge and basis of the new GRAAL model[END_REF]. The modeling of flow and transport phenomena in fractured media is a topical issue in Earth sciences. The nuclear energy industry has shown interest in this subject because fractures could eventually enable the release of radionuclides in the ground waters of deep geological waste repositories. In fact, multiple studies have modeled the migration of radionuclides in fractures of host rocks surrounding deep geological repositories for highlevel waste [START_REF] Andra | Dossier 2005 Argile -Tome Évolution Phénoménologique du Stockage Géologique[END_REF][START_REF] Zhang | Thermo-hydro-mechanical behavior of clay rock for deep geological disposal of high-level radioactive waste[END_REF][START_REF] Nykyri | Radionuclide release and transport RNT-2008[END_REF][START_REF] Cvetkovic | Stochastic simulation of radionuclide migration in discretely fractured rock near the Aspo Hard Rock Laboratory[END_REF].

By merging the lessons learned from both fields, it is possible to imagine quantitative simulations of the coupled chemical evolution of the glass and the migration of fluids and solutes within its fracture network at the block scale. To that end, good characterization of the fracture network is a key step before quantifying the long-term immobilization properties of the matrix. However, a discrete fracture model is considered to be too computationally demanding for conducting the reactive transport modeling (RTM) applied to a block of nuclear glass; instead, construction of an equivalent fracture network model and determination of its equivalent hydraulic and diffusive parameters are necessary.

In this article, a new methodology is introduced to characterize a fracture network of a glass package at scale 1 and to simulate multiple realizations of its equivalent representation to enable estimation of the uncertainties of the quantity of altered glass (QAG) resulting from the uncertainties of the fracture network-related parameters such as diffusion, permeability, and porosity. Beginning with highresolution photography of the cross-section of an inactive glass container, the methodology used here relies on the combination of mathematical morphology techniques and geostatistical analysis to build a model of the fracture network. The results from thermo-mechanical simulations of the cooling phase of the glass [START_REF] Barth | Sur la modélisation et la simulation du comportement mécanique endommageable de verres borosilicatés sous sollicitation thermique[END_REF] are also used, particularly the calculated arrival times of the solidification front.

The paper is organized in the following manner. Section 1, Material and methods, presents an overview of the steps involved in the characterization methodology of the image processing-based fracture network and includes (i) characterization of a glass block fracture network by using mathematical morphology approaches, (ii) numerical analysis of a physical parameter proved to be indicative of a glass block fracturing pattern, (iii) geostatistical exploratory analysis and modeling of fracture geometry-related parameters, and (iv) creation of multiple realizations of a glass block equivalent fracture network. Section 2, Results, details the implementation of the methodology on the high-resolution photography of a section of an inactive glass container. The fracture networks of the training block images are compared with the realizations of equivalent fracture networks. Section 3,

Material and methods

The complete workflow consisted of several steps, as shown in Figure 1: 1) fracture network image analysis, 2) creation of a proxy model via analysis of the map of arrival times of the solidification front, 3) geostatistical modeling, and 4) construction of multiple realizations of fracture network equivalent tessellations. Application of the model, including upscaling of diffusion, permeability, and parameters governing the kinetics of glass dissolution, and subsequent RTM are beyond the scope of this paper but will be covered in a forthcoming article. Although major verification of the feasibility of the model is addressed in the second part of the research, where the results of the RTM are compared with the experimental data, in the Results section it is shown that the distribution and the anisotropy of fracture network static parameters are honored and reproduced. The variation of the fracture density is respected and the distribution of the fracture apertures is preserved.

The workflow asks for input of the pre-prepared grayscale image and progresses through seven stages: (i) fracture network extraction, (ii) measurement of fracture apertures, (iii) image segmentation, (iv) anisotropy characterization of the segmented image, (v) calculation of average aperture for each zone (segment) of the segmented image, (vi) construction of a proxy model based on the analysis of the map of the arrival times of the solidification front [START_REF] Barth | Sur la modélisation et la simulation du comportement mécanique endommageable de verres borosilicatés sous sollicitation thermique[END_REF], (vii) geostatistical modeling of the position of mesh seeds and averaged aperture values, and (viii) simulation of equiprobable realizations of fracture network equivalent tessellations. Free libraries MAMBA 2 and Smil3 written in C and Python were adopted to conduct morphological image analysis. Analysis of the map of the arrival times of the solidification front was conducted in Python, R, and Paraview. For geostatistical exploratory analyses, simulations of equivalent tessellations were performed in R, specifically the RGeostats package 4 .

Glass block section

In our research, morphological mathematics methods were engaged to characterize the morphology of the fracture networks via extensive analysis of high-resolution photography. The images of the longitudinal cross-sections of specimens of non-radioactive packages of nuclear glass were used as a testbed for this study. Two packages were produced at the Vitrification Production Facilities in Marcoule in a cold crucible induction furnace: uranium-molybdenum 5 (UMo) and uranium oxides6 (UOx) [START_REF] Vernaz | Nuclear Waste Conditioning[END_REF]. Molten glass was poured in a steel canister 1.1 m high with an internal diameter of 0.42 m. The injection from the crucible occurred in two phases, with each casting being half of the nominal volume of 400 kg in total. The duration between the two castings was close to 6 h. After a waiting time of about 2 h below the furnace, the package was removed and was placed into a room in which the cooling progressed further with no implication of forced airflow. Soon after their fabrication, the canisters were cut in two parts along their long axis. Some precautions were taken to preserve the integrity of the intrinsic fracture network. Hence, viscous resin was poured on top of the UMo glass prior to the sawing operations. However, resin could not fill the entire network. In a similar manner, the integrity of block UOx was not fully assured during the cutting operations; the batch of grout poured on the top of the glass package could not penetrate the network.

Nevertheless, we assumed in this work that the resultant fracture networks separated from the pullout zones were representative and fully appropriate to serve as objects for further modeling.

A series of photographs of the two block halves, at 200 kg with dimensions of 430 mm × 1335 mm, were shot at different resolutions. The methodology developed in this study is presented by using the UMo image (Figure 2) and includes a complete process from extraction of the fracture network to simulation of its multiple equivalent realizations. The UOx image is presented subsequently when comparing the results of the equivalent network construction to the real fracture networks. The image parameters for both blocks are indicated in Table 1.

Thermo-mechanical simulation

The thermo-mechanical model used to complement the photography information is based on the results of the thermo-mechanical fracturing program [START_REF] Barth | Simulation of cooling and solidification of three-dimensional bulk borosilicate glass: effect of structural relaxations[END_REF][START_REF] Barth | Modeling and simulation of the cooling process of borosilicate glass[END_REF] of the VESTALE project initiated and financed by the French Atomic Energy Commission (CEA), National Agency for Radioactive Waste Management (ANDRA), and Orano (previously Areva). The objective of the program was to develop a tool for simulating the fracturing of nuclear borosilicate glass used for radionuclides confinement. For this purpose, the finite element method was used to model the viscoelastic behavior and the damage behavior, evolving in response to the structural relaxation of the glass. This program was established to build models from the thermo-mechanical parameters determined experimentally on SON68 glass, as given by CEA, and to provide simulations on the scale 1 package with the most representative scenarios of the industrial process. The 3D representation of the thermo-mechanical model is shown in Figure 3a.

Although the simulation results showed the evolution of the glass fracture network surface due to radiogenic heating, it is not expected to be significantly different from the fracture networks analyzed in this research [START_REF] Barth | Sur la modélisation et la simulation du comportement mécanique endommageable de verres borosilicatés sous sollicitation thermique[END_REF] The radioactivity has the effect of slightly reducing the fracturing of the blocks. The inactive glass used in this study can therefore be considered conservative in the sense that radioactive glass would develop a similar or reduced fracture network.

The position of the solidification front was obtained by analyzing the delay between the fictive temperature (T f ) and the glass temperature at each time step T(t). In fact, it is very difficult to determine the glass transition temperature (T g ) directly because its value depends on two fundamental phenomena that cannot be summed up in a single variable expressed by the temperature. On the one hand, T g is influenced by the thermal history T(t) of the glass network and on the other hand, the vitreous transition depends on the property of the glass under consideration (e.g. specific volume, thermal conductivity, specific heat capacity, viscoelasticity, refractive index, etc.). For this reason, Barth used T f generalized by Narayanaswamy and Moynihan in the model "TNM" (Tool-Narayanaswamy-Moynihan) [START_REF] Moynihan | Dependence of the Fictive Temperature of Glass on Cooling Rate[END_REF][START_REF] Narayanaswamy | A Model of Structural Relaxation in Glass[END_REF]. The state of structural relaxation, obtained through T f , quantifies the progress of structural relaxation over the thermal history. T f is calculated in every point of the volume of glass. The internal glass properties change in a different way according to the state of the glass, which can be liquid, solid or semiliquid / semisolid. In consideration of the T f and T evolution, it was possible to determine at each time step of the simulation the position of the solidification front (Figure 3b). The assembled map of the arrival times of the solidification front (Section 1.4) will be used to supplement pointwise data from the image analysis with global spatial distribution of the key parameters.

Figure 3 Results of thermo-mechanical study: a) three-dimensional (3D) mesh of the package for thermomechanical simulations showing longitudinal and transverse sections, obtained from [START_REF] Barth | Sur la modélisation et la simulation du comportement mécanique endommageable de verres borosilicatés sous sollicitation thermique[END_REF]; b) Variable expressing the delay between the fictive temperature and the temperature inside the glass (ΔT= T f -T), indicating the glass internal state (blue: liquid; red: solid). Evolution of the solidification front approximated by the isoline of the delay ΔT= 10.

Fracture network characterization by mathematical morphology

Mathematical morphology methods were used to characterize the fracture network of a block of nuclear glass. First developed by Matheron [START_REF] Matheron | Eléments pour une théorie des milieux poreux[END_REF] and Serra [START_REF] Serra | Introduction à la Morphologie Mathématique[END_REF] in 1960s, it is currently applied to many domains in which image processing in both two and three dimensions is requested. Its high demand in the geosciences is easy to explain because it offers a nondestructive approach of heterogeneous medium characterization by combining non-intrusive analysis such as micro-tomography [START_REF] Tokan-Lawal | Investigating flow properties of partially cemented fractures in Travis Peak Formation using image-based pore-scale modeling[END_REF], scanning electron microscopy [START_REF] Blunt | Multiphase flow predictions from carbonate pore space images using extracted network models[END_REF][START_REF] Jouini | Multiscale characterization of pore spaces using multifractals analysis of scanning electronic microscopy images of carbonates[END_REF], and ordinary digital imaging [START_REF] Chatterjee | Image-based quality monitoring system of limestone ore grades[END_REF]. Because mathematical morphology is based on the assumption that any image consists of structures that can be handled by the set theory; these sets are studied to determine important network attributes such as the porosity, connectivity, and size of objects that comprise the fracture network.

Fracture network extraction

Image processing includes several steps and aims to extract the fracture network to be studied afterward. In the first step, the chosen images are cropped so that the stainless steel container edges are removed and only the glass matrix remains. The images are then transformed into grayscale images of 8 bit/pixel. The gray images are filtered via opening by reconstruction, defined as that by erosion of the image using a hexagon as structuring element followed by a morphological reconstruction [START_REF] Vincent | Morphological Grayscale Reconstruction in Image Analysis[END_REF]. This operator enables elimination of light reflections that cannot contain the structuring element, leaving other image features unaltered.

The fracture network is then extracted by the following operations:

1) Black top hat transformation, defined as the difference between the closing of the image and the image itself, is performed by a hexagon structuring element of size 20 to remove slow trends and to perform contrast enhancement. This operator brings out fractures whose thickness is half the chosen size; as a consequence, it facilitates the separation of the glass matrix and the fracture network. 2) Thresholding based on the histogram of the image is conducted to separate the clear gray background of a glass matrix from the foreground dark gray fracture network together with the pullout zones. 3) Supremum 7 of linear openings by reconstruction is performed to preserve objects (fractures) greater in size than the chosen opening size in at least one direction. This operator ensures conservation of the elongated fractures, with the noise coming from the thresholding removed. All retained fractures are reconstructed to their real forms by geodesic reconstruction.

As mentioned previously, both fractures and pullout zones created during the sawing operation 8 are present in the foreground. The next step is devoted to their separation from the fracture network so that the artificially created pullout zones are removed. For this purpose, the operators used are closing, which closes fractures that are thinner than the pullout zones, and hole filling, which closes 7 In this article, "supremum" and "union" as well as "infimum" and "intersection" are used as interchangeable synonyms; supremum (or point-wise maximum) and infimum (or point-wise minimum) replace the union and the intersection set operators for grayscale images. 8 The bloc has been split in two longitudinal cross-sections to expose the existing fracture network.

the pullout zones. Because certain fractures cross the extracted pullout zones, their prolongation is assured by tracing the skeleton of influence zones. This operation does not add redundant connections but enables preservation of the connectivity of some fine fractures in the artificially damaged zones. Final slight filtering is performed by extraction of the fracture skeleton, followed by removal of the skeleton ramifications and the geodesic reconstruction of the resulting skeleton.

Measurement of fracture apertures

Fracture aperture is a principal network characteristic governing the intensity of the glass aqueous alteration. Three methods based on morphological mathematics are proposed for measuring the fracture apertures.

Discrete way of aperture measurement via bounding boxes

The first method is based on measuring the Feret diameters of each rotated bounding box containing a single, non-connected labeled component, or each simple branch, of the fracture network. The separation of the network into non-connected branches is performed by subtraction of the dilated multiple points9 from the extracted fracture network (Figure 4, Figure 5a). It should be noted that each bounding box, i.e. the rectangle that circumscribes its corresponding fracture, is aligned with the inertia tensor computed from the second-order moments of this fracture (Equation 4). This makes the values of minimum Feret diameter indicative of the average fracture aperture.

Continuous way of measurement via quench and distance functions

The second and third methods rely on the techniques of homotopic thinning. In the framework of these methods, the aforementioned problem of the fracture network separation is resolved by implementing several operations:

1) homotopic thinning with the L structuring element [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF];

2) removal of the skeleton's endpoints (pruning); 3) extraction of multiple points at the intersections of fracture branches; and 4) subtraction of multiple points from the resulting L-skeleton.

After implementing this sequence of operations, we obtained simple arcs that represent disconnected ramifications of the fracture network. For each ramification, its barycenter and aperture average will be found subsequently.

We used two techniques of aperture measurement: the quench function10 (Figure 6a), associated with the skeleton of maximal balls, and the distance function (Figure 7), recorded on the traces of simple arcs 11 . Both of these functions give the (n + 1) value, where

• for the first technique, n is the size of the maximal disc put at point i on the skeleton by maximal balls; and • for the second technique, n is equivalent to the level of the distance function.

These techniques of measurement thus provide information on the thickness of the fracture branches in a continuous manner. However, the averaged apertures of the fractures are not expected to be identical for several reasons. First, on the periphery of fractures, the traces of skeleton branches are expected in which the quench function values are equal to 1. This occurs because fracture walls are not ideal flat lines, and the skeletonization by maximal balls is highly sensitive to even small variations in its boundaries. These irrelevant branches of the skeleton (Figure 6b) should not be considered in the course of aperture average computation. This step can be performed by finding the intersection of the traces of simple arcs and the skeleton of maximal discs 12 . It should be noted that the aperture average value measured by this technique is not always available for all simple arcs. For some simple arcs, intersection of the trace of skeleton by simple arcs and the trace of the skeleton by maximal balls does not exist. This issue can be resolved by reading the values of distance function directly on the traces of the simple arcs without calling upon their intersections with the skeleton by maximal balls. The distance function of the fracture network image is calculated by using the MAMBA image library. The values of average aperture of each fracture network ramification were calculated by formula (1):

∑(2(F-1)+1)/nb = ∑(2F-1)/nb, ( 1
)
where F is the value of the distance function or quench function read on the trace of the corresponding simple arc, and nb is the number of points making up the trace. In formula (1) on the left-hand side, one comes from the fact that there is an influence zone (Figure 7b) that is considered when measuring fracture width.

Image segmentation

In mathematical terms, segmentation is the partition into disjoint nonempty sets known as segments [START_REF] Soille | Morphological Image Analysis: Principles and Applications[END_REF]. Usually, an image that has been segmented is represented as a label image that enables separate processing of all compound elements including disjoint sets named segments, as discussed below. The main objective of the segmentation is to define segments that represent parts of the glass matrix separated by the extracted early fracture network. This is realized by a sequence of operations based on watershed transformation [START_REF] Beucher | Watersheds of functions and picture segmentation[END_REF]. To avoid oversegmentation, we imposed markers represented by the labeled maxima of the reconstructed half distance function. The development of individual tiny segments positioned inside large segments is prevented by applying the supremum of the small size opening at the marker construction stage. The edges of the catchment basins delineate the glass segments (Figure 8).

Characterization of the segmented image

This section details the use of the segments obtained as a result of image segmentation. Introduction of the representation by segments facilitates geostatistical simulation of the block fracture network, which makes uncertainty estimation possible.

The concept behind the adoption of the representation by segments is that it enables replication of the block facture network by reproducing its segmented image. That is, our objective is to portray a previously extracted fracture network by reproducing parts of the matrix separated by this network via any type of tessellation that could reconstruct, or model, (i) the internal texture of the segmented image to identify the relationships between segments and (ii) the spatial structure and distribution of average fracture aperture. The internal structure of the segmented image is determined by the position of segment centers, the rotation angles of the segments, and the segment shape factor.

Here, we focus on the characterization of the segmented image by first-and second-order moments. This approach is widely used in image analysis [START_REF] Flusser | Moments and Moment Invariants in Pattern Recognition[END_REF][START_REF] Flusser | 2D and 3D Image Analysis by Moments[END_REF] and was adopted from the mechanics of rigid and deformable bodies [START_REF] Lubarda | Areal moments of inertia revisited: on the distinction between the principal directions[END_REF].

The characterization includes computation of the following parameters: lengths of the main axes (A, B) of each segment, segment rotation angles (theta), and positions of segment centers. The rotation angle is defined as the angle associated with the largest axis of a segment and a horizontal plane and is positive when the major axis is turned in a counterclockwise direction. Image moments are defined as weighted averages of the pixels' intensities (i, j). In our case, we are concerned with the segmented binary image in which the (p,q) 13 moment M p,q gives equation ( 2):

, = ∑ , ∈ . (2)
The first-order normalized moments give the coordinates of the barycenter in the horizontal and vertical directions, and the second-order central moments determine the equivalent ellipse, furnishing the orientation and the dimensions of the major and minor axes.

The three central moments of second-order central moments ( , , , , , ) form the inertial tensor of the rotation of the object about its gravity center, expressed by equation ( 3):

( ) = ! ",# ! $,$ ! $,$ ! #,"
%. (3) The main inertial axes corresponding to the semi-major and semi-minor axes of the ellipse are derived from the eigenvalues of this tensor. The orientation of the object is defined as the tilt between the x-axes and the axis, around which the object is rotated with the minimal inertia and is calculated by expression (4):

& = #. ( )* )+ "! $,$ ! ",# ,! #," . (4)
For all segments, the position of segment barycenter (X i , Y i ), axis length (A i , B i ), and rotation angle (theta i ) are computed and further exploited at the stages of the geostatistical modeling (Section 1.5.2) and the creation of fracture network equivalent tessellation (Section 1.5.1) 14 . 13 p,q denote the order of a moment 14 Voronoï anisotropic tessellation built in RGeostats package, Mines ParisTech

Calculation of aperture average inside segments of the segmented image

Despite the high accuracy of the continuous approach of aperture measuring and its judicious results, processing of the fracture network via simple arc representation is problematic. In fact, the aperture values of each fracture ramification kept at their centers are not usable, at least directly, in the framework of the segment representation approach. In fact, at the stage of fracture network modeling by the technique of mesh growth, we expect to determine the positions of the mesh boundaries in which the thickness values need to be provided at the centers of the meshes.

The aperture average of the fractures associated with each segment should therefore be calculated. As discussed in Section 1.3.2, we first separated the fracture network into its simple constituents, and we then calculated the aperture average for all ramifications via one of the proposed approaches such as bounding boxes, quench function, or distance function. To calculate the aperture average inside the segments of the segmented image, we used the results of the aperture measurements via distance function because this continuous method of fracture aperture measurement is considered to be more accurate, and the distance function approach is believed to be less biased in comparison to the quench function approach. Afterward, the fractures were grouped according to their position in the segments.

All fractures with simple arc centers positioned inside a segment are used to calculate the aperture average assigned to that segment barycenter position (Figure 9). The operation is repeated for all segments.

Analysis of the solidification front map

In geostatistical modeling, it is recommended to make use of all available soft data that can supplement the raw data, particularly if the raw data are scarce and are difficult to interpret. As previously mentioned, the raw data include the values of parameters obtained by image analysis of one block of nuclear glass such as the values of apertures and positions of segment centers. This database cannot be considered as exhaustive. Moreover, it should be noted that the analyzed block could be very likely damaged at the cutting stage. Thus, photography analysis is supplemented with the results of the thermo-mechanical simulations (Section 1.2) to provide a physically meaningful parameter that helps to constrain the spatial structure and distribution of the parameters of the fracture network.

For this purpose, an exhaustively known parameter has been identified whose transforms can supplement the raw data at the stage of the equivalent fracture network simulation. Of particular interest are the knowledge of the direction of the fractures and the aspect ratio of the cells. These data can be determined through analysis of the segmented image by first-and second-order moments (Section 1.3.4). However, these pointwise data known only at the segment centers are not sufficient for construction of an equivalent tessellation. Indeed, values of the rotation angles and the ratio of the scaling factors of the meshes should be provided in the form of a continuous map.

Numerous studies [START_REF] Barth | Modeling and simulation of the cooling process of borosilicate glass[END_REF][START_REF] Crevoisier | Semi-stochastic generator (FRAGMA) of 2D fractured media by mechanistic analogy. Application to reactive transport in a fractured package of vitrified nuclear waste[END_REF] have reported the fabrication process of a full-scale package specimen containing nuclear glass, which is a close analogue of the block analyzed in our study. According to the previous research, the map of the solidification front arrival times indicates that fissuring occurs first along the perimeter of the package and then propagates radially to the block center [START_REF] Crevoisier | Semi-stochastic generator (FRAGMA) of 2D fractured media by mechanistic analogy. Application to reactive transport in a fractured package of vitrified nuclear waste[END_REF]. The network density is highest on the periphery owing to the highest intensity of the temperature drop that initiates fracturing immediately after the bloc fabrication ends. In contrast, the fracture density is significantly lower close to the package core, where the temperature remains homogeneous. In addition, the upper part of the package is less fractured because it represents a stress-free surface, where dissipation of the external energy is not limited.

After analyzing several variables that describe the thermo-mechanical evolution of the internal state of glass after being poured in a rigid canister, we decide to select a variable referred to throughout this paper as the solidification front. This variable indicates the temperature delay between the fictive temperature 15 and the temperature inside the glass and serves to indicate the transition from liquid to solid state inside the glass.

First, we assembled a map of the arrival times of the solidification front by selecting front positions for all ranges of time values existing in the thermo-mechanical model. Later, we interpolated the arrival time values by applying kriging estimation with a linear model. The result is shown in Figure 10a. Second, we calculated the angles between the axis parallel to the block axis and all pairs of neighboring points of the solidification front sequentially for each arrival time. We used the same angle convention as mentioned earlier. Then, we interpolated the values of the rotation angle by applying kriging estimation with a linear model. The map of the fronts line rotation angles is presented in Figure 10b. Third, we proceeded with analysis of the lines of the glass solidification fronts by calculating their curvature using the formula of the Menger curvature (5) [START_REF] Pajot | Analytic Capacity, Rectifiability, Menger Curvature and the Cauchy Integral[END_REF] for each set of three points of the solidification front:

(-, ., /) = $ 0 = 1 * 3 |-,.||.,/||/,-| , ( 5
)
where S is the surface of the triangle spanned by the x, y, and z points At this stage, we followed the concept such that in the first approximation, the shapes of the solidification fronts determine the shapes of the glass matrix segments. That is, fractures divide the glass into segments such that along the long sides of the block, it is probable that elongated segments meet with A/B >> 1. However, in the zone of the re-liquefaction at the boundary between the first and the second castings and in the zone of the first casting in the middle of the block, the fractures have identical lengths in two dimensions (A/B ≈ 1). In fact, when estimating the A/B ratio of scaling factors, it is useful to employ observations that can be explained by the physics of glass structural relaxation. To be precise, we used the map of front line curvature, presented in Figure 10c, as an external drift in order to generate the map of the A/B16 ratio of the mesh scaling factors. The external drift kriging technique [START_REF] Hudson | Kriging Temperature in Scotland using the External Drift Method[END_REF] implies that the overall shape of the target variable is given by an external drift map, whereas residuals are responsible for the fluctuations around this drift.

Geometric and geostatistical modeling

As mentioned in Section 1.3.3, the objective of the segment representation is to enable replication of the block facture network by reproducing its segmented image, with the seeds of the cells and the thickness of borders provided by the geostatistical modeling. In the following section, geostatistical modeling of seed positions and average thickness of the cell borders will be discussed.

To better understand the meaning of these parameters, which will be discussed in Section 1.5.2, we here explain the technique intended to be used for constructing multiple realizations of fracture network equivalent tessellation.

The generation of multiple equivalent realizations of the block fracture network is based upon the geostatistical spatial analysis of the data coming from the image analysis and the information collected from the map of front of solidification.

Creation of fracture network equivalent tessellation

To model an equivalent block fracture network, we used the anisotropic version of Voronoï algorithm [START_REF] Anton | The Voronoi Diagram of Circles and Its Application to the Visualization of the Growth of Particles[END_REF], which represents a version of mesh growth models and delineates mesh boundaries considering the cell scaling factors and the cell rotation matrix.

It should be noted that in the remainder of this paper, the term "cell" is applied as an analogue of the term "segment" by referring to a separate unit of an equivalent representation of the segmented image. Moreover, the term "seed" or "centroid" will be used to replace the term "center" in reference to the origin of the corresponding cell growth.

According to the Voronoï algorithm, a cell C (P i ) generated by the growth of P i centroid (seed) is Mathematically, the cell is defined by formula (5): where υ is the constant radial rate of growth equal to 1.

First, the distance d (8, 7 5 ) between the datum and the target site is computed by considering the rotation matrix A BCDE * DFGE * -DFGE * BCDE * I, where & * is the rotation angle specified in all points M belonging to the cell.

Next, the scaling (defined by a M i , b M i ) is operated in the rotated system. The component is unchanged if the scaling factor is equal to 1, reduced if it is higher than 1, and increased if it is lower than 1.

Finally, the boundaries of a newly created equivalent tessellation are thickened by considering the results of the geostatistical modeling of the average aperture calculated inside each segment (Section 1.5.2).

The average values of the apertures are attached to the cell seed positions, whereas the rotation angle and scaling factors in all directions 17 are defined on a fine grid.

As presented in Sections 1.3.4 and 1.4, the maps of the rotation angle and the ratio of scaling factor had been obtained by analyzing the map of the arrival times of front of glass solidification, while modeling of seed positions and average aperture of cell boundaries will be covered in the next section.

Geostatistical analysis

It is apparent that the fracture network morphologies and inherent characteristics vary among blocks. The reasons for these variations are numerous and are related mainly to the procedure of block fabrication and some random fluctuations. It is evident that when dealing with only one specimen, it is impossible to evaluate the uncertainty of the altered glass quantity resulting from the variation in parameters related to the fracture network. Therefore, in the next step, we conduct exploratory spatial data analysis followed by geostatistical simulations of two variables. The first is the position of cell seeds. The second is the average aperture assigned to cell seeds that will determine the thickness of the boundaries between the adjacent cells, which in the context of our research represents the fracture average aperture. The results of the simulations will be used to construct multiple realizations of an equivalent fracture network based on the Voronoï algorithm explained in the preceding section.

Position of cell seeds

Calculation of the point density of the training image is based on the results of the image analysis where the segments centers have been located within the upper part of UMo segmented image (Figure 17a, Y > 1215); the lower part is obtained by a mirror symmetry (the reason why the lower part (Figure 17a, Y < 1215) is discarded is explained in Section 2). Calculation of the point density of the training image is based on the results of the image analysis, specifically the positions of segment centers of the upper part of the UMo segmented image reflected by mirror symmetry to obtain the lower part. This step includes (i) calculation of the number of centers divided by the surface 18 of counting window with windows of varying dimensions and (ii) regression of these density values against the surface values. This regression analysis is used to evaluate the impact of the moving window size on the density regionalization. Later, the values of point density, determined on the grid support, 19 are migrated to the position of the segment centers. Further, the geostatistical analysis and modeling is based on the values of density collected at 352 segment centers. In the same manner, the values of the solidification front arrival time are migrated to the positions of the segment centers. Next, the values of density are regressed against the values of the solidification front arrival time, and the procedure is iterated for all window sizes. The spatial variability of densities and their residuals are studied by using the experimental variograms and their best-fit theoretical models. Then, the parameters derived from the fitted models for density and residuals, and the values of correlation coefficients are examined to determine the impact of the moving window dimensions on the density values and to eventually determine the density values to be used as the input data for further modeling. After the residuals are normal-score transformed to correct for the slight skewness of the histogram, the non-conditional simulations (using the turning band method [START_REF] Matheron | The intrinsic random functions and their applications[END_REF]) are conducted to produce N realizations of the spatial distribution of residuals. After the back-transformation to the raw scale and the use of the coefficients of the linear correlation, N maps of the varying Poisson intensity are obtained. They are subsequently used to generate the spatial Poisson point process with varying intensity [START_REF] Lantuejoul | Geostatistical Simulation: Models and Algorithms[END_REF]. Moreover, a rejection step is used to ensure that the simulated seeds are separated by a minimum distance that is equivalent to the minimum distance between the centers of the training image. It should be mentioned that the N realizations of the intensity maps constructed in the previous step are calibrated to the number of the centers in the training image.

Averaged apertures assigned to cell seeds

After determining the position for cell seeds, the averaged apertures assigned to cell seeds were modeled. As input data, measurements of the averaged apertures inside segments of the segmented UMo image (Section 1.3.5) were used 20 .

Next, we verified the correlation between 1) the values of average apertures and the arrival times of the solidification front and 2) the values of average apertures and the values of Poisson density calculated in Section 1.5.2.1.

Because the correlation for both pairs of variables was inadequate, we defined a geostatistical model of averaged which does not involve any explanatory function.

The geostatistical modeling of average aperture on the cell seeds was performed in four steps: 1) normal score transform of the data, 2) calculation of the experimental variogram and determination of its theoretical best-fit model, 3) non-conditional simulation of the aperture average values in the normal space by using the turning band method, and 4) back-transform of simulated values in the original scale.

Results

In this section, we first show the results of the fracture network characterization of the UMo training image, and we present the maps obtained by the analysis of the map of the solidification front arrival times. Next, we demonstrate some realizations of the geostatistical modeling of seed positions and average aperture, and we shows the resultant realizations of the equivalent fracture network constructed by considering all previously demonstrated results.

Figure 11 illustrates the fracture network extracted separately from the pullout zones. To select the size of the structuring element necessary to perform the morphological operations, we calculated and analyzed the distribution curves with several types of opening, including a hexagonal opening and that by erosion-reconstruction (Figure 12). When selecting the size of the structuring element, we adhered to the following rules:

1) The chosen size of the opening by erosion-reconstruction is higher than the maximal size of the reflections. According to the distribution curves computed with morphological erosion and with the opening by erosion-reconstruction, light reflections are expected up to 5 pixels 2) The size of the black top hat transformation is chosen with respect to the range of the fracture thicknesses. According to the distribution curves, fractures are expected to be up to 20 pixels in width.

3) The threshold operator is applied to separate the matrix background and the identified fractures. Because the black top hat transformation puts the matrix values close to zero, we fix the threshold to 30, which is slightly higher than the maximal expected fracture width. By doing so, we assumed that the crack depths are greater that their widths. 4) The size of the morphological closing implemented to set apart fractures and pullout zones is equal to the size of the average fracture width, at 10 pixels. The histograms of apertures measured by applying the three proposed approaches are shown in Figure 13 and in Table 2. All approaches gave similar values of fracture aperture; however, the values obtained by the bounding boxes method were discrete more often. With respect to the results, there is no evidence that the bounding boxes approach overestimates the values of fracture thickness, which means that correction for shape curvature made at the stage of results post-processing is justified.

Regarding the continuous way of aperture measurement, the use of the distance function read along the traces of simple arcs is considered to be more judicious and gives more accurate results because it does not imply the use of the skeleton by maximal balls, which has technical limitations. The segmentation results are illustrated in Figure 14. The segmentation procedure is realized by marker controlled watershed (see Section 1.3.3) in order to avoid oversegmentation that could have occurred owing to the presence of the spurious minima. Undoubtedly, the choice of the marker image is important because it determines the output. Because our objective was to replace fracture representation with segment representation, it was critical that we obtain a segmented image that has no excessive segments inside the segments to ensure that the segments were separated exclusively by fractures or their prolongations in the absence of fractures. At this stage, we considered that a satisfactory result was achieved. Therefore, characterization of the segmented image by first-and second-order moment was next conducted.

As was stated in Section 1.3.4, the objective of the segmented image characterization is to determine attributes of the image that must be reproduced in order to mimic the anisotropic nature of fracture network of the glass blocks. The result of the calculation of segment rotation angles and the ratio of their axis lengths is demonstrated in Figure 15 in the form of the maps, where, the calculated parameters are used to fill the segmented image for enhanced visualization. Unfortunately, these data were not sufficient for use as input maps needed to effectuate coordinate transformation, including scaling and rotation of the fracture network equivalent tessellations in order to consider fracture network anisotropy. Therefore, we used the soft data presented by the continuous map of solidification front arrival times obtained from the thermo-mechanical simulation.

Figure 10 provides the results of the analysis of the thermo-mechanical parameter representing the transition from liquid to solid state inside the glass. As shown in Figure 10a, the transition from one arrival time to the following arrival time is rather smooth; however, zones were noted in which the transition was sharp such as the abscissa axis, from 3800 to 4000 and the ordinate axis, from 0 to 500 and from 2000 to 2430. These zones are associated with the complex phenomena of recasting at the boundary between the first and the second glass casting. The zones are even more noticeable in the map of the normalized curvature (Figure 10c).

Following calculation of the curvature, estimation of the ratio of scaling factors A/B was performed using the external drift kriging technique. It should be noted that the estimation was conducted inside a specified polygon for one half of the block. Indeed, because the external drift map is symmetric, as the whole model of the thermo-mechanical behavior of a glass block, and some doubt existed on the quality of the fracture network of the lower part of the UMo block 21 , we conducted the estimation based only on the A/B measurements for the upper half of the block (Y > 1215, Figure 17a). The lower part (Y < 1215) of the map was obtained by applying the mirror symmetry. Moreover, it should be noted that the polygon did not cover the entire half of the block. The analysis of the segmented image in which the segments are filled with the values of A/B (Figure 17a) and some observations of transversal and axial cross-sections of glass blocks [START_REF] Crevoisier | Semi-stochastic generator (FRAGMA) of 2D fractured media by mechanistic analogy. Application to reactive transport in a fractured package of vitrified nuclear waste[END_REF] led us to presume that in the center of the block and close to the upper part of the block, large rather homogeneous segments likely form in which the A/B values could be assigned by default to 1.

Figure 16 represents the scatter plot showing the regression between the values of the external drift (curvature of the solidification front lines) and the available data of the searched variable A/B obtained through image analysis. The scatter plot above shows a weak linear correlation (ρ = -0.15) between the values of the curvature and the values of segment length ratio of the block UMo. In fact, along the long sides of the container in the zones in which the elongated segments are located, some isotropic segments are also present.

Despite the weak correlation, we used the map of normalized curvature as an explanatory deterministic map for construction of the map of the A/B ratio of the scaling factors. It fact, we did not look for an accurate assessment of this ratio because this map served only as a background map while constructing fracture network equivalent tessellations. Nevertheless, for future research, we strongly recommend analysis of the regression of the scaling factor ratio against the values of the normalized curvature for a new block (having ensured integrity).

Figure 17b illustrates the result of the A/B estimation. This map 22 together with the map of the rotation angle, were used to consider the anisotropy. Although it is possible to precisely measure the apertures of fracture ramifications, these data are not generally used at the stage of modeling of multiple realizations of equivalent fracture network. Thus, in an attempt to preserve the accuracy of the aperture measurements and to diminish the amount of data to be processed in the course of the geostatistical analysis, fracture aperture averages were computed inside the segments of the upper half of the segmented image and were assigned to their centroids. The same values of the average facture apertures were assigned to the centers of the lower half. The results of the averaged aperture measurements are shown in Figure 18. The results of the geostatistical modeling of seed positions and aperture averages that provide information on the thickness of boundaries of equivalent tessellations are given subsequently.

As explained in Section 1.5.2.1, Poisson density modeling is based on the calculation of the segment center density with different moving windows. The results are summarized in Figure 19. From these data, it is clear that by increasing the size of the moving window, the regularization of the density is increased, reflected by a decrease in the sills of the variograms. In consideration of the pros and cons of regularization, specifically variogram sill decrease and better correlation with the explanatory variables versus loss of the variability, we selected the density values calculated with a moving window size of 13 × 13 blocks, given that the total size of the grid is 123 × 49 blocks.

A theoretical nested model with nugget, exponential and spherical structures fitted to the experimental variogram of residuals, shown by an asterisk in Figure 19c, was used to perform nonconditional simulations of residuals. The parameters of the regression between the density and the explanatory variable, shown by an asterisk in Figure 19b, were involved when computing maps of the Poisson intensity. For each realization of the Poisson intensity map, the seeds were generated so that the minimum distance between seeds is larger than 25 pixels and their total number is close to 350. Both parameters were fixed with reference to the results of the analysis of the training image.

Afterward, for each realization of seeds, aperture average modeling was conducted. The nonconditional simulations were based on the theoretical nested-type model with nugget, exponential and spherical structures fitted on the experimental variogram (Figure 18b) and the model of the normal-score transform (Figure 18c). Three realizations of seed positions generated according to a spatial Poisson point process with varying intensity and the histograms of average apertures simulated at the seed locations are given in Figure 20.

Fracture network extraction

We believe that by having implemented operations for contrast enhancement (black top hat) and noise removal (opening by erosion-reconstruction) with a carefully chosen size of the structuring element, we obtained satisfactory results. Nevertheless, we admit that the resulting fracture network is not a unique representation that could be obtained from the raw image. In fact, changing the size of the structuring element of any operation involved in fracture extraction would certainly give a different result. However, the proposed procedure of fracture network extraction is believed to be useful for future analysis of cross-sections of glass blocks, especially because it offers the possibility of generating multiple realizations of the fracture network model.

ii. Measurement of fracture apertures

Although the results of aperture measurement provided by the three approaches (Figure 13) are compatible, each has weak and strong points that are worth consideration:

Advantages of the discrete method of aperture measurement via the bounding boxes approach:

1. Implementation ease 2. Rapid calculation Disadvantages of the discrete method of aperture measurement via the bounding boxes approach:

1. The facture network must be disintegrated into non-connected branches at the first step of the measurements. 2. The aperture values are derived from the measurements of the shortest axes of the bounding boxes surrounding the fractures and represent the maximal apertures of these fractures. By applying this method, it is impossible to consider the variability of the fracture aperture. The calculated aperture values are integers. For fractures in which the shapes are significantly curved, the length of the shortest bounding box axis does not reflect the real aperture magnitude. Therefore, to indicate the average aperture, these values of Feret diameter need to be corrected by using the ratio of the area of the object over the length of the longest bounding box axis.

General advantages of the continuous method of aperture measurement:

1. The problem of fracture separation24 can be bypassed.

2. This method is appropriate for characterizing fractures of non-constant thickness.

Disadvantages of the continuous method of aperture measurement via quench function:

1. The apertures are average values of the quench function found at the intersection of the skeleton of maximal balls and the traces of simple arcs. This could jeopardize the accuracy of the calculations because not all relevant points of the skeleton of maximal balls lie on the same trace as the simple arcs (Figure 6b). 2. The quench function is not secured from parity bias; in certain cases, it is prone to underestimate the real fracture aperture. Moreover, the maximal balls skeleton is not connected and can have thickness of more than 1 pixel, which might increase the error of average aperture estimation [START_REF] Beucher | Digital skeletons in Euclidean and geodesic spaces[END_REF].

Disadvantages of the continuous method of aperture measurement via distance function:

1. The distance function is not secured from the parity bias, which means that values of average fracture aperture can be also slightly underestimated.

Considering the above-enumerated advantages and disadvantages of the aperture measurement techniques, for further research, the continuous method of aperture measurement via the distance function is recommended to be retained because has been proved to be the most accurate and least biased.

iii.

Analysis of the map of front of solidification

A comparison of the values of the rotation angles measured for segments of the training UMo image (Figure 15a) and the map of the front line rotation angles (Figure 10b) revealed their good agreement.

However, the value of the linear correlation coefficient between the values of the curvature and those of the segment length ratio calculated for the training image was rather unsatisfactory. There are several possible explanations for this cohabitation of segments of different levels of lengthening generated in similar thermo-mechanical conditions, according to the map of internal structural relaxation. One is that the block chosen for the image analysis was not fully representative and could have carried imprints of the artificial damage, particularly in fragile peripheral zones, which would have masked the real fracture pattern.

iv. Geostatistical and geometrical modeling

Figure 20a proves that spatial distribution of seeds generated by the Poisson process and the centers of the segmented image are comparable; however, for each realization, the number of seeds generated in the central part of the blocks was higher than in the UMo training image. These small discrepancies are not believed to bias the result, particularly when considering that for the UOx block, more segments were observed in the central part (Figure 23b).

As shown in Figures 23, the resulting realizations of equivalent fracture network tessellation exhibited strong similarities with the fracture networks of nuclear glass block. Moreover, they have the potential for modeling the fracture network at scale 1 and for conducting the reactive transport modeling. However, the demonstrated fracture network equivalent tessellations possess some features that make them different from the real fracture networks. First, on the 2D equivalent network representations, all cells of the tessellations were connected, whereas the analysis of the 2D crosssection of the UMo fracture network showed some dying-out fractures (Figure 11). Second, although the total number of seeds of equivalent was close to the total number of centers of the training image, the simulated segments were smaller in the central part of the equivalent tessellations than in the average in the training UMo image. On average, the porosity of the presented equivalent tessellations was higher than that of the UMo training image. However analysis of the fracture networks of both UMo and UOx blocks and comparison with the presented equivalent tessellations indicate that simulated tessellations represent effective reproduction of fracture networks at scale 1. Moreover, they enable consideration of network structure variability.

Conclusions and perspectives

A crucial issue in reducing the uncertainty associated with a fracture network model is to consider all data that can be of very distinct natures. In our case, a combination of hard and soft data was used. Hard data were presented by the results of the image analysis on one block of nuclear glass, whereas the soft data correspond to the map of the arrival times of the solidification front and some derivatives of this map.

The first part of the workflow, extraction and determination of the fracture characteristics, appeared to be reliable and offered important insight into the meticulous examination of fracture networks at scale of a glass block. However, it should be noted that the block used as a test bed of the workflow application is not fully indicative because its fracture network was damaged in the course of sawing. Thus, the resulting values of the apertures demonstrated here should not be considered as incontestable and characteristic of a real intrinsic fracture network of a block of nuclear glass. In pursuit of more relevant values of fracture apertures, further analysis using the same methodology on carefully prepared sections is advised.

The second part of the work dealing with the geostatistical modeling and construction of fracture network equivalent tessellations is believed to be a practical and rapid method of fracture network simulation allowing the variability of fracture distribution and aperture to be considered.

The presented research is methodological and strives to elaborate the workflow needed to conduct the reactive transport modeling at the scale of the fractured glass package. Nevertheless, in order to demonstrate the feasibility of the proposed procedure, we presented some results of the training image characterization and used them for further modeling. The findings, while preliminary, suggest that at relatively low expense, it is realistic to simulate a fracture network of blocks of nuclear glass by using the proposed anisotropic Voronoï tessellation. Equivalent tessellation is a good approximation and represents object modeling that considers both hard and soft input data. However, the final test of the quality of the model will come from applications: a comparison should be made of the results of reactive-transport simulation on the direct support (training images) and on the realizations of the equivalent tessellation.

Although the characterization and modeling efforts are important, it is critical to improve the quality of glass block preparation by ensuring the integrity of the fracture network or by using noninvasive techniques. At the same time, it could be useful to deepen the understanding of the relationships between thermal damage parameters and fracture distribution and geometries. In fact, the choice of the actual proxy map was based on the fact that the speed of the solidification front makes it possible to dissipate internal stress more or less easily, thus controlling the morphology of the fracture network. Other types of maps might be useful if they are justified by mechanical considerations.

Although this research focuses on a 2D model, by unifying information obtained by analyzing images of longitudinal and transversal cross-sections of glass canisters, it is realistic to construct a realistic 3D model. This is especially appealing, since almost all tools necessary to accomplish this exercise are already prepared. First, the knowledge of the fracture pattern of a transverse plane of a fracture package can be obtained by applying the proposed workflow. Second, the data required for geostatistical modelling are identified. Moreover, the maps of the solidification front arrival times necessary to support the image analysis results can be procured easily and the way in which they should be analyzed is already known. Besides, the algorithm used to generate the anisotropic Voronoї tessellation (realized in RGeostats package) should potentially be able to do it in 3D after some proper adaptations are made. These additional implementations concern the creation of the connected 3D fracture planes. However, certain caution must be applied because the fracture network is likely to vary along the z-axis owing to the presence of the different thermo-mechanical environments such as a zone of re-liquefaction or a stress-free surface. Nevertheless, it could be interesting to create realizations of 3D equivalent tessellation based on the data obtained by analyzing two transversal cross sections and one longitudinal cross section. The transversal cross sections should belong to two different castings, in order to capture the difference of the fracture network morphology related to the two-stage manufacturing procedure.

As a general conclusion, we emphasize that the workflow presented here, combining (i) image analysis, (ii) analysis of a physical parameter indicative of glass internal structural relaxation, (iii) geostatistical modeling, and (iv) reproduction of equivalent fracture network tessellations is a powerful tool that enables the combination of hard and soft data concerning nuclear glass fracturing at scale 1. As such, it provides the basis for future investigation of the impact of fracturing on the aqueous alteration of borosilicate nuclear glass.

The demonstrated workflow is part of a more complex workflow and data integration process aimed at characterization of fracture network on the scale of a nuclear glass block, upscaling of permeability, diffusion and kinetics governing parameters, and subsequent reactive transport modeling with consideration for variability in the fracture network. Future work includes the use of the model fracture network for reactive transport simulations at the block scale. These simulations will bring quantitative assessment of the degradation of glass and the release of radioanuclides under different scenarios of evolution, bridging the gap between laboratory scale (micrometers, years) and geological repository scale (meters for the block, tens of thousands of years).

Variogram Function describing the degree of spatial dependence of a spatial random field or stochastic process. It is calculated as a half the average squared difference between all points separated at distance h.

Figure 1

 1 Figure 1 Image processing-based fracture network characterization method. The parts contoured with dashed lines are not presented in this paper.

Figure 2

 2 Figure 2 Grayscale image of the UMo bloc of nuclear fractured glass. Dimensions: 1044 mm × 413 mm.

Figure 4

 4 Figure 4 Preparation of fracture network for aperture measurement via the bounding boxes approach a) Fragment of UMo image with dimensions of 318.9 mm × 122.2 mm; b) multiple points at the locations of fracture branching shown in yellow. For visibility in this image, the multiple points and the skeleton by homotopic thinning are dilated with a hexagon structuring element of size equals to three.

Figure 5 a

 5 Figure 5 a) Fragment of UMo image with dimensions of 68 mm × 19.9 mm. Patchwork palette was used to mark the fracture ramifications separated by dilated multiple points; b) Each ramification was circumscribed by a bounding box in which the minimum Feret diameter (min D F ) was used to measure the fracture aperture average. The shown ramification has min D F = 12 px (2.04 mm).

Figure 6

 6 Figure 6 Fragment UMo image with dimensions of 27 mm × 30.3 mm: a) sketch showing superposition of fractures and the skeleton by homotopic thinning with overlapped maximal balls in which the radii are recorded to evaluate the fracture aperture average; b) patchwork palette application to visualize the quench function that represents the locus of the skeletal points with associated values of the radii of maximal balls.

Figure 7 a

 7 Figure 7 a) Fragment of UMo image with dimensions of 127.2 mm × 43.5 mm. The visualized distance function used the values read on the traces of the simple arcs to calculate the fracture aperture average; b) sketch showing a hexagon structuring element and an influence zone taken into account when calculating quench and distance functions.

Figure 8

 8 Figure 8 Fragment of UMo image with dimensions of 314.8 mm × 98.6 mm): a) part of the fracture network and b) its corresponding watershed lines.

Figure 9

 9 Figure 9 Fragment of UMo image, where fractures belonging to the chosen segment are identified. a) Superposition of a chosen segment (in beige) and arcs (in red, dilated for better visibility with a hexagon structuring element of size equals to two) retained for measurement of the aperture average assigned to the segment barycenter. The fragment's dimensions are 138.6 mm × 97.4 mm b) Fractures in red are retained for measurement of the aperture average assigned to the segment barycenter. The fragment's dimensions are 134.0 mm × 96.6 mm.

Figure 10

 10 Figure 10 Maps used to introduce anisotropy when constructing fracture network equivalent tessellations. a) Map of arrival times of solidification front. The solidification fronts were analyzed for all arrival times, from 9 e03 to 54 e03 s after the start of glass being poured into a container, in order to produce maps of the rotation angles and curvature of the front lines. b) Map of the fronts line rotation angles. c) Map of the normalized curvature, calculated by using the formula of the Menger curvature for each triple of set of three points of solidification front.

  defined by X i , Y i ; a M i , b M i , & * 5 . It contains all points M which are closer to the seed P i than to any other seed, according to the parameters a M i , b M i , & * 5 . Here, X i , Y i are coordinates of a cell center, a M i , b M i are scaling factors, and & * 5 is the rotation angle defined in all points M.

Figure 11

 11 Figure 11 Result of the fracture network extraction (image UMo: dimensions are 1044 mm × 413 mm).

Figure 12

 12 Figure 12 Size distribution curves calculated with a) morphological opening and b) opening by erosionreconstruction for the 8 bit input UMo image.

Figure 13

 13 Figure 13 Histograms of apertures computed via distance and quench functions and Feret diameter calculation.

Figure 14

 14 Figure 14 Results of the UMo image segmentation, where the dimensions are 1044 mm × 413 mm. For better visibility, the watershed lines are eroded with size = 3.

Figure 15

 15 Figure 15 Segmented UMo image in which the segments are filled with the values of a) the rotation angles and b) the ratio of axis lengths (image dimensions are 1044 mm × 413 mm).

Figure 16

 16 Figure 16 Scatter plot between the values of front line curvature and those of segment length ratio.

Figure 17

 17 Figure 17 Analysis of A/B ratio of mesh scaling factors. a) Segmented image in which the segments are filled with the values of ratio axes lengths (UMo image). b) Upper half of the map of A/B ratio of mesh scaling factors obtained by kriging with external drift. The lower half ( Y< represents the mirror reflection of the upper part (Y> 1215).

Figure 18 a

 18 Figure 18 a) Histogram, b) omni-directional variogram, and c) normal score transform of aperture averages computed inside the segments of the segmented image.

Figure 19

 19 Figure 19 Results of the computation of point density for different sizes 23 of moving windows: a) regression of density against normalized surface (in number of blocks, each block covers 2500 px 2 ); b) regression of density against arrival time of solidification front ; c) omni-directional experimental variograms of point density (black) and residuals (blue) and their theoretical models; d) histograms of residuals; e) map of point density in which values read at the positions of the segment centers were used afterward as input data for the seed generation.

Figure 20

 20 Figure 20 Construction of equivalent fracture network tessellations: a) positions of seeds obtained by Poisson point process (blue) vs. positions of barycenters of UMo image (brown); b) histograms of aperture average.

Figure 21

 21 Figure21gives the results of the equivalent tessellations generated using the anisotropic Voronoï algorithm for three realizations of seeds and simulated apertures. In order to evaluate the representativeness of the obtained equivalent tessellations, they were compared with the fracture network of UMo and UOx training images presented in Figure23. Moreover, to discern the ability of equivalent tessellations in reproducing a 2D nuclear glass fracture network, the results in Figure22were analyzed.

Figure 21

 21 Figure 21 Equivalent fracture network by Voronoï anisotropic tessellation. Left columns including a), c), e), and g): equivalent fracture network tessellation (dimensions are 6145 px ×2430 px with 1 px corresponding to 0.17 mm. Right columns including b), d), f), and h): porosity (defined as the ratio of the surface occupied by fractures to the sum of the surfaces of the glass matrix and the fracture network) calculated for meshes of the 41 × 17 simulation grid.

Figure 22

 22 Figure 22 Reproduction of segment representation by equivalent tessellation: a) reproduction of the segmented UMo image by anisotropic Voronoï tessellation; b) UMo segmented image in which watershed lines are thickened in accordance with the values of the aperture averages computed inside each segment (dimensions are 6144 px × 2430 px with 1 px corresponding to 0.17 mm); c), d) their corresponding porosities.

Figure 23

 23 Figure 23 Effect of replacement of fracture network by segment representation and its reproduction by equivalent tessellations. Top left a): UMo; top right b): UOx; middle c), d): segmented images; bottom e), f): equivalent fracture network. Equivalent fracture network tessellations are expected to mimic the structures of segmented images, whereas fracture apertures are considered by thickening boundaries of equivalent tessellations

Table 1

 1 Characteristics of the examined images

Name of image Size (px) of input image Resolution (mm) Assembly of N images

	UMo	6144 × 2430	0.17	1
	UOx	6592 × 2856	0.1443	1

Table 2

 2 Statistical characteristics of aperture distribution obtained by discrete and continuous approaches.

	Measurement function /	Max		Min		Mean		Standard deviation
	object of application	# px	mm	# px	mm	# px	mm	# px	mm
	Feret diameter/
	bounding boxes	19	3.23	1	0.17	5.72	0.97	3.01	0.51
	Quench function/
	simple arcs	19	3.23	1	0.17	5.92	1.01	3.60	0.61
	Distance function/
	simple arcs	19	3.23	1	0.17	5.82	0.99	3.60	0.61

http://www.mamba-image.org.html

http://smil.cmm.mines-paristech.fr/doc/index.html

http://rgeostats.free.fr/

This is a vitrocrystalline matrix with 13 wt% molybdenum oxide developed to confine solutions used for UMo fuel treatment

This is ordinary industrial glass for confining fission products arising from UOx fuel treatment

The procedure of multiple point extraction is described in section 1.3.2.2

This function represents the locus of the skeletal points together with their minimal distance to the fracture walls.

The skeleton by homotopic thinning with removed multiple points of the skeleton gives a multitude of disconnected ramifications named simple arcs.

Even if the L homotopic skeleton and the maximal balls skeleton do not always occupy the same position, they still have overlapping points.

According toBarth [9], the state of structural relaxation obtained through fictive temperature quantifies the progress of structural relaxation as a function of thermal history.

A and B are lengths of the main axes of each mesh of equivalent fracture network tessellation.

Two scaling factors (or their ratio as in the example above) if 2D tessellation is build.

Surface of a moving window lying inside an output regular grid.

Procedure is repeated for all maps of density calculated for several sizes of moving window.

Similar to that used to determine the positions of centers, we retained the measurements of the averaged apertures assigned to the centers of segments located in the upper half of the image, and we assigned the same averaged apertures values to the centers of the lower half of the block obtained by mirror reflection.

The lower half of the UMo block was found to be unacceptably damaged. Indeed, several indications of artificial splitting in the fracture network parallel to the block's long wall have been noticed.

The final version of the map includes the reconstructed lower part and the replaced by 1 values of NA.

One moving window is composed of NxN blocks, with each block being 50 × 50 pixels. N is indicated in the header of the column of the figures.

To be precise, fractures are still separated into simple ramifications, but the measurement is applied to the network in which the integrity is preserved, i.e., when we performed aperture measurement by the continuous approach, we apply it to the initial network.

Acknowledgements

The authors wish to thank H. Beucher (MINES ParisTech) and S. Beucher (MINES ParisTech) for fruitful discussions at the early stage of the project and during its implementation. We thank the three anonymous reviewers for their constructive critics, suggestions and comments.

Research data availability

The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study. The data that support the findings of this study will be later available from the corresponding author, Frederic Bouyer, upon request.

Glossary

Binary image

The value of a pixel of this image is either 1 or 0 depending on whether the pixel belongs to the foreground or to the background. Black top-hat transformation Is defined as the difference between the image closing ϕλ(f) and the image itself f, it extracts black and narrow features. In this study, black top hat transformation was applied to grayscale images at the stage of image preprocessing to remove slow trends and to perform contrast enhancement. This operator brought out fractures whose thickness is half the chosen size of the SE and, as a consequence, it facilitated the separation of the glass matrix and the fracture network.

Bounding box

Smallest enclosing box for an object.

Cell

Applied as an analogue of the term "segment" by referring to a separate unit of an equivalent representation of the segmented image. Synonyme: mesh.

Center

Geodesic center of the segment (shard, crystal). It is obtained by the homotopic full thinning with the double structuring element D (Golay's alphabet). Synonym: barycenter.

Centroid

Origin of the corresponding cell growth. Synonym: seed.

Closing

The closing of an image f by a structuring element B is denoted by ϕB(f) and is defined as the dilation of f with a structuring element B followed by the erosion with the reflected SE. The idea behind the morphological closing is to build an operator able to recover the initial shape of the image structures that have been dilated. This is achieved by eroding the dilated image.

Dilation

The dilation of an image f by a structuring element B is denoted by δB(f) and is defined as the maximum of the translation of f by the vectors -b of B. Like for the erosion, the dilated set is the locus of points where the answer to the question is affirmative, however the question slightly changes and becomes: "Does the structuring element hit the set?"

Distance function Defined as the length of the shortest path drawn on the Euclidian grid between two points. In the digital version, each section at level i of the distance function of a set X corresponds to the erosion of size i-1 of X.

Equivalent tessellation

One realization of the equivalent fracture network constructed by the anisotropic Voronoï algorithm.

Erosion

The erosion of an image f by a structuring element B is denoted by εB(f) and is defined as the minimum of the translations of f by the vectors -b of B. Or to put it simply, the eroded set is the locus of points where the answer to the question "Does the structuring element fit the set?" is affirmative.

Feret diameters Distances between the parallel planes restricting the object. Fracture aperture Defined as the width of the fracture between its two opposite walls on the 2D image. For peripheral part of the longitudinal cross-section, it represents a true fracture width, since here fracture network is perpendicular to the cross-section plane. For central part of the longitudinal cross-section, it differs from the true fracture width by cosα, where α is the angle between the plane parallel to the longitudinal cross-section and the plane orthonormal to the fracture whose aperture is measured. Synonymes: aperture, fracture opening.

Geodesic reconstruction

Allows the reconstruction of all connected components of an image g from a marker image f. It is effectuated by iterating geodesic dilatations of f with respect to g until stability. The opening by reconstruction, in contrast to simple opening, preserves the contours of the objects that have not been removed by the erosion (binary image) and removes unconnected light features (grayscale image).

Grayscale image

The range of the values of the pixels of this image is extended from 0 to 2n -1 for pixel values coded on n bits. Image segmentation Its partition into different regions, each one having certain properties.

L-skeleton by homotopic thinning

Obtained by performing the sequential thinning with six rotations of L structuring element (Golay's alphabet) until stability. This skeleton is always connected, of unit thickness and is less prone to generating spurious branches due to irregularities of object frontiers.

Moving window

Rectangle defined by 4 points (x 1 ,y 1 ; x 2 ,y 1 ; x 2 ,y 2 ; x 1 ,y 2 ) which includes certain number of samples i so that their coordinates (x i , y i ) respect the conditions: x 1 ≤ x i ≤ x 2 and y 1 ≤ y i ≤ y 2 .

Multiple points

Used for detecting connections between one-pixel width curves such as those obtained after the homotopic skeletonisation. They were useful in separating fracture ramifications.

Opening

The opening γ of an image f by a structuring element B is denoted by γB(f) and is defined as the erosion of f by B followed by the dilation with the reflected SE.

The idea behind the morphological opening is to dilate the eroded image to recover as much as possible the original image.

Pullout zone

Zone of damage, created during sawing operations.

Quench function

Calculated together with the skeleton by maximal balls, associates the radius of the ball to the ball center.

Rotation matrix

Matrix used by the anisotropic version of Voronoï algorithm to perform a rotation in Euclidean space for defining the anisotropy on distances.

Scaling factor

Number by which the distance along the chosen axis (after possible rotation) is divided.

Seed

Origin of the corresponding cell growth. Synonym: centroid.

Segment

In mathematical terms, it represent a disjoint nonempty set and it is obtained by image segmentation. It the specific case of this research it represents a part of the glass matrix separated by the extracted early fracture network. Synonyms: shard, crystal.

Skeleton

Compact representation obtained by thinning of the object to a set of lines condensing the information of the original object while keeping its homotopy (topology preservation). The skeletal points and their distance to the boundary of the set enable the reconstruction of the set. There are several techniques to perform skeletonisation. In this research, we used a skeleton by maximal balls and a skeleton by homotopic thinning.

Skeleton by maximal balls

Defined by the centers of maximal balls, where a maximal ball B is considered maximal in X object if there exists no other ball included in X and containing B.

Solidification front Corresponds to the isoline of the delay ΔT= T f -T(t)=10 between the fictive temperature (T f ) and the glass temperature at each time step T(t).

Structuring element (SE)

Set of known shape used to probe the analyzed image. There are several types of structuring elements: hexagon, triangle, linear segment etc. The shape of the SE is usually chosen according to some preliminary knowledge about the geometry of the relevant and irrelevant image structures.


MD-PhD Francesco Macri
email: francesco.macri@yahoo.it

PhD b ; Joel Greffier

MD-MS Elina Khasanova

MD-PhD Pierre-Geraud Claret

MD-PhD Sophie Bastide

MD b ; Ahmed Larbi

MD c ; Xavier Bobbia

PhD Fabricio R Pereira

Jean-Emmanuel De

La Coussaye

Jean P Beregi

MD-PhD Jean De La Coussaye

; Jean

MD-PhD Paul Beregi

Minor Blunt Thoracic Trauma in the Emergency Department: Sensitivity and Specificity of Chest Ultralow-Dose Computed Tomography Compared With Conventional Radiography

Keywords: Minor blunt thoracic trauma, Ultra-low-dose chest CT, Chest X-ray, Emergency room, Radiation safety

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Background

Thoracic injuries are frequent, with significant morbidity and mortality.

Management of patients mostly depends on the mechanism of the injury, age of the patient and gravity of the symptoms [START_REF] Rodriguez | NEXUS Chest : validation of a decision instrument for selective chest imaging in blunt trauma[END_REF][START_REF] Langdorf | Prevalence and clinical import of thoracic injury identified by chest computed tomography but not chest radiography in blunt trauma : multicenter prospective cohort study[END_REF]. The Vittel criteria [START_REF] Babaud | Benefit of the Vittel criteria to determine the need for whole body scanning in a severe trauma patient[END_REF] are part of a triage algorithm used by first response crew in France for managing trauma patients prior to admission. The diagnosis of severe trauma is made based on the presence of one single criterion (other than predisposition), which results in the most critical patients being directed toward a specialized traumatology center. A "severe" trauma patient, according to the Vittel criteria, corresponds to a patient with at least one potentially life-threatening or disabling injury, or who sustained trauma whose mechanism or severity suggests that such lesions may exist. When no criterion is present, Chest X-ray (CXR) is the first-line imaging for suspicion of chest lesion in thoracic trauma. Although Computed Tomography (CT) has a higher sensitivity and specificity than CXR in the evaluation of these patients [START_REF] Langdorf | Prevalence and clinical import of thoracic injury identified by chest computed tomography but not chest radiography in blunt trauma : multicenter prospective cohort study[END_REF][START_REF] Exadaktylos | Do we really need routine computed tomographic scanning in the primary evaluation of blunt chest trauma in patients with "normal" chest radiograph?[END_REF], it is usually used in second-line because of its higher irradiation.

Importance

In the recent years, technological improvements such as iterative reconstruction (IR) have become available and have improved the radiation safety of the patients in CT. The IR reduces the noise of a CT image generated with a lowered dose [START_REF] Greffier | Dose reduction with iterative reconstruction: Optimization of CT protocols in clinical practice[END_REF]. The IR allows an ultra-low-dose-CT (ULD-CT), with a dose in the range of the CXR, while preserving a diagnostic image quality despite a noisier image [START_REF] Macri | Value of ultra-low-dose chest CT with iterative reconstruction for selected emergency room patients with acute dyspnea[END_REF]. Hence the ULD-CT can obviate excessive radiation dose exposure of a standard CT for the evaluation of minor-blunt-thoracic-trauma (MBTT) patients.

Objective

The primary objective of this study was to compare the sensitivity and specificity (globally and per lesion type) of chest ULD-CT against CXR for MBTT patients in the emergency department, the diagnostic performances being evaluated against a reference CT (Ref-CT).

Materials and Methods

Study design

This was a cross-sectional single-center diagnostic study. Patients were prospectively recruited from October 2016 to March 2017. The Institutional Review Board of our institution (Level I trauma center) approved this study (N°1886494), providing that total irradiation per patient (CXR + CTs) did not exceed the national diagnostic reference level (DRL) of a chest CT [START_REF] Government | Arrêté du 24 octobre 2011 relatif aux niveaux de référence diagnostiques en radiologie et en médecine nucléaire[END_REF].

Population

Participants were identified in our emergency and traumatology department and were included in the availability periods (convenience sample) of the ULD-CT scan (8AM to 6PM from Monday to Friday). Inclusion criteria were: hemodynamically stable patients older than 18 years that sustained a kinetic energy trauma (motor vehicle accident, fall down the stairs, etc.). The exclusion criteria were patients with unstable hemodynamic conditions (fall of the blood pressure >30% or systemic pressure <110mmHg), dyspnea, with loss of consciousness or with at least one Vittel criterion [START_REF] Babaud | Benefit of the Vittel criteria to determine the need for whole body scanning in a severe trauma patient[END_REF].

Settings

Patients systematically underwent the CXR as a first step, with rib series when rib fractures were suspected by emergency physicians or radiologists (125kVp and 70kVp, respectively; Digital Diagnostic X-ray room; Philips, Amsterdam, Netherlands), immediately (<30 minutes later) followed by a CT examination as a second step. Chest CT images were performed on a 64-Multidetector-computed-tomography (Definition AS+, Siemens Healthineers, Forcheim, Germany). Patients lay in supine position for a Ref-CT (100kVp, 60mAs, CareDose-4D activated, pitch 1.5) and an ULD-CT (100kVp, 10mAs fixed, pitch 1.5). Both acquisitions were acquired consecutively without intravenous injection of contrast media. The tube current parameters were fixed for ULD-CT to be below the national dose reference level for a CXR (Effective dose (E): 0.225 mSv) [START_REF] Government | Arrêté du 24 octobre 2011 relatif aux niveaux de référence diagnostiques en radiologie et en médecine nucléaire[END_REF]. ULD-CT images were reconstructed with IR (Sinogram-

Affirmed-Iterative-Reconstruction [SAFIRE]).

Data collections

Two independent study radiologists (SRs), with nine (F.M. The subjective image quality of the CXR and ULD-CT was evaluated using a 5point Likert scale: 1, excellent; 2, good; 3, fair; 4, poor; 5, unacceptable [START_REF] Macri | Value of ultra-low-dose chest CT with iterative reconstruction for selected emergency room patients with acute dyspnea[END_REF]. SRs expressed for each patient and for the two modalities a diagnostic confidence level (DCL) as certain or uncertain.

The radiation dose was evaluated by calculating the effective doses (E). Dose-Area-Product (DAP) for the CXR and the Dose-Length-Product (DLP) for the ULD-CT and the Ref-CT were recorded from the dose report of each patient and converted to E. E was calculated by multiplying the DLP and the DAP by normalised effective dose conversion coefficient for the chest: 0.014 mSv.mGy - 1 •cm -1 for the CT and 0.00018 mSv.mGy -1 .cm -2 for the radiography.

Statistical analyses

The sensitivity and specificity of ULD-CT and CXR were estimated, with 95% of confidence intervals [95% CI], against the Ref-CT (reference standard) globally and per lesion type. SR1 results were used for diagnostic performance estimations. Concordance between SRs (inter-observer agreement) was assessed with kappa and weighted kappa statistics [with their 95% CI] globally and per lesion type. Hodges-Lehmann estimator was computed to test the effective dose differences between ULD-CT versus CXR and ULD-CT versus Ref-CT. All analyses were performed using R software using a two-sided type 1 error rate of 5% as a threshold for statistical significance. The sample size was calculated on the expected difference between the sensitivities of ULD-CT and 8 CXR. With a hypothesis that 60% of the study population would present at least one lesion, 154 patients are necessary to show a sensitivity increase from 80% to 95% (with a two-sided type 1 error rate of 5% and a 90% power). To be conservative, 160 patients were planned to be included.

Results

Patients

A total of 160 patients were included and all evaluated by CXR and CTs (14 minutes average time lapse between the two examinations; range from 5 to 30) to detect chest lesion. Patients mean age was 65±21 years old with a mean body mass index (BMI) of 23.9±5.4 kg.m -2 . Of 160 patients, 98 (61%) were men (age 61±20 years; BMI 23.8±5.8 kg.m -2 ) whereas 62 (39%) were women (age 69±22 years; BMI 24.0±5.0 kg.m -2 ).

Image assessment

Among these 160 patients, Ref-CT showed pathological findings in 104 patients (65%) with 187 lesions. A single lesion was found in 52 (50%) of the patients. Two lesions were found in 32 patients (31%). A total of 16 patients (15%) presented three lesions and four patients (4%) had at least four lesions.

∆

Limitations of the study

The first limitation of the study concerns practices that differ between centers in terms of trauma management. Our population, older and with a low BMI, might not reflect the population from other centres. Definition of pleural or pericardial effusions could be assessed differently, depending on the size of the effusion and the impact for the management of the patient. Difference between small or extensive effusion, number of rib fractures, or clinical impact on patients was not recorded. ULD-CT cost is higher than a CXR. A comparison with standard chest CT (5 to 7mSv) was not done as for four years we have been routinely performing the low-dose chest CT (1 to 2 mSv) whose diagnostic reliability is promoted by a vast literature. In spite of the rigorous method of analysis between the two SRs, data were anonymised rather than de-identified, which might have influenced SRs reviewing.

Discussion

This study demonstrated that a chest ULD-CT with IR allows a more indepth assessment than a CXR for MBTT patients, delivering a radiation dose lower than the national CXR DRL. Moreover, the radiologist with a short training (SR2) in ULD-CT imaging was more confident on such CT imaging, despite the modification of the image quality due to the lowered dose and IR application, than on standard CXR.

The SRs found more injuries related to MBTT on ULD-CT than on CXR. Small pneumothorax and pneumomediastinum should not be underestimated because of the incidence of aggravation and potential broncho-digestive damage [START_REF] Rezende-Neto | Occult pneumomediastinum in blunt chest trauma: clinical significance[END_REF], respectively. ULD-CT, compared to CXR, had a very high sensitivity and specificity in detecting pulmonary contusions. The ULD-CT allowed a thorough The major issue of CT dose reduction is reduction in image quality due to increased noise. SAFIRE decreases the noise but when high strengths are applied, a smoothing effect can occur, which can initially confuse a nonaccustomed radiologist. SRs had a fair inter-observer agreement for ULD-CT. SR1 scored the majority of cases excellent for image quality whereas SR2 rated them good. SR2 was trained for ULD imaging specifically for this study, whereas SR1 had dealt with ULD for three years. It is noteworthy that SRs never considered the image quality poor or unacceptable for ULD-CT. SR2 and SR1 had a higher DCL on ULD-CT than on CXR (96% vs 80%). This result highlights that even though SR2 was new and less satisfied with the image quality of ULD-CT, finally proved more confident on ULD-CT than on standard radiographs.

In conclusion, ULD-CT with IR allows a thorough study of MBTT patients with lower radiation-induced cancer risks, suggesting a new approach for the MBTT work-up. Other studies should be conducted to define whether this initial major expense, which allows a more reliable and quick assessment, avoids further elevated health costs in the management of MBTT patients.

  [SR1]) and six (E.K. [SR2]) years of experience as board certified radiologists, independently assessed anonymized CXR and ULD-CT datasets on clinical workstations with PACS viewer system (Centricity, GE Healthcare, Wisconsin, USA) in a random fashion and with at least two weeks between each data reading session. SRs were not aware of clinical data of the patients. SR1 had three years of experience with ULD-CT images whereas SR2 had no experience in chest ULD-CT (only two-week training period in view of this study). Ref-CT images were read by the scheduled radiologists during the regular shifts (excluding SR1 and SR2) who were informed of the final diagnosis. SRs evaluated the imaging to detect the presence of pneumothorax, subcutaneous emphysema, pneumomediastinum, pulmonary contusion, fractures of the thoracic skeleton (thoracic spine, rib cage, clavicles and sternum), pleural effusion, pericardial (>2mm thickness) and pleuro-pericardial effusions. The primary outcome was the presence of at least one lesion. Number of rib fractures or severity of effusions were not evaluated.

Figure 1

 1 Figure 1 describes all the lesions and whether they were detected by each SR for

Finally

  the inter-observer agreement for DCL was certain for CXR and ULD-CT respectively (CXR: k=0.92 [0.84 -1.00]; ULD-CT: k=0.83 [0.59 -1.00]). Regarding the dose level, ULD-CT was significantly lower than Ref-CT ( ∆ =0.93; [0.86 -1.00]) but not significantly different from CXR ( =0.06; [0.10 -0.20]).

  study of the rib cage, thoracic spine and occasionally of the shoulder girdle. The versatility of 2D/3D CT images and their multiplanar reconstructions avoided the inconvenience of structure superimposition generated by plain radiographs, with misleading interpretations. Rib fractures are extremely painful during regular breathing and if overlooked they are a critical factor of mortality increase in elderly patients, especially with comorbidity [9]. Overall, SRs counted more pleural effusions on ULD-CT than on CXR. The pericardial effusions were seen only on ULD-CT due to low volumes making them non-visible with CXR. Pleural and pericardial effusions are frequent in this study due to conservative cut-off for effusion (definition >2mm) with unpredictable clinical course [10]. Effusions are not necessarily hemorrhages and could have been present prior to the trauma. Also, elderly individuals with low BMI form a particular population, more fragile to traumas involving minor mechanisms of injury. Sensitivity and specificity values of the ULD-CT were higher than those of CXR, in keeping with the published literature for the standard CT [4].

Figure legend Figure 1 .

 legend1 Figure legend

Table 1

 1

	presents the results of sensitivities

Table 2 .

 2 The inter-observer agreement for image quality was moderate (k=0.26 [0.17 -0.35]) for ULD-CT and good (k=0.85 [0.74 -0.96]) for CXR. The DCL on CXR images was rated certain in 81.8% [75.9 -87.9] and 79.3% [73.1 -85.7] of cases, respectively for SR1 and SR2. The DCL for ULD-CT was certain in 96.8% [94.2 -99.6] and 95.6% [92.4 -98.8] of cases for SR1 and SR2 respectively.

Table 2 .

 2 Subjective image quality of chest X-Ray and chest ultra-low-dose CT Study radiologist (SR) 1 and 2 with inter-observer agreement (weighted kappa estimations with 95% Confidence Interval). CXR: chest X ray; ULD-CT: ultra-lowdose computed-tomography.

			CXR			ULD-CT
		SR1 SR2	Kappa	SR1 SR2	Kappa
	Excellent	140 135		132 59
	Good	9	14		21 94
	Fair	5	6	0.85	7	7	0.26
	Poor	6	5 [0.74 -0.96] 0	0 [0.17 -0.35]
	Unacceptable 0	0		0	0

Acknowledgements

We thank Sarah Kabani for the English revision of the paper. The authors state that this work has not received any funding.

Author contributions.

-Francesco Macri; design and responsible for authorization and implementation of the study, one of the blind readers for images, writing