M R Hall

D Allinson

Transient numerical and physical modelling of temperature profile evolution in stabilised rammed earth walls

Three established stabilised rammed earth (SRE) mix types (433, 613, 703) were identified for analysis, in the form of 300mm thick test walls, by being subjected to different static air temperature and relative humidity differentials. The predictive numerical model outputs from WUFI Pro v4.1 hygrothermal simulations displayed good accuracy when validated against experimental data from physical modelling conducted using test walls in a climatic simulation chamber. The wall temperature profile evolution and resultant steady state gradients were very similar regardless of mix type indicating that the majority of the wall remained relatively dry. Unless liquid water is present, the thermal resistance and heat capacity of these materials does not change sufficiently to make significant differences to temperature profile evolution regardless of soil mix type. Little scope exists to intelligently modify the ability of SRE walls to absorb and store heat energy simply by manipulation of particle size distribution (PSD) and the resultant bulk density/void ratio relationships, under these conditions. Only the outer layers of the walls appear to interact with moisture in the air, and the predicted transient responses indicate that significant potential exists to

intelligently modify the ability of SRE walls to absorb, store and release moisture vapour from the surrounding air simply by manipulation of PSD and the resultant bulk density/void ratio relationships. an occupied building, industrial premises, or for the archiving/storage of sensitive media [START_REF] Hall | Analysis of the Hygrothermal Functional Properties of Stabilised Rammed Earth Materials[END_REF][START_REF] Houben | Earth construction: a comprehensive guide-Reissue edition[END_REF]. The passive buffering of relative humidity inside earth buildings has been reported [START_REF] Minke | The Earth Construction Handbook[END_REF] and the role of other absorbent building materials measured [START_REF] Padfield | The role of absorbent building materials in moderating changes of relative humidity[END_REF][START_REF] Olaken | Moisture buffering capacity of hygroscopic building materials: experimental facilities and energy impact[END_REF]. It is hypothesised that the functional properties controlling this behaviour are a product of controllable material characteristics in earth materials and that these can be optimised to suit the application. In order to begin to understand this requires a transient hygrothermal approach for predictive numerical modelling, which therefore requires validation by coupled physical modelling. One must consider both the climatic conditions and the material properties, both of which dynamically change as a function of time.

Keywords

The authors have previously shown that whilst the dry-state thermal conductivity for compacted earth material types appears to vary little, their moisture-dependent thermal conductivity varies significantly due to the wide range of possible internal pore structures. As a result, the moisture-dependent volumetric heat capacity can experience even greater fluctuations [START_REF] Hall | Analysis of the Hygrothermal Functional Properties of Stabilised Rammed Earth Materials[END_REF][START_REF] Hall | Assessing the moisture-content dependent parameters of stabilised earth materials using the cyclic-response admittance method[END_REF]. The authors have also previously found that (by analogy) the reverse is true for moisture, where the sorption isotherms between possible SRE types varies little, but that liquid water sorptivity/desorptivity, and also water vapour permeability, can vary over a significantly large range. Since temperature gradients induce mass transfer inside the material, and the thermal properties are moisture content-dependent, one can appreciate the fully transient nature of the coupled heat & mass transport phenomena associated with material behaviour in relation to air temperature and relative humidity buffering. The diagram in figure 1 illustrates hygrothermal behaviour in relation to the wall of a building [START_REF] Ashrae | Handbook of Fundamentals: Chapter 25 -heat, air and moisture control in building assemblies -fundamentals[END_REF]. In addition, phase change of moisture transporting within the material has to be considered in relation to the heat flux associated with its enthalpy of vaporisation.

This paper describes how a transient hygrothermal numerical model was used to simulate the behaviour of walls, made from three contrasting SRE material types, in response to changes in adjacent environmental conditions. These simulations were replicated as a full-scale physical model by constructing the said walls in a climatic simulation chamber and imposing realistic indoor/outdoor environments either side of the walls. The objective of the modelling exercises was to validate the numerical predictions for temperature profile evolution by assessing the hygrothermal behaviour of the materials.

Transient numerical modelling of heat & mass transport/storage

When considering simultaneous transient modelling of heat & mass transport/storage in porous media the thermodynamic behaviour is referred to as hygrothermal. Transient hygrothermal simulation tools include the following transport and storage phenomena:

• heat storage in dry, unsaturated and saturated porous media

• dry and moisture-dependent thermal conductivity (heat transport)

• latent heat transport as a result of vapour diffusion (including phase change due to vapour evaporation/condensation)

• water (mass) storage by water vapour sorption and liquid capillarity

• water vapour transport by diffusion

• liquid transport by surface diffusion and capillary conduction Numerous models exist for the definition of these phenomena (see [START_REF] Hens | Heat, Air and Moisture Transfer in Insulated Envelope Parts, Modelling, Common Exercises[END_REF] for a good summary).

Note that the benchmark standard BS EN 15026:2007 assumes that there is no total pressure difference across the component under consideration and air movement is ignored [START_REF] Bsi | Hygrothermal performance of building components and building elements-Assessment of moisture transfer by numerical simulation[END_REF]. The well-established WUFI ® model [START_REF] Künzel | Simultaneous Heat and Moisture Transport in Building Components -One-and two-dimensional calculation using simple parameters[END_REF][START_REF] Karagiozis | Advanced Hygrothermal Models and Design Models[END_REF] uses two coupled equations; one for heat transport (Eq 1) and the other for mass transport (Eq 2). They are solved simultaneously and numerically using an implicit finite volume method of discretisation with a tridiagonal matrix algorithm (known as the Thomas algorithm).

( )

( ) sat p p D t w ϕ δ ϕ ϕ ϕ ϕ ∇ + ∇ ⋅ ∇ = ∂ ∂ ⋅ ∂ ∂ (Eq. 1) ( ) ( ) ( ) sat p p T t T T H H vap ϕ δ λ ∇ ⋅ ∇ + ∇ ⋅ ∇ = ∂ ∂ ⋅ ∂ ∂ (Eq. 2)
In equation 1, the left side of the moisture balance terms represents water storage which is proportional to the derivative of the water retention curve ( w/ ), i.e. defined by the sorption isotherm. The divergence of liquid and vapour flow is computed using the transport terms on The hygrothermal material properties for each of the three SRE material types are taken from published data by the authors [START_REF] Hall | Analysis of the Hygrothermal Functional Properties of Stabilised Rammed Earth Materials[END_REF] and were used as the input parameters for the WUFI Pro simulations. As with most porous civil engineering materials (e.g. concrete, brick, soil) SRE is not homogenous because the pore size distribution and tortuosity vary greatly within the material. The conventional, accepted approach to measuring coefficients related to transport processes is to test a representative sample of the material and use this output as the volume averaged coefficient for that material, e.g. thermal conductivity. Obviously, no one individual part of the material has this property because it has effectively been statistically weighted, but across a representative area at full scale it is valid. The higher the sample population during testing, the greater the accuracy of the mean value. The required number of samples needed to achieve this, as detailed in the relevant BS EN ISO standards, was used for the testing conducted by the authors for this work [START_REF] Hall | Rammed earth sample production: context, recommendations and consistency[END_REF].

Both in the numerical model (WUFI Pro) and the physical model (climate chamber tests), the indoor and outdoor initial conditions for the simulation were taken as 21°C and 45% relative humidity. At time, t = 0 the outside climatic set points were adjusted to 8°C and 70% RH.

Further set point changes were introduced at 48 hours and 72 hours, when temperatures were reduced to 0°C and -8°C respectively whilst maintaining a constant 70% RH. The simulations ran for a total of 96 hours. The surface resistances for heat transfer and vapour diffusion were given typical values of 0.125m 2 K/W and 45.45x10 6 m 2 s Pa/kg respectively. These values were the same for both sides of the test walls to suit the conditions found in the climatic simulation chamber experiments.

Physical modelling in the climatic simulation chamber

The simulation chamber itself is composed of two self-contained halves called the 'design' side and the 'climate side' that join together on either side of the test specimen. The dimensions of each chamber are 4m long, 3m wide and 2.6m in height. The normal operating temperature range is +20°C to -15°C (± 5°C) depending on the internal conditions specified [START_REF] Taylor-Firth | Climatic Simulation and Environmental Monitoring -A Facility for Realistic Assessment of Construction Materials In-Service Performance[END_REF]. A schematic diagram illustrating the basic operation of the climatic simulation chamber is shown in Figure 2. Full-scale test components (e.g. walls) can be constructed inside the climatic simulation chamber following normal construction practices, and realistic climatic effects (on both sides of the walls) including various weather conditions can be accurately simulated and monitored under controlled laboratory conditions [START_REF] Laycock | Damp Towers -Understanding and Controlling the Ingress of Driven Rain through Exposed Solid Masonry Wall Structures[END_REF][START_REF] Hall | Assessing the Environmental Performance of Stabilised Rammed Earth Walls using a Climatic Simulation Chamber[END_REF]. The design side of the simulation chamber represents the interior conditions of a building and is typically configured to maintain standard indoor room conditions of 20°C ±1° and 45% RH ±5%. The climate side of the simulation chamber creates realistic sequences of different external weather conditions that can run in a fixed mode, sequence mode, or cyclic mode. The advantage of this is that the start and finish temperatures can be specified, as well as the rates of change of other weather components such as air pressure, relative humidity and rainfall.

The 300mm-thick SRE test walls were constructed inside the climatic simulation chamber on top of a concrete plinth, as described in previous work published by Hall [14]. The walls were made using the 433, 613 and 703 mix types detailed in previous studies [START_REF] Hall | Analysis of the Hygrothermal Functional Properties of Stabilised Rammed Earth Materials[END_REF][START_REF] Hall | Assessing the moisture-content dependent parameters of stabilised earth materials using the cyclic-response admittance method[END_REF][START_REF] Hall | Assessing the Environmental Performance of Stabilised Rammed Earth Walls using a Climatic Simulation Chamber[END_REF][START_REF] Hall | Rammed earth sample production: context, recommendations and consistency[END_REF]. The aggregates and silty clay, along with the mix types, were specified and controlled to be precisely the same as those used for the small representative test specimens. Prior to compaction the soils were mixed to their optimum moisture content and placed in layers 150mm thick inside oiled proprietary formwork and compacted using an Atlas Copco RAM30 pneumatic tamper. The test walls were separated from one another by a 100mm thick timberframe incorporating the formwork end boards. Following their completion, the test walls were allowed to cure for a minimum period of 28 days in laboratory conditions at 22°C ±1° and 40% RH ±5%. The wall edges were sealed with phenolic resin-coated plywood and caulked with silicone sealant. The remaining height above the test walls was completed with timber frame walling and covered with 1200g damp proof membrane sealed with black cloth tape.

The resultant dimensions for the exposed test face of each wall were 500mm wide by 900mm high (area = 0.45m 2 ).

An array of electronic sensors was retrofitted to each of the four rammed earth test walls by drilling a small hole from the inside of the wall to the required depth and then carefully inserting the sensor using metal rod marked at different depths. Once in place the sensor hole was capped and sealed using silicone caulking. In each wall, six type-T copper/constantan thermocouples (protected by hollow plastic tubing) were installed to depths of 25mm, 50mm, 75mm, 100mm, 125mm, and 150mm. In addition, armoured thermocouples were mounted onto both the internal and external wall faces.

The simulation parameters involved creating various fixed temperature & humidity gradients across the test walls. This was achieved in the climate chamber by maintaining a difference in temperature/humidity between the design side and the climate side of the chamber for a 24hour period. All of this testing was performed in the absence of any rainfall or pressure differential. The indoor (design side) conditions were maintained at 20°C ±1° and 40% RH ±5% to represent a comfortable indoor living environment. The outdoor conditions (climate side) were kept at a constantly high level of relative humidity: 75% ±5% to represent a damp, temperate climate. The test run lasted for a total period of four days and included three different outdoor temperature set points. On days 1 and 2 the 'outdoor' conditions were set to 8°C 75% RH, on day 3 they were set to 0°C 75% RH, and on day 4 they were set to -8°C 75% RH, during which snowfall began to occur. During all four days of the test run the sensor array was used for recording the performance of the test walls.

Analysis of simulated temperature and moisture profile evolution

The predictive numerical model outputs from WUFI Pro simulations were compared with experimental data from the physical model using the climatic simulation chamber. The apparent steady state temperature depth profiles appeared to be remarkably similar to one another for each soil mix type. This is consistent with previous findings for dry-state thermal conductivity on these soil mix types [START_REF] Hall | Analysis of the Hygrothermal Functional Properties of Stabilised Rammed Earth Materials[END_REF], and hence the thermal resistance of each wall is very similar indicating that the majority of the wall remained relatively dry. The presence of significant quantities capillary-domain water can cause large increases in the thermal conductivity of the material, which is significantly different between soil mix types according to their void characteristics.

The change across three external climate set points (8° C 75% RH; 0° C 75% RH; and -8° C 75% RH) presented a situation where the partial vapour pressure, and hence absolute moisture content, of the air on one side of the walls was significantly reduced in each case, as illustrated by the psychrometric chart in Figure 9. Figures 10,11 and 12 show the transient simulation results for relative humidity (of the internal pore structure) and moisture content of each test wall at 48, 72 and 96 hour increments. It is apparent that only the internal and external outer layers (approx. 5cm thick) were involved in the interaction with changing humidity over a daily period. One can appreciate the drying appears to occur on the external side when the temperature (and hence the air moisture content at 75% RH) is sufficiently low, and the absorption that occurs on the internal side as a result of the vapour pressure differential across the wall. Interestingly, there appears to be a large disparity between the three soil mix types in terms of i) the mass of absorbed water in the wall core, and ii) the rate of change in moisture content and relative humidity in the outer wall layers. This suggests that significant potential exists to intelligently modify the ability of SRE walls to absorb, store and release moisture vapour from the surrounding air simply by manipulation of PSD and the resultant bulk density/void ratio relationships. This merits further investigation into optimisation of the moisture buffering capacity of SRE materials and possible applications.

Some small discrepancies appeared between numerically and physically modelled temperature data, which could be explained by a number of specific issues. Due to a brief start up period prior to the beginning of the climate chamber simulations, a gradual decrease in the outside wall surface temperatures had already begun whereas for the computer simulation an instantaneous temperature drop from 21°C to 8°C was assumed. The climate side and design side air temperatures within the simulation chamber were not steady state throughout the test period, reflected by the small periodic fluctuation in the recorded wall surface temperatures (i.e. on/off cycling of the control plant). Also, some temperature stratification could reasonably be expected due to imperfect mixing of the air in the chamber and the induced temperature step changes did not occur instantaneously, as they did in the model.

Previous research found that the macro-scale properties of dry density and porosity of large scale SRE test walls could be matched very closely to small Proctor test specimens by accurately controlling PSD and optimum moisture content, regardless of varying compaction input energy. However, other scale effects such as surface finish, compaction layer thickness and degree of homogeneity could reasonably cause disparity between functional hygrothermal properties of the SRE test walls compared to the 1 litre test specimens used to determine the input parameters for the WUFI Pro simulations. Whilst the degree of accuracy in hygrothermally modelling temperature profiles in SRE materials is acceptably accurate, the opportunities for further improvements outlined above require much further research.

Conclusions

The numerically modelled data from WUFI Pro v4. The temperature profile evolution and resultant steady state gradients, in each test wall, were very similar regardless of soil mix type indicating that the majority of the wall remained relatively dry. This is further supported by the predicted liquid and relative humidity content and distribution within each test wall. From previous research [START_REF] Hall | Analysis of the Hygrothermal Functional Properties of Stabilised Rammed Earth Materials[END_REF][START_REF] Hall | Assessing the moisture-content dependent parameters of stabilised earth materials using the cyclic-response admittance method[END_REF] it is known that unless liquid water is present, the thermal resistance and volumetric heat capacity of these materials does not change sufficiently to make significant differences to temperature profile evolution regardless of soil mix type. Therefore, there appears to be little scope to intelligently modify the ability of SRE walls to absorb and store heat energy simply by manipulation of PSD and the resultant bulk density/void ratio relationships, under these conditions.

Only the outer layers of the walls appear to interact with moisture in the air. The predicted transient response of material moisture content indicate that significant potential exists to intelligently modify the ability of SRE walls to absorb, store and release moisture vapour from the surrounding air simply by manipulation of PSD and the resultant bulk density/void ratio relationships. The WUFI Pro model appears to show that the relationship between thermal and hygric (hygrothermal) functional properties must be simultaneously balanced in order to control the moisture sorption/desorption rates and storage capacity that govern air temperature and relative humidity buffering. Future research must attempt to scale up the simulation so that whole building modelling and real weather data can be included, and seasonal or annual timeframes can be modelled. Liquid water transport will also be studied so that rain fall absorption can be simulated, which will necessitate modelling over periods of many months. The potential of SRE walls to offer improved thermal comfort and operational energy savings from air conditioning inside buildings is still to be quantified. This work has identified the use of a hygrothermal approach to tackle this problem and has proven that data from numerical and physical models can be accurately validated under controlled climatic scenarios.

:

  transient thermal modelling; stabilised rammed earth; hygrothermal properties; climatic simulation chamber, temperature buffering; WUFI -based) materials have a good potential to buffer air temperature and relative humidity fluctuations in response to climatic variables in a controlled environment, e.g. inside

the right side of equation 1 .

 1 The partial vapour pressure differential is the motivational force for vapour diffusion, and is strongly temperature dependent since p sat exp(T). Liquid flow, on the other hand, is dominated by capillary forces which are chiefly a function of relative humidity (Kelvin equation). Note that the WUFI model ignores the influence of temperature on capillary-induced liquid flow. In equation 2, the sensible heat flow is simply represented by the divergence of heat conduction, whilst the divergence of vapour flow multiplied by the enthalpy of vaporisation represents the latent heat term on the right hand side. The derivative of total enthalpy, on the left hand side of equation 2, represents energy storage and contains the heat capacity of the dry material plus the enthalpy of absorbed water, which is itself phase-dependent (i.e. solid, liquid or vapour).

  1 has been validated to a good level accuracy against the physically modelled experimental data from the climatic simulation chamber. Experimentally determined hygrothermal properties for SRE materials can be inputted to WUFI Pro and used to accurately simulate transient behaviour of full sized test walls in response to changes in climatic conditions as a function of time.

Figure 1 : 4 ]Figure 2 :

 142 Figure 1: Hygrothermal behaviour illustrating the relevant fluxes acting on the wall of a

Figure 3 :

 3 Figure 3: Comparison between numerically and physically modelled temperature depth

Figure 4 :

 4 Figure 4: Comparison between numerically and physically modelled temperature depth

Figure 5 :Figure 6 :

 56 Figure 5: Comparison between numerically and physically modelled temperature depth

Figure 7 :

 7 Figure 7: Transient temperature profile evolution as a function of time for the 613 mix type

Figure 8 :

 8 Figure 8: Transient temperature profile evolution as a function of time for the 703 mix type

Figure 9 :

 9 Figure 9: A psychrometric chart to illustrate the change in environmental conditions in

Figure 10 :

 10 Figure 10: Transient liquid water and RH inside the pore structure of the 433 mix type in

Figure 11 :

 11 Figure 11: Transient liquid water and RH inside the pore structure of the 613 mix type in

Figure 12 :

 12 Figure 12: Transient liquid water and RH inside the pore structure of the 703 mix type in

  The temperature of the inside surface, mid-point and outside surface of each wall are plotted against elapsed time in Figures 3, 4 and 5, for each of the test walls. It can be seen that the simulation mainly predicts the temperatures to a good level of accuracy in relation to the recorded data from the climatic simulation tests.The temperature profile evolution, as a function of time, was calculated for the three walls and representative values at 0, 2, 6, 12, 24, 48, 72 and 96 hours were plotted on the graphs in Figures6, 7, and 8. It is notable that the temperature profiles do not become linear (i.e. apparent steady state conditions) until around 48 hours after the first temperature set points were initiated, and then at 72 and 96 hours respectively in response to the second and third temperature set point initiations. This enforces the important point that steady state heat loss assumptions are (in the absence of static temperature gradients) inappropriate for simulating the thermal behaviour under the varying environmental conditions that occur in nature. The small variations in volumetric heat capacity of the three soil mix types tested here appeared to make little difference to the time-dependent evolution of the apparent steady state temperature profiles. This suggests that little scope exists to intelligently modify the ability of SRE walls to absorb and store heat energy simply by manipulation of PSD and the resultant bulk density/void ratio relationships.

Figure 2

Acknowledgements

The authors wish to acknowledge the support of the Engineering and Physical Sciences Research Council for the research presented in this paper.


MD Salvatore Pellegrino

Vincenzo Villanacci

Naire Sansotta

Rosa Scarfì

Gabrio Bassotti

Giuseppe Vieni

Alessandro Princiotta

Concetta Sferlazzas

Giuseppe Magazzù

Giovanni Tuccari

Consolare Valeria

Redefining the intraepithelial lymphocytes threshold to diagnose gluten sensitivity in patients with architecturally normal duodenal histology

Keywords: X keyword = no topic, Coeliac disease < Disease-based, Diagnostic tests < Topics, Histopathology < Topics celiac disease, intraepithelial lymphocytes, sensitivity, specificity

INTRODUCTION

Celiac disease is a chronic inflammatory disease of the small intestine that is triggered by ingestion of gluten in genetically predisposed individuals [START_REF] Aga Institute | AGA Institute Medical Position Statement on the Diagnosis and Management of Celiac Disease[END_REF] . Virtually all patients with celiac disease have heterodimer HLA-DQ2 or HLA-DQ8 and the absence of these alleles provides a negative predictive value for the disease of close to 100% 2 . Thus, celiac disease is virtually excluded if the individuals lack these disease-associated alleles.

Although serologic testing for anti-tissue transglutaminase antibodies (anti-tTG) and antiendomysium antibodies (AEA) is highly sensitive and specific for the diagnosis of celiac disease, duodenal biopsy demonstrating villous atrophy (Marsh III) is considered the diagnostic gold standard [START_REF] Aga Institute | AGA Institute Medical Position Statement on the Diagnosis and Management of Celiac Disease[END_REF] . However, according to the Marsh-Oberhuber classification, celiac disease progresses through various stages and villous atrophy (Marsh III C) represents only the end stage of mucosal lesions, while the early changes are characterized by intraepithelial lymphocytosis (Marsh I) and subsequently hypertrophy of the crypts (Marsh II) [START_REF] Marsh | Gluten, major histocompatibility complex, and the small intestine. A molecular and immunobiologic approach to the spectrum of gluten sensitivity (celiac sprue)[END_REF][START_REF] Oberhuber | The histopathology of coeliac disease: time for a standardized report scheme for pathologists[END_REF] .

It is now recognized that even patients with mild enteropathy celiac disease in absence of villous atrophy, referred to as latent or potential celiac disease, are of sufficient clinical importance to be diagnosed and treated with a gluten-free diet as they may develop any of the manifestations encountered in celiac disease sufferers with villous atrophy. In fact, these patients may suffer from gastrointestinal symptoms and osteopenia or osteoporosis responding to gluten avoidance [START_REF] Kaukinen | Celiac disease without villous atrophy. Revision of criteria called for[END_REF][START_REF] Tursi | The symptomatic and histologic response to a gluten-free diet in patients with borderline enteropathy[END_REF] . Furthermore, they are at increased risk of autoimmune disease [START_REF] Hadjivassiliou | Does cryptic gluten sensitivity play a part in neurological illness?[END_REF][START_REF] Valentino | Markers of potential coeliac disease in patients with Hashimoto's thyroiditis[END_REF] and mortality [START_REF] Ludvigsson | Small-intestinal histopathology and mortality risk in celiac disease[END_REF] .

A recent randomized, controlled clinical study showed that patients with positive AEA and mild enteropathy (Marsh I-II) only benefit from a gluten-free diet, making long delay in diagnosis unacceptable, highlighting the importance of early recognition of mild enteropathy celiac disease [START_REF] Kurppa | Diagnosing mild enteropathy celiac disease: a randomized, controlled clinical study[END_REF] . Increased intraepithelial lymphocytes (IELs) along superficial epithelium of villi in duodenal mucosa has been considered a sensitive marker of the effect of gluten, but not specific, since this lesion can be found in several other diseases [START_REF] Chang | Pathological and clinical significance of increased intraepithelial lymphocytes (IELs) in small bowel mucosa[END_REF][START_REF] Brown | Intraepithelial lymphocytosis in architecturally preserved proximal small intestinal mucosa: anincreasing diagnostic problem with a wide differential diagnosis[END_REF] . However, until now, what the normal number of duodenal IELs should be has not been completely defined [START_REF] Dickson | Coeliac disease: an update for pathologists[END_REF][START_REF] Veress | Duodenal intraepithelial lymphocyte-count revisited[END_REF] .

Accordingly, several studies have focused on other early mucosal findings suggesting celiac disease. Counting IELs present at the tips of the villi might suggest mild enteropathy celiac disease in patients with an architecturally normal duodenal mucosal [START_REF] Goldstein | Morphologic features suggestive of gluten sensitivity in architecturally normal duodenal biopsy specimens[END_REF][START_REF] Mino | Role of lymphocytic immunophenotyping in the diagnosis of gluten-sensitive enteropathy with preserved villous architecture[END_REF][START_REF] Järvinen | Villous tip intraepithelial lymphocytes as markers of early-stage celiac disease[END_REF][START_REF] Biagi | Intraepithelial lymphocytes in the villous tip: do they indicate potential coeliac disease?[END_REF][START_REF] Salmi | Diagnostic methods beyond conventional histology in coeliac disease diagnosis[END_REF] .

Furthermore, IELs distribution pattern along villi may be more sensitive than formal IELs counts and it might also help in screening for mild enteropathy celiac disease [START_REF] Goldstein | Morphologic features suggestive of gluten sensitivity in architecturally normal duodenal biopsy specimens[END_REF] .

Nevertheless, since both counts and pattern distributions of IELs have a high rate of false positive cases, it has been suggested that they cannot be considered as a diagnostic marker of mild enteropathy celiac disease [START_REF] Chang | Pathological and clinical significance of increased intraepithelial lymphocytes (IELs) in small bowel mucosa[END_REF][START_REF] Dickson | Coeliac disease: an update for pathologists[END_REF] .

However, only one [START_REF] Järvinen | Villous tip intraepithelial lymphocytes as markers of early-stage celiac disease[END_REF] of five studies [START_REF] Goldstein | Morphologic features suggestive of gluten sensitivity in architecturally normal duodenal biopsy specimens[END_REF][START_REF] Mino | Role of lymphocytic immunophenotyping in the diagnosis of gluten-sensitive enteropathy with preserved villous architecture[END_REF][START_REF] Järvinen | Villous tip intraepithelial lymphocytes as markers of early-stage celiac disease[END_REF][START_REF] Biagi | Intraepithelial lymphocytes in the villous tip: do they indicate potential coeliac disease?[END_REF][START_REF] Salmi | Diagnostic methods beyond conventional histology in coeliac disease diagnosis[END_REF] Based on these inclusion criteria and considering the predictive value of celiac serology and HLA to diagnose celiac disease, we selected 3 study groups.

Pathology Assessment

All duodenal biopsies taken from the first and second portion of the duodenum were correctly oriented on acetate cellulose filters (Bio Optica Milan Italy), fixed in 10% neutral buffered formalin for 8-12 hrs at room temperature and then included in paraffin at 56°C.

Two serial 4 micron-thick sections were cut from each tissue block: one slide was stained with H&E and another was utilized for the IHC study. For this latter staining, sections were mounted on silane-coated slides and then subjected to antigen retrieval treatment by heating slides placed in 0.01 M citrate buffer pH 6.0 in a microwave oven for 3 cycles x 5 min each, before adding the primary antibody. The immunohistochemical procedure was given as the mean/100 epithelial cells, after counting of at least 300 epithelial cells along both sides of 5 villous bodies.

Statistical Analysis

Statistical analysis was performed using Statistica release 7 (StatSoft) and GraphPad Prism 5. D'Agostino-Pearson normality test was used to assess Gaussian distribution data.

Quantitative data were expressed as means, SD, ranges and 95% confidence intervals (CI). For categorical variables, Fisher's exact test was used. Comparison between multiple groups was carried out with one-way analysis of variance, followed by Bonferroni's Multiple comparison test; P value <0.001 was considered significant. The threshold of the IEL counts was defined as the mean + 2 standard deviations (SD) of non-celiac controls [START_REF] Hayat | Quantitation of intraepithelial lymphocytes in human duodenum: what is normal?[END_REF][START_REF] Mahadeva | Is a raised intraepithelial lymphocyte count with normal duodenal villous architecture clinically relevant?[END_REF] . For correlations between different IELs counts, the Pearson two-tailed t test was used.

RESULTS

The results of the IELs counts of the different stainings and study groups are summarized in Table 3.

IELs/20 enterocytes counts: the mean densities of tip IELs, both at H&E and CD3stainings (Figure 1), in the potential celiac disease group were higher than the control group (p<0.0000001), but not different from the latent group (p=0.12 and p=1 respectively). Comparing tip IELs counts of potential and latent celiac disease groups with the control group, there was no, or only slight, overlap following CD3-stainings, while there was slightly more overlap following H&E (Figure 2).

The upper limit of the normal range of tip IELs (mean + 2 SD) in the control group following H&E was 3.5 IELs, and for CD3-stainings 3.2 IELs per 20 enterocytes, respectively. The upper limit of 95%CI was 2.4 IELs/20 enterocytes at H&E and 2.8

IELs/20 enterocytes at IHC, respectively.

IELs/100 enterocytes counts: with regard to IELs counts per 100 enterocytes, the mean densities of tip IELs, both at H&E and CD3-stainings, in the potential celiac disease group were higher than the control group (p<0.0000001), but not different from the latent group (p=1 and p= 0.53, respectively). Comparing IELs counts of the potential and the latent groups with the control group, there was no overlap following CD3-stainings, while there was some overlap following H&E (Figure 3).

The cut-off value (mean +2 SD) for IELs in the control group was set at 18 IELs/100 enterocytes at H&E and 17 IELs/100 enterocytes at IHC. The upper limit of 95%CI was 13.6 IELs/100 enterocytes at H&E and 13.3 IELs/20 enterocytes at IHC, respectively.

Diagnostic accuracy of IELs counts: taking into account that there was no difference in

IELs counts between the potential and the latent celiac disease group, the accuracy of IELs to diagnose mild enteropathy celiac disease was calculated by merging both groups in a single patient population.

Considering the upper limit of the normal range of IELs (mean +2 SD) in control group, sensitivity and specificity of tip, and per 100 enterocytes counts following H&E and CD3stainings, are summarized in Table 4.

Correlation between IELs counts and between stainings: there was a high correlation between IEL counts per 20 and 100 enterocytes, both at H&E and CD3-staining (r=0.914 and r=0.929; p<0.0000001 respectively). Similarly, there was a high correlation between Considering the whole patient population, sensitivity and specificity of IELs pattern along villi at H&E are shown in Table 4. There was a difference comparing the whole patient population to control group regarding IELs pattern (p=0.00001).

DISCUSSION

Detecting early changes of celiac disease in patients with normal duodenal architecture has clinical relevance, but sometimes it is a challenge.

An increased number of IELs is the earliest pathological sign of mild enteropathy celiac disease [START_REF] Ferguson | Clinical and pathological spectrum of celiac disease: active, silent, latent, potential[END_REF] and several morphologic criteria that raise the possibility of gluten sensitivity in the context of an architecturally normal duodenal mucosa have been proposed, as recently reviewed by Dickson BC et al [START_REF] Dickson | Coeliac disease: an update for pathologists[END_REF] . Both sensitivity and specificity of a diagnostic test lie in the possibility of applying it to a control population explicitly not affected by the disease, and a study population known to be affected by the disease. Therefore, it is crucial, to avoid misleading results, to determine the upper normal value of duodenal IELs in control populations in which celiac disease has been excluded with certainty. In this setting, the ideal non-celiac disease control population should not only have negative serology, but, above all, the absence of any HLA allele associated with celiac disease, while the negative predictive value of HLA is virtually 100% an abnormal response to ingested gluten. Third, HLA haplotypes associated with celiac disease were seen in all patients belonging to the potential celiac disease group, although 3 out of 21 (14%) subjects had incomplete HLA-DQ2 heterodimer. However, 3 out of 8 patients with latent celiac disease had only one celiac disease-associated HLA allele, namely DQ2-DR7, and therefore we speculate that the presence of incomplete HLA heterodimer may be one factor that delays the development of intestinal villous atrophy.

Fourth, clinical and serological gluten dependence was present in all 9 patients of the potential celiac disease group who underwent a gluten-free diet.

With regard to villus tip IELs counts , the mean and range calculated, both at H&E and CD3 stainings, in the non-celiac control group, were comparable to the study of Jarvinen and coworkers [START_REF] Järvinen | Villous tip intraepithelial lymphocytes as markers of early-stage celiac disease[END_REF] , the only study that used a control group lacking HLA alleles associated with celiac disease. Accordingly, this method showed better diagnostic accuracy, especially regarding specificity, than previous studies [START_REF] Goldstein | Morphologic features suggestive of gluten sensitivity in architecturally normal duodenal biopsy specimens[END_REF][START_REF] Biagi | Intraepithelial lymphocytes in the villous tip: do they indicate potential coeliac disease?[END_REF][START_REF] Salmi | Diagnostic methods beyond conventional histology in coeliac disease diagnosis[END_REF] including the study of Jarvinen and coworkers [START_REF] Järvinen | Villous tip intraepithelial lymphocytes as markers of early-stage celiac disease[END_REF] who calculated diagnostic accuracy of tip IELs with a heterogeneous control population in which the HLA-negative group only represented a part.

With regard to villous IELs per 100 enterocytes counts, currently, the accepted normal limits are 30 or 25 IELs depending on classification schemes [START_REF] Dickson | Coeliac disease: an update for pathologists[END_REF][START_REF] Mahadeva | Is a raised intraepithelial lymphocyte count with normal duodenal villous architecture clinically relevant?[END_REF] In this study, comparing tip IELs against per 100 enterocytes counts , both at H&E and CD3-staining, classical counting on complete villous was slightly superior in terms of sensitivity in detecting mild enteropathy celiac disease regardless of staining. Several studies suggested higher diagnostic sensitivity of tip IELs than IELs counts on complete villus to disclose early-stage celiac disease [START_REF] Veress | Duodenal intraepithelial lymphocyte-count revisited[END_REF][START_REF] Mino | Role of lymphocytic immunophenotyping in the diagnosis of gluten-sensitive enteropathy with preserved villous architecture[END_REF][START_REF] Järvinen | Villous tip intraepithelial lymphocytes as markers of early-stage celiac disease[END_REF][START_REF] Biagi | Intraepithelial lymphocytes in the villous tip: do they indicate potential coeliac disease?[END_REF] . However, 2 of these studies only [START_REF] Mino | Role of lymphocytic immunophenotyping in the diagnosis of gluten-sensitive enteropathy with preserved villous architecture[END_REF][START_REF] Biagi | Intraepithelial lymphocytes in the villous tip: do they indicate potential coeliac disease?[END_REF] reported diagnostic accuracy data of both IELs analyses, where tip IELs counts , namely IELs per 20 enterocytes, were compared with complete villus counts expressed as the total IELs per mm of villus epithelium, while we performed IELs counts per 100 enterocytes.

With regard to histological stainings, it was suggested that the value of the IELs counts is higher following immunohistochemistry with CD3 antibody than in H&E-stained specimens, although CD3-staining is not currently routine application [START_REF] Dickson | Coeliac disease: an update for pathologists[END_REF] . In the present study overall diagnostic accuracy of IELs counts following CD3-staining and H&E were comparable, as demonstrated by the good correlation between H&E and CD3-staining, considering both villus tip and per 100 enterocytes counts, making IELs counting feasible following H&E where CD3-staining is not available in clinical practice.

Finally, considering the distribution of IELs along the villi, the loss of normal pattern of IELs, following H&E-staining, showed high sensitivity but unacceptable specificity in our study. However the interpretation of this Gestalt approach is highly observer-dependent and this may explain, regarding the different stainings used in these studies, the reported diagnostic variability of this method in previous studies [START_REF] Goldstein | Morphologic features suggestive of gluten sensitivity in architecturally normal duodenal biopsy specimens[END_REF][START_REF] Järvinen | Villous tip intraepithelial lymphocytes as markers of early-stage celiac disease[END_REF] . In this regard, Mino et al. [START_REF] Mino | Role of lymphocytic immunophenotyping in the diagnosis of gluten-sensitive enteropathy with preserved villous architecture[END_REF] , using a less observer-dependent method as such the IELs tip-to-base ratio, reported, following CD3-stainings, a sensitivity of 87.5%, and specificity 100% when mild In conclusion, the present study on diagnostic value of IELs findings was carried out using a control group certainly not affected by celiac disease and a population with potential celiac disease selected exclusively on the basis of multiple positivity of both AEA and anti-tTG.

Our data indicate that IELs counts both at tip and along villi are diagnostic feasible tools, in terms of both sensitivity and specificity, for disclosing mild enteropathy celiac disease.

The main result of our study is that the threshold of normal duodenal IELs may be lower than currently accepted and therefore, there is a need to update diagnostic criteria for celiac disease.

The application of a lower duodenal IELs cut-off, as suggested by the present study, may allow other centers to re-evaluate pending diagnosis of mild enteropathy celiac disease, thus avoiding complications of disease due to a delayed diagnosis. The number of participants satisfying the criteria for inclusion who did or did not undergo the index tests and/or the reference standard; describe why participants failed to undergo either test (a flow diagram is strongly recommended).

F

Not applicable

Test results

17 Time-interval between the index tests and the reference standard, and any treatment administered in between.

  Redefining the intraepithelial lymphocytes threshold to diagnose gluten sensitivity in patients with architecturally normal duodenal histology

  established normal values of these parameters in patients in whom celiac disease was excluded with certainty because of the absence of genetic predisposition, namely HLA-DQ2 and DQ8 alleles. Latent celiac disease represents a condition in which patients have positive celiac serology and genetic predisposition but normal duodenal mucosa architecture while on glutencontaining diet; subsequently they develop duodenal villous atrophy compatible with celiac disease. Therefore, latent celiac disease provides the opportunity to verify whether morphological features found in the absence of villous atrophy might predict the onset of overt celiac disease subsequently. The aim of this study was to evaluate the accuracy of different IELs counts and distribution pattern to diagnose mild enteropathy celiac disease in patients with architecturally normal duodenal mucosa. A potential celiac disease group and a latent celiac disease group represented the targeted patient populations. A non-celiac group lacking HLA celiac disease-associated alleles represented the control population.

  retrospective study designed to fulfil the STARD (Standards for Reporting of Diagnostic Accuracy) recommendations 20 . It was carried out at the Center for Celiac Disease of the University of Messina, which follows more than 2,500 patients and makes about 100 new diagnoses of celiac disease each year. We reviewed our computerized files of patients attending the Center during the years 2002-2009 and selected all cases corresponding to the following criteria: (1) histology showing architecturally normal duodenal mucosa. Each duodenal biopsy had to be well oriented with villous height/crypt depth ratio of at least 3:1; (2) availability of celiac serological data; (3) availability of genomic HLA-DQ2 and HLA-DQ8 typing.

  moist chamber with the following steps: 1) 0.1% H 2 O 2 in methanol to block intrinsic peroxidase activity; 2) normal sheep serum to prevent unspecific adherence of serum proteins; 3) with monoclonal primary antibody (mouse anti-human CD3, DakoCytomation, Denmark; working dilution 1:200, overnight at 4 °C); 4) anti-mouse immunoglobulin antiserum (Behring Institute; w.d. 1:25, 30 min at room temperature); 5) mouse anti-horseradish peroxidase-antiperoxidase complexes (DakoCytomation; w.d. 1:25, 30 min at room temperature). For the demonstration of peroxidase activity, sections were incubated in darkness for 10 min with 3-3' diaminobenzidine tetrahydrochloride (Sigma Chemical Co., Stl Louis, MO, USA), 100 mg in 200 ml 0.03% hydrogen peroxide in phosphate-buffered saline. Nuclear counterstaining was performed by Mayer's haemalum.To test the specificity of CD3 immunostaining in order to exclude the possibility of a nonspecific reaction, additional 4 micron thick serial sections of each specimen were tested by replacing the anti-mouse immunoglobulin antiserum by either phosphate buffered saline or normal rabbit serum: the results obtained were constantly negative. Moreover, the CD3 immunoreactivity demonstrated in lymphocytes outside of intestinal epithelial cells, but present in the lamina propria of duodenal biopsies, was considered as positive control.Finally, in order to test the inter-run variability of CD3 staining, the same CD3-positive lymphnode tissue section was utilized in every run. All routinely and immunohistochemically stained slides were blindly examined by two experienced pathologists; neither the final diagnosis nor the relationship of the biopsy to the patient's diet was known at the time of assessment. Histomorphological analysis was done on light microscopy using a x20 and x40 objective lenses and x10 eyepiece. Only peri-or supranuclear lymphocytes above the basal membrane of villi were assumed as IEL. IELs counts were initially performed on H&E stained sections and after, on the corresponding CD3 immunostained slides. The relative IELs counts were separately recorded by determining the quantification of the following parameters: a) existence of a pattern of sparse IEL infiltrate at H&E from the basal portion of villi to the middle and upper regions of villi, according to Goldstein 21 ; b) villus tip IEL score calculated by counting the number of IELs per 20 enterocytes from 5 randomly chosen villus tips and then computing the mean value of the 5 villus tip IEL counts; c) the number of IELs was

10 H&E

 10 and IHC staining considering IELs counts per 20 and 100 enterocytes (r=0.761 and r=0.848; p<0.0000001, respectively).

  patients with preserved villus architecture were compared to control subjects.

  celiac disease; IELs=intraepithelial lymphocytes; H&E= hematoxylin and eosin; SD= standard deviation; CI= confidence Interval a Reference value considered was the current lowest limit accepted of 25/100 enterocytes b Statistically significant difference compared to controls (p<0.0000001) c Statistically not significant difference compared to latent CD (p>0.05) d Statistically significant difference compared to controls (p= 0.00001)

a

  Potential CD group plus latent CD group b Non-CD patients without CD-associated HLA alleles served as controls to make comparisons and cut-off of IELs was defined as the mean + 2 standard deviations c Pattern of IELs along villi was evaluated only at H&E according to Goldstein (ref. 21).

Figure 1 .

 1 Figure 1. A,B,E,F: normal villi ( 3/1 villous-crypt ratio) with normal number of intraepithelial

Figure 2 .

 2 Figure 2. Villus tip intraepithelial lymphocytes (IELs) countings per 20 enterocytes at

Figure 3 .

 3 Figure 3. Villus tip intraepithelial lymphocytes (IELs) countings per 100 enterocytes at

  Figure 1.

  Figure 2.

F

Table 1

 1

	mean follow-up of these patients was 2.4 years (range 1-8 years). Nine of 21 potential
	celiac disease patients included in the study underwent a gluten free diet and clinical and
	serological gluten dependence was demonstrated in all cases
	The latent celiac disease group comprised 8 patients. One patient belonged initially to the
	potential celiac disease group, but during the study developed villous atrophy, and was
	then included in the celiac disease latent group. These patients were initially at the stage

summarizes demographic and clinical data in different groups. The potential celiac disease group comprised 21 patients at risk of celiac disease with positive celiac serology and HLA-DQ2 or DQ8 alleles. To exclude false positive results of celiac serology, only patients with both anti-tTG and AEA positive, in at least two measurements at different times were included. Indication for upper endoscopy in these patients was the suspect of celiac disease on the basis of positive celiac serology. The

  [START_REF] Karell | European Genetics Cluster on Celiac Disease. HLA types in celiac disease patients not carrying the DQA1*05-DQB1*02 (DQ2) heterodimer: results from the European Genetics Cluster on Celiac Disease[END_REF] . In our study, reference values for IELs were determined in a control population with low pre-test probability of celiac disease, negative celiac serology and a lack of HLA alleles associated with the disease, while, to our 11 subjects increases with age[START_REF] Savilahti | Density of gamma/delta+ T cells in the jejunal epithelium of patients with coeliac disease and dermatitis herpetiformis is increased with age[END_REF] . Therefore, we would have had an even lower IELs counts if we had enrolled younger controls As mentioned above, in evaluating a diagnostic test, it is important to have a group known to be affected by the disease, and there were several indicators that potential celiac disease patients in the study effectively had early-stage celiac disease. First, all patients have positive AEA and anti-tTG values, in at least two measurements at different times making false positive results highly improbable. Second, we included patients with latent celiac disease who initially presented positive celiac serology with a normal mucosa

	F o r
	P
	e
	e r
	R
	e
	v i e
	w

knowledge, only one study

[START_REF] Järvinen | Villous tip intraepithelial lymphocytes as markers of early-stage celiac disease[END_REF]

had a group lacking DQ2 or DQ8 haplotype as part of a nonceliac disease control group. Our control patients have an average age higher than the potential celiac disease group but it is noteworthy that intestinal density of IELs in healthy (potential celiac disease), but who later developed villous atrophy, to assess how the histological features of patients with potential celiac disease were comparable to those with latent celiac disease when they had a normal mucosa. In this regard, IELs count values of potential celiac disease group were comparable to those in latent patients, but different from the non-celiac disease group, meaning that changes in duodenal mucosa seen in our study group might predict future overt celiac disease or, in any case, indicate

Table 1 :

 1 Demographic and clinical data in study groups

	o r
	P
	e
	e r
	R
	e
	v i e
	w

Table 2 .

 2 Indicators of mild enteropathy celiac disease in potential and latent celiac disease groups with normal duodenal villous architecture

	CD=celiac disease

Table 3 .

 3 Mean values, SD, 95% confidence intervals, number of abnormal values of IELs counts in small bowel mucosal histological findings in study groups

		Study groups
	Potential CD	Latent CD	Non-CD
			Controls

Table 4 .

 4 Sensitivities and specificities with 95% confidence interval of IELs countings and distribution pattern in detecting mild enteropathy celiac disease in patient population a

		Sensitivity (%) (95%CI)	Specificity b (%) (95%CI)
		H&E	CD3-staining	H&E	CD3-staining
	IELs/20 enterocytes	89.6	89.6	100	100
		(72.6-97.8)	(72.6-97.8)	(76.8-100)	(76.8-100)
	IELs/100 enterocytes	93.1	93.1	100	100
		(77.2-99.2)	(77.2-99.2)	(76.8-100)	(76.8-100)
	IELs pattern c	96.5		78.5
		(82.2-99.9)		(49.2-95.3)
	CI= confidence interval; H&E= hematoxylin and eosin; IELs= intraepithelial lymphocytes;

  Technical specifications of material and methods involved including how and when measurements were taken, and/or cite references for index tests and reference standard. Definition of and rationale for the units, cut-offs and/or categories of the results of the index tests and the reference standard.7,8 10The number, training and expertise of the persons executing and reading the index tests and the reference standard. 7 11 Whether or not the readers of the index tests and reference standard were blind (masked) to the results of the other test and describe any other clinical information available to the readers.Statistical methods 12Methods for calculating or comparing measures of diagnostic accuracy, and the statistical methods used to quantify uncertainty (e.g. 95% confidence intervals).

	Section and Topic	Item		On page #
		#
	TITLE/ABSTRACT/	1	Identify the article as a study of diagnostic accuracy (recommend MeSH	1,2
	KEYWORDS		heading 'sensitivity and specificity').
	INTRODUCTION	2	State the research questions or study aims, such as estimating diagnostic	4
			accuracy or comparing accuracy between tests or across participant
			groups.
	METHODS
	Participants	3	The study population: The inclusion and exclusion criteria, setting and	5
			locations where data were collected.
		4	Participant recruitment: Was recruitment based on presenting symptoms,	5,6
			results from previous tests, or the fact that the participants had received
			the index tests or the reference standard?
		5	Participant sampling: Was the study population a consecutive series of	5
			participants defined by the selection criteria in item 3 and 4? If not,
			specify how participants were further selected.
		6	Data collection: Was data collection planned before the index test and	5
			reference standard were performed (prospective study) or after
			(retrospective study)?
	Test methods	7	The reference standard and its rationale.	5, 6
		8 9	o r F o r F r o	6-8
			P P P	7
			e e e	8
	RESULTS	13	e r e r e Methods for calculating test reproducibility, if done. r	Not done
	Participants	14 15	R R When study was performed, including beginning and end dates of recruitment. R Clinical and demographic characteristics of the study population (at least	5 5,6,19,20
		16	e e information on age, gender, spectrum of presenting symptoms). e
			v i e v i e v i e
			w w w

  of the results of the index tests (including indeterminate and missing results) by the results of the reference standard; for continuous results, the distribution of the test results by the results of the reference standard.

				5,6
		18	Distribution of severity of disease (define criteria) in those with the target	Not
			condition; other diagnoses in participants without the target condition.	applicable
		19	A cross tabulation Not done
		20	Any adverse events from performing the index tests or the reference	Not
			standard.	applicable
	Estimates	21	Estimates of diagnostic accuracy and measures of statistical uncertainty	8-10,22,23
			(e.g. 95% confidence intervals).
		22	How indeterminate results, missing data and outliers of the index tests	Not
			were handled.	applicable
		23	Estimates of variability of diagnostic accuracy between subgroups of	Not done
			participants, readers or centers, if done.
		24	Estimates of test reproducibility, if done.	Not done
	DISCUSSION	25	Discuss the clinical applicability of the study findings.	10-14

Specific author contributions:

All authors listed below approved the final draft submitted. Salvatore Pellegrino: planning and conducting the study, collecting and interpreting data, statistical analysis, drafting the manuscript; Vincenzo Villanacci: planning the study, interpreting data, drafting the manuscript; Naire Sansotta: collecting and drafting the manuscript; Rosa Scarfì: collecting data; Gabrio Bassotti: interpreting data and critical revision of the manuscript for important intellectual content; Giuseppe Vieni: collecting and interpreting data; Alessandro Princiotta: collecting data; Concetta Sferlazzas: planning the study and interpreting data ; Giuseppe Magazzù: planning the study, interpreting data, drafting the manuscript, critical revision of the manuscript for important intellectual content;

Giovanni Tuccari: conducting the study, interpreting data, drafting the manuscript, critical revision of the manuscript for important intellectual content.

Financial support: none

Potential competing interests: none

Acknowledgements: We thank Trays MacDonnell for revising the English of the article.


Philippe Angot
email: angot@cmi.univ-mrs.fr

Jean-Paul Caltagirone

Pierre Fabrie
email: pierre.fabrie@math.u-bordeaux1.fr

A Spectacular Vector Penalty-Projection Method for Darcy and Navier-Stokes Problems

Keywords: Vector penalty-projection, Penalty method, Splitting method, Multiphase Navier-Stokes/Brinkman, Anisotropic Darcy problem, Incompressible flows MSC 2010: 35Q30, 35Q35, 65M12, 65M85, 65N12, 65N85, 74F10, 76D05, 76D45, 76M25, 76R10, 76S05, 76T10

, to efficiently compute the solution of unsteady Navier-Stokes/Brinkman problems governing incompressible multiphase viscous flows. The method is also efficient to solve anisotropic Darcy problems. The key idea of the method is to compute at each time step an accurate and curl-free approximation of the pressure gradient increment in time. This method performs a two-step approximate divergence-free vector projection yielding a velocity divergence vanishing as O(ε δt), δt being the time step, with a penalty parameter ε as small as desired until the machine precision, e.g. ε = 10 -14 , whereas the solution algorithm can be extremely fast and cheap. The method is numerically validated on a benchmark problem for two-phase bubble dynamics where we compare it to the Uzawa augmented Lagrangian (UAL) and scalar incremental projection (SIP) methods. Moreover, a new test case for fluid-structure interaction problems is also investigated. That results in a robust method running faster than usual methods and being able to efficiently compute accurate solutions to sharp test cases whatever the density, viscosity or anisotropic permeability jumps, whereas other methods crash.

Introduction to model incompressible multiphase flows

Let Ω ⊂ R d (d =2 or 3 in practice) be an open bounded and connected domain with a Lipschitz continuous boundary Γ = ∂ Ω and n be the outward unit normal vector on Γ . For T > 0, we consider the following unsteady Navier-Stokes/Brinkman problem [START_REF] Khadra | Fictitious domain approach for numerical modelling of Navier-Stokes equations[END_REF] governing incompressible non-homogeneous or multiphase flows where Dirichlet boundary conditions for the velocity v |Γ = 0 on Γ , the volumic force f and initial data v(t = 0) = v 0 , ϕ(t = 0) = ϕ 0 ∈ L ∞ (Ω ) with ϕ 0 ≥ 0 a.e. in Ω , are given. For sake of briefness here, we just focus on the model problem [START_REF] Ph | Analysis of singular perturbations on the Brinkman problem for fictitious domain models of viscous flows[END_REF][START_REF] Ph | A penalization method to take into account obstacles in incompressible viscous flows[END_REF][START_REF] Ph | Vector penalty-projection methods for the solution of unsteady incompressible flows[END_REF] where d(v) = (∇v + (∇v) T )/2, as a part of more complex fluid mechanics problems.

ρ (∂ t v + (v• ∇)v) -2 ∇• (µ d(v)) + µ K -1 v + ∇p = f in Ω × (0, T ) (1)
∇• v = 0 in Ω × (0, T ) (2) ∂ t ϕ + v• ∇ϕ = 0 in Ω × (0, T ). ( 3
)
The permeability tensor K in the Darcy term is supposed to be symmetric, uniformly positive definite and bounded in Ω . We refer to [START_REF] Ph | Analysis of singular perturbations on the Brinkman problem for fictitious domain models of viscous flows[END_REF][START_REF] Khadra | Fictitious domain approach for numerical modelling of Navier-Stokes equations[END_REF] for the modeling of flows inside complex fluid-porous-solid heterogeneous systems with the Navier-Stokes/Brinkman or Darcy equations. The equation [START_REF] Ph | Vector penalty-projection methods for the solution of unsteady incompressible flows[END_REF] for the positive phase function ϕ governs the transport by the flow of the interface between two phases, either fluid or solid, respectively in the case of two-phase fluid flows or fluid-structure interaction problems. The force f may include some volumic forces like the gravity force ρ g as well as the surface tension force to describe the capillarity effects at the phase interfaces Σ . The advection-diffusion equation for the temperature T is not precised here and we assume some given state laws: ρ = ρ(ϕ, T ) and µ = µ(ϕ, T ) for each phase, where the functions are continuous and positive.

2 The fast vector-penalty projection method (VPP ε )

The (VPP ε ) method for multiphase Navier-Stokes/Brinkman

We describe hereafter the two-step vector penalty-projection (VPP ε ) method with a penalty parameter 0 < ε ≪ 1; see more details in [START_REF] Ph | A fast vector penalty-projection method for incompressible non-homogeneous or multiphase Navier-Stokes problems[END_REF]. For ϕ 0 with ϕ 0 ≥ 0 a.e. in Ω , v 0 and p 0 ∈ L 2 0 (Ω ) given, the method reads as below with usual notations for the semi-discrete setting in time, δt > 0 being the time step. For all n ∈ N such that

(n + 1) δt ≤ T , find ṽn+1 , v n+1 , p n+1 ∈ L 2 0 (Ω ), ϕ n+1 ∈ L ∞ (Ω ), such that: ρ n ṽn+1 -v n δt + (v n • ∇)ṽ n+1 -2 ∇• µ n d(ṽ n+1 ) + µ n K -1 ṽn+1 + ∇p n = f n (4) ε δt ρ n vn+1 -∇ ∇• vn+1 = ∇ ∇• ṽn+1 (5) v n+1 = ṽn+1 + vn+1 , and ∇(p n+1 -p n ) = - ρ n δt vn+1 (6)
p n+1 = p n + φ n+1 with φ n+1 reconstructed from ∇φ n+1 = - ρ n δt vn+1 (7) ϕ n+1 -ϕ n δt + v n+1 • ∇ϕ n = 0 (8)
with: ṽn+1 |Γ = 0, or for non homogeneous Dirichlet conditions: ṽn+1 |Γ = v n+1 D , and vn+1 • n |Γ = 0. Here v n , p n are desired to be first-order approximations of the exact velocity and pressure solutions v(t n ), p(t n ) at time t n = n δt. Since the end-ofstep velocity divergence is not exactly zero, the additional spherical part λ ∇• v I of the Newtonian stress tensor is included within the dynamical pressure gradient ∇p.

Once the equations (4-8) have been solved, the advection-diffusion equation of temperature can be solved too for T n+1 and we can find: ρ n+1 = ρ(ϕ n+1 , T n+1 ) and µ n+1 = µ(ϕ n+1 , T n+1 ).

The key feature of our method is to calculate an accurate and curl-free approximation of the momentum vector correction ρ n vn+1 in (5). Indeed (5-6) ensures that ρ n vn+1 is exactly a gradient which justifies the choice for ∇φ n+1 = ∇(p n+1 -p n ) since we have:

ρ n vn+1 = δt ε ∇ ∇• v n+1 ⇒ ∇(p n+1 -p n ) = - ρ n δt vn+1 = - 1 ε ∇ ∇• v n+1 . (9)
The (VPP ε ) method effectively takes advantage of the splitting method proposed in [START_REF] Ph | A new fast method to compute saddlepoints in constrained optimization and applications[END_REF] for augmented Lagrangian systems or general saddle-point computations to get a very fast solution of (5); see Theorem 1. When we need the pressure field itself, e.g. to compute stress vectors, it is calculated in an incremental way as an auxiliary step. We propose to reconstruct φ n+1 = p n+1 -p n from its gradient ∇φ n+1 given in (6) with the following method.

Reconstruction of φ n+1 = p n+1 -p n from its gradient. By circulating on a suitable path starting at a point on the border where φ n+1 = 0 is fixed and going through all the pressure nodes in the mesh, we get with the gradient formula between two neighbour points A and B using the mid-point quadrature: [START_REF] Sarthou | Eulerian-Lagrangian grid coupling and penalty methods for the simulation of multiphase flows interacting with complex objects[END_REF] with h AB = distance (A, B). The field φ n+1 is calculated point by point from the boundary and then passing successively by all the pressure nodes. This fast algorithm is performed at each time step to get the pressure field p n+1 from the known field p n . We refer to [START_REF] Ph | A fast vector penalty-projection method for incompressible non-homogeneous or multiphase Navier-Stokes problems[END_REF] for more details and validations on the present method.

φ n+1 (B) -φ n+1 (A) = B A ∇φ n+1 • dl = - B A ρ n δt vn+1 • dl ≈ - ρ n δt |v n+1 | h AB

The (VPP ε ) method for anisotropic Darcy problems

We present below the fast solution to incompressible Darcy flow problems in porous media with the (VPP ε ) method. The model problem reads in dimensionless form:

s ∂ t v + µ K -1 v + ∇p = f in Ω × (0, T ) (11) ∇• v = 0 in Ω × (0, T ) (12) v• n = 0 on Γ × (0, T ) (13)
where the viscosity µ > 0 is constant and the permeability tensor K is supposed to be symmetric, bounded in Ω and uniformly positive definite. The dimensionless stationarity parameter s > 0 includes the Darcy number: Da = K re f /L 2 re f and thus we have s ≪ 1 for most practical problems or even s = 0 for the steady anisotropic Darcy problem. The equations [START_REF] Vincent | Augmented Lagrangian and penalty methods for the simulation of two-phase flows interacting with moving solids. Application to hydroplaning flows interacting with real tire tread patterns[END_REF][START_REF] Youngs | Time-dependent multimaterial flow with large fluid distortion[END_REF](13)) also model flows inside heterogeneous porous-solid systems by letting the permeability tend to zero inside the impermeable media; see also [START_REF] Ph | Analysis of singular perturbations on the Brinkman problem for fictitious domain models of viscous flows[END_REF][START_REF] Khadra | Fictitious domain approach for numerical modelling of Navier-Stokes equations[END_REF] for the analysis and validations of the so-called L 2 volume penalty method.

The (VPP ε ) method with r = O(ε) > 0 and 0 < ε ≪ 1 to solve (11-13) reads as follows. For all n ∈ N such that (n + 1) δt ≤ T , find ṽn+1 , v n+1 and p n+1 such that:

s ṽn+1 -v n δt + µ K -1 ṽn+1 -r ∇ ∇• ṽn+1 + ∇p n = f n (14) ε s δt + µ K -1 vn+1 -∇ ∇• vn+1 = ∇ ∇• ṽn+1 (15) v n+1 = ṽn+1 + vn+1 , and ∇(p n+1 -p n ) = - s δt + µ K -1 vn+1 -r ∇ ∇• ṽn+1 (16)
p n+1 = p n + φ n+1 with φ n+1 reconstructed from its gradient ∇φ n+1 (17) with the boundary conditions: ṽn+1

• n |Γ = 0 and vn+1 • n |Γ = 0 on Γ . The space discrete solution to the prediction step ( 14) is explicit for s and r sufficiently small to invert a perturbation of the Identity matrix with a Neumann asymptotic expansion.

On the fast discrete solution to the (VPP ε ) method

The great interest for solving [START_REF] Ph | A fast vector penalty-projection method for incompressible non-homogeneous or multiphase Navier-Stokes problems[END_REF] or (15) instead of a usual augmented Lagrangian problem lies in the following result issued from [START_REF] Ph | A new fast method to compute saddlepoints in constrained optimization and applications[END_REF] which shows that the method can be ultra-fast and very cheap if η = ε/δt is sufficiently small. Let us now consider any space discretization of our problem. We denote by B = -div h the m × n matrix corresponding to the discrete divergence operator, B T = grad h the n × m matrix corresponding to the discrete gradient operator, whereas I denotes the n × n identity matrix with n > m and D the n × n diagonal nonsingular matrix containing all the discrete density values of ρ n > 0 a.e. in Ω . Here n is the number of velocity unknowns whereas m is the number of pressure unknowns. Then, the discrete vector penalty-projection problem corresponding to [START_REF] Ph | A fast vector penalty-projection method for incompressible non-homogeneous or multiphase Navier-Stokes problems[END_REF] with ε = η δt reads:

D + 1 η B T B vη = - 1 η B T B ṽ, with v η = ṽ + vη . (18)
We proved in [START_REF] Ph | A new fast method to compute saddlepoints in constrained optimization and applications[END_REF] the crucial result below due to the adapted right-hand side in the correction step (18) which lies in the range of the limit operator B T B. Indeed, (18) can be viewed as a singular perturbation problem with well-suited data in the righthand side. More precisely, we give in Theorem 1 the zero-order term of the solution vη to (18):

vη = - 1 η D + 1 η B T B -1 B T B ṽ ( 19
)
when the penalty parameter η is chosen sufficiently small; see the asymptotic expansion of vη and the proof in [4, Theorem 1.1 and Corollary 1.3].

Theorem 1 (Fast solution of the discrete vector penalty-projection).

Let D be an n × n positive definite diagonal matrix, I the n × n identity matrix and B an m × n matrix. If the rows of B are linearly independent, rank(B) = m, then for all η small enough, 0 < η < 1/ S -1 where S = BD -1 B T , there exists an n × n matrix C 1 bounded independently on η such that the solution of the correction step (19) writes for any vector ṽ ∈ R n :

vη = C 0 ṽ + η C 1 ṽ with C 0 = -D -1 B T S -1 B = -D -1 B T (BD -1 B T ) -1 B. ( 20
)
If rank(B) = p < m, there exists a surjective p × n matrix T such that B T B = T T T and a similar result holds replacing B by T . Hence, for a constant density ρ > 0 and choosing now η = ρ ε/δt, we have: D = I, S = BB T and C 0 = -B T S -1 B = -B T (BB T ) -1 B. Moreover, if rank(B) = p ≤ m ≤ n, the zero-order solution v = C 0 ṽ in (20) is the solution of minimal Euclidean norm in R n to the linear system: B v = -B ṽ by the least-squares method, and the matrix B † = B T (BB T ) -1 is the Moore-Penrose pseudo-inverse of B such that C 0 = -B † B. Indeed, a singular value decomposition (SVD) or a QR factorization of B yields: C 0 = -I 0 where I 0 is the n × n diagonal matrix having only 1 or 0 coefficients, the zero entries in the diagonal being the np null eigenvalues of the operator B T B.

Hence, for η small enough, the computational effort required to solve (18) amounts to approximate the matrix C 0 which includes both D and D -1 inside non commutative products. Thus, we always use the diagonal preconditioning in the case of a variable density which makes the effective condition number quasi-independent on the density or permeability jumps. We also use the Jacobi preconditioner in the prediction step (4) to cope with the viscosity or permeability jumps as performed in [START_REF] Khadra | Fictitious domain approach for numerical modelling of Navier-Stokes equations[END_REF]. However, for a constant density when D = I, we get C 0 = -I 0 . This explains why the solution can be obtained with only one iteration of a suitable preconditioned Krylov solver whatever the size of the mesh step or the dimension n; see the numerical results in [START_REF] Ph | A new fast method to compute saddlepoints in constrained optimization and applications[END_REF].

Numerical validations with discrete operator calculus

The (VPP ε ) method has been implemented with discrete exterior calculus (DEC) methods, see the recent review in [START_REF] Blair Perot | Discrete conservation properties of unstructured mesh schemes[END_REF], for the space discretization of the Navier-Stokes equations on unstructured staggered meshes. The (DEC) methods ensure primary and secondary discrete conservation properties. In particular, the space discretization satisfies for the discrete operators: ∇ h ×(∇ h φ ) = 0 and ∇ h • (∇ h ×ψ) = 0, which is not usually verified by other methods; see [START_REF] Blair Perot | Discrete conservation properties of unstructured mesh schemes[END_REF]. Hence, the (VPP ε ) method is now validated on unstructured meshes both in 2-D or 3-D.

The structure and solver of the computational code are issued from previous works, originally implemented with a Navier-Stokes finite volume solver on the staggered MAC mesh and using the Uzawa augmented Lagrangian (UAL) method to deal with the divergence-free constraint; see [START_REF] Khadra | Fictitious domain approach for numerical modelling of Navier-Stokes equations[END_REF]. We refer to [START_REF] Ph | Analysis of singular perturbations on the Brinkman problem for fictitious domain models of viscous flows[END_REF][START_REF] Ph | A penalization method to take into account obstacles in incompressible viscous flows[END_REF][START_REF] Khadra | Fictitious domain approach for numerical modelling of Navier-Stokes equations[END_REF] and the references therein for the analysis and numerical validations of the fictitious domain model using the so-called L 2 or H 1 -penalty methods to take account of obstacles in flow problems with the Navier-Stokes/Brinkman equations. Hence, our approach is essentially Eulerian with a Lagrangian front-tracking of the sharp interfaces accurately reconstructed on the fixed Eulerian mesh, see e.g. [START_REF] Sarthou | Eulerian-Lagrangian grid coupling and penalty methods for the simulation of multiphase flows interacting with complex objects[END_REF][START_REF] Vincent | Augmented Lagrangian and penalty methods for the simulation of two-phase flows interacting with moving solids. Application to hydroplaning flows interacting with real tire tread patterns[END_REF] and the references therein. Thus we use no Arbitrary Lagrangian-Eulerian (ALE) method, no global remeshing nor moving mesh method.

Fig. 1 Benchmark for 2-D bubble dynamics with (VPP ε ) method, ε = 10 -8 : motion of a circular bubble with surface tension at time t = 3 and Re = 35 -bubble initial diameter ⊘ = 0.05, ρ 1 /ρ 2 = 1000/100 = 10, µ 1 /µ 2 = 10/1 = 10, domain 0.1 × 0.2, mesh size 128 × 256, δt = 0.007143, circular bubble initially with no motion at height y = 0.05. LEFT: isobars and isoline ϕ = 0.5 of the phase function at interface. RIGHT: superposition of isoline ϕ = 0.5 at interface for (UAL), (SIP), (VPP) and vertical velocity field (in absolute referential).

Multiphase flows: dispersed two-phase bubble dynamics

The (VPP ε ) method is numerically validated for multiphase incompressible flows by performing with the three methods (UAL), (SIP) and (VPP), the benchmark problem studied in [START_REF] Hysing | Quantitative benchmark computations of two-dimensional bubble dynamics[END_REF] for 2-D bubble dynamics. In that problem, we compute the first test case which considers an initial circular bubble of diameter 0.05 m with density and viscosity ratios equal to 10 which undergoes moderate shape deformation. In this case, the bubble is driven up by the external gravity force f = ρ g, whereas the surface tension effect on the interface Σ between the two fluid phases is taken into account through the following force balance at the interface Σ :

[[v]] Σ = 0 and [[ -p I + µ ∇v + (∇v) T • n]] Σ = σ κ n |Σ , or f st = σ κ n |Σ δ Σ
where σ = 24.5 is the surface tension coefficient, κ the local curvature of the interface, n |Σ the outward unit normal to the interface and δ Σ the Dirac measure supported by the interface Σ . The solution of the phase transport (3) is carried out by the so-called VOF-PLIC method, i.e. the famous VOF method using a piecewise linear interface construction proposed in [START_REF] Youngs | Time-dependent multimaterial flow with large fluid distortion[END_REF] to precisely reconstruct the sharp interface Σ at the isoline ϕ = 0.5, with ϕ 0 = 0 in Ω 1 and ϕ 0 = 1 in Ω 2 ; see [10, 11].

The results of the three methods (UAL), (SIP) and (VPP) after 420 time iterations are presented in Figure 1 by superposing the different fields to get a more precise comparison. We observe an excellent agreement both between the three methods and the reference solution in [START_REF] Hysing | Quantitative benchmark computations of two-dimensional bubble dynamics[END_REF]. However, the (VPP) method runs faster. Fig. 2 ACF11-ball with (VPP ε ) method, ε = 10 -6 : free fall of a heavy solid body in air at time t = 0.15 and Re = 7358 -Cylinder diameter ⊘ = 0.05, ρ s = 10 6 , ρ f = 1, µ s = 10 12 , µ f = 10 -5 , domain 0.1 × 0.2, mesh size 256 × 512, δt = 0.0002, cylinder initially with no motion at height y = 0.15. LEFT: isobars and isoline ϕ = 0.5 of the phase function at interface. RIGHT: vertical velocity field and horizontal velocity isolines.

A test case for fluid-structure interaction problems

To evaluate the robustness of the (VPP ε ) method with respect to large density or viscosity ratios, we compute the motion of an heavy solid body which freely falls vertically in air with the gravity force f = ρ s g. The rigid behaviour of the body is obtained by letting the viscosity µ s tend to infinity inside the ball in order to penalize the tensor of deformation rate d(v). This fictitious domain method using a penalty was studied in [START_REF] Ph | Analysis of singular perturbations on the Brinkman problem for fictitious domain models of viscous flows[END_REF] (see the references therein) to design a numerical wind-tunnel, then numerically validated in several works, e.g. [START_REF] Vincent | Augmented Lagrangian and penalty methods for the simulation of two-phase flows interacting with moving solids. Application to hydroplaning flows interacting with real tire tread patterns[END_REF], and also analyzed theoretically in [START_REF] Ph | Analysis of singular perturbations on the Brinkman problem for fictitious domain models of viscous flows[END_REF][START_REF] Ph | A penalization method to take into account obstacles in incompressible viscous flows[END_REF] where optimal global error estimates are proved for the H 1 penalty method. Moreover, this fictitious domain method allows us to easily compute the forces applied on the obstacle, see [START_REF] Khadra | Fictitious domain approach for numerical modelling of Navier-Stokes equations[END_REF]; the error estimate being proved in [START_REF] Ph | Analysis of singular perturbations on the Brinkman problem for fictitious domain models of viscous flows[END_REF] when the nonlinear convection term is neglected inside the solid obstacle.

The results obtained by the (VPP ε ) method are presented in Figure 2 at time t = 0.15 s after 750 time iterations when the ball velocity reaches: V b = gt = 1.4715 m/s. The computation shows that the strain rate tensor inside the ball Ω s vanishes as d(v) L 2 (Ω s ) = O(µ f /µ s ), i.e. of the order of the machine precision. Hence, the (VPP ε ) method efficiently ensures both the rigidity of the solid body and a velocity divergence vanishing as O(ε δt) [START_REF] Ph | A fast vector penalty-projection method for incompressible non-homogeneous or multiphase Navier-Stokes problems[END_REF], whereas it avoids the blocking effect observed with other methods; see e.g. [START_REF] Vincent | Augmented Lagrangian and penalty methods for the simulation of two-phase flows interacting with moving solids. Application to hydroplaning flows interacting with real tire tread patterns[END_REF]. The (SIP) method crashes after a few time iterations. The (UAL) method is still able to compute the flow with a larger velocity divergence and the computation is far more expensive than with the (VPP ε ) method.


Philippe Angot
email: angot@cmi.univ-mrs.fr

Well-posed Stokes/Brinkman and Stokes/Darcy problems for coupled fluid-porous viscous flows 1

Keywords: Transmission problems, jump embedded boundary conditions, Stokes/Brinkman problem, Stokes/Darcy problem, fluid/porous coupled flows, well-posedness analysis, asymptotic analysis, vanishing viscosity, singular perturbation PACS: 02.30.Jr, 02.60.Lj, 45.10.Db, 47.10.Ad, 47.56.+r -MSC (2000): 35J25, 35J55, 35Q35, 65J20, 76D03, 76D07 76M45, 76S05

    Well-posed Stokes/Brinkman and Stokes/Darcy problems for coupled fluid-porous viscous flows

COUPLED FLUID-POROUS VISCOUS FLOWS

Notations. Let the domain Ω ⊂ R d (d = 2 or 3 in practice) be an open bounded and Lipschitz continuous domain. Let an interface Σ ⊂ R d-1 , Lipschitz continuous, separate Ω into two disjoint connected subdomains: the fluid domain Ω f and the porous one Ω p such that Ω = Ω f ∪ Σ ∪ Ω p ; see Fig. 1. For any quantity ψ defined all over Ω, the restrictions on Ω f and Ω p are denoted by ψ f and ψ p respectively. For a function ψ in H 1 (Ω f ∪ Ω p ), let ψ -and ψ + be the traces of ψ |Ω p and ψ |Ω f on each side of Σ respectively, ψ |Σ = (ψ + + ψ -)/2 the arithmetic mean of traces of ψ, and [[ψ]] Σ = (ψ + -ψ -) the jump of traces of ψ on Σ oriented by n. There exist in the literature different models with physically relevant stress or velocity jump boundary conditions for the tangential momentum transport at the fluid-porous interface Σ, see e.g. [20,14,18]. When the homogeneous porous flow is to be governed by the Brinkman equation, cf. [1,2,8,10,[START_REF] Hornung | Homogenization and porous media[END_REF], the interface condition below linking the jump of shear stress with a continuous velocity was derived with volume averaging techniques by Ochoa-Tapia and Whitaker [19] instead of the usual stress and velocity continuity boundary conditions at the interface [2]:

µ∇v f • n - µ φ ∇v p • n Σ • τ = µ β otw √ K v Σ • τ and v f = v p = v Σ on Σ, (1)
where the dimensionless parameter β otw is of order one; see [14,12,22] for its characterization. We show below, as a by-product of our general framework [5] recalled in the next Section that stress jump interface conditions of this type yield a well-posed fluid-porous Stokes/Brinkman problem whatever the dimensionless parameter β otw ≥ 0; see [6] for more details. This was not already stated up to our knowledge.

When the porous flow is governed by the Darcy equation, see e.g. [START_REF] Hornung | Homogenization and porous media[END_REF], the well-known Beavers and Joseph interface condition [9] must be used. It links the shear stress at the interface with the jump of tangential velocity:

(µ∇v f • n) |Σ • τ = µ α b j √ K v f -v p Σ • τ and v f • n = v p • n = v• n Σ on Σ, (2)
where the dimensionless parameter α b j = O( 1 √ φ ) depends on the porosity φ and may vary between 0.1 and 4. The approximate Saffman interface condition [21], derived by homogenization techniques in [16], is also written when the porous filtration velocity can be neglected with respect to the fluid velocity at the interface:

|v p Σ • τ| ≪ |v f Σ • τ|, i.
e. for a permeability value K or Darcy number Da = K/H 2 sufficiently small. The global solvability of the Stokes/Darcy problem with the Saffman condition for v p Σ • τ ≈ 0 is proved with a mixed hybrid formulation in [17] whatever the dimensionless parameter α b j ≥ 0, and then by many others with various formulations, see e.g. the recent review [13]. The only result of well-posedness for the full form of Beavers and Joseph condition is recently established in [11] for α 2 b j sufficiently small. We show further by a singular perturbation in our general framework with a vanishing viscosity that the above Beavers and Joseph interface conditions yield a well-posed Stokes/Darcy problem whatever the parameter α b j ≥ 0; see [6] for more details. Here, the main difficulty lies in how to give a sense to the tangential trace of the porous velocity on the interface with minimal regularity assumptions. This is particularly relevant for thin fluid layers as for conducting fractures in porous media flows [7,11].

WELL-POSED MODELS FOR COUPLED FLUID-POROUS VISCOUS FLOWS

We first describe the general framework with jump embedded boundary conditions proposed and studied in [5] which is also useful for fictitious domain methods in a similar way as in [4]. It is derived by a generalization to vector elliptic systems of a previous model stated for scalar problems [3,4]. (∇v + ∇v t ) being the strain rate tensor. We consider the following Stokes/Brinkman problem including jump embedded boundary conditions (J.E.B.C.) [5] on the interface Σ which link the trace jumps of both the stress vector σ (v, p)• n and the velocity vector v through the interface Σ:

Stokes/Brinkman transmission problem with jump embedded boundary conditions

-∇• σ (v, p) = f in Ω f , (3)
-∇• σ (v, p) + µ K -1 v = f in Ω p , (4)
∇• v = 0 in Ω f ∪ Ω p , (5)
v = 0 on Γ f ∪ Γ p , (6)
[[σ (v, p)• n]] Σ = M v |Σ on Σ, (7) σ (v, p)• n |Σ = S [[v]] Σ on Σ. (8)
Here, the viscosity coefficient µ and effective viscosity μ in the porous medium are bounded positive functions such that µ 0 = min(µ, μ) > 0, the symmetric permeability tensor K ≡ (K i j ) 1≤i, j≤d is uniformly positive definite, and the transfer matrices S, M on Σ are measurable, bounded and uniformly semi-positive matrices. With usual notations for Sobolev spaces, we now define the Hilbert spaces:

H 1 0Γ f (Ω f ) d ≡ w ∈ H 1 (Ω f ) d ; w |Γ f = 0 on Γ f , H 1 0Γ p (Ω p ) d ≡ w ∈ H 1 (Ω p ) d ; w |Γ p = 0 on Γ p , W ≡ w ∈ L 2 (Ω) d , w |Ω f ∈ H 1 0Γ f (Ω f ) d and w |Ω p ∈ H 1 0Γ p (Ω p ) d ; ∇• w = 0 in Ω f ∪ Ω p
equipped with the natural inner product and associated norm in

H 1 (Ω f ∪ Ω p ) d .
The model with the J.E.B.C. (7-8) also allows a possible pressure jump [[p]] Σ = 0 in H -1 2 (Σ) with additional regularity assumptions. Then, as a consequence of the general framework stated in [5], the problem (3-8) satisfies in Ω the nice weak formulation below:

Find v ∈ W such that ∀w ∈ W, a(v, w) = l(w) with a(v, w) = 2 Ω f µ d(v) : d(w) dx + 2 Ω p μ d(v) : d(w) dx + Ω p µ K -1 v• w dx + Σ M v |Σ • w |Σ ds + Σ S [[v]] Σ • [[w]] Σ ds l(w) = Ω f• w dx. ( 9
)
Then, we proved in [5, Theorem 1.1] that the problem (3-8) with f ∈ L 2 (Ω) d has a unique solution (v, p) ∈ W × L 2 (Ω) satisfying the weak form ( 9) for all w ∈ W.

Stokes/Brinkman problem with Ochoa-Tapia & Whitaker interface conditions

We now consider that μ = µ/φ , where φ ∈]0, 1] is the porosity of the porous medium, and stress jump interface conditions of Ochoa-Tapia & Whitaker's type [19] like in (1), the original ones reading with β τ = β otw and β n = 0:

[[σ (v, p)• n]] Σ = M v with M j j = µ β τ √ K τ , j = 1, • • • , d -1, M dd = µ β n √ K n and [[v]] Σ = 0 on Σ, ( 10
)
where M is a positive diagonal matrix with β τ , β n ≥ 0 a.e. on Σ and K τ , K n permeability coefficients. Then, as a consequence of the general framework stated in [5], the problem (3-6,10) satisfies in Ω the weak formulation below:

Find v ∈ V = {u ∈ H 1 0 (Ω) d ; ∇• u = 0} such that, 2 Ω f µ d(v) : d(w) dx + 2 Ω p µ φ d(v) : d(w) dx + Ω p µ K -1 v• w dx + Σ M v• w ds = Ω f• w dx, ∀w ∈ V. ( 11
)
Theorem 2.1 (Global solvability of Stokes/Brinkman problem with OT-W). If usual ellipticity assumptions hold, the problem (3)(4)(5)(6)10) with f ∈ L 2 (Ω) d has a unique solution (v, p) ∈ V × L 2 (Ω) satisfying the weak form (11) for all w ∈ V and such that p f = p f 0 +C 1 /2 and p p = p p 0 -C 1 /2 with p 0 ∈ L 2 0 (Ω) and the constant C 1 defined by:

C 1 = 1 |Σ| [[σ (v, p 0 )• n]] Σ -M v , n -1 2 ,Σ .
We can also interpret this solution as the limit solution of the problem (3-8) with penalized velocity jumps on Σ when the penalty parameter ε > 0 tends to zero; see [6] for the details.

Stokes/Darcy problem with Beavers & Joseph interface conditions

We consider the problem (3-8) with the Dirichlet boundary condition (6) on Γ p replaced by the stress boundary condition of Neumann where ν is the outward unit normal vector on Γ p and q ∈ H -1 2 (Γ p ) d given, e.g. q = -p e ν:

v = 0 on Γ f and σ (v p , p p )• ν = -p p ν + μ∇v p • ν = q on Γ p . (12)
Let us define the Hilbert space W N equipped with the natural inner product and norm in H 1 (Ω f ∪ Ω p ) d :

W N ≡ {w ∈ L 2 (Ω) d , w |Ω f ∈ H 1 0Γ f (Ω f ) d and w |Ω p ∈ H 1 (Ω p ) d ; ∇• w = 0 in Ω f ∪ Ω p }.
Then, as a corollary of [5, Theorem 1.1], see also [5,Theorem 2.1], the problem (3)(4)(5)(7)(8)12) has a unique solution (v, p) ∈ W N × L 2 (Ω).

For any ε > 0, let us now consider the solution (v ε , p ε ) ∈ W N × L 2 (Ω) of the problem (3)(4)(5)(7)(8)12) with a vanishing viscosity μ = ε for the Brinkman problem in Ω p . The condition (12) avoids the creation of a spurious boundary layer along Γ p for the Darcy problem when ε → 0. The J.E.B.C. (7)(8) are also calibrated as follows to obtain interface conditions of Beavers & Joseph's type [9] with a jump of tangential velocity (2) allowing a possible pressure jump:

[[σ (v, p)• n]] Σ = M v |Σ with M j j = 0, j = 1, • • • , d -1, M dd = µ β n √ K n on Σ, (13)
σ (v, p)• n |Σ = S [[v]] Σ with S j j = µ α τ √ K τ , j = 1, • • • , d -1, S dd = 1 ε on Σ, (14)
where M, S are positive diagonal matrices with α τ = α b j , β n ≥ 0 a.e. on Σ and K τ , K n permeability coefficients.

Let us define the Hilbert spaces

W S/D ≡ w ∈ L 2 (Ω) d , w |Ω f ∈ H 1 0Γ f (Ω f ) d , w |Ω p ∈ L 2 (Ω p ) d ; ∇• w = 0 in Ω f ∪ Ω p
equipped with the natural inner product and norm in

H 1 (Ω f ) d × L 2 (Ω p ) d and W S-D ≡ w ∈ W S/D ; ∇• w ∈ L 2 (Ω), [[w]] Σ ∈ L 2 (Σ) d , [[w• n]] Σ = 0
equipped with the norm defined by: w

2 W S-D = w 2 1,Ω f + w 2 0,Ω p + ∇• w 2 0,Ω + [[w]] Σ 2
0,Σ . We prove in [6] the following convergence result which also ensures the well-posedness of the Stokes/Darcy problem with Beavers & Joseph's type interface conditions (2,13) whatever the coefficients α τ , β n ≥ 0 a.e. on Σ. Theorem 2.2 (Convergence to Stokes/Darcy problem with B-J). With the data f ∈ L 2 (Ω) d and q = 0, the solution (v ε , p ε ) in W N × L 2 (Ω) for any ε > 0 of the problem (3)(4)(5)12,13,14) with a vanishing viscosity μ = ε weakly converges to the solution (v, p) in W S/D × L 2 (Ω) of the Stokes/Darcy problem with the interface conditions (2,13) on Σ when ε → 0. Indeed, in the porous domain Ω p , v p and p p satisfy the Darcy equation, i.e. Eq. ( 4) with μ = 0, and p p belongs to H 1 (Ω p ) such that p p = 0 on Γ p . With additional regularity assumptions such that v p ∈ H 1 (Ω p ) d , then v ∈ W S-D ∩ W N and we have the global error estimate with C > 0 depending on the data, ∇v 0,Ω p , ψ 0,Σ and ψ defined as the weak limit of 1 ε [[v ε • n]] Σ in L 2 (Σ):

||v ε -v|| 1,Ω f + √ ε ||v ε -v|| 1,Ω p + ||v ε -v|| 0,Ω p + ||p 0ε -p 0 || 0,Ω + 1 √ ε [[v ε • n]] Σ 0,Σ ≤ C ψ 0,Σ √ ε.

FIGURE 1 .

 1 FIGURE 1. Configuration for fluid-porous flows inside the domain Ω = Ω f ∪ Σ ∪ Ω p .

  Let σ (v, p) ≡ -p I + 2 μ d(v) denote the Newtonian stress tensor defined with the effective viscosity μ in the porous domain Ω p , with μ = µ in the fluid domain Ω f and d(v) ≡ 1 2

©2010 American Institute of Physics (AIP) Conference Proceedings -8th ICNAAM

2010, Rhodes (Greece), 19-25 sept. 2010, AIP-CP 1281, pp. 2208-2211, 2010.