J Fruth
email: fruth@statistik.tu-dortmund.de

O Roustant
email: roustant@emse.fr

T Muehlenstaedt

The fanovaGraph Package: Visualization of Interaction Structures and Construction of Block-additive Kriging Models

Keywords: Sensitivity analysis, FANOVA decomposition, Sobol indices, High-order interactions, Superset importance, Additive structure

FANOVA graphs have a broad range of application in model interpretation, sensitivity analysis and analysis of computer experiments. They give a clear and easy to understand visualization of function structures, reveal blockadditive decompositions and also effectively improve Kriging model predictions by kernel adaption. These several fields are combined in a new released R package fanovaGraph. This tutorial presents the implemented methods step by step together with theoretical background and illustrations. Additionally it includes methods in simulation and the crucial question of threshold decision, where graphical solutions are suggested.

Introduction

Many computer experiments behave like black box functions in the sense that relations between input and output factors are not obvious from the function. If the experiment is expensive, usually a meta model like a Kriging emulator (but also other models like Neural Networks or Random Forests) is estimated from a low number of runs as a replacement for the experiment. Then again, the meta model provides a function which needs to be explored. To analyse and interpret those black box functions, common approaches are scatterplots or sensitivity analysis methods, which become problematic in the presence of interactions (see for instance [START_REF] Saltelli | Sensitivity analysis. Wiley series in probability and statistics[END_REF] for a review). This problem has recently been approached by using graphical displays of the interaction structure (see [START_REF] Hooker | Discovering additive structure in black box functions[END_REF] and [START_REF] Muehlenstaedt | Datadriven Kriging models based on FANOVA-decomposition[END_REF]).

The R package fanovaGraph [START_REF] Fruth | fanovaGraph: Building Kriging models from FANOVA graphs[END_REF] is an efficient tool for a clear and easy to understand graphical visualization of the global interaction structure of a black box function. Interactions are depicted by a so-called FANOVA graph, whose edge thicknesses indicate the interaction strength. Estimators for these interactions, based on the superset importance estimation in [START_REF] Liu | Estimating mean dimensionality of analysis of variance decompositions[END_REF], have been shown to be asymptotical normal and efficient in [START_REF] Fruth | Total interaction index: A variance-based sensitivity index for interaction screening[END_REF].

Following the procedure suggested in [START_REF] Muehlenstaedt | Datadriven Kriging models based on FANOVA-decomposition[END_REF], the detected interaction structure (block-additive decomposition) by the FANOVA graph can further be used to construct block-additive Kriging models. The package fanovaGraph therefore contains in a second part all methods for block-additive Kriging analysis based on the package DiceKriging [START_REF] Roustant | DiceKriging: Kriging methods for computer experiments[END_REF]. This includes estimation, prediction and simulation of those models.

The paper is organised as follows. At first the main procedure is demonstrated at an example function allowing an easy access and short overview. Section 3 then gives detailed information about the procedure's steps and explains the statistical background. New confidence intervals for the estimators are presented. In Section 4, the question of choosing the right threshold value for the graph edges is adressed and a graphical decision plot for the threshold is suggested. Finally Section 5 shows how further models, here Neural Networks and Support Vector Machines can be analyzed by the package.

Throughout this paper, we use the following six dimensional function as a toy example for the "unknown" black box function.

f (x) = cos((1, x 1 , x 3 , x 5 )β T ) + sin((1, x 2 , x 4 , x 6 )γ T )
with β = (-0.8, -1.1, 1, 1.1) and γ = (-0.5, 1, 0.9, -1.1) and x ∈ [-1, 1] 6 . The function has the block-additive form

f (x) = g(x 1 , x 3 , x 5 ) + h(x 2 , x 4 , x 6 ). > d <-6 > domain <-c(-1, 1) > fun <-function(x) { + beta <-c(-0.8, -1.1, 1, 1.1) + gamma <-c(-0.5, 1, 0.9, -1.1) + g <-cos(cbind(1, x[, c(1, 3, 5)]) %*% beta) + h <-sin(cbind(1, x[, c(2, 4, 6)]) %*% gamma) + return(g + h) + }

A demo

Before detailing the functions of fanovaGraph, it is worth giving an overview of their usage in a practical case study. We presume basic knowledge about Kriging and sensitivity analysis and refer to [START_REF] Saltelli | Sensitivity analysis. Wiley series in probability and statistics[END_REF], [START_REF] Fang | Design and modeling for computer experiments[END_REF].

There are two main goals of fanovaGraph: 1. To provide a description of the interaction structure of a black box function, say f , visualized in a so-called FANOVA graph, which at the same time gives a block-additive decomposition of the function. 2. To use the block-additive decomposition given by the FANOVA graph to build an advanced Kriging type metamodel f app .

Of course the steps 1 and 2 can be done independently: When only the interaction structure of the function f is of interest (step 1) or when external information is available to construct a block-additive kernel (step 2). This section shows the frequent situation where both steps 1 and 2 are done sequentially: the modeling of a costly black box function, e.g. a finite element computer experiment. This situation is presented carefully in fanovaGraph by the way of two demos. We detail it here by means of the 6D example presented in Section 1.

When dealing with a costly function, even the first step requires a metamodel of f . Such functional approximation can be built in different ways, by using linear models, neural networks, polynomial chaos expansion, Kriging models (or equivalently Gaussian process regression), splines, etc. In the following, the metamodel is based on a Kriging model built from the package DiceKriging based on a maximin design sampled with the package lhs [START_REF] Carnell | lhs: Latin Hypercube Samples[END_REF], but other solutions might have been chosen as well. The code corresponding to the initialization step is shown below.

> library(lhs) > L01 <-maximinLHS(100, d) > x <-L01 * (domain [START_REF] Fang | Design and modeling for computer experiments[END_REF] -domain [START_REF] Carnell | lhs: Latin Hypercube Samples[END_REF]) + domain [START_REF] Carnell | lhs: Latin Hypercube Samples[END_REF] > y <-fun(x) > KM <-km(~1, design = data.frame(x), response = y)

Obviously, the accuracy of the first metamodel should be verified before continuing further. This step is not shown here, for the sake of brevity. The model's prediction function now serves as a metamodel. It is in fanova-Graph implemented by the kmPredictWrapper function, based on predict from DiceKriging. The model (here KM) must be given as an extra argument.

The usual first step of a sensitivity analysis, estimation of first order Sobol indices, can then be done with the package sensitivity [START_REF] Pujol | sensitivity: Sensitivity Analysis[END_REF]. The result, visualized in Figure 1, shows that all 6 variables are influent, and indicates the presence of interactions. > i1 <-fast99(model = kmPredictWrapper, factors = d, + n = 3000, q = "qunif", q.arg = list(min = domain [START_REF] Carnell | lhs: Latin Hypercube Samples[END_REF], + max = domain [START_REF] Fang | Design and modeling for computer experiments[END_REF]), km.object = KM) > plot(i1) However, this analysis does not show the way how variables interact and the functions two additive blocks are not revealed. To do so, we estimate a FANOVA graph. This is done in fanovaGraph by the estimateGraph function. The obtained graph object g can then be plotted: > g <-estimateGraph(f.mat = kmPredictWrapper, d = d, n.tot = 30000, + q.arg = list(min = domain [START_REF] Carnell | lhs: Latin Hypercube Samples[END_REF], max = domain [START_REF] Fang | Design and modeling for computer experiments[END_REF]), km.object = KM) > plot(g, plot.i1=FALSE)

The corresponding figure can be seen in Figure 2, on the left. Clearly, two groups of variables are visible, and can be identified by choosing an appropriate threshold. How to choose the threshold value is detailed in Section 4. Here δ = 0.01 seems to be a suitable value, as shown on the same figure, on the right. We can use the information for the second part to build up a blockadditive Kriging model. Mathematically, the two unconnected (additive) parts are the two cliques of the graph g

> Cliques <-g.cut$cliques > print(Cliques) [[1]] [1] 1 3 5 [[2]] [1] 2 4 6
The clique structure suggests that a good kernel for this function, among kernels of stationary processes, has the block-additive structure:

k(h) = k 135 (h 135 ) + k 246 (h 246 )
with the set notation h I = {h i , i ∈ I} for all set I, and where k 135 , k 246 are kernels defined over 3-dimensional subspaces. Several kernel types are available, and we choose for both the Gaussian tensor-product kernel, assuming that f is a smooth function. The estimation of the corresponding Kriging model can now be done with the MLoptimConstrained function of fanovaGraph, and predictions computed with yhat. Finally, the results are compared to the 'standard' Kriging model, based on a 6D Gaussian tensor-product kernel, and confirm the clear improvement in prediction of the new Kriging model (Figure 3). Note that at no point we made any priori assumptions about the block-additive form of the black-box function.

> parameter <-kmAdditive(x, y, cl = Cliques) > xpred <-matrix(runif(d * 1000, domain [START_REF] Carnell | lhs: Latin Hypercube Samples[END_REF], domain [START_REF] Fang | Design and modeling for computer experiments[END_REF]), ncol = d) > y_new <-predictAdditive(xpred, x, y, parameter, cl = Cliques) > y_old <-kmPredictWrapper(xpred, km.object = KM) > y_exact <-fun(xpred) > par(mfrow = c(1, 2)) > plot(y_exact, y_old, asp = 1, xlab="y, exact", + ylab="y, predicted", main="Standard Kernel") > abline(0, 1) > plot(y_exact, y_new[, 1], asp = 1, xlab="y, exact", + ylab="y, predicted", main="Modified Kernel") > abline(0, 1)

-2 -1 0 1 2 -2 -1 0 1 2 Standard Kernel y, exact y, predicted -2 -1 0 1 2 -2 -1 0 1 2

Modified Kernel

y, exact y, predicted

Detailed description of functions

Let us now generally introduce the procedure and its implementation in fanovaGraph. Say that the function under investigation f is defined over d factors x 1 , . . . , x d , which can be described as independent random variables with probability measure ν. Usually, if no special distribution is predetermined, a uniform distribution in the factor domain is chosen for ν.

Estimation of indices

The key point in building the FANOVA graph is the estimation of the edge thicknesses. An edge shall represent the interaction strength of the two factors that are represented by the vertices it combines. Therefore a Sobolbased index is chosen, that measures the overall influence of the two factors, called total interaction index (TII) in [START_REF] Fruth | Total interaction index: A variance-based sensitivity index for interaction screening[END_REF]. For two factors x i and x j it is defined by the sum of all Sobol indices that contain both indices:

D ij := I ⊇{i,j} D I , (1)
where D I is the unscaled Sobol index of the group of factors (X k ) k∈I . [START_REF] Liu | Estimating mean dimensionality of analysis of variance decompositions[END_REF] show, that this index can be computationally obtained by

D ij = 1 4 E ∆ i,j (z i,j , x i,j , x -{i,j} ) 2 = 1 4 ∆ i,j (z i,j , x i,j , x -{i,j} ) 2 dν i,j (z i,j )dν(x) (2)
with ∆ i,j (z i,j , x i,j , x -{i,j} )) =

f (x i , x j , x -{i,j} ) -f (x i , z j , x -{i,j} ) -f (z i , x j , x -{i,j} ) + f (z i , z j , x -{i,j} )
and where z i (resp. z j ) is an independent copy of x i (resp. x j ). Since f is assumed to be a black box function, the total interaction index can not easily be computed by [START_REF] Fang | Design and modeling for computer experiments[END_REF]. One way to estimate it is by partly fixed Monte Carlo integration of the integrals. Denote by x k and z k , k = 1, . . . , N two independent samples of (sufficiently large) size N drawn from ν. Then the total interaction index is estimated by

D ij = 1 4 × 1 N N k=1 f (x k i , x k j , x k -{i,j} ) -f (x k i , z k j , x k -{i,j} ) -f (z k i , x k j , x k -{i,j} ) + f (z k i , z k j , x k -{i,j} ) 2 . ( 3
)
In fanovaGraph the TII are estimated by estimateGraph and can be extracted from the output by $ tii or respectively $ tii.scaled for the TII scaled by the overall variance. To get insight about the certainty of the estimates, asymptotic 95% confidence intervals are computed additionally exploiting the estimator's asymptotic normal distribution of the Liu and Owen estimator in (3). Here an example code for the previously defined prediction function of the 6D example. Interactions of variables of the same block (e. g. X1 and X3) are clearly higher than from different blocks: ,d = d,+ n.lo = 1000, q.arg = list(min = domain [START_REF] Carnell | lhs: Latin Hypercube Samples[END_REF], + max = domain [START_REF] Fang | Design and modeling for computer experiments[END_REF]), method = "FixLO", km.object = KM) > head(g$tii) totalInt Std.Error lower upper X1*X2 0.0007741586 6.330495e-05 0.0006500832 0.0008982341 X1*X3 0.0411745532 3.110340e-03 0.0350783986 0.0472707078 X1*X4 0.0010549226 1.050348e-04 0.0008490582 0.0012607871 X1*X5 0.0553378302 4.649570e-03 0.0462248410 0.0644508194 X1*X6 0.0013846099 1.167610e-04 0.0011557626 0.0016134572 X2*X3 0.0009118153 1.316350e-04 0.0006538154 0.0011698151 Due to the consistency of the estimator, the higher the sample number, the more exact is the estimation. In estimateGraph the number of Monte Carlo samples can be specified directly by n.lo or in terms of overall function evaluations by n.tot. Further parameters are the dimension size d, the input distribution q (defaults to uniform distribution) with parameters q.arg and the estimation method method. Beside the one from Liu and Owen, three further estimation methods, based on FAST, RBD-FAST and Sobol's pick and freeze method, are implemented. See [START_REF] Fruth | Total interaction index: A variance-based sensitivity index for interaction screening[END_REF] for a detailed description.

> g <-estimateGraph(f.mat = kmPredictWrapper

Plots

Total interaction indices are easily visualized in the FANOVA graph. Input variables are represented by the graph's vertices and are connected by an edge if the total interaction index is positive. The thickness of an edge is relative to the size of the corresponding total interaction index. With the plotiGraph function, graph objects obtained from estimateGraph can directly be plotted. There are several modifications for the graph plotting. The position of the indices in the plot is determined by the layout argument, adopted from igraph0, which searches by default for the best way. To compare graphs, the layout can be fixed. Standard main effect Sobol indices can be represented additionally by the vertex thicknesses. Moreover, pure second order Sobol indices (D ij ) can be added as additional inner edges since D ij < D ij by definition. They can give insight about the amount of active interactions of orders higher then two.

> plot(g, ploti1 = FALSE, layout = layout.circle) > plot(g, ploti1 = TRUE, layout = layout.circle) > plot(g, ploti1 = FALSE, i2 = c(0, 0.03222, 0, 0.04093, + 0, 0, 0.0135, 0, 0.023, 0, 0.03222, 0, 0, 0.01753, 0), + layout = layout.circle)

Block additive Kriging

As shown in Section 2, the information brought by the FANOVA graph about total interactions can be used to build a suited Kriging model. More precisely, the clique decomposition of the graph provides a block-additive kernel of the form (see [START_REF] Muehlenstaedt | Datadriven Kriging models based on FANOVA-decomposition[END_REF] for details):

k(h) = L i=1 k C i (h C i ) ( 4
)
where C 1 , . . . , C L are the cliques, and k C 1 , . . . , k C L the corresponding kernels, associated to independent stationary centered Gaussian processes Z C 1 , . . . , Z C L . Note, that the cliques can have common vertices. In fanova-Graph, each kernel k C i has a tensor-product structure:

k C i (h C i ) = σ 2 C i j∈C i k 1 (h j ; θ C i ,j ) (5)
where σ 2 C i is the variance of Z C i , k 1 is a 1-dimensional kernel, and the θ C i ,j 's (j ∈ C i ) denote the other covariance parameters. The implementation of k C i relies on DiceKriging. In particular, k 1 can be chosen among a list of stationary kernels (at present time: Exponential, power-exponential, Matérn 3/2, Matérn 5/2 or Gaussian (Default)), according to the smoothness of the function of interest. We refer to [START_REF] Roustant | DiceOptim: Two R packages for the analysis of computer experiments by krigingbased metamodelling and optimization[END_REF] for more details. All covariance parameters (including the variance one) are estimated by maximum likelihood (see [START_REF] Muehlenstaedt | Datadriven Kriging models based on FANOVA-decomposition[END_REF]). However, the number of parameters equals L i=1 (1 + |C i |), and can be large. Hence, a useful option is to constrain the parameters in one clique to have a common value. We call such cliques isotropic. It is often a good idea in practice to consider isotropic cliques, and several examples are provided in [START_REF] Muehlenstaedt | Datadriven Kriging models based on FANOVA-decomposition[END_REF]. To illustrate its usage in fanovaGraph, let us consider again the 6-dimensional function f of Section 2. Looking at the estimated FANOVA graph, it seems reasonable to consider an isotropic kernel for the clique C 1 = {1, 3, 5}, since both total effects and total interaction indices are similar in between the variables of the clique. Hence, we add the constraint θ C 1 ,1 = θ C 1 ,3 = θ C 1 ,5 in equation 5, which gives only 6 parameters to estimate instead of 8. This change is done in kmAdditive by setting the argument iso, which is a vector of Boolean whose TRUE coordinates indicates the isotropic cliques. Looking again at the cliques order (see Section 2):

> print(Cliques) [[1]] [1] 1 3 5 [[2]] [1] 2 4 6
we see that C 1 is indeed the first clique in the list. Hence, to make C 1 isotropic and leave C 2 unchanged, we must choose:

> iso <-c(TRUE, FALSE)
Parameter estimation is then done by writing: > param <-kmAdditive(x, y, cl = Cliques, covtype = "gauss", + iso = iso)

The accuracy of the corresponding Kriging metamodel can be compared to the one based on an anisotropic block-additive kernel (see Section 2). The RMSE value shows that the predictive power looks similar, with the benefit of reducing the number of parameters.

[1] 0.006851003 Actually, for this example, the predictive power is also very good when the two cliques are assumed to be isotropic, which reduces the number of parameters to 4 only : > iso <-c(TRUE, TRUE) > param <-kmAdditive(x, y, cl = Cliques, covtype = "gauss", + iso = iso) > y_block_iso <-predictAdditive(xpred, x, y, param, + cl = Cliques, iso = iso) > sqrt(mean((y_block_iso[, 1] -y_exact)^2)) [START_REF] Carnell | lhs: Latin Hypercube Samples[END_REF] 0.006463344 Finally generating samples from block-additive kernels for simulation purposes can be done by the function simAdd. Here for a given kernel the Kriging covariance matrix R at desired data points x of size n sim is computed together with its Cholesky decomposition R = C T C. A simulation with mean µ and covariance matrix R is then obtained by

y = µ + C T ε with ε i ∼ N (0, 1), i = 1, . . . , n sim independent random numbers.
In the following code two simulations are generated to compare graphically a standard tensor-product kernel and an additive kernel (see Figure 5). > x1 <-x2 <-seq(-1, 1, , 20) > x.grid <-expand.grid(x1, x2) > parameter <-list(list(alpha = 1, theta = c(0.5, 0.5))) > y.sim <-simAdditive(x.grid, mu = 0, parameter, covtype = "matern5_2", + cl = list(1:2)) > persp(x1, x2, matrix(y.sim, 20), theta = -40, col = "lightgreen", + zlab = "y.sim") > parameter <-list(list(alpha = 0.5, theta = 0.5), + list(alpha = 0.5, theta = 0.5)) > y.sim <-simAdditive(x.grid, mu = 0, parameter, covtype = "matern5_2", + cl = list(1, 2)) > persp(x1, x2, matrix(y.sim, 20), theta = -40, col = "lightgreen", + zlab = "y.sim")

x 1

x 2 y.sim (a) Tensor-product kernel

Thresholding decision

The problem

Errors in the initial metamodel and in the index estimation usually pertubate the FANOVA graph, so that inactive indices are not necessarily zero. Therefore a thresholding cut is necessary to separate high from low indices and thus find inactive interactions. This separation can easily be done within the package by applying the threshold function to the graph object with a given cut value delta. See below for a graphical support to find a suitable value. By default, the cut is performed at the scaled indices. This can be changed through the argument scaled=FALSE. Here we perform a threshold cut at δ = 0.01, which sets all indices that are smaller then δ to zero. > g.cut <-threshold(g, delta = 0.01) > head(g.cut$tii) totalInt Std.Error lower upper X1*X2 0.00000000 0.00000000 0.00000000 0.00000000 X1*X3 0.04117455 0.00311034 0.03507840 0.04727071 X1*X4 0.00000000 0.00000000 0.00000000 0.00000000 X1*X5 0.05533783 0.00464957 0.04622484 0.06445082 X1*X6 0.00000000 0.00000000 0.00000000 0.00000000 X2*X3 0.00000000 0.00000000 0.00000000 0.00000000

Decision plot

Decision plots are a popular way to help the user find data dependent values. We present a decision plot called delta jump plot that helps finding a suitable threshold value δ (see Figure 6).

> plotDeltaJumps(g) On the x-axis it shows all relevant threshold values, i.e. the ordered values of indices in the FANOVA graph in equal distances. In the bottom part of the graph, the real distances of the indices are shown, so that high jumps point to big differences between successive indices. In a moderately pertubated graph, the difference between inactive and active indices should be high and thus be revealed by a jump. High jumps are additionally highlighted by shading of the relevant interval, the darker the shading the higher the jump. Values of delta that lie inside dark intervals indicate good choices for the threshold.

0
Simultaneously the upper part of the plot shows the number of cliques a threshold would lead to. A small number of cliques might be preferable in order to obtain a clearer structure and less parameters to estimate. Figure 6 reveals that in this case a suitable δ lies between 0.0021 and 0.025, because the jump to right of this interval is the highest and the number of cliques is very low.

fanovaGraph with further initial metamodels

As mentioned in Section 2 the metamodel is not restricted to Kriging but can be any prediction model, chosen suitable for the underlying problem. Here we show how further fitted models can easily be analyzed by fanova-Graph. Using the same data as before, we exemplary fit a Neural Network model and a Support Vector Machine (SVM) regression model, using standard R functions. First, for Neural Networks the output has to be scaled to lie between 0 and 1.

> data <-data.frame(y,x) > data$y <-(data$y-min(data$y))/(max(data$y)-min(data$y))

Now we fit the models and construct a prediction function wrapper wrapp, that can be used for both functions.

> library(nnet) > library(e1071) > nnet.mod <-nnet(y ~., data=data, size=2) > svm.mod = svm(y ~., data=data) > wrapp <-function(newdata, object) { + colnames(newdata) <-c("X1","X2","X3","X4","X5","X6") + as.numeric(predict(object, newdata = newdata)) + } From the prediction functions the FANOVA graph can be estimated and plotted as before (see Figure 7). Both metamodels discovered the interaction structure of function f very good, though not as good as Kriging. But of course the choice of the metamodel always depends on the situation. > nnet.g <-estimateGraph(f.mat = wrapp, d = d, n.tot = 30000, + q.arg = list(min = domain [START_REF] Carnell | lhs: Latin Hypercube Samples[END_REF], max = domain [START_REF] Fang | Design and modeling for computer experiments[END_REF]), object = nnet.mod) > svm.g <-estimateGraph(f.mat = wrapp, d = d, n.tot = 30000, + q.arg = list(min = domain [START_REF] Carnell | lhs: Latin Hypercube Samples[END_REF], max = domain [START_REF] Fang | Design and modeling for computer experiments[END_REF]), object = svm.mod) > plot(nnet.g, plot.i1=FALSE) > plot(svm.g, plot.i1=FALSE)

Conclusion

The FANOVA graph is a useful tool for visualization of interaction structures of various functions, and can especially be used for Kriging model improvement via block-additive decomposition. The R package fanovaGraph provides a user friendly implementation of FANOVA graphs and their application including estimation, plotting, block-additive Kriging analysis and threshold decision graphs.

One point for further improvement of the package is the metamodel error. Kriging models with second order polynomial trend have proven to better catch the interaction structure than standard Kriging models. Other Kriging kernels can be taken into consideration like FANOVA kernels, where each interaction is modelled explicitely. Also iterative repetition of model fitting and graph fitting could be a possible alternative. Furthermore, other fields of application of the block-additive structure can be implemented like dimension reduction and optimization problems, where special clique structures are necessary.

Figure 1 :

 1 Figure 1: Main effects (white box) and total effects (large box) for the 6D-function f , obtained with the package sensitivity ([8]).

Figure 2 :

 2 Figure 2: The estimated FANOVA graph for the 6D function f .

Figure 3 :

 3 Figure 3: Predictions over a 1000-point uniform test set of two Kriging models: A standard one, based on a 6D Gaussian tensor-product kernel, and a new one, based on a sum of two 3D Gaussian tensor-product kernels k 135 + k 246 .

Figure 4 :

 4 Figure 4: Different modifications of FANOVA graph plots.

Figure 5 :

 5 Figure 5: Simulations from Kriging kernels.

Figure 6 :

 6 Figure 6: Delta jump plot.

Figure 7 :

 7 Figure 7: Estimated FANOVA graph for two other metamodels on the 6D function f .

Acknowledgements

This paper is based on investigations of the collaborative research centre SFB 708, project C3. The authors would like to thank the Deutsche Forschungsgemeinschaft (DFG) for funding. We also thank Sonja Kuhnt, Simon Neumaerker and Stefan Hess for their helpful support.


Simon Rouchier

Hans Janssen

Carsten Rode

Monika Woloszyn

Geneviève Foray

J.-J Roux

Jean-Jacques Roux

Characterization of fracture patterns and hygric properties for moisture flow modelling in cracked concrete

Keywords: concrete, fracture, digital image correlation, x-ray radiography, moisture transfer modelling

à la diffusion de documents scientifiques 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

The hygrothermal performance of porous building materials is now an important issue of building physics. Moisture transfer has a critical influence on the durability and sustainability of built structures, as well as on the health and comfort of their occupants. Any analysis of moisture transfer requires dependable moisture transfer properties, for which various experimental and numerical methods are currently available. However, the understanding of how these properties evolve during the lifespan of building materials is far from complete. The pore structure of the material itself may change over time, or be altered by cracks and defects. Several years after their installation, building materials such as concrete present signs of ageing in the form of fractures covering a wide range of sizes, 1 from microscopic to macroscopic cracks. All sizes of fractures can have a strong impact on heat and moisture flow in the building envelope, and their influence is to be accounted for in any performance assessment, not only of building and building components, but of any built structure in general. The infiltration rate of moisture is known to have an important role in the degradation of buildings [START_REF] Torres | The influence of the thickness of the walls and their properties on the treatment of rising damp in historic buildings[END_REF].

Fractures in building materials are caused by a large variety of external factors: mechanical loading, temperature and moisture content fluctuations, infiltration of chemicals such as carbon dioxide causing carbonatation in concrete [START_REF] Černý | Transport processes in concrete[END_REF] or chloride ions causing corrosion. An accurate prediction of all ranges of damage in a building component, from microscopic to macroscopic cracks, would either require to include all these effects in a simulation or to assume important simplification hypotheses. This level of complexity being difficult to model, an other approach might be considered when studying the hygrothermal performance of an ageing building component. Fluid flow modelling in a porous medium including a discrete fracture network can be performed following two main methodologies: either by considering an already existing fracture pattern, or by coupling a mechanical model to the fluid flow model. The latter is a coupled mechanical and hygrothermal modelling methodology which consists in predicting the propagation of localized [START_REF] Roels | A coupled discrete-continuum approach to simulate moisture effects on damage processes in porous materials[END_REF][START_REF] Moonen | Continuous-discontinuous modelling of hygrothermal damage processes in porous media[END_REF] or smeared [START_REF] Ožbolt | Modelling the effect of damage on transport processes in concrete[END_REF] cracks during fluid ingress in the porous medium eventually including poromechanical effects. This allows the study of the mechanical behaviour of the material under hygrothermal influences, but is difficult to achieve in the case of building materials due to the complexity of the crack-inducing phenomena. The other option is to run moisture transfer simulations in fractures of known geometry. Setting the crack patterns prior to fluid flow modelling can allow to estimate equivalent properties of the fractured porous network at a higher scale [START_REF] Réthoré | A two-scale approach for fluid flow in fractured porous media[END_REF], to compare the performance of several modelling methodologies [START_REF] Alfaiate | On the use of strong discontinuity formulations for the modeling of preferential moisture uptake in fractured porous media[END_REF] or study the propagation of the crack under hygric influences. The present paper is based on an extension of this approach: the fracture geometries in concrete samples were experimentally characterized before moisture transfer modelling was applied. This method provides the exact distribution of the fracture network without the need for an extensive mechanical modelling and allows estimating the hygrothermal performance of a damaged building component.

The present work serves two main purposes: the first one is to provide experimental measurements of the fracture network inside samples of concrete in order to assemble a complete validation set for 2D models of moisture flow in a fractured heterogeneous porous material. The second purpose is to estimate the ability of non destructive observation techniques to provide the geometrical characteristics of the fractures which may have an influence on the moisture flow.

Digital image correlation (DIC) was performed during the fracturing of concrete samples, in which moisture uptake was then monitored with X-ray radiography. Concrete is a composite material made of porous, permeable cement paste binding aggregates of various sizes with very low permeability. Therefore, in such a material, fractures caused by mechanical loading have a tortuous shape and variable aperture due to the presence of aggregates or initially distributed shrinkage cracks. DIC allows to quantify these geometrical characteristics by providing a complete mapping of all displacements at the surface of the sample. The samples are then put in contact with water in order to monitor its infiltration by X-ray radiography inside the characterized fracture network. The resulting moisture concentration profiles serve as a basis for the validation of a finite-element model, which is set to automatically detect the fracture positions from the optical measurements and to simulate moisture flow through the fractured porous material.

The paper is organized as follows: the complete methodology is first entirely outlined starting from the elaboration of the specimen. The experimental and numerical procedure for the quantification of the fracture network, the observation of water infiltration and the moisture transfer modelling in cracked concrete are then described. The outcome of the whole methodology is then displayed and the simulation results are compared to measured moisture profiles.

Methodology

A formulation of steel fibre reinforced concrete was prepared, based on the proportions shown in table 1, using Aalborg portland cement with a water to cement ratio of 0.5. It follows the proportions from previous studies [START_REF] Pease | Influence of concrete cracking on ingress and reinforcement corrosion[END_REF] performed at the Technical University of Denmark. The steel fibres have a length of 12.5 mm and a diameter of 0.4 mm. Their purpose is to allow the interruption of the mechanical loading before complete failure of the specimen, and to ensure their post-crack cohesion in order to use them for subsequent moisture uptake tests.

Water/cement ratio 0. The mixture is cast in prismatic forms of dimensions 350×100×100 mm, and first wrapped in plastic sheets during 48h after casting as to initially maintain a high relative humidity (RH), then stored during 28 days in a stable environment of 50% RH. After maturation, the edges of the beams were removed since their cement/aggregate ratio is not representative of the average volumetric distributions of the material. Smaller samples of dimensions 140×50×15 mm were cut out from the beams and notched in their middle as depicted in Fig. 1. The purpose of the notch is to ensure stress concentration and fracture initiation in a concentrated area on which the observation can be focused. The specimen are loaded by uniaxial tension with a 10 kN force cell imposing a constant longitudinal displacement speed of 0.5 mm/min in the direction shown in Fig. 1. A reference picture of the surface is taken with the camera before loading. The tests were then manually interrupted after initiation of a macroscopic fracture. A total of 15 samples was tested and 6 were kept for the moisture uptake measurements. 2 shows the response to tensile loading of two specimen of steel fibre reinforced concrete, one of which was fully fractured whereas the other was interrupted before complete failure. The profiles are normalized over their maximal stress and peak displacement. The material shows a quasi-brittle behaviour until its tensile strength is reached. After an important drop in the stress profile indicating fracture initiation, cohesive forces still apply between both sides of the main crack where its aperture is still low [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF]. In order to dispose of partially fractured specimen, the tests are interrupted during this stage of the loading, and a second picture of the sample surface is taken with the camera after unloading. This picture is to be used for the image correlation algorithm by comparison with the reference image taken before loading, as explained below. The fractured specimen are then kept in a 45 • C oven until their water content stabilises at a nearly entirely dry state. Prior to a water uptake test, the sides of each sample were insulated with impermeable tape in order to avoid 3-dimensional hygric effects such as evaporation. Their lower side was put in contact with a water surface as shown on Fig. 1 and the upper side was left open as to avoid entrapping of the air initially contained by the porous network as the water front rises inside the sample. The water content distribution was monitored during the uptake by 2-dimensional X-ray radiography applied on an observation zone centered on the main fracture originating from the notch.

A finite-element model for moisture transfer is then applied to recreate the moisture distributions measured by X-ray radiography, using the deformations of the specimen measured by image correlation. For the purpose of the modelling procedure, the hygric properties of undamaged concrete have been characterized as well, using dessiccators and pressure plates for the equilibrium moisture content and one-dimensional water uptake monitoring for the permeability of the raw material.

Experimental and numerical procedures

Digital Image Correlation

Several techniques are available for observing fracture processes in materials such as concrete, each of which has advantages and limitations and offers different opportunities for the interpretations of the results [START_REF] Shah | Nondestructive techniques for studying fracture processes in concrete[END_REF]: interferometry, and tomography can reveal existing crack patterns, whereas acoustic emission measurements can be used as a mean to quantify existing or developing damage [START_REF] Ohtsu | Quantitative damage estimation of concrete by acoustic emission[END_REF].

Optical techniques such as Digital Image Correlation (DIC) [START_REF] Sutton | Determination of displacements using an improved digital correlation method[END_REF] are non-destructive and therefore do not disturb eventual further testing on the samples. DIC provides a full map of the deformations at the surface of a specimen and allows following the fracture development without restriction of number or size of the cracks. It also presents the advantage of being easily implemented into most experimental setups, requiring no strict operating conditions or time consuming preparation, nor does it require gauges in contact with the specimen which might interfere with the experiment. It is now commonly used for the determination of constitutive laws [START_REF] Hild | Digital image correlation: from displacement measurement to identification of elastic properties -a review[END_REF] and has also been applied to the study of the fracture behaviour of concrete [START_REF] Choi | Measurement of deformations on concrete subjected to compression using image correlation[END_REF][START_REF] Corr | Digital image correlation analysis of interfacial debonding properties and fracture behavior in concrete[END_REF]. The principle of DIC is briefly summed up below. For a more complete view of its theory and applications, one can refer to a recent book on the subject [START_REF] Sutton | Image correlation for shape, motion and deformation measurements: basic concepts, theory and applications[END_REF].

Principle

The principle of digital image correlation is the conservation of the optical flow between two pictures of a specimen, taken at different stages of its deformation. The planar displacement u (x) at each point x of a surface is defined as:

g (x + u) = f (x)
where f is a reference image, typically taken at the start of the mechanical loading, and g a deformed image. The target is to determine the best possible approximation for u by minimizing the error between u and an approximation v. This is done by searching the minimum value for the functional φ defined as:

φ (v) = g (v + x) -f (x) 2 dx (1)
Assuming a certain smoothness of the investigated displacement field, g is replaced by its first-order Taylor expansion in Eq. 1 and v is decomposed into a linear combination of functions

v (x) = N i (x) v i (2) ∇g (v + x) ⊗ ∇g (v + x) : N j ⊗ N i dx v i = f (x) -g (x) ∇g (x) N i (x) dx (3)
Since f , g and N i (x) are known, Eq. 3 is a linear system of unknown v i , from which the approximate displacement field v can be obtained. Understandably, the choice of the basis functions N i (x) is important for a good approximation of the displacement field, and several options have been presented in the literature for different types of mechanical tests. In the present case, the computation of the displacements is achieved with the Icasoft software 1 , which uses bilinear finite element shape functions. The displacement field is analysed by dividing the image into grid elements, each of which is assigned two displacement components in the decomposition of Eq. 2.

The result of Eq. 3 is a rough estimate of the real solution, since it only computes an integer number of pixels for the displacement fields. Subpixel interpolation algorithms are then implemented for finer results.

The result of the algorithm is an approximation of the displacement field u of each point of the surface. While gradients of u define the strain fields, a sharp difference in the displacement of two adjacent points indicate a crack.

Settings

The samples are observed during tensile loading with a camera focused on the surrounding of the notch as shown in Fig. 1. The setup uses a 50 mm fixed focal macro lens with a 2.8 focal ratio. Because of the low depth of field of the camera in this configuration, the focus is manually fixed as to prevent an automatic adjustment of the lens during the tests. Depending on the distance between the specimen and the camera, the spatial resolution varies from 10 to 20 µm/pixel. The observation zone is extracted from larger images of 3057×4592 pixels, in order to cover for a possible translation of the sample during loading. In order to facilitate the procedure of the correlation algorithm, a speckle pattern was applied on the surface of the specimen with a black paint spray. A grid of 8×8 pixel elements is used for the calculation.

X-ray radiography

The observation of the two-dimensional moisture concentration distributions in porous building materials is now a common procedure in the characterization of their hygric properties. This can be done using gravimetric methods [START_REF] Drchalová | A simple gravimetric method for determining the moisture diffusivity of building materials[END_REF] or techniques such as nuclear magnetic resonance, γ-ray attenuation, capacitance method or X-ray projection [START_REF] Černý | Transport processes in concrete[END_REF]. Several purposes can be served by these methods: extrapolation of the liquid conductivity from transient one-dimensional moisture uptake tests [START_REF] Carmeliet | Determination of the liquid water diffusivity from transient moisture transfer experiments[END_REF], observation of local variations of the flow speed due to fractures [START_REF] Roels | Analysis of moisture flow in porous materials using microfocus x-ray radiography[END_REF] or of the infiltration of chemicals, studies of reinforcement corrosion in concrete [START_REF] Pease | Influence of concrete cracking on ingress and reinforcement corrosion[END_REF]. . . In the recent years, these techniques have been increasingly used in material sciences, and new operating methods were presented for their optimization and for a more efficient interpretation of their results.

Principle

The principle of X-ray radiography is to measure the attenuation of an X-ray photon beam through a specimen. A camera placed on the other side of the specimen detects the amount of transmitted photons, which in the case of building materials allows measuring the concentration of water or other chemicals. If the intensity of the source is denoted I 0 , the intensity I transmitted by a sample of thickness d is expressed by Beer's law:

I = I 0 exp -µd ( 4
)
where µ is the attenuation coefficient of the material. This coefficient indicates the probability of a photon to interact with the material (and to not reach the camera) per unit length. It is a function of the photon energy, the density and effective atomic number of the material. In case of a moisture uptake experiment, the local water concentration inside a specimen can be considered as an additional layer of thickness d w and attenuation coefficient µ w . The measured intensity from a wet sample I wet can then be expressed as separately attenuated by the dry material and this water layer.

I wet = I dry exp -µ w d w = I 0 exp -µd -µ w d w (5)
The moisture content w is related to this fictitious moisture thickness by:

w = ρ w d w d (6)
The water concentration inside the porous medium can therefore be obtained from the intensities measured through the material before and after its contact with the water:

∆w = - ρ w µ w d ln I wet I dry ( 7
)
where ∆w is the increase in water content from the initial state, which is not necessarily entirely dry. The linear attenuation coefficient of water is considered unknown at first and will be determined by coinciding the X-ray intensity transmitted by a water-saturated concrete sample with the measured capillary moisture content. The procedure of the experiment is as follows: the sample is initially placed in a container and a first image is taken to measure the dry intensity I dry . Water is then poured in the container as to be in contact with the bottom surface of the sample. Measurements of intensity transmitted by the wet material I wet follow regularly as the moisture front rises inside the material. According to Eq. 7, the result w shown at one location of the 2-dimensional plane of observation is the average of the water mass per volume over the thickness of the sample. If the material is considered homogeneous in comparison with its dimensions and the size of the observation zone, w will reach a homogeneous distribution at saturation. However, concrete is a highly heterogeneous material, and the specimen of the study have a thickness of only a few times the size of the largest aggregates. The assumption of homogeneity can therefore not be made here. Aggregates have a much lower equilibrium moisture content than cement paste, and the observation of the local water concentration does not depict the actual wetting state of a specimen correctly. This is corrected for by considering the non-dimensional water content s, or saturation degree of the liquid phase: s = ln I wet /I dry ln I sat /I dry [START_REF] Drchalová | A simple gravimetric method for determining the moisture diffusivity of building materials[END_REF] where I sat is the recorded intensity transmitted by a capillary saturated sample. This definition implies that a series of measurements must be carried on until saturation of the specimen in order to normalize the water content between an initial and a final state.

Settings

The moisture profile tests were performed with the GNI X-Ray system at the Department of Civil Engineering of the Technical University of Denmark. The intensities are recorded as 16-bit images showing up to 65535 greyscale levels. The size of each image is 252×256 pixels covering an area of 12×12 mm, i.e. with a pixel size of 50 µm. The X-ray source and the camera can be displaced following a pre-programmed path, which allows the user to visualise larger areas by assembling images taken at several adjacent locations. In this study, the entire observation zone for each series of measurements covers a rectangle area of 3×3 images (or 3×4 images for observation of uncracked samples, as explained below), positioned so that the main fracture of the specimen is slightly away from the center (see Fig. 1). This setup allows including the entire height of the sample between the water front and the notch, and to cover a large enough area in order to capture the evolution of the water front within several hours of measurements.

Since non-negligible fluctuations in the intensity of the X-ray source are expected during the measurements, each image of the assembling is the averaged intensity measured over 40 pictures of the same location, each of which taken with a camera integration time of 1500 ms. The entire acquisition of the observation area, including the time of displacement of the camera and source between the locations, lasts around 10 minutes. Such a time interval between the first and the last images of a series is not problematic here because of the low permeability of concrete, and all images can be considered taken at the same time. The ImageJ software2 is used for the assembling of all images.

The energy of the X-ray source is set as to ensure a good resolution of the water content profiles. According to [START_REF] Roels | Analysis of moisture flow in porous materials using microfocus x-ray radiography[END_REF], the resolution of moisture content R w is defined as one greyscale step of an image and is mainly a function of the sample thickness for each type of material.

R w = - ρ w µ w d ln I dry -1 I dry (9)
The drawback of a high energy is a possible instability of the source which as a consequence would make longer integration times necessary. A fair compromise was found at E = 85 keV, which coincides with the energy suggested by [START_REF] Roels | Analysis of moisture flow in porous materials using microfocus x-ray radiography[END_REF] for ceramic brick samples of similar 15 mm thickness. This results in an average resolution of R w < 0.5 kg•m -3 .

Moisture transfer modelling in concrete

After the observation of moisture concentration profiles within the fractured samples, a first attempt to implement DIC results into a code for moisture transfer simulations was made. The first step into this procedure is the hygric characterization of concrete. This data was then implemented into a finite element model which was programmed to automatically detect fractures and to run 2-dimensional transient moisture transfer simulations on this domain.

Basic concepts

The 2-dimensional conservation equation for isothermal moisture transfer in building materials is written with the capillary pressure as driving potential:

∂w ∂p c ∂p c ∂t = ∇ • K w ∇p c ( 10
)
The right side of the equation is often split in two terms, one describing liquid moisture transfer driven by capillary forces and one describing water vapour diffusion in the gaseous phase. In the present case however, both phases are described by a total permeability of the material K w . The left side of the equation indicates the storage of moisture in the pore structure and is described by the moisture retention curve, i.e. the equilibrium moisture content of the material w versus relative humidity or capillary pressure p c . The full hygric characterization of porous building materials includes measurements of these storage and transport properties, in both hygroscopic and over-hygroscopic ranges. In the hygroscopic range, moisture transfer mostly occurs through diffusion or advection phenomena of water in the gaseous phase. The over-hygroscopic range represents higher concentrations of water in the material, where transfer in the liquid phase prevails. The transport and equilibrium properties of concrete are to be measured during wetting, as opposed to previous studies concerning the drying of several formulations of concretes and cement pastes [START_REF] Baroghel-Bouny | Characterization and identification of equilibrium and transfer moisture properties for ordinary and high-performance cementitious materials[END_REF].

The procedures for experimental characterization of storage and transport properties are explained below: the moisture retention curve was characterized by weighing samples placed in stable humidity conditions. The permeability was measured following the Boltzmann technique [START_REF] Carmeliet | Determination of the liquid water diffusivity from transient moisture transfer experiments[END_REF], which requires the monitoring of moisture concentration distributions in the material during capillary uptake.

Hygric characterization

The complete profile of moisture content w versus capillary pressure p c was measured using two separate methods in order to cover the entire humidity range: dessiccators were used in the hygroscopic range between 11% and 97% relative humidity (p c > 4 × 10 6 Pa) and pressure plates were used in the over-hygroscopic range (p c < 4 × 10 6 Pa). The samples used for the dessiccator method are initially oven-dried, whereas the pressure plate method proceeds by extracting water from initially vacuum-saturated samples. Therefore, the former method provides the sorption curve for the hygroscopic range, while the latter provides the desorption curve for the over-hygroscopic range. This was corrected for by subtracting the recorded difference between capillary and vacuum saturation contents to the desorption (extractors) measurements. This difference being relatively low, the assumption could be made that the result is a full sorption curve.

Aside from the observations of water infiltration in cracks, the X-ray equipment was also used to monitor moisture uptake in uncracked samples of concrete. Transport properties can either be characterized by means of a water absorption coefficient estimated by weighing samples during free water uptake experiments [START_REF] Gummerson | Water movement in porous building materials-ii. hydraulic suction and sorptivity of brick and other masonry materials[END_REF][START_REF] Vejmelková | Free water intake as means of material characterization[END_REF], or by monitoring concentration profiles in the samples. The latter approach allows characterizing the entire profile of unsaturated moisture permeability. In order to estimate the moisture transport properties of the material, let us consider a simplified formulation of Eq. 10 for the case of one-dimensional moisture transfer with w as the driving potential (note: in the present case, moisture transfer occurs in the vertical direction y):

∂w ∂t = ∂ ∂y D w ∂w ∂y ( 11
)
where D w is the moisture diffusivity. As explained in [START_REF] Carmeliet | Determination of the liquid water diffusivity from transient moisture transfer experiments[END_REF], D w can be derived from transient 1D moisture profile measurements through a Boltzmann transform λ = y t -1/2

D w = - 1 2 w 0 λdw ∂w/∂λ ( 12
)
The one-dimensional moisture profiles were monitored in a non-damaged concrete sample by X-ray radiography. Several consecutive measurements of the vertical distribution of the moisture content w were performed during a water uptake test and superimposed as to obtain a single humidity profile. The moisture diffusivity is calculated from the resulting mean moisture concentration profile using Eq. 12. The permeability was then derived following its definition:

K w = D w ∂w ∂p c ( 13
)
where ∂w/∂p c is the slope of the moisture retention curve. The calculation of K w therefore requires the knowledge of both the equilibrium moisture content and the moisture diffusivity. One of the main advantages of the Boltzmann method to derive the transport properties of a material is the fact that the term K w includes moisture transport in both liquid and gaseous phases.

Modelling settings

Simulations were run using a new model currently under development, which aims at including experimental observations of the fracture network geometries into heat and moisture simulations. The present work aims at recreating isothermal moisture flow from the main fracture through the surrounding porous material. Eq. 10 was discretised over a triangular finite-element mesh using the Galerkin weighted-residual method, implicit temporal discretisation and a Newton-Raphson iterative scheme. The complete procedure for numerical implementation of Eq. 10 is explained in [START_REF] Janssen | Conservative modelling of the moisture and heat transfer in building components under atmospheric excitation[END_REF]. As a first use of the model, the fractures in contact with the water front will be considered initially water-filled before starting the simulations. This is consistent with observations and previous simulations [START_REF] Roels | Measuring and simulating moisture uptake in a fractured porous medium[END_REF] of these fractures filling within seconds after contact of the sample with water, and with the fact that concrete has a low permeability in comparison to other building materials such as brick or sandstone. This hypothesis implies that the fracture is considered a hydraulic continuity, as opposed to other existing contact configurations [START_REF] De Freitas | Moisture migration in building walls -analysis of the interface phenomena[END_REF]. The model will be extended in upcoming works as to account for fracture flow and imperfect hydraulic contact between crack sides.

A first non-structured mesh was generated on the domain where the surface displacement fields have been measured by DIC. The correlation results were then interpolated at the node coordinates of the finite element mesh. If two adjacent nodes show an important difference in displacement, their support is considered crossed

Results and discussion

The results of the experimental and numerical procedure are now displayed and discussed. As a first step, the results of the fracture characterization by digital image correlation are displayed as to determine whether this technique can provide a satisfactory observation of the crack patterns within the material. The resulting cartographies of deformations are then superimposed to the moisture profiles obtained by X-ray radiography in order to estimate the reliability of non-destructive fracture characterization as a mean to predict moisture infiltration. The characterized transport and storage moisture properties of concrete are then displayed as a basis for the upcoming transfer modelling. Finally, the simulation results are shown and confronted to moisture ingress measurements by X-ray in the fractured porous material. All data presented in the present section concerns the same sample of concrete, whereas numerical data is provided in appendix A concerning the crack geometry of a second sample.

Observation of fracture patterns

The DIC algorithm works by comparing a reference picture of the surface of a specimen taken before loading, and a deformed picture taken after unloading (see Fig. 3(a)). The result of the calculation is the displacement field in the direction of the loading. It is shown in the form of contour lines on Fig. 3(b), each one representing a difference of 5 µm of the displacement field in the x (horizontal) direction. A local concentration of the contour lines indicates a steep local difference in the displacement, and therefore locates a fracture.

The mechanical loading of the sample was interrupted before its entire failure. The main fracture can be seen on Fig. 3(a) originating from the notch in the middle of the upper boundary, and propagating downwards. This crack presents a tortuous shape and secondary discontinuities can be seen propagating from it, mostly along the cement/aggregate interfaces. The horizontal contour lines on both sides of the main crack do not indicate the presence of fractures, but are caused by the notch on the upper side of the sample: as the crack expands, both sides of the sample tend to slightly rotate around the third axis in opposite directions, causing a gradient of horizontal displacement along the vertical coordinate. With further loading, the main fracture (shown by the higher concentration of contour lines) expands and tends to cause the surrounding cracks to close. It can therefore be ex- pected that the moisture content distributions observed during water uptake tests will be influenced by the state of the sample when mechanical loading is interrupted.

From the results shown on Fig. 3(b), it can be seen that with the present settings, digital image correlation is able to detect cracks of down to 10 µm aperture on a relatively large area of observation. In terms of spatial resolution, other techniques such as electronic speckle pattern interferometry or micro-tomography have been found to detect cracks of under 1 µm opening [START_REF] Shah | Nondestructive techniques for studying fracture processes in concrete[END_REF]. The spatial resolution of DIC is mostly a matter of the physical pixel size, i.e. of the settings of the camera (spatial position, focal length). A compromise must be found between searching for smaller crack apertures and covering a large enough observation area.

The technique also provides the shape and tortuosity of all fractures, as well as the local variations of their aperture. All of these parameters are likely to have an influence on the moisture flow inside the fractures [START_REF] Vandersteen | A network modeling approach to derive unsaturated hydraulic properties of a rough-walled structure[END_REF] and in the surrounding porous material. DIC therefore has the potential to be a significant improvement to the moisture flow prediction in comparison to an arbitrary choice of the fracture geometry.

DIC has two main downsides when applied to the localization of fractures for flow simulations: first, since it is a mapping of the surface of the samples, the discontinuities appearing within the volume of the material do not show on the results. As a consequence, subsequent water infiltration may follow paths which are not captured by the technique. This issue can be corrected by considering samples of relatively low thickness, in which most discontinuities are expected to show on the sample surface. This will be discussed in the next section. The second downside of the technique is the fact that DIC must be applied while deformations occur, since a reference picture of the sample is required for the calculations. This prevents its applicability to in situ characterization of building components. This matter will be discussed in the conclusion of the paper.

Superposition of DIC and X-ray measurements

As explained above, each moisture content distribution provided by the equipment is an assembling of pictures taken at several locations, each of which is the average intensity measured during 60 seconds of integration. During the earliest stage of wetting, water movement near the water front and inside the fracture is fast in comparison with the time needed for measuring. Therefore, in order to be able to consider an assembling of pictures as coming from one single time coordinate, the first moisture content distribution within each specimen was measured around 8 minutes after contact with water. The result is shown on Fig. 4 for the sample shown above on Fig. 3.

The bottom of the image show the water meniscus in contact with the specimen. The x-axis was centered on the position of the notch. The water content is depicted by the color mapping, while the fractures and defor-mations are shown by the contour lines previously depicted. The following observations deduced from the example shown on Fig. 4: first, the moisture concentration profiles on the early stages of wetting follow the shape of the main cracks observed with digital image correlation, which are visible as the highest concentration of contour lines. Fractures under 100 µm aperture also contribute significantly to the flow, as shown by the smaller fractures which are quickly filled with water as well. However, in some cases, water can take additional paths which have not been captured by DIC. As mentioned above, this observation technique only captures deformation processes appearing on the surface of the sample, and not all cracks opening within the volume of the material are visible.

These observations underline the applicability of DIC as a mean to predict preferential moisture paths, but also show its limitations. The next part of the results handles the observation and modelling of the moisture flow outside of the fracture.

Hygric characterization

The sorption isotherm was measured with dessiccators and pressure plates as described above. In order to interpolate the profile to the entire range of capillary pressure, the results are approximated with a van Genuchten function [START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[END_REF]:

w p c = w sat 1 + α|p c | n -m (14)
where the parameters m, n, α are determined by approximation of the measured equilibrium moisture contents. The results of this process are shown on Fig. 5 and table 2. w sat (kg.m -3 ) 0.2191 1.2566 6.165×10 -7 101.44

Table 2: Coefficients of the sorption isotherm approximation for concrete X-ray radiography was then used to observe capillary moisture uptake in an undamaged concrete sample. Fig. 6(a) was obtained by stacking 6 consecutive moisture profile measurements, performed between 100 and 400 minutes after wetting of the sample. It can be seen that the data is quite scattered due to the heterogeneity of the material, but the distributions globally coincide. All profiles are then assembled as one average distribution, which is then smoothed and approached by a cubic spline interpolation.

As stated in [START_REF] Carmeliet | Determination of the liquid water diffusivity from transient moisture transfer experiments[END_REF], a boundary moisture content w B must be assigned at the inflow surface in order to satisfy the Boltzmann boundary conditions. This value can be considered equal to the capillary moisture content of the material, i.e. the maximal moisture concentration w sat reachable by water uptake under atmospheric pressure. Under this assumption, the moisture diffusivity D w was calculated as follows: first, D w was computed using Eq. 12 from one-dimensional concentration profiles based on an arbitrary value of the photon attenuation coefficient of water µ w . The diffusivity profile showed a maximum value where the boundary moisture content w B is supposed to be reached. The attenuation coefficient was then arranged so that w B matches the known value of the capillary saturation content w sat , which results in a value of µ w = 17.12 m -1 .

The resulting diffusivity profile shown on Fig. 6(b) is used to calculate the permeability through Eq. 13, which was then approximated by a multi-modal Durner function [START_REF] Durner | Hydraulic conductivity estimation for soils with heterogeneous pore structure[END_REF]:

K w = K sat s τ   i l i α i 1 -1 -s 1/m i i m i i l i α i   2 (15)
where K sat is the saturated permeability. In this expression, the retention function s is divided into subcurves s i , each one representing a pore size subsystem:

s = i l i s i = i l i 1 + α i |p c | n i -m i (16)
The parameters l i , α i , m i , n i estimated by approximation of expression 15 with the experimental results do not necessarily have the same values as those obtained from the retention curve (see table 2). In the present study, a better fit of the permeability profile was reached with the parameters displayed in table 3, as shown on Fig. 7.

Moisture transfer modelling and monitoring

As stated above, the depiction of the wetting state of the material by the water content w is not always appropriate due to the high heterogeneity of the material. Fig. 4 shows local fluctuations of the water content which are not only due to the presence of cracks, but also to the distribution of aggregates and cement paste. In order to reveal a moisture front and attempt to measure its advance, the non-dimensional moisture content s is now considered (see Eq. 8).

Due to the low permeability of concrete, the water concentration distributions were both measured and simulated over several hours after contact of the sample with the water front. The are based on a grid of 0.7 mm elements, with time steps of 2 minutes. In order to account for the initially nearly dry state of the material, the capillary pressure at t = 0 is set to p c = -5 • 10 8 Pa.

The comparative results of the moisture concentration monitoring and modelling are shown on Fig. 8 as two consecutive measurements of the saturation distributions for the same sample. The displacement field at the surface of the sample is shown as contour lines on the upper left image. The simulations show a fair accuracy in the shape and evolution of the moisture fronts from the bottom water surface and from the fractures. This confirms the relevance of the hypothesis of initially water-filled fractures.

At earlier stages of wetting, it is visible on both simulations and experiments that the shape and tortuosity of the cracks has a noticeable impact on the shape of the moisture fronts (Fig. 8(a) and 8(b)). The crack width shows no noticeable impact on the moisture flow in the neighbouring porous material: because concrete has a low permeability, the flow inside the fracture occurs at a much smaller time scale than through the porous network. Therefore, the moisture saturation of a water-filled crack does not appear affected by the migration of moisture through its walls to the porous network. This effect is likely to not be valid for other building materials with higher permeability such as brick or gypsum boards, making it necessary to include a fracture flow model in the simulations.

Under atmospheric conditions, the eventual contact of a liquid water phase with a side of a building component is caused by wind-driven rain [START_REF] Hall | Water movement in porous building materials-v. absorption and shedding of rain by building surfaces[END_REF][START_REF] Blocken | A review of wind-driven rain research in building science[END_REF]. The water flow at the surface is then no longer comparable to a continuous water front and the fracture is likely to present a variable moisture content. Under these conditions, a fracture flow model is necessary as well, for which the complete fracture geometry provided by DIC can be an important input as it includes the distributions and variable apertures of the cracks. This underlines the importance of the present experimental and numerical methodology, and shows the necessity of extending the current modelling procedure to include flow within the fractures.

Conclusion

Cracked samples of concrete have been submitted to moisture uptake measurements by X-ray radiography, after their fracture patterns were observed by digital image correlation. The target was to evaluate the ability of DIC to provide a reliable characterization of the fracture geometries in order to predict moisture uptake through a damaged porous material.

It was established that DIC is a reliable non-destructive technique which can provide accurate quantitative measurements of the fracture patterns: microcracks of down to 10 µm aperture were observed on a relatively large area. The overall shape, tortuosity and variable aperture of the surface fractures can be accurately captured. The technique is however limited to a 2-dimensional observation of the surface of a thin sample, which prevents its implementation to the case of relatively larger building components, in which fractures might propagate without appearing on the monitored surface. Moisture uptake was monitored as to confirm the ability of computer vision to predict the moisture uptake paths in the earlier stages of wetting. The water infiltration from the fractures to the surrounding porous material was then observed. These moisture content profiles were then approached by finite-element moisture uptake modelling based on the measurements of deformation fields at the surface of the samples. The 2-dimensional infiltration of water from a water-filled fractured could be recreated similarly to the experimental observations. The in situ hygrothermal characterization of damaged building materials can follow the same procedure, by replacing the DIC results by an other similar non-destructive technique for fracture localization. Further work is to be performed as to investigate the ability of such techniques to provide 3-dimensional crack geometries in existing building components. The moisture flow modelling is to be improved as well, as to include the possibility of simulations under climatic conditions of fractured porous media.

A Fracture geometry

In order to provide a "real" 2-dimensional fracture geometry for fracture flow simulations, numerical data is provided here. This data concerns a sample showing a single tortuous macroscopic crack and no microscopic side cracks. A fair approximation of its trajectory can be reached by cubic spline interpolation between the points listed in Table 4. The table also lists the fracture aperture at each of these points.

x (mm) y (mm) a (×10

Figure 1 :

 1 Figure 1: Setup of the samples during mechanical loading (above) and capillary uptake experiments (below). The arrows respectively show the directions of the tensile loading and moisture uptake, the dotted lines respectively show the DIC and X-ray observation zones

Figure 2 :

 2 Figure 2: Complete (dotted line) and interrupted (continuous line) tensile loading profile of steel fibre reinforced concrete

Figure 3 :

 3 Figure 3: 3(a) Picture of a specimen surface after mechanical loading, 3(b) surface displacement mapping of a concrete sample from Digital Image Correlation (5 µm/line)

Figure 4 :

 4 Figure 4: Water content distribution 10 minutes after wetting (kg•m -3 ) and sample deformation contour lines (5×10 -6 m/line)

Figure 5 :

 5 Figure 5: Moisture retention curve measurements versus capillary pressure (circles) and approximation by van Genuchten model

Figure 6 :

 6 Figure 6: (a) Averaged w/λ distribution and standard deviation over 6 consecutive measurements, and spline interpolation (continuous line), (b) Moisture diffusivity extrapolated from transient measurements

Figure 7 :

 7 Figure 7: Moisture permeability from measurements (dots) and Durner model (line) versus capillary pressure

Figure 8 :

 8 Figure 8: Measured (left) and simulated (right) distributions of the saturation degree s after 10 minutes (above) and 40 minutes (below) of wetting. The first image shows the position of the cracks captured by DIC.

Table 1 :

 1 Concrete proportions

		5
	Cement	3.465 kg
	Sand 0-4 mm	8.043 kg
	Aggregates 4-8 mm	11.580 kg
	Steel fibres	0.2 kg

Table 3 :

 3 Coefficients of the multi-modal approximation of the moisture permeability profile

	-1 )	l	K sat (s)	τ
	0.6148 2.5963 5.5383×10 -7 0.5062 2.2182×10 -13 -4.6975 0.1913 1.2366 2.2493×10 -8 0.4938

Table 4 :

 4 Position and aperture of the main crack in the second reference sample

			-4 m)
	-3.3042	0	1.4889
	-3.6836	1.2236	2.7256
	-3.4596	2.4471	2.8292
	-3.5335	3.6707	2.7422
	-3.2970	4.8943	2.9217
	-3.1608	6.1179	3.0143
	-3.0388	7.3414	2.8012
	-2.1352	8.5650	2.7121
	-2.1257	9.7886	2.6523
	-1.0648 11.0122	2.4571
	-0.9687 12.2357	1.7370
	0.6810	13.4593	3.0582
	1.2109	14.6829	2.8954
	0.8548	15.9065	2.8565
	1.6427	17.1300	2.7664
	2.0858	18.3536	2.9611
	2.2049	19.5772	2.8261
	1.9665	20.8008	2.8060
	1.2050	22.0243	2.8996
	0.9245	23.2479	3.0533
	0.9347	24.4715	2.9321
	1.4069	25.6951	2.9470
	0.9609	26.9186	2.8931
	0.0112	28.1422	2.9744
	-0.1474 29.3658	3.0566
	-0.2061 30.5894	3.0740
	0.1037	31.8129	1.9138
	-0.0044 33.0365	2.9545
	1.1235	34.2601	3.4714
	0.5347	35.4837	2.5649

Icasoft Digital Image Correlation Software -http://icasoft.insa-lyon.fr/

ImageJ image processing and analysis in Java, www.rsbweb.nih.gov/ij/

Acknowledgements

The results presented in this paper have been obtained during the Renovbat project for the study of flow in fractured building components. Funding was provided by the Energy research cluster of the Rhône-Alpes region, which the authors would like to thank for their financial support through the project.

The authors would also like to thank K.K. Hansen and the researchers of the Construction Materials section of the Civil Engineering Department at the Technical University of Denmark for their help in the experimental processes.