Vagelis G Papadakis

Costas G Vayenas

Michael Fardis

Physical and Chemical Characteristics Affecting the Durability of Concrete

Keywords:

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.

The durability of reinforced concrete is influenced by those physical characteristics of concrete that control the di f of gases, such as CO, and 0,, or of liquids (mainly water) through its pores, and the diffusion of io. ns, such as c1-, dissolved in the pore water. These physical characteristics depend on the composition of concrete, the chemical composition and type of cement, and the re/atiVe humidity and temperature of the environment. In the present paper, these characteristics of concrete are determined analytically and/or experi mentally in terms of the composition parameters and environmental conditions: the molar concentration of those constituents that are susceptible to carbonation; the porosity o.. nd pore-size distribution; the degree of saturation of the pores; and the effective di f of gases through the concrete.

• Keywords: concrete durability; diffusion; diffusivity; hydration; porosi ty; pozzolans; reinforced concret � s; saturation; specifi c surface.

•

The good performance of concrete in service, includ ing its durability, is a major factor for its success as a / construction material. This is particularly true of i t s application, in the form o f reinforced concrete; as a structural material, since its major competito ; s, steel and timber, have inherent durability limitations. How ever, the last two decades have seen a ilisconcerting in crease in examples of the unsatisfactory durability of concrete structures, especially reinforced concrete ones.

Since the state of the art of concrete and reinforced concrete durability has enKToyed a relatively short re search history compared .iO strength, durability has re placed strength as the µumber one issue concerning the engineering community today. Deterioration of concrete in service may be the result of a variety of physical and chemicr.! processes, 1•2 such as attack by acids, sulfates, or alkalis, alkali-aggregate reactions, freeze-thaw cycles, e.�c. In reinforced con crete, the inost serious deterioration mechanisms are those leading to corrosion of the reinforcement, result ing in a reduction in the effective cross-sectional area of reinforcing bars, and in ultimate disruption of the con crete due to spalling of the concrete cover. Reinforce ment corrosion occurs only after depassivation due to carbonation of the surrounding concrete, penetration of chloride ions, or a combination of both. 1 • 2 All proc esses causing deterioration of the concrete itself or car� rosion of the reinforcement involve transport phenom ena through the pores of the concrete. These transport mechanisms inclu4e diffusion of such gases as 0 2 , C0 2 , or S02 in the gaseous phase of the pores, from the en vironment to interior regions where the gas concentra tion is low due to gas consumption (by carbonation for C02 and corrosion of the steel for 0 2 ). They also in clude diff u sion of aggressive ions dissolved in the pore water,, / again from the external surfaces to interior regions of consumption and low concentration. These dis s olved aggressive ions can be c1-, H ... , or SO�-; o r iginating from dissociation of dissolved acids, alkali ions, Na ... , etc. Flow of water through the pores is also possible due to capillary suction or hydraulic pressure gradients. Finally, it is in the pore water that most of the chemical reactions between the externally supplied aggressive substances and the constituents of cement and concrete take place. Therefore, to study, describe quantitatively, and ultimately limit or control these transport phenomena and the deterioration of con crete, knowledge of the structure of its pores and dis tribution of the pore volume into gaseous and aqueous phases is absolutely necessary. In addition, one needs to know the evolution in time of the concentrations of those constituents of cement, hydrated or not, that are expected to react with the externally supplied aggres sive substances.

Significant work on hydration reactions and mecha nisms has already been done by • previous workers. The literature is also rich in concrete porosi• y and poresize distribution measurements, s-20 studies of water be havior in hardened cement,

RESEARCH SIGNIFICANCE

Quantitative prediction and control of the evolution in time of physicochernical processes preceding corro sion of reinforcing bars, such as the reduction in alka linity due to carbonation of concrete or the diffusion of c1in the pore water, and processes that take place si multaneously with corrosion, such as diffusion of 02 through the pores, require knowledge of some physical characteristics of concrete, such as porosity, pore-size distribution, degree of saturation of the pores, and gas diffusivity, as well as the amount of carbonatable ma terials in concrete, etc. Computed and measured values of these characteristics or parameters are presented in this paper for portland cement "in terms of the compo sition parameters of concrete, chemical composition of cement, and environmental parameters, such as ambi ent relative humidity. In a separate paper, the effect of pozzolans on these characteristics or parameters will be discussed.

MOLAR CONCENTRATION OF CARBONATABLE CONSTITUENTS OF HARDENED CEMENT

PASTE

The major constituents of hardened cement paste • subject to carbonation in i.he presence of moisture are Ca(OH)2 and calcit!'ll silicate hydrate (CSH). Prior to their hydration, calcium silicates C3S and C2S are also subject to carbonation. For the other hydrated or un hydrated constituents of hardened cement paste, car bonation seems to be limited to a surface zone, without affecting the bulk of the crystaHites. s Therefore, it is the molar concentrations of Ca(OH)i, CSH, and the yet unhydrated amount� of C3S and C2S that are of interest as far as carbonation is concerned. These constituents of cement are produced or consumed by the chemical reactions of hydration. Using the notation of cement technology,* the hydration reactions• of ordinary port land cement (OPC) in the presence of gypsum are 3

•1 f H,C3S 2C3S + 6H -C3S2H3 + 3CH (l) rH,C2S 2CiS + 4H -C3S2H3 + CH (2) C4AF + 2CH + 2CSH2 (3)
rH,C4AF

+ I8H -C8AFS2H24 rH,C3A C3A + CSH2 + lOH -C4ASH 1 2 (4)
rH,C4AF

C4AF + 4CH + 22H -C8AFH26 (5)
rH.C3A

C3A + CH + 12H -C4AH13 (6)
The previous hydration reactions take place at molar rates per unit concrete volume rH, i (i = C;S, C2S, C4AF, C3A) in mol/m 3

• sec. The presence of gypsum has only a minor effect on the rates of hydration of c;A and C4AF, so that Reactions (3) and (4) can be taken to have the same rate as Reactions [START_REF] Bensted | Hydration of Portland Cement[END_REF] and (6), respec tively ,6 i.e., rtt,c¢ and rH , C J A • Expressions for the rates rH,; can be obtained from measurements of the fraction F,{t) of compound i, which bas been hydrated at time t (in sec) after mixing. Such measurements have been presented in the past by Brunauer and Copeland 3 and more recently by Taylor ,7 and point to expression of the form

fH,I := -d(i )/df = kH, i [i]"i/[i] 81-I (7) F/J ) = 1-[1 1/(tlo = 1-(l -kH,I t(l -n1))11C1-n1) (8)
in which [i] and [i]0 are the current and initial (at t = 0) molar concentrations of compound

•

In the presence of gypsum, Reactions (3) and (4) dominate over Reactions [START_REF] Bensted | Hydration of Portland Cement[END_REF] and (E),s.r. so Reactions [START_REF] Bensted | Hydration of Portland Cement[END_REF] and (6) take place only after all the gypsum has been consumed. This happens at a time t*, which can be de termined from the condition The value oft* can be found numerically, from an ex ponential in t* equation, which results from substitu tion of Fc3A (t*) and Fc.iAF (t*) from Eq. ( 8) in terms of t*. The second term in Eq. (9a) can be considered neg� ligible, because SO�ions react preferentially with C3A rather than with C4A F ,6 so the value of i* can be com puted in approximation (9b)

The molar concentrations of the c::arbonatable con stituents at time t after mixing are, therefore, equal to

3 1 [Ca(OH)2J = 2 [C3SJ0Fc 3s + 2 [C2Slo Fc2s -2[C4A F]0Fc4AF• 0 � t � t* (lOa) � 2. 0 0. 4 n E OPC WIC=O. 65 a/C=3 ' Ca(OH)2 � 1. 5 >- o. 3 I- H V1 n 0 b a::: C-+ 0 a.
)( 1. 0 0.2 [

Ca(OH)2] = 2 [C3SJ0Fc3s + 2 [C2SJ0Fc2s • -4[C4AF]0Fc 4A F -[C3A]o F�A + [CSH2]o, t ;;;i: t* (10b) (12) (13)
The Fig. 2 shows the evolution in time of the molar con centrations of the carbonatable constituents of con-crete, according to Eq. (10) through [START_REF] Hansen | Phys ical Structure of Hardened Cement Paste. A Classical Approach[END_REF] for ordinary portland cement, and for w/c = 0.65 and ale = 3. Fig. 3 presents the final concentrations of Ca(OH) 2 and CSH for various mix proportions of concrete and for fully hydrat ed OPC.

POROSITY OF CONCRETE

The porosity of concrete E, defined as the ratio of pore volume to the total volume of concrete, decreases with time due to the processes of cement hydration and of carbonation [START_REF] Gunter | Effect of Curing and Type of Cement on. the Resistance of Concrete to Freezing in Deicing Salt Solutions[END_REF] in which Eo is the porosity of fresh concrete and AEH(t), Aec correspond to the reduction in porosity due to hy dration and carbonation, respectively. The initial value of porosity Eo is the sum of the volume fractions of mixing water• and entrapped (or entrained plus en trapped, in air-entrained concrete) air. The volume fraction of entrapped (or entrained) air, denoted by Ea;,, is in the range of a few percent, depending on the max imum aggregate size,1.6 whereas the volume fraction of water in fresh concrete equals the mass of water per unit volume of concrete divided by the mass density of wat er. In terms of the composition parameters of con crete, i.e., of the water-cement ratio wlc. and the ag gregate-cement ratio ale, the initial porosity equals

w Pc ) --(1 -fai r C Pw Eo = --------+ fai r ( W Pc a Pc) I+--+ - c Pw C Pa (17)
The reduction in porosity due to hydration Ae,).,t) is due to the fact that the molar volume of the solid products of hydration, in the right-hand side of Eq. (1) through ( 6), exceeds that of the solid reactants, in the left-hand side of these equations. So, if A Vc3s. A Vc 2s. A Vc4AF.s• A i7c3A,s• A Vc4AF• and A Yc3A are the differ ences in molar volumes between solid reaction products and solid reactants in the hydration reactions, Eq. (1) through [START_REF] Fri Gione | Gypsum in Cement[END_REF], respectively (giiven in Table I, in m1 per re acting mole of compound i), then ---------- 0:: 0. 5 1--r--- Fig. 4 presents the porosity e p (t) of the hardened ce ment paste, which is related to the total porosity E(t), as follows

Ae ,).,t) = [C3S]0Fc3sA Vc3s + [CzS]o Fc2sAVc2s + [C3AJ0Fc3A.6 Vi:3A,s + [CAF]0.Fc4AF AVc4AF,S• O�t��* (18a) .6e,).,t) = [C 3 SJ0Fc3sA Vc3s + [C2S]0F ciSAVc2s + [C3A]0Fc3A (t*)Li Vc3A.s + [C3A]0(Fc3A -Fc 3 A(t*))A V c 3 A + [C,AF]0Fc4AF(t*)AV c4AF.s + [C4AF]0(Fc4AF -Fc4AF(t*))AVc4Af>t � t* (18b)
___________ 7 _ _ _ _ z w u z OPC --Co (OH) 2 --CSH 8
W Pc 1 + -- ep(t) = e(t) (l + � ) c Pw ( 20
)
The dependence of ey(t) on w/c is shown separately for noncarbonated and carbonated concrete (or paste), the former for t = 1, 10, 100, and 1000 days, and the latter for t = 100 and 1000 days only. All the quanti ties in Fig. 4 (i.e., e0, Ae,At), Aec, and .ep) refer to unit volume of cement paste (plus entrapped or entrained air). Measured values of porosity are also shown in Fig. The pore-size distribution {PSD) of concrete or hardened cement paste affect" its effective diffusivity and controls the degree of saturation of the pores (for given relative humidity) and their specific surface area.

Two methods mainly used today in research and by in dustry for the determination of the pore-size distribu-. tion of porous solids are: mercury porosimetry and nitrogen desorption according to the BET theory. Al though some cement chemists1 1 • 14 maintain that one of these methods is applicable to hardened cement pa ste preted as a strong indication that the pores of hardened cement paste with size below IO to 15 A are not acces sible, not only to mercury but to nitrogen as well. This is not very surprising, since a similar behavior has been reported in the literature for coal, where pores below IO to 15 A have been found to be accessible to C02 but not to nitrogen.

The PSD of concrete has been determined herein by a combination of nitrogen desorption for small pore diameters with mercury porosimetry for large ones. The measurements from these two techniques were com bined so that the same PSD was obtained on the aver age over their common range of application, i.e., be tween -60 and -500 A.

Some measured pore-size distributions are presented! herein in the form of the probability density function f d of the volume of pores with diameter equal to d (in µm), which is defined so that the fraction of the total porosity of concrete, due to pores with the logarithm of d ranging from logd to logd + dOogd), is equal to fd (logd) • d(logd). Pore-size distributions shown in Fig. 5 refer to essentially fully hydrated cement. The chemical composition and physical properties of cr;ment are pre sented in Table 2. The aggregate used was a graded sand of Type A gradation (given in Tabk 2). Fig. 5(a) describes the effect of the water-cement ratio on the PSD of noncarbonated OPC mortar, with an aggre gate-cement ratio equal to 3, whereas Fig. 5(b) shows the effect of carbonation on the PSD of hardened OPC paste with a water-cement ratio equal to 0.50. As is well known from previous studies, 11 •14•1 8 the effect of reduc-ing the wlc is to shift the PSD to the left. As shown in Fig. 5(a), the shifting refers only to the macropores. On the contrary, the effect of carbonation is to shift the PSD to the left, over the entire range of pore diame ters. Fig. 5{b) shows also the effect of the presence of aggregates .

The carbonation reactions (and other chemical reac tions related to durability) take place at the interface of the pores and the solid volume of concrete, at a rate which is proportional to the specific surface area of the pores a,, i.e., the total pore surface area per unit vol ume of concrete. The value of a, can be computed from the PSD and the porosity of concrete, or alternatively from the total pore surface area per unit weight of con crete S (measured by water vapor or by the nitrogen adsorption method), the density of the solid phase ip,, and the porosity of concrete, as follows "' 4

a , = i: J d f ,fl O og d) (21a) a, = SKS,(I -KR) {21b)

DEGREE OF SATURATION OF PORES

The degree of saturation of the pores with water f is the fraction of the pore volume filled with water.

Knowledge off is important for durability, because the fraction of the pore volume available for diffusion of gases, such as C02 for carbonation and 0 2 for corro sion, equals 1f, whereas it is in the remaining frac tion of the pore volume, i.e., inf, that the diffusion of ions that dissolve the passivating layer of reinforcing bars, such as c1-, as well as the carbonation reaction of Ca(OH)2 and• C02 (with both of them dissolved in wa ter), take place.

The degree of saturation can be computed if the PSD and the solid surface chemistry are known, assuming hygrothermal equilibrium between the pores and the environment. Under this condltion, all pores with di ameter less than the Kelvin diameter dK are completely filled with water, while the walls of all pores with di ameters d exceeding dK will be covered by a continuous film of water, of thickness w.

f K = J f d d(logd) . . ( 24
)
whereas the aqueous film of the pores with diameter d > dK occupies the following fraction of the total vol ume of the pores

KL __ 2 Cxd,. r 1 "' [l + (C-l)x](l -x) J d to¢x ( l -K|12dw) J "'1 (logd)
The degree of saturation is the sum of fx and/...,

(26) The dependence off on the water-cement ratio is shown in Fig. 6(a). For the same value of RH, f de creases as the water-cement ratio increases, due to the reduction of the proportion of small-diameter pores in the PSD. Pig 5(1J) and (c) show that the aggregate-ce ment ratio and the aggregate gradation have a negligi ble effect on f (as shown in Table 2, aggregate A was a graded sand, and Aggregate B a mixture of sand and gravel). r;.' his is due to the fact that the presence and gradation of aggregates affect very little the proportion and the sizes of small-diameter pores. No theoretically obtained values off are shown in Fig. 6(b) for aggre gate-cement ratios greater than 3 .0 because, for such ratios, reliable measurements of the PSD by mercury porosimetry and by the BET technique were diff icult to obtain. Io the same figure, the contributions of fK and /,.,to the total value off are separated. Finally, Fig. 6(d) shows that, although carbonation increases the relative proportion of small-diar.1etP.r pores, it causes a reduc tion in the value off. This reduction is due to the chemistry of the poie surface, which corresponds to lower values of A and C in Eq. ( 23) and [START_REF] Hurling | Oxygen Permeability of Concrete[END_REF] and, hence, to less water in the pores of concrete after car bonation.

I (c) � o . . o� -- • ---� ---�---...._ __.

EFFECTIVE GAS DIFFUSIVITY OF CONCRETE

The molecules of a gas A, such as C02, 02 or water vapor, move within a medium such as air from regions where the concentration (and hence the partial pres sure) of gas A is high to those where it js low. This transport mechanism of gases, called diffusion, follows Eq. [START_REF] Papadakis | Reac tion Engineering Approach to the Problem of Concrete Carbona tion[END_REF].

From the point of view of durability of reinforced concrete, we are interested in the diffusion of C02 and 02 through the gas phase of the pores, as it is the for mer that carbonates the Ca(OH)2 of the hardened ce ment paste, therefore reducing the alkalinity of con crete, anct ii. is the latter that oxidizes the already de passivated (because of the drop in alkalinity or chloride penetration or both) reinforcing bars. The effective diffusivity of a gas A, De. A is inversely proportional to the square root of its molar weight MWA since, in con crete Knudsen diffusion controls,30 as the mean free path t.. of gases is of the same order as the pore diame ter. So, the effective diffusivity of a gas B, D.,8 can be obtained from that of A by multiplying it by (MWA/ MW8)v.. This concept can be applied to find the eff ec tive diffusivity of C02, De,co2 from that of N2, which is Similarly, the effective diffusivity of 0 2 can be com puted from that of C0 2 by multiplying it by 1.17 and cement paste (open circles) specimens, as long as the water-cement ratio and the type of cement are the same. In other words, the presence of aggregates does not seem to affect significantly the value of effective diffusivity, apparently because diffusion along the ag gregate-paste interface makes up for the lack of diffu sion through the aggregate particles themselves. There fore, the value of effective diffusivity is considered herein as a function of the porosity of hardened ce ment paste ep(t), given by Eq. ( 20). Carbonation is found to reduce diffusivity, as it reduces the value of e P '

The significant effect of relative hurnidty can be taken into account by considereing the value of D, as a func tion off, which also includes the effect of pore-size distribution. Alternatively, the effective diffusivity can be considered directly as a function of (1 -RH/100).

Regression analysis of the experimental data yielded the following expression for the effective diffusivity of C0 2 (in m2/sec)

D 1.64 • 10 -6 Ep 1 . 8 (1 -RH/ 100) 2 • 2 •. co2 =- (29a)
The predictions of Eq. (29a) are also presented in Fig. 9, and are in good agreement with the corresponding experimental data. The effective diffusivity of 02 is given by D •. 02 = l.92• 10 -6 e�•8 (l -RH/ 100)2• 2 (29b) It is worth mentioning that the latest draft of the CEB 1990 Model Code4 0 states in Chapter 2.1.9 (''Transport of Liquids and Gases in Hardened Con crete") that the effective diffusivity of C02 and 02 in concrete is in the range 0.5 • lQ-8 to 6 • rn-8 m 2 /sec, and proposes the following expression for concrete "sheltered from rain" (meaning relative humidity around 60 percent) D, ,co i (m 2 /sec) = 10 -1 . 10 -o. 0 2 s fck [START_REF] Smith | Chemical Engineering Kinetics[END_REF] in which !c k is the characteristic compressive strength of the concrete in MPa. For concrete exposed to rain (meaning relative humidity around 70 percent, it is proposed to reduce the values given by Eq. ( 30) by 50 percent. Eq. (29a) is in good overall agreement with these code proposals, but gives a more complete quan titative description of the dependence of the effective diffusivity on the controlling variables.

CONCLUSIONS AND SUMMARY

In this paper, quantitative information is given on some physical and chemical characteristics of concrete on which its durability depends. This information can be used as input data to quantitative models of the processes that lead to deterioration of concrete and re inforced concrete, such as carbonation. able constituents of hardened and noncarbonated ce ment paste, Ca(OH)2, CSH, C3S 1 and C2S, are given by Eq. ( 10) through [START_REF] Hansen | Phys ical Structure of Hardened Cement Paste. A Classical Approach[END_REF] in terms of their initial concen trations given by Eq. ( 14) and ( 16) and the evolution of hydration described by Eq. ( 8) and Table 1. The evo lution of the porosity of concrete with time, as affected by the evolution of hydration and the occurrence of carbonation, is also computed analytically [Eq. ( 16) through [START_REF] Day | Measurement of Po rosity in Blended Cement Pastes[END_REF]]. Pore-size distributions as measured by a combination of two techniques, so that the entire spec trum of the distribution can be obtained, are presented as functions of the composition of concrete, etc. (Fig. 5). From the computed porosity and the measured pore-size distributions, the specific surface area of the pores can be obtained from Eq. ( 21). The maximum diameter of concrete pores that are completely filled with water, and the thickness of the aqueous film cov ering the walls of the rest , are given by Eq. ( 22) and [START_REF] Lawrence | Transport of Oxygen through Concrete[END_REF], respectively, in terms. of the ambient relative hu midity. The volume of water in these two groups of pores are computed by Eq. ( 24) and ( 25), giving the de gree of saturntion of the pores [Eq. ( 26)). Finally, the effective diffusivity of gases in concrete, such as C02 or 02, is measured for various compositions of concrete and values of the ambient relative humidity. Empirical expressions, Eq. ( 29a) and (29b), fitted to the test re sults, give the effective diffusivity in terms of the po rosity of the paste and the relative humidity of the am bient air.

  cal Characteristics Affecting the Durability •of Concrete by Vagelis G. Papadakis, Costas G. Vayenas, and Michael N. Fardis

Fig. 1 -

 1 Fig. 1-Speed of hydration of the four major constitu ents of OPC, fill ed points by Brunauer and Copeland1; open points by Tay/or7

Fig. 2 -

 2 Fig. 2-Evolution with age of the molar concentration of the carbonatab/e constituents and porosity of OPC•

  initial molar concentrations of compound i (i = C 3 S, C2S, C4AF, C3A) and of gypsum, denoted by (110 and [CSH2](), respectively, in OPC concrete can be computed• from: a) the weight fraction of clinker and gypsum in OPC, denoted by mc1 and mgy = I -mc1, re spectively; b) the weight fraction of compound i in the clinker m; which can be calculated from oxide analy sis;4 c) the mix proportions of concrete, i.e., the water cement ratio w/c, the aggregate-cement ratio ale, and the volume fraction of concrete in entrapped or en trained air fair m;m c1 Pc (1fa;,) [ 1lo = -------- MW; (1 +�Pc+ �Pc� c p,..

  in which p0 p,.., and Pa (in kg/m 3 ) are the densities of cement, water, and aggregates, respectively, and MW; , MW gy are the molar weights of compound. i and CSH2, respectively, •given in Table I in kg/mo I.

  .......... 2. 0 .----.-.--.--r--.--r--.--r---.---, ,.., o r--------=-=-=-�;�•=-==---�

4 0Fig. 3 -

 43 Fig.3-Effect of water-cement and aggregate-cement ratio on the molar concentration of carbonatable ma terials for completely hydrated OPC

4Fig. 4 -Fig. 5 -

 45 Fig. 4-Dependence of cement paste porosity on age and water-cement ratio and comparison with measure ments at t = JOO days

  8• 2 8• 29 The values of w and dx can be determined if the absolute temperature T (K) and the relative humidity RH (percent) are known (both in µm) Cxdw d w = -----"-----(Ix 12t1w ) [l + (Cl)v] (1 •• Y) A dk = 2w + T In (llx) in which x = RH/100, the constant A is equal to 0.6323 for noncarbonated material and to 0.2968 for carbonated, C is the BET-constant for the particular temperature and water vapor-solid surface system, and d,. is the molecular diameter of water (equal to 3 • 10-4 µm). C-values for water vapor adsorption on some solid surfaces are given elsewhere. 8• 2 1• 29 From the pres ent work it was found that for noncarbonated cements C is approximately equal to 100 and for carbonated ce ments C equals 1. This large difference in the values of C and A is due to the fact that the pore surface of non carbonated concrete is hydrophilic, due to the presence of Ca(OH)2, whereas that of carbonated concrete is not, due to its absence. The completely filled pores, i.e., those with diameter d < dK, take up the following fraction of the total pore volume loga'K

Fig. 6

 6 Fig. 6 shows experimentally observed values off as a function of the ambient RH for several values of wlc, ale, and fully carbonated or not cement by using the static method.21 According to this method, the degree of saturation is computed from the steady-state weight of a concrete sample of known total volume, main tained in an environment of controlled relative humid ity. Open symbols in Fig. 6 ref er to data obtained by desorption (i.e., for reduction of relative humidity from 100 percent to 0), whereas filled ones refer to data ob tained by adsorption of water. For the same value of RH, the value off obtained during desorption exceeds that obtained by adsorption. This hysteresis phenome non is typical of drying processes in porous solids, and can be explained either by the extremely slow rate at which ultimate hygrothermal equilibrium is established (which implies that the actual state of equilibrium is in between those obtained by adsorption or desorption), or by the ink-bottle shape of the pores, resulting • in wa ter being trapped in the chambers along the pore dur ing desorption or blocked from entering C1e chambers during adsorption.

Fig. 6

 6 Fig. 6 also shows the results of the analytical rela tions Eq. (22) through (26), using in the latter the ex perimentally measured pore-size distributions. These "theoretical" values off lie close to the adsorption data for low relative humidities and the desorption ones for high, and are in between for intermediate relative hu midities. This type of behavior is consistent with an ink-bottle shape of the pores.

  0 ._ -_ -_-.._. _____ ..__ _ ....._ _ ..__ _ ....._ ___.

Fig. 6 -

 6 Fig. 6-Effect of relative humidity on degree of saturation of the pores: (a) effect of w!c; (b) effect of ale and ag gregates gradation; (c) effect of aggregates' presence; (d) effect of carbonation (C = carbonated, N = noncarbon ated)

Fick

  in which NA is the rate of transfer of gas A (in moles of A per sec, per m 2 of the medium normal to the direc tion of the flux), d[A]/dx is the concentration gradient of A (in moles of A per m, per m3 of the medium), and DA is the diffusion coefficient or diffusivity of A within the medium in question (in m 2 /sec). In porous media, like concrete, diffusion of gases essentially takes place only within the gaseous phase of the pores. Neverthe less, . it is more convenient to use in Eq. (27) "effec tive" values of NA and DA, which refer to the total cross-sectional area of the porous medium normal to the direction of the flux and not just to that of the gas eous phase of the pores. The subscript e is used to de note these effective values, so the effective diffusivity of gas A is denoted by D •. A' Diffusivity should not be confused with permeabil ity, which is the coeficient KA in Darcy's law consid ered to govern the convective flow of a fluid (gas or liquid) A through a porous medium, due to a gradient in total pressure dP!dx and not in partial pressure, the latter being proportional to d[A]ldx in Eq. (27). With the exception of concrete structures containing liquids or gases under pressure, total pressure gradients are very uncommon in concrete structures, so, the flow of gases or water through the pores of concrete is not the most important transport mechanism for the sub stances affecting durability (capillary flow due to cap illary suction, which may serve as a transport mecha nism for ions �r others substances, is not considered herein, as it is a transient phenomenon, contributing to the establishment of the hygrothermal equilibrium in the pores assumed in this work). The dominant trans port mechanism is diffusion of molecules or ions in the gaseous or aqueous phase of the pores, according to

Fig. 7 -

 7 Fig. 7-Schematic diagram of the effective diffusivity measurement apparatus (steady-state operation; dashed lines: dynamic operation) easier to measure exp . erimentally, as N2 is inert with re spect to concrete, by multiplying the latter by 0.80.

  .

8 .Fig. 8 -Fig. 9 -

 889 Fig. 8-Dijfusion cell

  2 7* Specifically, the evolution in time of the concentration of carbonat-•Papadakis, V. G.; Vaycnas, C. G.; and Fardis, M. N., "Fundamental Modeling and Experimental Investigation of Concrete Carbonation," submit ted to the AC/ Materials Journal, 1989.

  2

  are listed in Table 1, and the resulting F,(t)-curves are compared in Fig. I to measured values. Table I and Fig. l can be considered to apply to OPC of normal fineness, corresponding to Type I cement.

	i, respectively (in 3 ). Fitted values of the exponents n; and the co efficients kH,; mol/m

Table 1 -

 1 Parameters of major constituents of ordinary portland cement

	Compound Exponent n,, Eq. (8) Coefficient kH,i (20 C) x 10' (s-•), Eq. (8) Molar weight MW, x 10', kg/mo! • Molar volumes differences c,s c,s C,AF 2.65 3.10 3.81 1.17 0.16 l.00 228.30 172.22 485.96 270.18 172.17 C,A CSH, 2.41 -2.46 -
	Hydration reaction	(1)	(2)	(3)	(4)	(5)	(6)
	AV x 106, m'/mol	53.28	39.35 -220 155.86 -230 149.82

Table 2 -

 2 Composition and physical properties of portland cement-Gradation and physical properties of aggregates Chemical composition of clinker; percent by weight for the entire spectrum of pore sizes, it is generally ac cepted today that for hardened cement paste, mercury porosimetry is not reliable for diameters below 60 to 100 A, whereas nitrogen desorption is ap plicable for pore diameters only up to 300 to 500 A. Moreover, the

	Cao	SiO,	Al,O,	Fe,O,	MgO	SO,	Residue
	65.28	23.55	6.12	2.51	1.30	0.47	0.77
		OPC: 94.9 percent clinker, 5.1 percent gypsum
	Specific gravity: 3.16; fineness, according to Blaine test, 300 m'/k
			Aggregate gradation, percent by weight
	Particle size I-Va in. Va-J,i in. � in.-4 4-8 8-16 16-30 30-50	50-
	Aggregate A	-	-	-	-20.5 40.4 19.0	20.1
	Aggregate B	28.5	7. 5	18.0	13.0 10.5	7.5	6.0	9.0
				Specific gravity: 2.61

total surface area and the total pore volume of hard ened cement paste, measured by either mercury porosi metry or nitrogen desorption, are lower than the ones measured gravimetrically on dried samples for pore volume or by X-ray scattering for surface area. 1 0• 19 This result, which has been verified by the porosimetry measurements of the present study, has been inter

ACKNOWLEDGMENTS

The Hellenic Geoeral Secretariat for Research aod Technology, and lhe Institute of Chemical Engineering and High Temperature Chem ical Processes have provided financial support for this work. Thanks are also due to Dr. E. Galanoulis of the TITAN Inc. cement com pany for helpful wscussions and technical support provided.


Andrea Saltelli

Making best use of model evaluations to compute sensitivity indices

Keywords: Sensitivity analysis, Sensitivity measures, Sensitivity indices, Importance measures

This paper deals with computations of sensitivity indices in sensitivity analysis. Given a mathematical or computational model y = f(x 1 ,x 2 ,...,x k ), where the input factors x i 's are uncorrelated with one another, one can see y as the realization of a stochastic process obtained by sampling each of the x i from its marginal distribution. The sensitivity indices are related to the decomposition of the variance of y into terms either due to each x i taken singularly (first order indices), as well as into terms due to the cooperative effects of more than one x i . In this paper we assume that one has computed the full set of first order sensitivity indices as well as the full set of total-order sensitivity indices (a fairly common strategy in sensitivity analysis), and show that in this case the same set of model evaluations can be used to compute double estimates of:

• the total effect of two factors taken together, for all such k 2 couples, where k is the dimensionality of the model; • the total effect of k -2 factors taken together, for all k 2 such (k -2) ples.

We further introduce a new strategy for the computation of the full sets of first plus total order sensitivity indices that is about 50% cheaper in terms of model evaluations with respect to previously published works.

We discuss separately the case where the input factors x i 's are not independent from each other.

Introduction

Global sensitivity analysis aims to quantify the relative importance of input variables or factors in determining the value of an assigned output variable y. A recent review of applications of this discipline can be found in [START_REF]Sensitivity Analysis, Probability and Statistics Series[END_REF][START_REF] Saltelli | Sensitivity analysis as an ingredient of modelling[END_REF]. If the input to the model y = f(x 1 ,x 2 ,...,x k ) is composed of independent random variables, the joint probability density function of the input is: E-mail address: andrea.saltelli@jrc.it (A. Saltelli). URL address: http://www.jrc.cec.eu.int/uasa. 1 P(x 1 ,x 2 ,...,x k ) = k i=1 p i (x i ).

(1)

Mean and variance of y can be computed as:

E(y) = ••• f(x 1 ,x 2 ,...,x k ) k i=1 p i (x i ) dx i , (2)
V(y)= ••• f(x 1 ,x 2 ,...,x k ) -E(y) 2 k i=1 p i (x i ) dx i = ••• f 2 (x 1 ,x 2 ,...,x k ) k i=1 p i (x i ) dx i -E 2 (y). (3)
If one of the input factors x j is fixed to a generic value xj , the resulting variance of y will be equal to:

V y | x j =x j = ••• f(x 1 ,x 2 ,..., xj ..., x k ) -E(y | x j =x j ) 2 k i=1 i =j p i (x i ) dx i = ••• f 2 (x 1 ,x 2 ,..., xj ..., x k ) k i=1 i =j p i (x i ) dx i -E 2 y | x j =x j . (4)
For the purpose of sensitivity analysis one is interested in eliminating the dependence upon the value xj by integrating V(y| x j =x j ) over the probability density function of xj , obtaining E V(y| x j ) = ••• f 2 (x i ,x 2 ,...,x j ..., x k ) k i=1 p i (x i ) dx i -E 2 y | x j =x j p j ( xj ) d xj .

(5)

Note that we have dropped the dependence xj from the left-hand side, as it disappears due to the integration. Subtracting Eq. ( 5) from Eq. (3) one obtains:

V(y)-E V(y| x j ) = E 2 y | x j =x j p j ( xj ) d xj -E 2 (y). (6)
The left-hand side of Eq. ( 6) is also equal to V(E(y | x j )), and is a good measure of the sensitivity of y with respect to factor x j . If one divides it by the unconditional variance, one obtains the so-called first order sensitivity index S j = V(E(y | x j ))/V (y).T heS i 's are nicely scaled in [0, 1]. Eq. ( 6) is computationally impractical. In a Monte Carlo frame, it implies a double loop: the inner one to compute E 2 (y | x j =x j ), and the outer to compute the integral over d xj . For this reason the integral in [START_REF] Iman | A distribution free approach to inducing rank correlation among input variables[END_REF] has been rewritten by Ishigami and Homma [START_REF] Ishigami | An importance quantification technique in uncertainty analysis for computer models[END_REF] as:

E 2 y | x j =x j p j ( xj ) d xi = ••• f x 1 ,
p i (x ′ i ) dx ′ i . (7)
The expedient of using the additional integration variable primed, allows us to realize that the integral in [START_REF] Ishigami | An importance quantification technique in uncertainty analysis for computer models[END_REF] is the expectation value of the function F of a set of (2k -1) factors:

F x 1 ,x 2 ,...,x j ,...x k ,x ′ 1 ,x ′ 2 ,...,x ′ j -1 ,x ′ j +1 ,...x ′ k = f(x 1 ,x 2 ,...x k )f x ′ 1 ,x ′ 2 ,...,x ′ j -1 ,x j ,x ′ j +1 ,...x ′ k . (8)
The integral ( 7) can be hence computed using a single Monte Carlo loop. The Monte Carlo procedure that follows was proposed by Saltelli et al. [START_REF] Saltelli | Some new techniques in sensitivity analysis of model output[END_REF]. Two input sample matrices M 1 and M 2 are generated:

M 1 =    x 11 x 12 ... x 1k x 21 x 22 ... x 2k ... x n1 x n2 ... x nk    , M 2 =    x ′ 11 x ′ 12 ... x ′ 1k x ′ 21 x ′ 22 ... x ′ 2k ... x ′ n1 x ′ n2 ... x ′ nk    , ( 9
)
where n is the sample size used for the Monte Carlo estimate. In order to estimate the sensitivity measure for a generic factor x j ,i.e.

S j = V(E(y| x j )) V(y) = (U j -E 2 (y)) V(y) , U j = E 2 y | x j =x j p j ( xj ) d xj (10)
we need an estimate for both E(y) and U j . The former can be either obtained from values of y computed on the sample in M 1 or M 2 . U j can be obtained from values of y computed on matrices M 1 and N j , the latter being defined as:

N j =    x ′ 11 x ′ 12 ... x 1j ... x ′ 1k x ′ 21 x ′ 22 ... x 2j ... x ′ 2k ... ... ... ... ... ... x ′ n1 x ′ n2 ... x nj ... x ′ nk    , (11)
i.e. by:

U j = 1 n -1 n r=1 f(x r1 ,x r2 ,...,x rk )f x ′ r1 ,x ′ r2 ,...,x ′ r(j-1) ,x rj ,x ′ r(j+1) ,...,x ′ rk . (12)
If one thinks of matrix M 1 as the "sample" matrix, and of M 2 as the "re-sample" matrix, then U j is obtained from products of values of f computed from the sample matrix times values of f computed from N j , i.e. a matrix where all factors except x j are re-sampled. In this way the computational cost associated with a full set of first order indices S i is n(k + 1).Onesetofn evaluations of f is needed to compute E(y),andk sets of n evaluations of f are needed for the second term in the product (12). The very same procedure for the computation of sensitivity indices was proposed by Sobol' [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]. The problem setting of Sobol' was that of identifying a subset of the k factors that could account for most of the variance of y. Imagine that the factors have been partitioned into a trial set u = (x i 1 ,x i 2 ,...,x i m ) and the remaining set v = (x l 1 ,x l 2 ,...,x l k-m ). Then according to Sobol' an idea of the overall effect of the subset u on the variance of the output can be estimated from:

U v = ••• f(u, v)f (u ′ , v) du du ′ dv, (13)
V E(y | v) = U v -E 2 (y), (14)
V E(y | u) + V E(y | u, v) = V(y)-V E(y | v) . (15)
In Eq. ( 15), V(E(y| u)) is the first order effect of the set u, while V(E(y| u, v)) is the interaction term between the sets u and v.IfV(y) ∼ = V(E(y| v)),thenu is non-influent, and all factors in u can be fixed in a subsequent analysis of the model. Formula [START_REF] Saltelli | Some new techniques in sensitivity analysis of model output[END_REF] shows the same expedient of the additional integration variables already described. The Monte Carlo estimate of U v is:

U v = 1 n -1 n r=1 f(x ri 1 ,x ri 2 ,...,x ri m ,x rl 1 ,x rl 2 ,...,x rl k-m )f x ′ ri 1 ,x ′ ri 2 ,...,x ′ ri m ,x rl 1 ,x rl 2 ,...,x rl k-m , (16)
i.e. to estimate the total effect of set u one must now re-sample the variables in the set u. One can easily see that ( 12) is a particular case of [START_REF] Saltelli | On the relative importance of input factor in mathematical models[END_REF]. Error estimates for U j 's are discussed in the original reference of Sobol'. A bootstrap based alternative is discussed in [START_REF] Archer | Sensitivity measures, ANOVA like techniques and the use of bootstrap[END_REF]. Eq. ( 15) is a particular case of a general variance decomposition scheme proposed by Sobol', whereby the total unconditional variance can be decomposed as:

V(y)= i V i + i j>i V ij +•••+V 12...k , (17)
where

V i = V E(Y | x i ) , V ij = V E(Y | x i ,x j ) -V i -V j
and so on. The development in [START_REF] Saltelli | Sensitivity analysis as an ingredient of modelling[END_REF] contains k terms of the first order V i ,k(k -1)/2 terms of the second order V ij and so on, till the last term of order k, for a total of 2 k -1t e r m s .T h eV ij terms are the second order (or two-way) terms, analogous to the second order effects described in experimental design textbooks (see, e.g., [2]). The V ij terms capture that part of the effect of x i and x j that is not described by the first order terms. Formula (17) has a long history, and various authors have proposed different versions of it. A discussion can be found in [START_REF] Archer | Sensitivity measures, ANOVA like techniques and the use of bootstrap[END_REF], as well as in [START_REF] Rabitz | Managing the Tyranny of Parameters in Mathematical Modelling of Physical Systems[END_REF]. Sobol's version of formula ( 17) is based on a decomposition of the function f itself into terms of increasing dimensionality, i.e.:

f(x 1 ,x 2 ,...,x k ) = f 0 + i f i + i j>i f ij +•••+f 12...k , (18)
where each term is function only of the factors in its index, i.e. f i = f i (x i ), f ij = f ij (x i ,x j ) and so on. Decompositions [START_REF] Saltelli | Sensitivity analysis as an ingredient of modelling[END_REF], [START_REF] Saltelli | A quantitative, model independent method for global sensitivity analysis of model output[END_REF] are unique provided that the input factors are independent and that the individual terms f i 1 i 2 ...i s in [START_REF] Saltelli | A quantitative, model independent method for global sensitivity analysis of model output[END_REF] are square integrable and have zero mean over the domain of existence. One important aspect of Sobol' development is that similar decompositions can be written by taking the factors into subsets, as shown by Eq. [START_REF] Saltelli | About the use of rank transformation in sensitivity analysis of model output[END_REF]. This prompted Homma and Saltelli [START_REF] Homma | Importance measures in global sensitivity analysis of model output[END_REF] to introduce the new estimate U -j :

U -j = 1 n -1 n r=1
f(x r1 ,x r2 ,...,x rj ,...,x rk )f x r1 ,x r2 ,...,x r(j-1) ,x ′ rj ,x r(j+1) ,...,x rk .

As before:

V E(y | x -j ) = U -j -E 2 (y), (19)
where V(E(y| x -j )) is the total contribution to the variance of y due to non-x j . This implies that the difference V(y)-V(E(y | x -j )) is equal to the sum of all terms in the variance decomposition (15) that include x j .W e illustrate this for the case k = 3:

S T 1 = V(y)-V(E(y| x -1 )) V(y) = E(V (y | x -1 )) V(y) = S 1 + S 12 + S 13 + S 123 , (20)
where, e.g., S 1 = V(E(y| x 1 ))/V (y). Analogous expressions can be written for S T 2 ,S T 3 . We have called the S T j 's "total effect" terms. The total effects are useful for the purpose of sensitivity analysis, as discussed in [START_REF] Saltelli | A quantitative, model independent method for global sensitivity analysis of model output[END_REF], as they give information on the non additive part of the model. It may be useful to observe here that for a purely additive model, k i=1 S i = 1, while for a given factor x j an important difference between S T j and S j flags an important role of interactions for that factor in y. The same information could be obtained by computing all terms in [START_REF] Saltelli | Sensitivity analysis as an ingredient of modelling[END_REF], but these are as many as 2 k -1. This problem has been referred to by Rabitz et al. [START_REF] Rabitz | Efficient input-output representations[END_REF] as "the curse of dimensionality". The computational cost of estimating all effects in [START_REF] Saltelli | Sensitivity analysis as an ingredient of modelling[END_REF] is in fact as high as n2 k , where again n is the sample size used to estimate one individual effect. 1 For these reasons we customarily tend to compute the set of all S i plus the set of all S T i , that gives a fairly good description of the model sensitivities. This would normally entail a computational cost of n(2k + 1) model evaluations, i.e. nk for the first order terms, again nk for the total effect terms, plus n for E(y). In fact we have found in Homma and Saltelli [START_REF] Homma | Importance measures in global sensitivity analysis of model output[END_REF], that better estimates for the first order terms are obtained if the E 2 (y) term in ( 10) is estimated as

E 2 = 1 n n r=1 f(x r1 ,x r2 ,...,x rk )f x ′ r1 ,x ′ r2 ,...,x ′ rk (21)
rather than from

E 2 = 1 n n r=1 f(x r1 ,x r2 ,...,x rk ) 2 , (22)
i.e. using sample estimates from both M 1 and M 2 matrices rather than from M 1 alone. Eq. ( 21) is a legitimate estimate of the squared sample mean given the independence of the two sample vectors. It is clear from [START_REF] Rabitz | Managing the Tyranny of Parameters in Mathematical Modelling of Physical Systems[END_REF] and ( 12) that the estimate of S j goes more naturally to zero for a non-influential factor x j when (21) is used, as can be seen from:

U j -E 2 (y) = 1 n -1 n r=1 (f (x r1 ,x r2 ,...,x rk )f x ′ r1 ,x ′ r2 ,...,x ′ r(j-1) ,x rj ,x ′ r(j+1) ,...,x ′ rk - 1 n n r=1 f(x r1 ,x r2 ,...,x rk )f x ′ r1 ,x ′ r2 ,...,x ′ rk . (23)
On the same ground one can see that the computation of the total effect sensitivity indices is better achieved using (22). V(y)is computed from M 1 for all indices. In conclusion, the standard computational strategy so far employed to compute the full set of total and first order indices entailed a total of n(2k + 2) model evaluations, two samples beingusedtoestimate E 2 . Many applications of this strategy to different models can be found in various chapters of Saltelli et al. [START_REF]Sensitivity Analysis, Probability and Statistics Series[END_REF].

An important economy in model evaluation, that is described in [START_REF] Saltelli | A quantitative, model independent method for global sensitivity analysis of model output[END_REF], is that the S T j and S j terms can also be estimated using an extended version of the Fourier Amplitude Sensitivity Test (FAST). When using extended FAST, the same set of n model evaluations that was used to estimate a given S T j can also be used to compute S j , that makes the entire analysis feasible at the cost of nk model evaluations. For this reason the extended FAST method has been considered so far as the most efficient way to compute the full set of S T j and S j indices. In the present work we introduce an extended version of Sobol' method that out-performs FAST.

Extension of the method

We imagine that a model has been characterised using the Sobol' method, computing all S T j 's and S j 's estimates at the cost of n(2k +2) model evaluations. Can the same coefficient be estimated with a smaller cost? Can additional estimates be obtained with the same sets used to compute the S T j 's and S j 's? Before we proceed we need to introduce some new notation.

Let us call V c i 1 i 2 ...i s a sensitivity measure that is closed within a subset of factors, i.e. V c i 1 i 2 ...i s is the sum of all

V i 1 i 2 ...i s terms in (17) that is closed in the indices i 1 ,i 2 ,...,i s : V c 1 = V 1 , V c ij = V i + V j + V ij
, and so on. Likewise V c

-i 1 i 2 ...i s will indicate the sum of all indices that are closed within the complementary set of i 1 ,i 2 ,...,i s ,i . e .

V c

-i 1 i 2 ...i s = V c l 1 l 2 ...l k-s ,wherei p = l q for all p ∈[1, 2,...,s],q ∈[1, 2,...,ks]. Let a i 1 i 2 ...i s denote the vector of length n containing model evaluations corresponding to the rows of the input factor matrix N i 1 i 2 ...i s . As in Eq. ( 11) above, the matrix N i 1 i 2 ...i s is obtained from matrix M 1 by substituting all columns except i 1 ,i 2 ,...,i s by the corresponding columns of matrix M 2 . a 0 will hence denote a set of model evaluations corresponding entirely to matrix M 2 , while a i 1 i 2 ...i k will indicate the vector of model evaluations corresponding entirely to matrix M 1 .

A few equalities are repeated below for reader's convenience:

V c i 1 i 2 ...i s = V E(Y | x i 1 x i 2 ...x i s ) = U i 1 i 2 ...i s -E 2 (y), (24)
U i 1 i 2 ...i s = 1 n -1 n r=1 f(x r1 ,x r2 ,...,x rk )f x ri 1 ,x ri 2 ,...,x ri s ,x ′ rl 1 ,x ′ rl 2 ,...,x ′ rl k-s , (25)
U -i 1 i 2 ...i s = 1 n -1 n r=1 f(x r1 ,x r2 ,...,x rk )f x ′ ri 1 ,x ′ ri 2 ,...,x ′ ri s ,x rl 1 ,x rl 2 ,...,x rl k-s (26)
with the special cases

S j = ( U j -E 2 (y)) V(y) , (27)
S T j = 1 - ( U -j -E 2 (y)) V(y) . ( 28
)
We are now ready to submit the following theorem:

Theorem 1. G i v e nam o d e ly = f(x 1 ,x 2 ,...,x k )
, it is possible to compute at the cost of n(k + 2) model evaluations:

(1) One estimate for each of the k indices of the first order S j .

(2) One estimate for each of the k total effect indices S T j .

(3) One estimate for each of the k 2 V c -ij closed effect indices.

a 0 V(y) a 1 S T 1 V(y) a 2 S T 2 V c -12
V(y)

a 3 S T 3 V c -13 V c
-23

V(y)

a 4 S T 4 V c -14 V c -24 V c
-34

V(y)

a 5 S T 5 V c -15 V c -25 V c -35 V c
-45

V(y)

a 2345 S 1 E 2 (y) V c 12 V c 13 V c 14 V c 15 V(y) a 1345 S 2 V c 12 E 2 (y) V c 23 V c 24 V c 25 V c -12
V(y)

a 1245 S 3 V c 13 V c 23 E 2 (y) V c 34 V c 35 V c -13 V c
-23

V(y)

a 1235 S 4 V c 14 V c 24 V c 34 E 2 (y) V c 45 V c -14 V c -24 V c
-34

V(y)

a 1234 S 5 V c 15 V c 25 V c 35 V c 45 E 2 (y) V c -15 V c -25 V c -35 V c
-45

V(y)

a 12345 E 2 (y) S 1 S 2 S 3 S 4 S 5 S T 1 S T 2 S T 3 S T 4 S T 5 V(y)
An additional theorem is the following:

Theorem 2. If we modify the setting of Theorem 1 by allowing for n(2k + 2) model evaluations (i.e. as many as in the procedure of Section 1), we can obtain:

(1) Double rather than single estimates for each of the S T j 's and S j 's indices. (2) Double estimates for each of the k 2 V c ij terms. (3) Double rather than single estimates for each of the k 2 V c -ij terms.

Theorems 1, 2 constitute the promised extension of Sobol' method. Let us illustrate the new procedures for the case k = 5. We have to use this value as lower values of k are special cases and will be treated afterwards. Table 1 gives for each cell what term can be computed by the corresponding a i 1 i 2 ...i s vectors.

Note that:

(1) Table 1 can be interpreted by referring to Eqs. ( 24)-(28) above. E.g., we have labelled the entry corresponding to the intersection a 0 and a 1 as S T 1 ,asa 0 • a 1 yields U -1 that in turn can be used to compute S T 1 (Eq. ( 28)) and so on for the other terms.

(2) The diagonal has been labelled as providing an estimate of V(y), as this is what can be obtained by the scalar product a 2 i 1 i 2 ...i s . In fact each of the 2k + 2 vectors a i 1 i 2 ...i s can yield an estimate of E(y).Known E(y) each a i 1 i 2 ...i s can again be used to estimate V(y).

(3) The row labelled a 12345 illustrates the same procedure as in Section 1 for the computation of the first order indices and total order indices, i.e. S T 4 is obtained from V(y) and V T 4 , this latter being computed from a 12345 , a 1235 . (4) Looking at the column a 0 , one sees that the same set of indices (first order plus total) can be computed from a 0 , a 12345 , and either of the sets {a 1 , a 2 , a 3 , a 4 , a 5 } or {a -1 , a -2 , a -3 , a -4 , a -5 }≡{a 2345 , a 1345 , a 1245 , a 1235 , a 1234 }, bringing the computational cost from n(2k + 2) down to n(k + 2), with a reduction in computational cost that tends to 50% at increasing k values.

(5) All indices in rows other than a 12345 and columns other than a 0 are novel, in the sense that they were overlooked in the procedure of Section 1. The alternative arranging of the a i 1 i 2 ...i s terms shows the additional terms that can be computed. [START_REF] Iman | A distribution free approach to inducing rank correlation among input variables[END_REF] The intersection of vectors a 1 and a 2345 has been labelled as an estimate of E 2 (y), as all columns in the two sampling matrices are different and the scalar product a i 1 i 2 ...i s , a j 1 j 2 ...j r provides an estimate of the square of E(y), as in Eq. ( 21) above. [START_REF] Ishigami | An importance quantification technique in uncertainty analysis for computer models[END_REF] The two vectors a 2 and a 2345 allow the computation of V c 12 as columns 1 and 2 are identical in the two sampling matrices. [START_REF] Mckay | Variance-based methods for assessing uncertainty importance in NUREG-1150 analyses[END_REF] The two vectors a 2345 and a 1345 allow the computation of V c 345 = V c -12 as columns 3,4,5 are identical in the two sampling matrices.

Looking at this table, for the setting of Theorem 1 (grey cells in Table 1), it is easy to see that we have produced the S T j 's and S j indices, with j ∈[1, 2, 3, 4, 5] at the reduced cost of n(k + 2) = 7n model evaluations, instead of n(2k + 2) = 12n, with a reduction of 42% in computational cost. Furthermore, we have produced one estimate for each of the 5 3 = 10 indices complementary to the second order ones, that for k = 5happentobeclosedindices of the third order. Note that for, e.g., k = 6 we would have obtained one estimate for each of the 6 4 = 15 closed indices of the fourth order and so on for larger values of k, and so on based on the known property that k j = k k-j for k j .

For the setting of Theorem 2, it is easy to see that double, rather than single, estimates for each of the S T j 's, S j and V c -ij have been produced. We have additionally obtained double estimates for each of the 5 2 = 10 closed indices of the second order. Additional estimates of E 2 (y) are also available as discussed.

The reader might wonder which among the various estimates of E 2 (y), V(y)should be used in Eqs. ( 24)-(28) to obtain, e.g., the V c i 1 i 2 ...i s from the U i 1 i 2 ...i s . In [START_REF] Homma | Importance measures in global sensitivity analysis of model output[END_REF] we suggest that the estimate of E 2 (y) in ( 21) obtained from a 0 , a 12345 should be used for the first order indices and that from (22) based on a 0 alone for the total effect ones.

In the context of the extended procedures (Theorems 1 and 2) presented here, the following approach was taken:

(1) Theorem 1 setting. Eq. ( 21) is used for the first order indices. This means that for computing, any of the S j , a 0 and a 12345 vectors are used to estimate E 2 (y) and a 12345 to compute V(y). Eq. ( 22) is used for the total effect, i.e. for any of the S T j , a 0 alone is used to estimate E 2 (y). V(y)is also computed from a 0 for the total effect indices. For the closed effects of order k -2 Eq. ( 22) is used for E 2 (y), and the vector to be used in ( 22) is selected as one of the two that concur in the estimation of the effect. E.g., for V c -15 (in the grey table area), computed from a 1 and a 5 ,the E 2 (y) is computed from a 1 alone (or identically from a 5 alone). V(y)is computed from the same vector used for E 2 (y) (either a 1 or a 5 ).

(2) Theorem 2 setting. Same procedure as Theorem 1 for all double estimates of (i) first order, (ii) total order indices and (iii) order k -2 closed indices. For the closed indices of the second order Eq. ( 21) is used, where E 2 (y) is computed using one of the vector that concur in the estimation of the index. E.g., for that estimate of V c 12 , that is computed from a 1345 and a 1 , E 2 (y) is computed from a 1345 and a 2 (or identically from a 1 and a 2345 ). The variance is correspondingly computed from a 1345 or a 1 . These arrangements can be easily understood by inspecting equations like (23) above.

As we said, k = 4 is a special case (Table 2). For this value of k we obtain for the setting of Theorem 2:

(1) Double estimates for each of the 4 S i and each of the 4 S T i .

(2) Quadruple estimates of the 4 2 = 6 second order terms V c ij .

All 4 estimates of each term V c ij are independent.

a 0 V(y) a 1 S T 1 V(y) a 2 S T 2 V c 34 V(y) a 3 S T 3 V c 24 V c 14 V(y) a 4 S T 4 V c 23 V c 13 V c 12 V(y) a 234 S 1 E 2 (y) V c 12 V c 13 V c 14 V(y) a 134 S 2 V c 12 E 2 (y) V c 23 V c 24 V c 34 V(y)
a 124 S 3 V c 13 V c 23 E 2 (y) V c 34 V c 24 V c 14 V(y) a 123 S 4 V c 14 V c 24 V c 34 E 2 (y) V c 23 V c 13 V c 12 V(y) a 1234 E 2 (y) S 1 S 2 S 3 S 4 S T 1 S T 2 S T 3 S T 4 V(y)
a 0 V(y) a 1 S T 1 V(y) a 2 S T 2 S 3 V(y) a 3 S T 3 S 2 S 1 V(y) a 23 S 1 E 2 (y) V c 12 V c 13 V(y) a 13 S 2 V c 12 E 2 (y) V c 23 S 3 V(y) a 12 S 3 V c 13 V c 23 E 2 (y) S 2 S 1 V(y) a 123 E 2 (y) S 1 S 2 S 3 S T 1 S T 2 S T 3 V(y)
For k = 3 we obtain (Table 3):

(1) Double estimates for each of the 3 S i and each of the 3 S T i .

(2) Double estimates of the 3 2 = 3 second order terms V c ij .

(3) Two more estimates for each of the 3 S i .

The case k = 2 is non-relevant, as V c 12 = V(y).

Discussion of the methodological advantages

What benefit does the new computational procedure offer? The main advantage of the new method is that, given the computational effort already made to compute a full set of S i , S T i estimates, one can also obtain additional estimates.

Imagine, for the case of k 5, that the reduced procedure of Theorem 1 has been adopted, and that a 0 , a i 1 ,i 2 ,...,i k , and either of the sets {a 1 , a 2 ,...,a k } or {a -1 , a -2 ,...,a -k } has been computed (at the cost of n(k + 2) model evaluations). Beside S i , S T i , we now also dispose of the S c -ij indices. These can be very useful. If S c -ij ≈ 1, it means that neither x j nor x i contribute appreciably to the variance of y, neither by themselves nor in cooperation with other factors. These factors could hence be fixed in a subsequent analysis. Note that the condition S c -ij ≈ 1is equivalent to S T i ≈ 0, S T j ≈ 0. If we use instead the extended procedure of Theorem 2, at the cost of n(2k + 2) model evaluations, we obtain double estimates of S i , S T i , double estimates of the S c -ij indices and finally double estimates of V c ij for the closed effects of couples of factors, i.e. V c ij = V i + V j + V ij . These can be used as such or converted into second order sensitivity coefficients S ij = ( V c ij -V i -V j )/ V(y). A full set of second order coefficients is likely to tell us much of what we need to know about a model sensitivity, also because interaction of higher orders should in general be less frequent, as discussed in [START_REF] Rabitz | Managing the Tyranny of Parameters in Mathematical Modelling of Physical Systems[END_REF].We know from the value of S j if a factor is influent at the first order, and from S T j whether it has important interactions. If this is the case, inspection of the S ij for all i = j will allow us to identify which factor x j interacts with.

The case of the correlated input

Sensitivity analysis for correlated input is discussed in [START_REF] Mckay | Variance-based methods for assessing uncertainty importance in NUREG-1150 analyses[END_REF][START_REF] Saltelli | On the relative importance of input factor in mathematical models[END_REF]. For this setting, the important computational simplifications described in Section 1 are not applicable, and Eq. ( 17) loses its uniqueness. In these cases there is no alternative to the computation of the double loop needed to estimate conditional variances such as V(E(y| x j )), Eq. ( 6). For the purpose of Monte Carlo simulations, correlated input can be generated using Markov Chain Monte Carlo (MCMC), or procedures based on Cholesky decompositions (see, e.g., [START_REF] Iman | A distribution free approach to inducing rank correlation among input variables[END_REF]) or on Latin Hypercube Sampling (LHS, [START_REF] Mckay | Variance-based methods for assessing uncertainty importance in NUREG-1150 analyses[END_REF]). The problem with correlated sample, in brief, is that the reduction in variance that can be achieved by fixing one factor depends on whether or not other factors have been fixed, and the incremental reduction in variance for each factor depends on the order in which factors are fixed.

We have discussed in [START_REF] Saltelli | On the relative importance of input factor in mathematical models[END_REF] two general settings for sensitivity analysis. Each setting is based on a bet posed on the model y = f(x 1 ,x 2 ,...,x k ), for the general case where the input can be correlated. In the first bet, one has to make a rational guess on which parameter would induce the largest reduction in variance if it were fixed to its "true" value. Because such true value is in general unknown, the bet can be rationally placed by computing the estimates V(E(y| x j )), whether or not the input is correlated.

For the second setting, of relevance in risk analysis and control theory, the bet is on the identification of the smallest subset of x capable of inducing a target reduction in the unconditional variance V(y), as in the work of Sobol' [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF], discussed in Section 1. For the uncorrelated case, a rational selection strategy for the subset of interest is based on the computation of the full sets of S j and S T j . This strategy is meant to fight the curse of dimensionality, as attempting all combinations of factors, in a brute-force search for the smallest subset of x that gives the desired reduction in V(y), would be computationally prohibitive; one would have to compute all 2 k -1 terms in Eq. ( 17) to start with. As described in [START_REF] Saltelli | On the relative importance of input factor in mathematical models[END_REF], an iterative procedure can be adopted for the uncorrelated case, that includes as a step the computation of the full set of S j and S T j . For the correlated case, one might still engage in a brute force search computing all possible closed terms V c i 1 i 2 ...i s . Note that for the correlated case the V c i 1 i 2 ...i s can no longer be decomposed meaningfully into a sum of lower dimensionality terms, but would still allow a perfectly informed choice, as would the V i 1 i 2 ...i s in the uncorrelated case. Also for the correlated case, we suggest in [START_REF] Saltelli | On the relative importance of input factor in mathematical models[END_REF] a cheaper, albeit approximate, alternative that involves the computation of the S j and S T j for the non-correlated problem. Hence, in the context covered by these problem settings, the procedure proposed in Section 2 can still be usefully applied. Some of the uncertain input factors in PMOD are time-dependent physical quantities: the rock's total organic carbon content (TOC), the rock's porosity, permeability and thickness. In the analysis, the values at the initial time point, that corresponds to 30 million years before present (mybp), have been considered, neglecting the time dependency.

One of the inputs, the PHI stress variable, describes the mechanic behaviour of the rock. Due to lack of data on the specific site, three different curves relative to similar sedimentary basins for other areas have been used, and the model selects randomly which of the three to use at runtime. Similarly for the so-called "KEM/FIZ" files, that describe the stoichiometrics and kinetics of the chemical system considered. Eight alternative such descriptions were generated by the experts, so that sampling from these might be considered as representative of the system chemistry's variability. Also in this case the model selects at runtime one of the eight files. Finally the model needs as input 4 highly correlated time series (temperature, pressure, effective stress and hydrostatic pressure). Thirty-two such multivariate series have been generated by the experts and the model selects one set at random for each execution of the PMOD model. A summary of the 7 input factors is given in Table 4. All these factors are considered independent from each other.

A sensitivity analysis for this model had already been performed before the algorithms presented in the present paper were developed ( [START_REF] Tarantola | Global sensitivity analysis techniques for the analysis of the oil potential of sedimentary basins[END_REF]), where all coefficients of the first and total order had been computed. We have hence repeated the analysis using the setting of Theorem 2. The model output is composed of cumulative expelled amounts of oil, gas (CH4) and wet gas (CHx) at selected time points (i.e. at 30, 11.5, 8.5, 4.8, 1.9, 0 mybp). Figs. 4(a)-4(c) show the first order sensitivity indices obtained with the two approaches at different time points for the ouput CH4, as almost identical results hold for CHx and oil. Similarly for the total order indices, Figs. 5(a)-5(c).

In Figs. 6(a)-6(c) the coefficients of the second orders have been also computed using the set of Theorem 2, and compared with estimates obtained previously using an independent sample of size n for each index (i.e. 21n additional runs, from [START_REF] Tarantola | Global sensitivity analysis techniques for the analysis of the oil potential of sedimentary basins[END_REF]).

For all these Figs. 5-6 there is a general agreement between the two methods, and the confidence bounds, computed exactly as in the previous case study, are lower for the estimate from Theorem 2, as expected.

Finally the closed indices of order 5 are given in Fig. 7, as computed with Theorem 2. The results for Theorem 1 are very similar, as expected, and not shown here.

With the new procedures, we have been able to compute at no extra (or at a reduced) computational cost the coefficients of order (k -2), that allow us to identify the non influential factors. Given the highly non linear and non additive nature of PMOD, there are time points, such as t = 8.5 mybp, where none of the factors is non-influent, apart from source thickness and porosity (Fig. 7(b)). At time t = 0 mybp the coefficients of order (k -2) help us to rule out as non-influent all factors except TOC.

Conclusions

The present paper has suggested some efficient procedures for numerical experiments aimed at sensitivity analysis of model output. We have focused here on the computation of sensitivity indices that are based on decomposing the variance of the target function in a quantitative fashion. The approach presented here opens a road to fight the so-called "curse of dimensionality", that hinders the use of quantitative sensitivity analysis for computationally expensive models. The analyst willing to use such methods disposes now of two approaches to tackle the system. One is a parsimonious procedure (Theorem 1) that gives all effects of the first and total order, plus all those of order k -2, at the cost of n(k + 2) simulations. We have thus both reduced the computational cost for the standard procedure of some 50% and extended it to compute the indices of order k -2.

A second possible approach is the more expensive procedure (Theorem 2, cost = n(2k + 2)) that gives more robust estimates of the index of the first and total orders, plus estimates of all indices of order 2 and k -2.

Even assuming for n the value of 1000, the two procedures appear affordable for models whose cost per run is in the range from milliseconds or lower to some minutes. For models whose execution is in the tenths of minutes to a day range, quantitative methods are not applicable and efficient qualitative methods such as that of Morris [START_REF] Morris | Factorial sampling plans for preliminary computational experiments[END_REF] should be used (see [START_REF] Campolongo | Design of experiment, Sensitivity Analysis, Probability and Statistics Series[END_REF] for a review).

Fig. 1 .Fig. 2 .Fig. 3 .

 123 Fig. 1. First order (a) and total effect (b) indices for the g function with bootstrap-estimated 95% confidence bounds using the standard procedure of Section 1.

Fig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .

 4567 Fig. 4. Main effect sensitivity indices, with bootstrap-estimated 95% confidence bounds, for CH4 at three different time points: 11.5 (a), 8.5 (b) and 0 (c) million years before present.

  f x 1 ,x 2 ,..., xj ..., x k f x ′ 1 ,x ′ 2 ,..., xj ..., x ′

		k	k
		k	p i (x i ) dx i	p i (x ′ i ) dx ′ i p j ( xj ) d xj
		i=1 i =j	i=1 i =j
		k	k
	=	••• f(x 1 ,x 2 ,...,x j ..., x k )f x ′ 1 ,x ′ 2 ,...,x j ..., x ′ k	p i (x i ) dx i
		i=1	i=1 i =j

x 2 ,..., xj ,...,x k k i=1 i =j p i (x i ) dx i 2 p j ( xj ) d xj = •••

Table 1

 1 Terms that can be estimated given the corresponding vectors of model evaluations, k = 5

	a 0	a 1	a 2	a 3	a 4	a 5	a 2345	a 1345	a 1245	a 1235	a 1234	a 12345

Table 2

 2 Terms that can be estimated given the corresponding vectors of model evaluations, k = 4

	a 0	a 1	a 2	a 3	a 4	a 234	a 134	a 124	a 123	a 1234

Table 3

 3 Terms that can be estimated given the corresponding vectors of model evaluations, k = 3

	a 0	a 1	a 2	a 3	a 23	a 13	a 12	a 123

Table 4

 4 List of the uncertain input factors and their stochastic properties

	Factor's	Factor's	Type	Range of	PDF
	name	number		values
	"KEM/FIZ" files	1	Discrete	1,...,8	Uniform
	"phi-stress" curves	2	Discrete	1, 2, 3	Uniform
	TOC	3	Continuous	0.005-0.05	Uniform
	Porosity	4	Continuous	Min = 0.04	Triangular
				Mode = 0.05
				Max = 0.09
	Permeability	5	Continuous	1.e-9-1.e-6	Log-uniform
	Source thickness	6	Continuous	Min = 907	Triangular
				Mode = 1814
				Max = 2721
	Time-series	7	Discrete	1, 2,...,32	Uniform

n(2 k -1) would be needed to compute all effects, and n more to compute E(y), V(y).

In Russian.

Acknowledgements

The author is grateful to Marco Ratto for MATLAB coding of the algorithm, and Stefano Tarantola, Michaela Saisana, for reviewing this manuscript. Paolo Ruffo and Anna Corradi of ENI (AGIP Division) kindly produced the PMOD evaluations and allowed their use in the present work. ENI-AGIP funded part of the this study.

Test cases

We illustrate the algorithm on non-correlated test cases. We have selected an analytic function due to Sobol' and known as "Sobol' g function", and a more complex numeric calculus test case originating from modelling of petroleum generation in sedimentary basins. The cost of computing the former can be assumed as zero, while the computation of a single output time series for the latter takes about 0.05 s on a 8-nodes Linux cluster, 16 CPU's with a Pentium 3 processor running under Linux RedHat.

The g function is a strongly non-monotonic, non-additive function of k factors x i , assumed identically and uniformly distributed in the unit cube I k ={x | 0 x i 1; i = 1, 2,...,k}.

with

For each of the g i (x i ) functions 1 0 g i (x i ) dx i = 1, and for

The importance of each factor x i is driven by its associated coefficient a i .F o ra i = 0, the factor is important (0 g i (x i ) 2).F or ,e . g. ,a i = 9 the factor in non-important (0.9 g i (x i ) 1.1), while for a i = 99 the factor can be considered as non-influent (0.99 g i (x i ) 1.01). Analytical expressions are available for the sensitivity indices ( [START_REF] Archer | Sensitivity measures, ANOVA like techniques and the use of bootstrap[END_REF][START_REF] Saltelli | About the use of rank transformation in sensitivity analysis of model output[END_REF]):

In Figs. 123we have computed the sensitivity indices for a 6-dimensional g-function with a ={0, 0.5, 3, 9, 99, 99} using first the standard procedure of Section 1, at the cost of n(2k + 2) model evaluations, then with the restricted procedure of Theorem 1 at the cost of n(k + 2) model evaluations. Finally we have used the setting of Theorem 2, at the cost of n(2k + 2) model evaluations. For all experiments n = 1024, and the standard error associated with the computation of the sensitivity indices was computed using bootstrap, as described in [START_REF] Archer | Sensitivity measures, ANOVA like techniques and the use of bootstrap[END_REF], with a bootstrap sample dimension of 10,000.

Comparing Figs. 1(a)-1(b) (S j and S T j by the standard procedure) with 2(a)-2(b) (S j and S T j by the procedure of Theorem 1), we can see that the quality of the estimates is the same. Fig. 2(c) shows the advantage brought by the term of the 4th order, especially to identify couples of non-influent factors ij ={45, 46, 56}.

Moving to the procedure of Theorem 2, Figs. 3(a)-3(b), we see that the confidence bound of the estimates is lower (each estimate is the average of 2). Similarly for Fig. 3(c), to be compared with 2(c). Additional insight into the structure of the model is offered by the new Fig. 3(d), with the second order indices.

The PMOD model computes the generation and expulsion of hydrocarbons from a host rock, and is part of a suite of computer models used to estimate the oil potential of sedimentary basins. PMOD has been originally developed at Lawrence Livermore Laboratory ( [START_REF] Braun | User's manual for PMOD, A pyrolysis and primary migration model[END_REF]) and adapted at ENI-AGIP for its basin modelling environment ( [12]).


Anne Marie Dupuy

Nils Kuster

Anne-Sophie Bargnoux

Sylvain Aguilhon

Fabien Huet

Florence Leclercq

Jean Luc Pasquié

François Roubille

Jean-Paul Cristol

Long term pronostic value of suPAR in chronic heart failure: reclassification of patients with low MAGGIC score

à 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.


Anne Marie Dupuy

Luna Ruffel

Anne-Sophie Bargnoux

Stéphanie Badiou

Jean-Paul Cristol

Analytical evaluation of the performances of a new procalcitonin immunoassay