Taotao Zhang

Tao Wang
email: twang@itpcas.ac.cn

Yingying Feng

Xichen Li

Gerhard Krinner

Keywords: Eurasian snow cover, Eastern China summer rainfall, summer North Atlantic Oscillation, emerging impact

Eurasian spring snow cover is widely considered as an important predictor of Asian summer monsoon rainfall, but its possible role in the formation of the north-south dipole structure of rainfall anomalies (NSDR)-a major mode of the eastern China summer rainfall variability-remains elusive. Here, we show that, there is a close connection between the western Eurasian spring snow cover (WESS) and NSDR during our research period 1967-2018, with less WESS tends to be accompanied by a wetter south-drier north pattern over eastern China, and vice versa. However, this relationship was not significant before the late 1990s, but has since become significant. Further analyses demonstrate that the shift in the WESS-NSDR relationship could be attributed to the modulation of summer North Atlantic Oscillation (SNAO). After the late 1990s, the WESS-related anomalous atmospheric circulations during summer are largely reinforced by the constructive superposition of those with same signs induced by SNAO, which in turn would intensify the impact of WESS and hence lead to a strong WESS-NSDR connection. In contrast, the influences of WESS are counteracted by those with opposite signs associated with SNAO before the late 1990s and thereby result in a weak snow-rainfall relationship. Our findings, along with the decline in Eurasian spring snow cover, provide a potential explanation for the recent 'South Flood-North Drought' pattern observed over eastern China.

Introduction

Over one billion people live in eastern China, where the climate regime has been renowned for being distinctly wet in the south and relatively dry in the north. This marked contrast hastened the construction of the huge South-to-North water diversion project to address northern water resource shortages. Moreover, an emergent meridional dipole structure of summer rainfall anomalies, known as the 'South Flood-North Drought' pattern (Ding et al 2008, Huang et al 2013[START_REF] Zhu | Contribution of the phase transition of Pacific Decadal Oscillation to the late 1990s' shift in East China summer rainfall[END_REF], He et al 2017), has further increased the uneven distribution of water resources over eastern China. An improved understanding of the formation of this dipole pattern is crucial for flood and drought disaster reduction and socioeconomic development in this densely populated region.

Eurasian spring snow cover offers some potential for prediction of Asian summer monsoon rainfall via snow-atmosphere coupling (Barnett et al 1989, Cohen and Rind 1991, Turner and Slingo 2011, Halder and Dirmeyer 2017). It has been suggested that a negative anomaly of spring snow cover over Eurasia can induce increased rainfall in southern China, and a rainfall deficit in northern and northeastern China, i.e. the 'South Flood-North Drought' pattern (Wu et al 2009). In this view, relatively excessive Eurasian spring snow cover before the late 1980s would have favored an opposite 'South Drought-North Flood' pattern, but such a dipole pattern was not seen during the 1970s and 1980s [START_REF] Zhu | Contribution of the phase transition of Pacific Decadal Oscillation to the late 1990s' shift in East China summer rainfall[END_REF]Wang 2015, He et al 2017). This discrepancy implies that the nature of the snow-rainfall relationship over eastern China is not yet well understood. Furthermore, there is clear evidence to show that reduced snow cover due to climate warming, could weaken the relationship between snow cover and Indian summer rainfall (Kripalani et al 2003, Wu et al 2014[START_REF] Zhang | The weakening relationship between Eurasian spring snow cover and Indian summer monsoon rainfall[END_REF]. It remains unclear whether such a weakening is currently occurring over eastern China.

In this study, we aim to decipher the temporal evolution of the impact of Eurasian spring snow cover on eastern China summer rainfall, by using the longest available satellite-based records of snow cover data. The main purpose of this study is to clarify the reliability of using Eurasian spring snow cover to forecast seasonal eastern China summer rainfall, and give some clues as to the cause of the reported 'South Flood-North Drought' phenomenon over eastern China.

Data and methods

Datasets

The snow cover data used in this study is the widely used Northern Hemisphere EASE-Grid 2.0 weekly snow cover and sea ice extent (version 4), provided by the National Snow and Ice Data Center (NSIDC) (Brodzik and Armstrong 2013). The NSIDC snow cover data is the currently longest satellite derived snow data covering the period from 1966 to the present, which has been approved to appropriate for the large-scale climate research (Brodzik and Armstrong 2013). Although NSIDC data has some uncertainties near the boundary of the continental snowpack, there is large-scale agreement between it and other products for the most of the snow season and most regions (Frei et al 2012, Bormann et al 2018). The original snow cover data are stored in flat binary files indicating the presence or absence of snow in a specific week. The raw data are provided in the polar stereographic projection at a spatial resolution of 25 km. To facilitate comparison with other data, we converted the data to a regular 1 • × 1 • longitudelatitude grid. Here, we used the snow cover fraction (SCF) to represent the variability of snow cover. The monthly SCF is calculated by summing the number of weeks in that month for which snow was present for a pixel and expressing this number as a fraction of the total number of weeks in that month (Bamzai and Shukla 1999). In this study, the average of March, April and May is referred to as spring and the average of June, July and August as summer.

The monthly mean precipitation data used in this study was derived from the Climate Research Unit, University of East Anglia (CRU TS v4.03) which has a spatial resolution of 0.5 2003). This index has been recognized as the optimal representation of the spatial-temporal variability associated with the NAO on all timescales in a comparison with other indices (Li and Wang 2003). Following Chen and Zhai (2015), the East Asia/Pacific (EAP) teleconnection pattern (also known as the Pacific-Japan teleconnection) index was defined based on three key anomaly center points, the Sea of Okhotsk (OK, 60 • N, 130 • E), East Asia (EA, 37.5 • N, 120 • E), and the western Pacific (WP, 20 • N, 120 • E) (Chen and Zhai 2015, Li et al 2018):

EAP = 1 3 H OK - 1 3 H EA + 1 3 H WP (1)
where H OK , H EA , and H WP denote the normalized 500 hPa geopotential heights of OK, EA, and WP, respectively.

Calculation of water vapor flux and wave activity flux

To diagnose the moisture transport, the vertically integrated water vapor flux (Q) and its divergence (Q div ) are calculated based on the monthly NCEP-NCAR reanalysis data as follows:

Q = - 1 g pt ps qVdp
(2)

Q div = ∇ • Q = - 1 g pt ps ∇ • (Vq) dp (3)
where g, q, p, V denote the gravitational acceleration, specific humidity, pressure, and horizontal wind components (u, v), respectively. The water vapor flux was integrated from the surface (p s ), 1000 hPa, to the top layer (p t ), 300 hPa, assuming negligible water vapor above 300 hPa.

The wave activity flux (WAF) proposed by Takaya and Nakamura (1997,2001) was used to indicate the propagation of stationary Rossby waves. The horizontal components of WAF were defined as follow:

WAF = p cos ϕ 2 |U| ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ U α 2 cos 2 ϕ ∂ψ ′ ∂λ 2 -ψ ′ ∂ 2 ψ ′ ∂λ 2 + V α 2 cos ϕ ∂ψ ′ ∂λ ∂ψ ′ ∂ϕ -ψ ′ ∂ 2 ψ ′ ∂λ∂ϕ U α 2 cos ϕ ∂ψ ′ ∂λ ∂ψ ′ ∂ϕ -ψ ′ ∂ 2 ψ ′ ∂λ∂ϕ + V α 2 ∂ψ ′ ∂ϕ 2 -ψ ′ ∂ 2 ψ ′ ∂ϕ 2 ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ (4)
where |U| represents the climatological magnitude of winds, α is the mean radius of the earth, p denotes the pressure scaled by 1000 hPa, ϕ and λ are latitude and longitude, U and V denote the climatological zonal and meridional wind, respectively, and ψ ′ is the perturbed geostrophic streamfunction, which was defined as the anomalies obtained by regression onto the involved time series.

Statistical analysis

We performed an empirical orthogonal function (EOF) analysis on the normalized rainfall anomalies to derive the spatiotemporal variability of eastern China summer rainfall. Here eastern China was defined as the region east of 105 • E. Since we only focus on interannual variability, a 9 year highpass Lanczos filter (Duchon 1979) was applied to all original variables before analyses, to avoid the potential contamination of interannual relationships by decadal signals. In addition, linear correlation and regression analyses are used to detect the statistical connections between variables. The significance of all the analyses was determined by the standard twotailed Student's t-test.

Results

Emerging Eurasian spring snow impact on rainfall dipole in Eastern China

Eastern China summer rainfall features a complex spatiotemporal variability. We performed an EOF analysis to extract the dominant mode of the interannual variation of summer (June-August) rainfall in eastern China for the period 1967-2018. 1

(b), 10 • E-45 • E, 45 • N-65 • N)
, a region including eastern Europe and the area around the Baltic Sea, where the most significant positive correlation occurs. We found that western Eurasia spring snow cover (WESS) variability, defined as the time series of spring SCF averaged over this region, could explain a significant fraction of the NSDR variance (R WESS-NSDR = 0.44, P < 0.01, figure 1(c)) during the period 1967-2018. Specifically, in years with below-normal spring snow cover over western Eurasia, northern China is anomalously dry but southern China is anomalously wet during the summer, and vice versa.

To examine the temporal stability of the WESS-NSDR correlation, we further analyzed the decadal change in R WESS-NSDR calculated with a 19 year moving window. The results show that R WESS-NSDR is small and not statistically significant until after the late 1990s (figure 2(a)). Thus, we suggest a strengthened impact of WESS on NSDR since the late 1990s. This result is robust to the choice of different lengths for the moving window (15, 17, 21 and 23 years) (figure S3). To facilitate further inspection, we selected two separate 19 year epochs: 1973-1991 and 2000-2018, which represent the weakest (R WESS-NSDR = 0.11, P > 0.1) and the strongest (R WESS-NSDR = 0.69, P < 0.01) periods of the relationship between WESS and NSDR, respectively. The shift in the snowrainfall relationship is also clearly illustrated in the spatial pattern, with an evident north-south dipole structure of summer rainfall anomalies over eastern China in the later epoch (figures 2(c) and S4(b)),

Dynamic interpretation of the shift in the snow-rainfall relationship

Numerous observational and modeling studies suggested that the delayed snow hydrological effect acts as a crucial role in linking the spring snow cover and summer atmospheric circulations through regulating the soil moisture (e.g. [START_REF] Xu | Snow-atmosphere coupling strength. Part II: albedo effect versus hydrological effect[END_REF], Halder and Dirmeyer 2017, Zhang et al 2017a[START_REF] Zhang | The weakening relationship between Eurasian spring snow cover and Indian summer monsoon rainfall[END_REF]. The anomaly of soil moisture induced by the spring snowmelt could persist into the following season, providing a lagged influence on the summer atmospheric circulations through its effects on the surface energy balance and the resultant atmospheric thermal condition. Here, we show that following the positive anomaly of spring snow cover over the western Eurasia, the significantly increased soil moisture are observed over this region during both spring and subsequent summer (figures S5(a) and (b)). The wetter soil consumes more energy to evaporate the excessive moisture and results in less energy (figures S5(c) and (d)) being used to heat the atmosphere, thereby leading to the decrease of air temperature at surface and mid troposphere (figures S5(e)-(h)), which then would induce the anomalous atmospheric circulations.

We interpreted the emerging influence of WESS on NSDR from above perspective, through statistical regression analyses of large-scale circulation and moisture fields onto WESS in the two epochs (figure 3). Following positive anomalies of WESS in the later epoch, there is an obvious wave train pattern with alternate positive and negative geopotential height anomalies at 500 hPa over the northern mid-high latitudes (figure 3(b)). The WAFs show this wave train extending eastwards from eastern Europe to northern Siberia, then southeastwards to the Lake Baikal region and to the Yangtze River valley (figure 3(b)), indicating that the snow cover variability signal over western Eurasia can propagate to downstream areas through this teleconnection. As a result, at 850 hPa, there is a prominent anomalous anticyclone centered over the southeast coast of China and an anomalous cyclone to the east of Lake Baikal (figure 3(d)), indicating a quasi-barotropic structure. On the one hand, the anomalous anticyclone over southern China leads to an anomalous descending motion (figure 3(d)) and moisture divergence (figure 3(f)), which is responsible for the lower-than-normal rainfall to the south of Yangtze River (figure 2(c)). On the other hand, the anomalous southwesterly flow, on the northwest flank of the anticyclone, transports abundant moisture from lower latitudes to north and northeast China, which then converges with the anomalous westerly flow on the southeastern flank of the cyclone centered to the east of Lake Baikal (figures 3(d) and (f)). Such convergence results in an anomalous ascending motion and moisture convergence that ultimately contributes to the higher-than-normal rainfall over north and northeast China (figure 2(b)).

In contrast, in the earlier epoch with nonsignificant R WESS-NSDR , there is no clear wave train pattern in response to the WESS anomalies (figure 3(a)). The signal of snow cover anomalies over western Eurasia is not able to propagate to downstream areas, resulting in no significant anomalies of circulations and rainfall over eastern China (figures 3(c) and (e)). We suggest that the shift in the snow-rainfall relationship may result from changes in WESS-related anomalous circulations.

Mechanisms responsible for the shift in the snow-rainfall relationship

A crucial issue that needs to be addressed is what contributed to the changes in WESS-related circulation anomalies before and after the late 1990s. We found that the anomalous atmospheric circulations responsible for the formation of NSDR, exhibit a wave-train like pattern originating from the North Atlantic, where the dipole structure of geopotential height anomalies closely resembles the spatial pattern of the summer NAO (SNAO) (Hurrell et al 2003, Folland et al 2009) (figure S6). The SNAO could exert a notable impact on the eastern Asian climate through a wave train teleconnection over Eurasia (Linderholm et al 2011, Sun andWang 2012).

Here, we find that the correlation coefficient between SNAO and NSDR is 0.47 (significant at the 99% confidence level), suggesting a possible contribution of SNAO to the formation of NSDR during the period 1967-2018. Interestingly, a partial correlation analysis indicates that R WESS-NSDR becomes statistically that, if SNAO occurs in-phase with WESS, it will amplify the impact of WESS on NSDR, while, if out-of-phase, it will counteract the impact. We confirmed this hypothesis by analyzing the relationship between WESS and NSDR in two groups of years when SNAO occurred in-phase and out-of-phase with WESS. After applying a 9 year high-pass filter to the original WESS and SNAO time series, the in-phase (out-of-the-phase) group was defined to include the years having the same (opposite) sign of filtered WESS and SNAO (table S1). For the in-phase group, R WESS-NSDR can reach as high as 0.71 (significant at the 99% confidence level) (figure 4(a)). In this group, following the positive WESS anomalies, there is a wave train propagating from the northern Siberia to the vicinity of Lake Baikal and then to the Yangtze River basin (figure S8(a)). Correspondingly, an obvious anomalous anticyclone over eastern China and a cyclone over the Lake Baikal region are observed (figure 4(e)), which then form a pronounced meridional dipole structure of rainfall anomalies over eastern China (figures 4(c) and S4(c)). For the out-ofphase group, the R WESS-NSDR is only 0.17 (P = 0.40) (figure 4(b)), and there is no clear rainfall anomaly pattern (figures 4(d) and S4(d)) because of weak circulation anomalies (figures 4(f) and S8(b)).

Using a bootstrap sampling method, we estimated that, in the earlier epoch, 77% of years have the outof-phase relationship between the WESS and SNAO, with the remaining 23% having the in-phase relationship (figure S9). For the later epoch, the relevant values were 20% for out-of-phase and 80% for inphase. Therefore, in the later epoch, the WESS-related anomalous circulations would be largely reinforced by the constructive superposition of the frequently occurring same-sign circulation anomalies associated with SNAO. A strengthened anomalous circulation is then formed, which favors a pronounced meridional dipole structure of rainfall anomalies over eastern China, leading to a statistically significant snowrainfall relationship. In contrast, in the earlier epoch, the WESS associated circulation anomalies would be largely counteracted by the destructive superposition of the opposite-sign anomalous circulation induced by SNAO. The WESS-related dipole pattern of rainfall anomalies would then be cancelled out by the opposite sign dipole mode caused by SNAO, thus leading to a weak snow-rainfall relationship. Moreover, an empirical statistical model indicates that the prediction skill would substantially improve when both the WESS and SNAO indices are used to predict the NSDR time series in comparison with that using only the WESS index (figure S10), further suggesting the important role of the SNAO in the connection between WESS and NSDR.

Discussion and conclusion

We also explored other possible mechanisms that might contribute to the shift in R WESS-NSDR . First, the ENSO, an important source of inter-annual climate variability over east Asia (Wang et al 2000[START_REF] Xie | Indian Ocean capacitor effect on Indo-western Pacific climate during the summer following El Niño[END_REF], Zhang et al 2017b). To examine the possible effect of ENSO on R WESS-NSDR , we conducted a partial correlation analysis between WESS and NSDR, whilst statistically controlling for the effect of ENSO. Our results show that the observed shift in R WESS-NSDR still persists after removing the ENSO signal (table S2). Similar results were found when the IOD, EAP teleconnection pattern, and AO were considered (table S2). Second, the Tibetan Plateau provides an elevated heat source for the atmosphere, and its thermal forcing plays a crucial role in summer precipitation variation over eastern China (Wu et al 2007(Wu et al , 2015)).

There is evidence to show that recent reductions in the intensity of the Tibetan Plateau heat source could contribute to the 'South Flood-North Drought' pattern over eastern China (Duan andWu 2008, Duan et al 2013). However, our partial correlation analysis shows that the shift in R WESS-NSDR is still observed after excluding the effect of spring thermal forcing (table S2). Third, one particularly conspicuous land surface change in eastern China over the past three decades is human land-use management induced vegetation greening (Chen et al 2019). Such greening could affect local and downwind precipitation through enhanced evapotranspiration and changes in atmospheric circulation (Piao et al 2020). The possible impact of greening on the snow-precipitation relationship requires investigation, particularly as a prominent greening trend has emerged over South China in the last two decades (Chen et al 2019). We used the first two leading EOFs of LAI over eastern China to represent different modes of the vegetation greening pattern. After removing the covarying effect of each mode with rainfall through a partial correlation analysis, changes in the vegetation greening pattern did not explain the observed shift in R WESS-NSDR (table S2).

In addition, some inter-decadal oscillation of SST anomalies such as Pacific Decadal Oscillation and Atlantic Multidecadal Oscillation experienced a phase transition around the 1990s and contributed to the decadal change in the eastern China summer rainfall (Zhu et al 2015, Si andDing 2016). However, the relatively short time period of observations hinders any attempt to distinguish the possible impact of such multi-decadal variations on the shift in the snowrainfall relationship. Although these issues could be principally addressed through perturbed or sensitivity simulations, a better representation of snowatmosphere coupling is required than that in the current generation of coupled models which are fraught with snow-related deficiencies.

In summary, we have combined observational and reanalysis data sets to reveal a strengthened connection between western Eurasian spring snow cover (WESS) and the meridional dipole structure of the summer rainfall anomalies over eastern China since the late 1990s. Such an amplified snow-rainfall relationship, along with the decline in Eurasian spring snow cover, constitutes a possible reason behind the 'South Flood-North Drought' pattern observed in eastern China since the 1990s. Our finding is exactly the opposite of a previously reported result on the Indian summer monsoon, which showed a weakened snow-rainfall relationship since the 1990s (Peings and Douville 2010, Wu et al 2014[START_REF] Zhang | The weakening relationship between Eurasian spring snow cover and Indian summer monsoon rainfall[END_REF]. We thus suggest a different paradigm for understanding the impact of Eurasian spring snow cover on the Asian summer monsoon, with south and eastern Asian subsystems.

Figure 1 .

 1 Figure 1. Relationship between Eurasian spring snow cover and eastern China summer rainfall during 1967-2018. (a) The first EOF mode of eastern China summer (June-August) rainfall, which shows an NSDR. (b) Spatial distribution of the correlation coefficients (R) of Eurasian spring (March-May) SCF with the NSDR (represented by the PC1 of the EOF1 mode). Hatched areas indicate where R is significant at the 95% confidence level. The black box denotes western Eurasia (10 • E-45 • E, 45 • N-65 • N). (c) Time series of the normalized WESS and NSDR.

Figure 2 .

 2 Figure 2. Shift in the relationship between WESS and eastern China summer rainfall. (a) 19 year sliding correlation (R) between NSDR and WESS. The gray dashed line denotes the 95% confidence level. (b) and (c) Show the anomalies of eastern China summer rainfall obtained by regression onto WESS for 1973-1991 and 2000-2018, respectively. Significant rainfall anomalies at the 95% confidence level are hatched.

Figure 3 .

 3 Figure 3. WESS-related circulation anomalies for 1973-1991 (left panel) and 2000-2018 (right panel). (a) and (b) Anomalies of 500 hPa geopotential height (shading) and the associated WAF (vectors; m 2 s -2 ) in summer (June-August) obtained by regression against the WESS. Geopotential height anomalies that reached the 95% confidence level are outlined by dashes. WAFs with magnitude less than 0.05 m 2 s -2 are not shown. (c) and (d) Same as (a) and (b) except for the 850 hPa wind (vectors; m s -1 ) and 500 hPa vertical velocity (shading). (e) and (f) Same as (a) and (b) but for vertical integrated water vapor flux (vectors; kg m -1 s -1 ) and its divergence (shading). Green dots and black vectors indicate values significant at the 95% confidence level.

Figure 4 .

 4 Figure 4. Snow-rainfall relationship in the years when WESS and SNAO are in-phase (left panel) and out-of-phase (right panel). (a) and (b) Scatterplots of WESS with NSDR. Two asterisks indicate that R is significant at the 99% confidence level. (c) and (d) Summer rainfall anomalies obtained by regression against the WESS. Significant rainfall anomalies at the 95% confidence level are hatched. (e) and (f) Anomalies of 850 hPa wind (vectors; m s -1 ) and 500 hPa vertical velocity (shading) in summer obtained by regression onto WESS. Green dots and black vectors indicate values significant at the 95% confidence level.

  • × 0.5 • and is available for 1901 to the present (Harris et al 2014). The monthly atmospheric variables, such as geopotential height, horizontal wind, vertical velocity, specific humidity, and sensible heat flux were obtained from the National Centers for Environmental Prediction-

	National Center for Atmospheric Research (NCEP-
	NCAR) reanalysis 1 (Kalnay et al 1996) for 1948 to the
	present, at a horizontal resolution of 2.5 • × 2.5 • . The leaf area index (LAI) data derived from the LAI3gV1
	(Zhu et al 2013), which has a 1/12 • spatial resolution
	and spans the period 1982-2018.
	The monthly indices of Arctic Oscillation (AO)
	and El Niño Southern Oscillation (ENSO, Niño 3.4
	region) were provided by the US National Oceanic
	and Atmospheric Administration Climate Predic-
	tion Center. The indices of El Niño Modoki and
	Indian Ocean Dipole (IOD) were obtained from the
	Japan Agency for Marine Science and Technology.
	The North Atlantic Oscillation (NAO) index was
	defined as the difference in normalized zonally aver-
	aged sea level pressure between 35 • N and 65 • N
	over the region ranging from 80 • W to 30 • E (Li and
	Wang

  Positive correlations occur over the western Eurasia, central Siberia, and the Tibetan Plateau, while sparse patches of negative correlations are found on the north coast of eastern Siberia. Next, we considered spring snow cover only over western Eurasia (the black box in figure

	over southern China. The corresponding principal
	component of EOF1 is used as an index for the north-
	south dipole structure of rainfall anomalies (NSDR)
	to elucidate the snow-rainfall relationship.
	Figure 1(b) shows the spatial distribution of the
	correlation coefficients between NSDR and gridded
	spring SCF over Eurasia during the period 1967-
	2018.
	The first
	EOF mode (EOF1), which explains 19.3% of the total
	variances, characterizes the south-north dipole struc-
	ture, with rainfall anomalies of opposite sign over the

south of the Yangtze River and in northern and northeastern China (figures 1(a) and S1 (available online at stacks.iop.org/ERL/16/054012/mmedia)). Similar results are also found when using the station rainfall dataset (figure

S2

). The positive phase of EOF1 represents greater than normal rainfall over northern and northeastern China and lower than normal rainfall

© 2021 The Author(s). Published by IOP Publishing Ltd

Acknowledgments

This study was supported by Key Research and Development Programs for Global Change and Adaptation (2017YFA0603604 and Y7Jk011001), Preliminary research on Three Poles Environment and Climate Change (2019YFC1509103), and the National Natural Science Foundation of China (41871104, 41530528).

Data availability statement

All data that support the findings of this study are included within the article (and any supplementary files).

Bamzai A S and Shukla J 1999 Relation between Eurasian snow cover, snow depth and the Indian summer monsoon: an observational study J. Clim. 12 3117-32 Barnett T P, Dumenil L, Schlese U, Roekler E and Latif M 1989 The effect of Eurasian snow cover on regional and global climate variations J. Atmos


Mathieu Mory

Nicolas Gratiot

A Manning

Hervé Michallet

M Mory8

N Gratioth

A J Manningc

H Michalleth

CBS layers in a diffusive turbulence grid oscillation experiment

Keywords: Fluid mud, turbulence, fl oculation, cohesive sediment suspensions

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

Special attention has recently focused on the existence of Concentrated Benthic Suspensions (CBS) of cohesive sediments in numerous natural sites. The concentration is large in these suspensions which are subjected to a high level of turbulence. The concept of a saturation concentration for cohesive sediment was introduced in 1996 by Winterwerp (see [START_REF] Winterwerp | On the dynamics of high-concentrated mud suspensions[END_REF]. For a turbulent shear fl ow, the latter obtained a relationship between the thickness ofa CBS layer, the concentration, and the hydrodynamic properties of the fl ow. The concept of saturation is based on the principle (Turner, 1973) that a turbulent shear flow collapses when the flux Richardson number Ri 1 reaches a value of about 0.15.

We investigate in this paper the concept of saturation concentration in an oscillating grid experiment. Natural mud sediment was introduced in a tank where turbulence was produced by an oscillating grid. Grid tank experiments were used before to study the resuspension of cohesive sediments by turbulence [START_REF] Tsai | A portable device for measuring sediment resuspension[END_REF]Wolanski et al ., 1989). In a grid tank experiment similar to ours, E and [START_REF] Hopfinger | Stratification by solid particle suspensions[END_REF] observed that the layer containing non cohesive particles resuspended by turbulence reached an equilibrium thickness, an observation which is in a qualitative agreement with the concept of saturation. Although oscillating grid experiments are far from producing flow conditions achieved in the field, they allow to study turbulence conditions with no mean flow. Turbulence is transported by turbulent diffusion. The isolated contribution of turbulence to mixing can be studied in oscillating grid devices.

A specific feature of the present study is that experiments were performed in a grid tank facility using fl uid mud mixtures, implying that the settling of sediments is modified by the flocculation of particles. Because the sediment concentration is sufficiently high, it is assumed in this study that the effect of particles on turbulence is solely related to density-stratification. The focus of the paper is on the interaction between the turbulence and density-stratification. By performing settling velocity measurements, the effect of turbulence on the settling velocity is however accounted for. Turbulence properties, sediment concentration and settling velocities were all measured in the laboratory. The flux Richardson number below the interface was determined from these measurements in order to investigate a relationship between the thickness of the CBS layer and the sediment and turbulence properties. The paper is divided as follows. A review of the effect of stratification on turbulence is presented in the next section. The fl ux Richardson number is defined and its relation to the Rouse number and the Richardson number as usually defined in oscillating grid turbulence is given. The different measurements are presented in sections 3 to 5. The results are used in section 6 to determine the flux Richardson number at the interface.

THE EFFECT OF STRATIFICATION ON GRID TURBULENCE IN SEDIMENT LADEN FLOW

The evolution of the turbulent kinetic energy in a sediment laden flow is subjected to four mechanisms: the production of turbulent kinetic energy by energy exchange with the mean flow, the diffusion, the dissipation, and gravity effects. The balance equation for the turbulent kinetic energy k is often written in a symbolic form (see [START_REF] Winterwerp | On the dynamics of high-concentrated mud suspensions[END_REF],

dk -=P+D-s-B dt (1)
in order to distinguish the four processes. B

' 1 =-P+D (3) (4) (5)
quantifies the relative importance of gravity in the balance equation of the turbulent kinetic energy in steady conditions. For turbulent shear flows the dif fusion term is often omitted [START_REF] Villaret | Effect of stratif i cation by suspended sediments on turbulent shear fl ows[END_REF][START_REF] Winterwerp | On the dynamics of high-concentrated mud suspensions[END_REF] because the production of turbulence by the mean flow gradient is usually the dominant term, except in the vicinity of the free surface.

In the absence of mean flow, such as in grid turbulence, the production term vanishes and the diffusion of turbulence is the only process by which the turbulent kinetic energy is transported. For a steady turbulence with no mean shear fl ow, the flux Richardson is simply rewritten as

R . B 1 1 =-B+& (6)
Equation ( 6) is more convenient than (5) since the dissipation term is simply expressed in terms of the turbulent rms velocity u and turbulent lengthscale £ u 3 c�A £

where A is a constant of order 0(1). Using equation (4), we then deduce

(7) (8)
The flux Richardson number is a local quantity which varies over the vertical as a function of the turbulent rms velocity u(z), the turbulent lengthscale Jl(z) and the mean concentration C(z). For the conditions considered in this paper, our experimental observations will show that the concentration in the CBS layer is uniform (see section 4). We therefore writeC(z) = l'J.C in the CBS layer, /'J.C being the sediment concentration jump across the interface. The flux Richardson number Ri 1 (z) is therefore expressed as a function of the Richardson number Ri( z) , as usually defined for grid turbulence experiments [START_REF] Hopfinger | Spatially decaying turbulence and its relation to mixing across density interfaces[END_REF],

Ri(z) = tg l'J.
� Jl(z) P wU (z) and the Rouse number

w Rou(z) = _s_ u(z) (9) (10)
Considering steady conditions, the integration of equation (1) over a layer having its lower boundary h 0 below the lutocline and its upper boundary h1 above the lutocline gives

F(z)=A' u 3 (z), ( 12
)
where A' is a constant of order 0(1). Because turbulence has vanished above the lutocline, the fl ux of turbulent kinetic energy at h 0 balances in eq. ( 11) the rate of increase of potential energy and the total dissipation within the layer. Equation (11) is of little use for determining the position of the lutocline unless the vertical variation of the dissipation rate is known. However it explains the existence of a maximum position of the lutocline for given conditions since the buoyancy term in eq. (11) increases linearly with the layer thickness. By comparing in eq. (11) the buoyancy term to the fl ux of turbulent kinetic energy, we observe that the parameter driving the lutocline occurrence is the product of the Richardson number by the Rouse number Ri(z) Rou(z) as already anticipated in the expression of the flux Richardson number (eq. 8). When Ri(z).Rou(z) << 1, gravity effects are negligible and the evolution of u and £ with z is not modified by the presence of sediments. Gravity effects will be active when Ri(z).Rou(z) reaches a value of order 0(1) .

Different authors have considered before the significance of the flux Richardson number. Turner (1973) and [START_REF] Ivey | On the nature of turbulence in stratified fluids; Part I : the energetics of mixing[END_REF] have shown that turbulence collapses when the flux Richardson number Ri 1 reaches a value of the order of 0.15 to 0.25. This experimental evidence was taken by [START_REF] Winterwerp | On the dynamics of high-concentrated mud suspensions[END_REF] as a basis for establishing the concept of a saturation concentration for cohesive sediments. The critical value of Ri f in the range 0.15 to 0.25 was obtained by Ivey and Imberger (1999) in view of experimental data for salt-stratified or temperature-stratified shear flows. The applicability of this result to fluid mud mixtures in a grid tank turbulence is questionable. On the one hand, the settling velocity of particles tends asymptotically to zero for salt stratification. In oscillating grid tank experiments, no steady state is obtained and the interface is continuously displaced upward for this type of stratification. On the other hand, diffusing turbulence and shear produced turbulence are different. To our knowledge, estimates of Ri f in purely diffusive flows are still lacking.

Previous experiments [START_REF] Hopfinger | Stratification by solid particle suspensions[END_REF][START_REF] Huppert | Sedimentation and entrainment in dense layers of suspended particles stirred by an oscillating grid[END_REF] show the existence of a saturation concentration for sediments maintained in suspension by a diffusing turbulence. The above considerations provide a physical understanding of the process, but the quantitative determination of the position of the lutocline is not available*. In our laboratory experiments, measurements of the vertical variations of the turbulent rms velocity, the integral lengthscale, the mean concentration and the settling velocity have allowed to estimate the value of Ri 1 at the lutocline for various CBS experimental conditions.

THE EXPERIMENTAL SET-UP, PROCEDURES AND OBSERVATIONS

The experiment was performed in a square tank, 53 cm wide. For all experiments, the water depth was 40 cm. The horizontal grid, made of 1.5 cm square bars with a mesh M = 7.5cm, oscillated vertically. The mean position of the grid was 5 cm above the bottom of the tank. The grid oscillation conditions were determined by the stroke S = 4.5cm (twice the amplitude) and the frequency F of oscillation. Three different values F = 3Hz , F = 4Hz and F = 6Hz were considered in the study. The tank and the grid geometry were similar to these used previously by [START_REF] Hopfinger | Spatially decaying turbulence and its relation to mixing across density interfaces[END_REF]. The latter authors have shown that these experimental conditions ensure that there is no mean flow at distances from the grid larger than 2M . The turbulence is transported away from the grid by a diffusive process. The turbulence properties are IDV, i.e. they vary only with the distance z from the grid. The position z=O is here the mean position of the grid and the Oz axis is oriented upward. Figure 1 is a sketch of the experimental set-up. [START_REF] Hopfinger | Stratification by solid particle suspensions[END_REF][START_REF] Huppert | Sedimentation and entrainment in dense layers of suspended particles stirred by an oscillating grid[END_REF] have determined from an equation similar to (11) a relationship between the layer thickness, the sediment concentration and the turbulence properties, but they assume that turbulence is not modified by gravity effects.

• E &

ADV-------� OBS--�

--•------...

;

- A lutocline appeared that moved downward until it became stable. Oscillations of the lutocline were observed but the mean level of the lutocline did not appear to change anymore. When the amount of sediment introduced inside the tank exceeded a certain level, the turbulence produced in the CBS layer was also not sufficient to maintain the whole amount of sediment in suspension. Sediment settled on the bottom to form a deposited bed. The equilibrium state consisted of three superposed layers: the deposited bed, the CBS layer and the clear water layer. This is called the saturated regime. When only a limited quantity of sediments was introduced in the tank a quasi-uniform suspension with no lutocline could be observed, and eventually no deposited bed was seen. This is the non saturated regime.

)I -- �-==-=bottle s amp ling ��- syringe 1' � m = L5c: �rrr. r --• '-------] - M=7 .5cm
Experiments were also carried out by starting the grid oscillation with the sediment initially deposited on the bottom and left to consolidate. For the same grid oscillation conditions, a smaller amount of sediment was maintained in suspension by turbulence when the sediment was initially on the bottom of the tank than when the sediment was introduced in suspension while the grid oscillated. We attribute this observation to consolidation effects in the deposited bed.

The tank was made of transparent Perpex walls which enabled the visualisation of the different sediment layers. Figure 2 shows a typical picture of the three layers stratification inside the tank taken in saturated steady conditions. Internal wave oscillations are visible at the lutocline. Three different natural mud have been used from the Gironde estuary (organic and inorganic) and from the Tamar estuary (organic). The inorganic mud from Gironde was chemically treated with potassium permanganate and passed through a 100 µm sieve. Sediment analysis performed by de [START_REF] De Croutte | Flume measurements of mud processes over a flat bottom under steady and unsteady currents[END_REF] revealed that the sediment's median diameter was about D50=12 µm and the sand content (63-100 µm) was below 3%. The clay fraction was comprised of illite (35%), kaolinite (27%), smectite (24%) and chlorite (14%).

CONCENTRATION AND TURBULENCE MEASUREMENTS FOR SATURATED CONDITIONS

Vertical profiles of the mean sediment concentration C(z) and of the rms turbulent velocity u(z) were measured for different saturated conditions. The integral lengthscale £(z) was also determined from velocity measurements for selected cases. Saturated conditions were obtained when a sufficient amount of sediment C was initially introduced in the tank. Table 1 summarises the different experimental conditions. Only the conditions that vary are given. Table 1 also Vertical profiles of the mean sediment concentration were measured from samples analysis and using an Optical Backscatter System (OBS) developed at the Institute of Marine Studies [START_REF] Christie | Optical backscatter turbidity sensor[END_REF]. The calibration curve was determined for each mud. The voltage output of the OBS increases monotonically with increasing concentration when the concentration is below 40 g/l while it decreases with increasing concentration above 45 g/I. The OBS was therefore not used for concentrations in the range 35 -50 g/l. Sample analysis were made to determine the appropriate part of the calibration curve.

Time series measurements of the velocity were made using an Acoustic Doppler Velocimeter (ADV) developed at Laboratoire des Ecoulements Geophysiques et Industriels. This system determines the velocity from the changes in the time shift of pulse-to-pulse back scattered signals. The principle of operation is not new in itself, but a specific study (Gratiot et al., 2000) was conducted to check its ability to measure the turbulent velocity in concentrated fluid mud mixtures. The adjustment of parameters of the ADV enabled to measure instantaneous velocities with a data rate of the order of 70 data/s for a concentration of 50 g/I. The data rate decreases to about 40 data/s for a concentration of 100 g/I. The rms turbulent velocity u was computed from time series of about 4000 measurements. The ADV system could measure simultaneously the velocity component parallel to the beam axis at 16 positions. The distance between two successive measurement volumes was 3 mm. The integral lengthscale £ was deduced from the computation of the correlation function of the velocity measurements made simultaneously at the different locations. Experimental conditions corresponding to the symbols are given in Table I.

Figure 4 shows different vertical profiles of the mean concentration. The distance z from the grid is made dimensionless using the distance z1 of the lutocline from the grid, so that all profiles are centred around the lutocline position. We observe that the concentration is homogenous inside the CBS layer. The homogeneity of the concentration is apparently even better for the highest concentration ( CcBs = 202g I I); the interface thickness is smaller. This is partly because internal wave oscillations of the lutocline have a larger amplitude for smaller values of CcBs , contributing to make the interface look thicker in Figure 4 when CcBs is small.

The decay of the rms turbulent velocity u with increasing distance from the grid is shown in Figure 5 for different saturated experimental conditions. The concentration in the CBS layer considered for the plots in Figure 5 ranges between 3.8 g/l and 202 g/l. Also superimposed are two vertical profiles of the turbulent rms velocity made in clear water. It is observed that the decay of u is similar for all concentrations in the CBS layer considered. The uncertainty in measured data is rather large, as indicated by error bars, but no general trend indicating that u decreases for increasing values of CcBs is observed. Actually the measurements made for the largest values of CcBs collapse on the decay plots obtained in clear water, and the highest difference is between measurements in clear water and for CCBs =51g/ l. For analysis ofthe data (section 6) the decay laws were interpolated by the power law function ( )-14

u(z) =0.092 _!__ . FM M ( 13
)
in the range 1.2 .$ z IM .$ 3.2 (M is the grid mesh) and for z<z1• Although the exponent of the power is of the order of those obtained by different authors for oscillating grid turbulence, no particular attention should be paid to the value of the exponent, considering the narrow range of values of z IM . Because of the lutocline, turbulent velocity measurement could not be made at large distances from the grid. Measurements of the integral lengthscale evolution with the distance from the grid are less common in the literature than measurements of the rms turbulent velocity, but there is a general agreement that the integral lengthscale varies linearly with the distance from the grid. This experimental result is supported by k-& modelling [START_REF] Sonin | Calibration of the k-E turbulent model for the diffusion of turbulence[END_REF] in clear water, when the decay of the turbulent kinetic energy follows a power law, anyhow the value of the exponent is. Measurements of the integral lengthscale were made for selected cases. The results are shown in Figure 6a together with the results obtained by [START_REF] Hopfinger | Spatially decaying turbulence and its relation to mixing across density interfaces[END_REF]. A correlation function R is also shown in Figure 6b. A number of correlation functions that were computed did not display such a nice and smooth behaviour as the case shown in Figure 6b. This explains the large scatter in Figure 6a. Nevertheless, the data indicate a clear increase of the integral lengthscale with increasing distance from the grid. Our measurements are not inconsistent with the results of Hopfinger and Toly. In view of this, the evolution of the integral lengthscale with the distance from the grid was taken for subsequent use of the data as

f(z) = 0.2z (14) 5. SETTLING VELOCITY MEASUREMENTS
The settling velocity of cohesive sediments varies very much with the concentration. When the concentration in the CBS layer was below 15 g/l, the settling velocity was determined using the INSSEV system video camera [START_REF] Fennessy | INSSEV: an instrument to measure the size and settling velocity of fl oes in situ[END_REF]. This apparatus was developed at the Institute of Marine Studies for in situ measurements of the settling velocity of individual floes. The size and the settling velocity of floes are measured from the analysis of images of settling floes in a column recorded by a video camera [START_REF] Manning | INSSEV (In Situ Settling Velocity instrument) -1.3: Operator Manual[END_REF]. Whereas the system is immersed in the water column for in situ measurements, samples were taken from the grid tank using a syringe of length 450 mm and diameter 8 mm, and introduced in the settling column for image recording. By comparing the floe properties measured in situ and measured after sampling, [START_REF] Manning | A study of the effect of turbulence on the properties of flocculated mud[END_REF] concluded that the procedure of sampling did not modify the floes. Figure 7 shows a typical diagram relating the floe diameter D1, to the settling velocity w •. The diameter D1 of the floes measured by INSSEV is the square root of the length normal to settling direction times that measured in the settling direction. An essential feature is the dispersion in floe diameter and in settling velocity, indicating that the fl oe density is as well subjected to a significant variability. For further analysis, the settling velocity of each fl uid mud suspension was determined by computing the mass flux of the suspension using the settling velocities measured for individual floes. Error bars in the estimates were determined from the scatter in measurements.

For concentrations in the CBS layer larger than 15 g/l, the sedimentation of the suspension inside a column was recorded and the settling velocity was determined as a function of the concentration using the Kynch method (1952). [START_REF] Toorman | Sedimentation of self weight consolidation: constitutive equations and numerical modelling[END_REF] recently reported the relevance of this method applied to highly concentrated suspensions of cohesive particles. Figure 8 shows the variation of the settling velocity with the concentration determined for the different mud considered in the study. Results obtained using INSSEV appear using symbols, whereas results obtained :from sedimentation experiments are displayed as curves. Because INSSEV was available only during a limited period of the experimental study, the settling velocity has not been measured for all fl uid mud suspensions concentration below 15 g/l. The variation of the settling velocity with the concentration in the CBS layer was interpolated :from Figure 8. In spite of the large scatter in settling velocity estimates, Figure 8 clearly shows that the settling velocity decreases by two orders of magnitude when the concentration increases from 4 g/l to 200 g/l. The techniques used to measure the settling velocity in low concentrated suspensions (INSSEV measures the settling velocity of individual floes) and in high concentrated suspensions (analysis of sedimentation experiments) are very different in their principles. It is however expected that the settling velocity decreases significantly when the sediment concentration is larger than 15 g/l, because of hindered settling effects.

THE FLUX RICHARDSON NUMBER AT THE INTERFACE

For all experiments listed in Table 1, the rms turbulent velocity u and the integral lengthscale R were determined at the lutocline from eqs. (13) and ( 14), respectively. Their values, as well as the mud concentration Cc s s and the settling velocity W8, are given in Table 2. These data were used to estimate the fl ux Richardson number at the lutocline ( eq. 8). The variations of Ri1 with the concentration in the CBS layer are shown in Figure 9. The error bars result mainly from the uncertainties in the measurements of the settling velocity and of the rms turbulent velocity u (the dissipation rate varies as u 3 ). The value of A= 1 for estimating the turbulent kinetic energy dissipation rate (eq. 7) was used for the plot in Figure 9. Similar curves were obtained for A= 0.5 and A= 2. They are not shown here, as the conclusions of the paper do not depend on this value. The major result of Figure 9 is that the fl ux Richardson number decreases from about 1 at the interface when the concentration in the CBS layer is low (a few g/l) to 0.1 when the concentration in the CBS layer increases to 200 g/l. This behaviour is due to the strong decrease in the settling velocity when the concentration in the CBS layer is larger than 50.

The estimates of the Rouse and Richardson numbers at the lutocline are included in Table 2.

The Rouse number is indeed very small for the large concentration cases, which correspond simultaneously to the highest values of the Richardson number (i.e. the strongest stratification).

In spite of the limited number of experiments carried out for large concentration conditions, their interpretation is less subjected to experimental uncertainty. For such conditions, the Richardson number is high (see Table 2). Because the amplitude of oscillations at the lutocline was small, velocity measurements could be performed close to the interface as shown in Figure 5. The rrns turbulent velocity in the CBS layer displays no signifi cant decrease as compared to measurements in clear water. Table 2 indicates that the turbulent Reynolds number based on u and f, is of the order of 200 for all conditions. The estimates of turbulence properties at the lutocline are based on measurements made close to the interface, enabling a confi dent estimate of the flux Richardson number at the interface. On the other hand, turbulence measurements could not be made close to the interface for the low concentration cases (less than 12 g/I) because of the oscillations of the lutocline. The estimates of u and f at the lutocline given in Table 2 for the low concentration cases were extrapolated from eqs. (13) and (14). Since the product of the Rouse and the Richardson number is high, the fl ux Richardson number is about one at the interface. We cannot exclude that a smaller value would be obtained using measurements closer to the lutocline.

Nevertheless, using measured values of turbulence, we estimate that Rou.Ri == 1 at a distance from the grid ranging between 11.5 cm and 13.8 cm for the low concentration range (3.8g//5'Ccss :::;1I.6g//) and at a distance of 14.5 cm for Ccss ==5lg//. The fl ux Richardson number is therefore larger than 0.5 at the lutocline for the low concentration cases, a value which is clearly different from that obtained for the large concentration cases. The observation of CBS layers of thickness independent of time is not new. We have considered if the saturation effect that fixes the lutocline position can be characterised in terms of the criterion involving the fl ux Richardson number. The fl ux Richardson number at the interface (Figure 9) has been computed from measurements of turbulence and concentration inside the CBS layer and estimations of the settling velocity. Although our experiments have not been successful for determining a criterion for the lutocline position, a clear difference was found between the high concentration range ( C ens > 100 g I I) for which Ri 1 "" 0.1 and the low concentration cases ( Cc B s < 50g II) for which Ri 1 > 0.5 . A 14 numerical study is currently under investigation using a lDV k-i: model in order to compute the vertical variation of Ri 1 below the lutocline for the different cases considered in this study.

A significant part of the study was the measurement of the settling velocity. A very strong decrease of the settling velocity with increasing concentration was measured in our fluid mud mixtures. The theoretical analysis of the effect of stratification in a diffusing turbulence (section 2) shows that the driving parameter is the product of the Rouse number by the Richardson number, and therefore the settling velocity is indeed the most sensitive quantity when large concentration are considered.

  P and D denote the production of turbulent kinetic energy by the mean fl ow and the diffusion of turbulent kinetic energy, respectively. They are simply expressed in the framework of the k-& model of turbulence. & is the dissipation rate of the turbulent kinetic energy. The buoyancy term is (2) (Sheng and Villaret, 1989), with r = (Ps -P w )/ P s . Ps and P w are the densities of solid particles and water, respectively. g is the acceleration of gravity. B is the rate of increase of potential energy required to maintain the sediment in suspension. Particles in a sediment laden fluid have a non-zero settling velocity W8• Steady conditions are obtained when the turbulent flux of sediment is equal to the settling flux of sediment as deduced from the mass balance equation , i.e. The buoyancy term B has therefore a simple dependence with the mean concentration C: Following Ivey and Imberger ( 1991 ), the flux Richardson number R .

F

  (z = h 0 ) = wctgl'J.C(h1 -h 0 )+ jc(z).dz ho (11) F(z) denotes the fl ux of turbulent kinetic energy at the position z. The flux of turbulent kinetic energy is written

Figure 1 .

 1 Figure 1. Experimental set-up.An homogenous premixed fluid mud mixture was introduced in the tank while the grid oscillated. The vertical distribution of sediments was redistributed during a transient evolution in which settling of sediments and turbulent transport of sediments both occurred. Settling occurred in the upper part of the tank where the turbulence level was not sufficient to maintain the sediment in suspension. A lutocline appeared that moved downward until it became stable. Oscillations of the lutocline were observed but the mean level of the lutocline did not appear to change anymore. When the amount of sediment introduced inside the tank

Figure 3

 3 Figure3illustrates the process that leads from non-saturated to saturated conditions. Vertical profiles of concentration, in steady conditions, are shown for four different initial concentrations C. C is defined as the mass of dry sediment divided by the total volume of fluid inside the tank. When C is less than lg 11 the concentration is clearly uniform over the vertical. A small gradient of concentration is noticed in the upper part for C == 1.3g I I and the lutocline is clearly identified for C == 2. 7 g I I . A further increase of C enhances the gradient across the lutocline which is simultaneously displaced downward.

Figure 2 .Figure 3 .

 23 Figure 2. Photograph of the three layers observed for steady saturated conditions. The concentration in the CBS layer is Cc s s ""'80g I I. 35 30 25 �20 N 15 10

  contains the mean concentration Cc8s measured in the CBS layer in steady conditions and the distance z1 of the lutocline from the mean position of the grid. A wide range of concentration CCBs was investigated.

•

  Table 1 Experimental conditions of saturated concentration experiments.mud Sa (ppt) F (Hz) c (g/l) CcBS (g/l) z, (cm) symbol Gironde Sa: salinity, F : frequency of grid oscillation, C: averaged sediment concentration in the tank, Cc B s: mean sediment concentration in the CBS layer, z1: distance of the lutocline from the grid. The symbols are those used inFigures 4, 5, 6 and 9.

Figure 4 .

 4 Figure 4. Vertical profiles of concentration measured in saturated steady conditions. Experimental conditions corresponding to the symbols are given in TableI.

Figure 5 .

 5 Figure5. Vertical variations of the nns turbulent velocity in the CBS layer for saturated steady conditions. Experimental conditions corresponding to the symbols are given in Table1. The solid light line and dashed dotted line are measurements in clear water for frequencies of grid oscillation F = 6Hz and F = 3Hz, respectively. The heavy solid line is eq. (13).

Figure 6 .

 6 Figure 6. 6a Integral lengthscale variations versus the distance from the grid in saturated steady conditions. Experimental conditions corresponding to the symbols are given in Table 1. The dashed line and the dashed dotted line are the resuhs obtained by Hopfinger and Toly ( 1976) for S = 4cm while M = 1 Ocm and M = 5cm, respectively. 6b. Correlation function measured at z = l3cm for C cs s = 5.lg I/.

Figure 7 .Figure 8 .

 78 Figure 7. Variation of the settling velocity with the size of floe measured using INSSEV in a fluid mud suspension: Tamar mud. C = I g/l.

  have been generated in a grid tank experiment using fl uid mud mixtures of concentration in the range 4 g/I to 200 g/I. The concentration inside the CBS layer was nearly homogeneous. An important result of the study was that the decay of turbulence inside the CBS layer does not vary with the concentration, excepted in the vicinity of the lutocline. The decay of turbulence are quantitatively similar in a CBS layer of concentration as high as 200 g/I and in clear water.

  carried as part of the COSINUS Program, which is funded by the European Commission (contract MAS3-CT97-0082). We express our gratitude to E. Hopfinger, who made the grid tank facility available, and J.M. Barnoud for his continuous technical assistance. A close collaboration was established with the Institute of Marine Studies (University of Plymouth) regarding measurements techniques. K. Dyer is thanked for his support.

Table 2

 2 Turbulence properties at the lutocline and settling velocity for the CBS conditions given inTable 1

	Cc s s g/l z, cm u mmls R cm ws mm/s Rou Ri	202 10.5 12.9 2.1 0.007 0.016 149 9.8 14.2 2.0 0.0005 0.0011 0.014 51 17.0 13.l 3.4 0.18 159 90.1 62.6	4.9 18.0 8.1 3.6 3.3 0.41 16.8	11.6 17.0 8.8 3.4 1.5 0.17 32.0	8.4 19.0 7.5 3.8 6.8 0.91 35.4	4.7 21.5 6.3 4.3 3.2 0.51 31.7	8.1 20.0 7.0 4.0 7.4 1.1 41.5	5.1 23.5 5.6 4.7 9.8 1.8 48.2	3.8 23.5 5.6 4.7 3.3 0.59 35.9
		0.8
		0.6
		0.4 0.2 --��t�!___________________ _

0 Figure 9. Variation of the flux Richardson number at the lutocline with the concentration in the CBS layer for steady conditions. Experimental conditions corresponding to the symbols are given in Table

1

.


Mahmoud Ali

Thomas Sayet

Alain Gasser

Eric Blond

NONLINEAR THERMOMECHANICAL MODELING OF REFRACTORY MASONRY LININGS WITH DRY JOINTS

In many high temperature industries, including iron and steel industry, refractory masonry with dry joints is extensively used for lining the internal surfaces of several large-scale components such as steel ladles, rotary kilns, and furnaces. The normal operating conditions of these linings are high thermal gradients, thermal shock, slag attack, and high thermomechanical stresses. The design and optimization of these linings is still an engineering challenge due to the coupled interactions between thermal, mechanical, and chemical fields. Thus, advanced simulation techniques are necessary especially when experimental work is still a challenging task do to the harsh working conditions. Conventional simple macro models of these linings are in many cases not sufficient to predict the appropriate mechanical response. Micro modeling techniques must be resorted to, but these are known to be very expensive from the computational point of view. In this contribution, we propose a more efficient multi scale approach for the design and optimization of these linings. Nonlinear effects caused by progressive closure of dry joints and viscoplastic behavior of refractories at high temperature are accurately captured. Comparisons between experimental and numerical results of refractory masonry walls subjected to different loading conditions at room and high temperature are presented. Good agreements between the experimental and numerical results are obtained.

INTRODUCTION

In the steel-making industry, steel ladles are used as refining units and transportation vessels of liquid steel. They are composed of different refractory layers and steel construction components. The different layers include a working lining (layer in contact with liquid steel), a safety lining (also called permanent lining), insulation layers, and a steel shell. The working lining is usually made of refractory masonry with dry joints (in some cases it could be castable). While, the safety lining could be built up with refractory masonry with mortar joints (in some cases it could be castable). The refractory linings and the insulation layers insulate the steel shell from the liquid steel, reduce heat losses from the steel shell and, therefore, maintain the liquid steel temperature within the required limit. Among the different layers in the steel ladle, the working lining is the most critical layer due to being in direct contact with liquid steel. A well-designed working lining offers longer ladle's life time, better steel quality and productivity.

Finite Element (FE) based numerical models represent an efficient way to study, design and optimize large scale masonry structures with dry joints. Existing FE models could be classified as detailed micro modeling and multi scale models. Micro modeling approach provides accurate results as the detailed geometry of the bricks and the joints are separately modeled and, almost, all the information about the masonry components are taken into account. Within the framework of micro modeling, different models were developed such as those reported in [START_REF] Andreev | Compressive behaviour of dry joints in refractory ceramic masonry[END_REF][START_REF] Lourenc ¸o | Dry joint stone masonry walls subjected to in-plane combined loading[END_REF][START_REF] Nguyen | Mechanical homogenisation of masonry wall without mortar[END_REF]. They differ by the constitutive laws used to describe the behavior of dry joints and the physical models of joints. The majority of these models assumed rigid or purely elastic bricks and all material non-linearity is concentrated in the joints. The main limitations of the micro modeling approach are related to huge computational time and solution convergence problems.

Multi scale modeling approach has a fair compromise between accuracy and computational cost, because the structural problem is usually divided into two scales: a homogeneous equivalent medium studied at the macroscopic scale with homogeneous properties determined by homogenizing the microscopic stress and strain fields in a carefully selected Representative Volume Element (RVE, microscopic scale). The homogenized properties are usually obtained using mechanical homogenization technique. Some examples of previous multi scale models of refractory masonry linings are reported in [START_REF] Nguyen | Mechanical homogenisation of masonry wall without mortar[END_REF][START_REF] Ali | Computational homogenization of elasticviscoplastic refractory masonry with dry joints[END_REF][START_REF] Ali | A Multiscale Model for Numerical Modelling of Homogenized Elastic-Viscoplastic Behavior of Mortarless Refractory Masonry Structures[END_REF].

The present paper reports on large scale experimental characterization of refractory masonry with dry joints and computationally efficient multi scale models developed at University of Orleans -France. These validated models could help in predicting and understanding the complex thermomechanical behavior of refractory masonry and assist in their design and optimization.

EXPERIMENTAL CHARACTERIZATION

Several large-scale uniaxial and biaxial compression creep tests as well as biaxial relaxation tests of alumina spinel refractory masonry walls were performed. For all tests, the dimensions of the walls were 1125 x 1100 x 140 mm 3 and were built up with half and full cuboid bricks. The dimensions of the full bricks were 150× 100× 140 mm 3 (length, height, and depth) and the dimensions of the half bricks were 75× 100×140 mm 3 . The dimensional tolerances of the bricks were ± 2 mm in the pressing direction (150±2 mm) and ± 1 mm in the directions perpendicular to the pressing direction (100±1 mm and 140±1 mm). These tolerances lead to non-uniform joint thickness in the wall and limited initial contact.

A schematic of the experimental setup used to perform the biaxial compression tests of alumina spinel refractory masonry walls at room and high temperatures is shown in Figure 5.

It is composed of a) a steel reaction frame with 3× 3×1.6 m 3 size. b) Two hydraulic jacks with maximum stroke of 90 mm. c) Four steel plungers, two of them are fixed and the other two are moving. The two moving plungers are connected to the two hydraulic cylinders and equipped with lateral steel guides to ensure their alignment during load application. d) Plunger linings, built up with magnesia chromite bricks, installed inside the steel plungers. The alumina spinel walls were built in the test field and were surrounded by the moving and fixed plungers. For measuring the displacement, four alumina tubes and Linear Variable Differential Transformers (LVDTs) were used. In addition, Digital Image Correlation (DIC) was used to obtain the full displacement fields (at room temperature). In the high temperature tests, the temperatures of the cold and hot faces of the wall were measured thanks to several thermocouples. An example of refractory masonry specimen (before heating and testing) subjected to biaxial compression load at high temperature is given in Fig. 2. The homogenized elastic viscoplastic behavior of each joint pattern is obtained using FE-base nonlinear homogenization technique. A Representative Volume Element (RVE) is selected from the periodic joint pattern. Next, combinations of uniaxial and simple shear loading simulations are performed on it. Then, the average stresses and strains are calculated by averaging the local stress and strain fields in the RVE. Finally, the homogenized elastic and viscoplastic tensors are obtained. Further details on the nonlinear homogenization analysis, identification of the homogenized elastic viscoplastic behavior of the four patterns and the transition criteria from one joint pattern to another are described in detail in references [START_REF] Ali | Computational homogenization of elasticviscoplastic refractory masonry with dry joints[END_REF][START_REF] Ali | A Multiscale Model for Numerical Modelling of Homogenized Elastic-Viscoplastic Behavior of Mortarless Refractory Masonry Structures[END_REF][START_REF] Ali | Thermomechanical modelling of refractory mortarless masonry wall subjected to biaxial compression[END_REF][START_REF] Ali | Transient thermo-mechanical analysis of steel ladle refractory linings using mechanical homogenization approach[END_REF][START_REF] Ali | Nonlinear thermomechanical modelling of refractory masonry Linings[END_REF].

The homogenized elastic behavior of the four joint patterns can be described using the orthotropic forms of Hooke's law as [START_REF] Ali | Computational homogenization of elasticviscoplastic refractory masonry with dry joints[END_REF][START_REF] Ali | Transient thermo-mechanical analysis of steel ladle refractory linings using mechanical homogenization approach[END_REF]:

Σ = C e : E e (1)
Where Σ and E e are the second order macroscopic stress and elastic strain tensors, while C e is the fourth order macroscopic elastic stiffness tensor. C e of the four patterns is obtained using FE-based homogenization technique. The orthotropic homogenized viscoplastic behavior of each pattern is obtained through averaging the local stress (σ ), viscoplastic strain rate (ε vp ) and equivalent stress (σ eq ) fields over the volume of the RVE and using a localization tensor (N) according to:

Ėvp = 1 2 A Σ eq n-1 N : Σ (2)
Where A and n are creep parameters of the base material of the bricks, Σ eq is the macroscopic equivalent stress. Tensor N accounts for the orthotropy of the structure and enables using the creep parameters of the constitutive material (i.e. bridge between the micro and macro scales).

To consider the gradual increase in the effective stiffness with the gradual closure of joints and reproduce the strain stiffening behavior, the joints were reduced to an interface, at the RVE level -micro modeling (see Fig. 4), with small thickness (t i ) where the constitutive normal and shear behaviors of the joints have been implemented. The behavior law of the interface is expressed in an incremental form as:

∆σ = k(u) ∆u (3)
or

∆    σ n τ x τ y    =   k n (u) 0 0 0 k sx (u) 0 0 0 k sy (u)   ∆    u n u x u y    (4)
Where σ n is the normal stress, τ x and τ y are the transverse stress in the x and y directions, respectively (see the coordinate system in Fig. 4). k n , k sx and k sy are the joint stiffnesses in the normal (n) and shear directions (x, y), respectively. u n is the normal deformation of the joint or joint closure. u x and u y are the shear deformations of the joint.

The variations of k n with the joints closure can be obtained from joints closure tests as reported in reference [START_REF] Allaoui | Digital image correlation measurements of mortarless joint closure in refractory masonries[END_REF], while k sx and k sy are expressed as [START_REF] Lourenc ¸o | Dry joint stone masonry walls subjected to in-plane combined loading[END_REF]: The constitutive material model presented in Fig. 3 is implemented into Abaqus through a user material subroutine (UMAT) and, then, used to simulate refractory masonry walls subjected to different loading conditions at room and high temperatures.

k sx (u) = k sy (u) = k n (u n ) 2(1 + ν) (5)

RESULTS AND DISCUSSION

Biaxial loading and unloading at room temperature

In this test series, two refractory masonry walls with 1100 × 1125 × 140 mm 3 dimensions were tested at room temperature. A 6 MPa biaxial compression load/unload was applied to the directions normal to bed and head joints. The FE model of the wall is presented in Fig. 5. The interactions between the contact surfaces of the wall, the plungers and the ground were considered in the FE model using frictional contact. The boundary conditions of the ground and the two fixed rigid plates are set to fully fixed.

Comparisons between the experimental and numerical force -displacement diagrams in the directions normal to bed and head joints during loading and unloading are presented in Fig. 6. It can be seen that the present numerical model can reproduce with reasonable accuracy the orthotropic displacement stiffening mechanical behavior of the wall. The reaction force increases with the increase in the displacement due to the increase in material stiffness with gradual closure of joints. The maximum displacement in the direction normal to head joints is smaller as compared to that in the direction normal to bed joints, because the number of head joints is less than the number of bed joints.

Biaxial creep behavior

In this test series, two refractory masonry walls with 1100 × 1125 × 140 mm 3 dimensions were tested at high temperature. The walls were heated up to 1500°C, then a constant biaxial compression load of 6 MPa was applied in the directions normal to bed and head joints. The temperature variations (with time) of the hot face (HF) and the cold face (CF) of the wall during heating and testing are given in Fig. 7. The five solid black lines represent the temperatures measured by the thermocouples in contact with the CF. The temperature fields of the wall as well as the deformed shape, by the end of the heating step, due to thermal expansion effects are given in Fig. 8. The modeling technique of this test series are as following: firstly, transient heat transfer analysis have been carried out to compute the temperature fields during heating and testing. Then, the temperature fields are employed as thermal loads in the transient thermomechanical analysis. The physical model of the thermomechanical analysis is the same as that reported in Fig. 5. During heating, load application, load holding and unloading steps, the boundary conditions of two fixed rigid plates and the ground are set to fully fixed. During heating, the two moving rigid plates were free to move then, during load application, load holding and unloading steps, the measured experimental reaction forces, in both directions, were applied as concentrated loads to the moving rigid plates.

Comparisons between the experimental and numerical displacement -time diagrams in the directions normal to bed and head joints during loading, holding and unloading steps are shown in Fig. 9. Good agreements between the experimental and numerical results can be observed. During load application step, and for both directions, an increase in the displacement due to the gradual closure of joints can be observed. The maximum displacement in the direction normal to bed joints is higher as compared to that in the direction normal to head joints due to the difference between the number of bed and head joints in the wall.

An example of the time variations of the stresses (from the end of heating step to the end of unloading step) in the CF and HF in the direction normal to bed and head joints (σ bed and σ head ) are reported in Fig. 10. By the start of the loading step, the HF carries almost all the applied load due to the wedged deformed shape of the wall (shown in Fig. 8) caused by the thermal expansion and the temperature gradient through the thickness of the wall. With the increase of the applied load and the deformation of the zone near the HF, the contact area between the moving plungers and the sides of the wall started to increase and, therefore, the stresses in the CF as well. After reaching the maximum load level and due to the lower creep rate in the CF (due to the lower temperature) as compared to the HF, a decay in σ bed and σ head HF was observed due to stress relaxation. On the other hand, an increase in the σ bed and σ head CF was observed and, then, remained almost constant during the holding time.

Biaxial relaxation behavior

The goal of this series was to investigate the relaxation behavior of the wall (constant strain loading conditions). The modeling technique is similar to that of biaxial creep behavior (transient heat transfer analysis then, transient thermomechanical analysis). The physical model of the thermomechanical analysis is the same as that reported in Fig. 5. During heating, the mechanical boundary conditions of the plungers and the ground are similar to those of biaxial creep behavior. During loading, the boundary conditions are almost the same as biaxial creep tests with only one difference, displacement boundary conditions were applied to the moving plungers and kept constant during the holding time.

Comparisons between the experimental and numerical time variations of the resulting reaction forces, in the directions normal to bed and head joints, during the two loading cycles are reported in Fig. 11. During loading (1 st cycle), the resulting reaction forces increased gradually to reach around 600 kN. Then, when the position of the plungers is locked, a decay in the resulting reaction forces was observed due to the relaxation behavior of the wall. In the beginning of the holding stage, a significant decrease in the reaction forces can be observed. Similar behavior was noticed for both loading cycles. The joints patterns in the wall by the end of the holding step and after unloading of the two loading cycles are presented in Fig. 12. It can be seen that some joints reopen after unloading. In addition, the percentage of the closed joints in the second loading cycle (Fig. 12-c) is higher as compared to that of the first loading cycle (Fig. 12-a). This can be attributed to the higher load applied to the wall in the second loading cycle as compared to that of the first loading cycle.

CONCLUSION

This paper reports on large scale experimental characterization of refractory masonry with dry joints and computationally efficient multi scale models developed at University of Orleans -France. The multi scale numerical model of masonry with dry joints takes into account the highly nonlinear orthotropic elastic viscoplastic behavior of refractory linings with dry joints and accounts for the gradual closure and reopening of joints. The model has been used to predict the mechanical behavior of alumina spinel refractory masonry subjected to a wide range of loading conditions including: in-plane biaxial compression loading/unloading at room temperature, biaxial compression creep at high temperature and biaxial relaxation at high temperature. It has been shown that the developed numerical model can reproduce with reasonable accuracy the nonlinear orthotropic elastic viscoplastic behavior of masonry with dry joints and accounts for gradual closure of joints and joints reopening. This model may help in the design and optimization of large scale refractory masonry linings used in many high temperature applications such as steel ladles, rotary kilns and furnaces.

Figure 1 :

 1 Figure 1: Top view of the biaxial compression test field.

Figure 2 :

 2 Figure 2: An example of refractory masonry wall subjected to biaxial compression at high temperature showing the arrangements of the LVDTs, moving and fixed plungers.

Figure 3 :

 3 Figure 3: Schematics of all possible joint patterns of refractory masonry with dry joints and joints closure and reopening criteria.

Figure 4 :

 4 Figure 4: Schematic of RVE in pattern O showing the physical model, brick and interface.

Figure 5 :

 5 Figure 5: FE model of masonry wall.

Figure 6 :

 6 Figure 6: Comparisons between experimental (2 tests) and numerical force -displacement diagrams of the masonry walls subjected to biaxial compression loading/unloading at room temperature.

Figure 7 :

 7 Figure 7: Time variations of the cold and hot faces temperatures during heating and mechanical testing.

Figure 8 :

 8 Figure 8: Temperature distributions of the masonry wall by the end of heating step showing the deformation of the wall due to thermal expansion.

Figure 9 :

 9 Figure 9: Time variations of the displacement in the masonry walls subjected to biaxial creep load at 1500 °C: experiments and simulations.

Figure 10 :

 10 Figure 10: Time variations of the stresses in the masonry walls subjected to biaxial creep load at 1500 °C: simulation results.

Figure 11 :

 11 Figure 11: Time variations of reaction forces, in the directions normal to bed and head joints, in a masonry wall subjected to a biaxial relaxation loading conditions at 1500 °C: experiments (2 tests and two load/unload cycles) and simulations.

FUNDING

  This work was supported by the funding scheme of the European Commission, Marie Skłodowska-Curie Actions Innovative Training Networks in the frame of the project ATHOR -Advanced THermomechanical multi-scale mOdelling Refractory linings 764987 Grant.

Figure 12 :

 12 Figure 12: Joints patterns in a masonry wall subjected to biaxial relaxation: (a) after first loading, (b) after first unloading, (c) after second loading and (d) after second unloading.

ACKNOWLEDGMENT

M. Ali is thankful to R. Oliveira (University of Coimbra) and U. Marschall (RHI Magnesita -Leoben) for their contribution to the experimental tests.