Mouna Ben Jaber

Benoit Anet

Abdeltif Amrane
email: abdeltif.amrane@univ-rennes1.fr

Catherine Couriol

Thomas Lendormi

Pierre Le Cloirec

Gilles Cogny

Romain Fillieres

Benoît Anet

Romain Fillières

Impact of nutrients supply and pH changes on the elimination of hydrogen sulfide, dimethyl disulfide and ethanethiol by biofiltration

Keywords: Biofiltration, Sulfur compounds, Nutrients, pH, Elimination efficiency

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

Impact of nutrients supply and pH changes on the elimination of hydrogen sulfide, dimethyl disulfide and

ethanethiol by biofiltration.

Introduction

Odor nuisances related to industrialization is an emerging problem attracting the attention of public authorities who implement an increasingly strict regulation concerning odor emissions [START_REF]Pollutions olfactives : origine, législation, analyse, traitement[END_REF]. Perception thresholds of pollutants, especially those of sulfur compounds, which are very low, force to achieve high abatement rates; otherwise the residual concentrations may cause significant impacts on the welfare of local populations [START_REF] Anet | Characterization of gaseous odorous emissions from a rendering plant by GC/MS and treatment by biofiltration[END_REF][START_REF] Devos | Standardized human olfactory thresholds[END_REF][START_REF] Nagata | Measurement of odor threshold by triangle odor bag method, Odor Measurement Review, Office of Odor, Noise and Vibration Environmental Management Bureau[END_REF].

Various technologies are available for the treatment of these gaseous emissions; among them biological processes show some favorable advantages, since they are especially suitable for high gaseous flow rates (in the order of a hundred thousand of cubic meter of air per hour) containing low pollutant concentrations (below one gram per cubic meter) [START_REF] Le Cloirec | Biological Waste Gas Treatments[END_REF] [START_REF] Moletta | Management of environmental problems of agro-food industries[END_REF][START_REF] Verguet | What can be expected from an odour control system at a wastewater treatment plant (in French)[END_REF]. In addition, biological processes in general and biofiltration in particular, show limited costs and an ease of use justifying their use for the deodorization of industrial emissions [START_REF] Kennes | Review: Waste gas biotreatment technology[END_REF][START_REF] Devinny | Biofiltration for air pollution control[END_REF][START_REF] Delhoménie | Biofiltration of Air: A Review[END_REF][START_REF] Shareefdeen | Biotechnology for odor and pollution control[END_REF].

A biofilter consists of a porous organic or inorganic bed, through which passes a humid polluted gaseous stream. In a biofilter, the contaminants are absorbed from the gas phase to the biofilm, consisting of a thin layer of water and microorganisms, where the biological reactions occur [START_REF] Cloirec | Introduction au traitement de l'air[END_REF]. Pollutant transfer and biodegradation kinetics rely on material characteristics such as: pH, water holding capacity, organic matter and nutrients content, and structural behaviour [START_REF] Kennes | Review: Waste gas biotreatment technology[END_REF][START_REF] Mudliar | Bioreactors for treatment of VOCs and odours -A review[END_REF]. Biofiltration is therefore a complex process, involving various phenomena, such as mass transfer, adsorption, absorption and microbial activity [START_REF] Malhautier | Biofiltration of volatile organic compounds[END_REF]; they are closely linked to the biofilm and hence the packing material [START_REF] Deshusses | Biological waste air treatment in biofilters[END_REF], which represents therefore the cornerstone of the process [START_REF] Mudliar | Bioreactors for treatment of VOCs and odours -A review[END_REF][START_REF] Elias | Evaluation of a packing material for the biodegradation of H 2 S and product analysis[END_REF]. Intrinsic packing material's properties induce the establishment of a more or less conducive environment for the development of an effective microbial consortium and homogeneous gas distribution through the bed.

For this purpose, the chemical, physical and microbial properties, as well as the economical impact, of 11 organic and inorganic packing materials potentially suitable for biofiltration have been previously investigated in order to select the most relevant for the treatment of rendering gaseous emissions.. Among them and according to its appropriate pH, water holding capacity and its highest nutrients content and colonization at the biofilter start up, composted would mulch showed the best odour removal efficiency. However, the pressure drop of the compost bed decreased suggesting structural changes which may impact the performances on long term. On the contrary, pine bark showed almost constant height and pressure drop during the operation time suggesting that its structure remained stable over time, while its performances remained significant [START_REF] Anet | Characterization and selection of packing materials for biofiltration of rendering odourous emissions[END_REF]. In addition, the use of pine bark has been recommended by other authors in odor treatment [START_REF] Luo | The use of pine bark and natural zeolite as biofilter media to remove animal rendering process odours[END_REF], while total H 2 S removal was observed in a biofilter containing compost-based bed [START_REF] Hort | Evaluation of sewage sludge and yard waste compost as a biofilter media for the removal of ammonia and volatile organic sulfur compounds (VOSCs)[END_REF].

Therefore, the packing material considered in the present work was constituted of a stratum of pine bark covered by a stratum of composted wood mulch to combine the physico-chemical benefits of the former and the efficiency of the latter. The packing medium should contain nutrients in a sufficient amount to allow microbial activity; however it is not always the case.

Consequently, the impact of nutrients supply on treatment performance was examined.

Most of the laboratory studies focus on the treatment of gaseous effluent containing only one pollutant; while industrial effluents usually contain a complex mixture of pollutants, and hence removal efficiency differs significantly from that observed in the case of a single pollutant [START_REF] Hwang | Biofiltration of waste gases containing both ethyl acetate and toluene using different combinations of bacterial cultures[END_REF][START_REF] Khammar | Link between spatial structure of microbial communities and degradation of a complex mixture of volatile organic compounds in peat biofilters[END_REF][START_REF] Lee | Degradation characteristics of methyl ethyl ketone by Pseudomonas sp. KT-3 in liquid culture and biofilter[END_REF][START_REF] Chan | Compounds interaction on the biodegradation of acetone and methyl ethyl ketone mixture in a composite bead biofilter[END_REF]. Removal capacities of biofilters differ according to the considered pollutant and biodegradability decreases as follows: alcohols > esters > ketones > aromatic compounds > sulfur compounds, showing the recalcitrance of sulfur compounds [START_REF] Anet | Characterization of gaseous odorous emissions from a rendering plant by GC/MS and treatment by biofiltration[END_REF][START_REF] Smet | Abatement of high concentrated ammonia loaded waste gases in compost biofilters[END_REF], which have been therefore selected for this study. They also show very low perception thresholds associated with high olfactory contributions [START_REF] Anet | Characterization of gaseous odorous emissions from a rendering plant by GC/MS and treatment by biofiltration[END_REF] and can be nearly always found in industrial emissions.

Among sulfur compounds, biodegradability decreases as follows: H 2 S > methyl mercaptan > dimethyldisulfur (DMDS) > dimethylsulfur (DMS) [START_REF] Myung Cha | Removal of organo-sulphur odour compounds by Thiobacillus novellus SRM, sulphur-oxidizing microorganisms[END_REF], and total H 2 S removal was usually the case [START_REF] Iranpour | Literature review of air pollution control biofilters and biotrickling filters for odor and volatile organic compound removal[END_REF][START_REF] Ramirez | Biofiltration of reduced sulphur compounds and community analysis of sulphur-oxidizing bacteria[END_REF]. Among them, H 2 S and DMDS were selected for this study, as well as ethyl mercaptan (EtSH) which is also often reported in gaseous emissions.

Biological oxidation was realized by sulfato-oxidative bacteria, leading to H 2 SO 4 and CO 2 [START_REF] Kelly | Organic sulfur compounds in the environment: biogeochemistry, microbiology, and ecological aspects[END_REF] and hence affecting significantly the pH, which has been therefore monitored during the course of the experiments.

Materials and methods

Experimental Set Up

The biofiltration unit consisted of three identical PVC columns with an internal diameter of 150.4 mm. For clarity reasons, only one column is shown in Figure 1.

Each column was composed of two stratums: a 50 cm pine bark layer chosen for its good physical proprieties supports, a 50 cm composted wood mulch stratum selected according to its propensity to bring diversified microbial communities and its low cost.

The compressed air was humidified at room temperature and fed in the bottom of the biofilter using a bubbler thus achieving a relative humidity close to 99%.

The reserve of liquid pollutants (DMDS and EtSH supplied by Acros Organics with respective purity of 98 and 99% ), supplying the two mass regulators disposed in series, was a 250 mL steel 316 Swagelok cylinder pressurized under 2 bar of compressed air. The operating principle was based on the evaporation of a liquid, diluted in a carrier gas stream. Rates of these two streams were respectively controlled by regulators μ-FLOW mass liquid and gaseous EL-FLOW (Bronkhorst ®, Montigny-Les-Cormeilles, France). The volatilization chamber was flushed by a controlled flow rate of a compressed air stream supplied under 5 bars (Bronkhorst ®, Montigny-Les-Cormeilles, France). High concentrated flow of hydrogen sulfide was provided to the system by a gas tank concentrated at 2500 ppm in nitrogen (Linde Gas France S.A., Noyal sur Vilaine, France). The feed rate of each column was measured by a flowmeter stainless steel ball (Brooks® GT1355, Hatfield, PA, USA).

Biofilter watering was ensured by a peristaltic pump for selecting a nutrient solution or deionized water at neutral pH and at a flow rate of 0.344 L day -1 . The dispersion was carried on by stainless steel full cone nozzle GG-SS1 FullJet supplied by Spraying Systems Emani® (130 ml min -1 , spray angle: 60°, at 0.5 bar) (Arcueil, France). This solution was used to maintain the growth of microorganisms in the packing material in order to improve pollutant removal in the biofilter.

Two pH measuring points were located respectively at 25 and 75 cm from the biofilter bottom. Eight points for gas sampling, separated by 12.5 cm each, were distributed over the entire height of each column.

The packing material

The characterization of packing materials properties was realized according to standards methods [START_REF] Anet | Characterization and selection of packing materials for biofiltration of rendering odourous emissions[END_REF][START_REF] Dorado | A comparative study based on physical characteristics of suitable packing materials in biofiltration[END_REF]; these characteristics are reported in Table 1.

The used pine barks were a mixture of black pines and maritime pines and supplied by Hobiofilter S.A.R.L. (Mazères, France). The composted wood mulch was supplied by Europe Environnement (Aspach Le Haut, France).

Inoculation

At the biofilter start-up, the three different columns were inoculated with 500 mL of diluted activated sludge (initial dry weight 1.2 g L -1 ) recovered from wastewater treatment plant (Beaurade, Rennes, France).

Analytical methods

Inlet and outlet H 2 S, EtSH and DMDS concentrations were measured by a TRS MEDOR® analyzer (Chromatotec, France). The separation was performed on a capillary column swept by reconstituted air under 230 mbar, followed by an electrochemical detection in a cell filled with CrO 3 at 10 g L -1 . Samples of 400 µL of gas were used for the measurements. The retention times were 70, 110 and 290 s for dimethyl disulfide (DMDS), hydrogen sulfide (H 2 S) and ethanethiol (EtSH), respectively. The detection limits ranged between 35 and 45 µg m -3 for these pollutants. The external calibration was managed with a mixture containing 20 mg m -3 (± 2 mg m -3 ) for each component supplied by Linde Gas France S.A. (Noyal sur Vilaine, France) and internal calibration was managed by a dimethylsulfide permeation furnace (107.4 ng min -1 at 45°C).

The pH of the leachates and in the packing materials was measured using a CyberScan pH510. In order to limit the gaseous distribution through the packing, the pH was measured with a penetration probe (Mettler Toledo, Greifensee, Switzerland).

Operating conditions and start-up

A 3.2 m 3 h -1 flow rate was applied leading to 60 s empty bed retention time (EBRT). During this study, an increase in pollutant concentrations was applied in order to compare the removal efficiencies of the three biofilters and to follow the influence of nutrients supply facing concentrations increase. The operating conditions are reported in Table 2. It should be observed that the considered pollutant concentrations were selected according to the range of values reported in the literature [START_REF] Anet | Characterization of gaseous odorous emissions from a rendering plant by GC/MS and treatment by biofiltration[END_REF][START_REF] Kastner | Comparison of chemical wet scrubbers and biofiltration for control of volatile organic compounds using GC/MS techniques and kinetic analysis[END_REF][START_REF] Legrand | Treatment of recalcitrant organic sulphur compounds by biofiltration: optimization of operating conditions for an industrial application[END_REF][START_REF] Soupramanien | Biological treatment of a mixture of gaseous sulphur reduced compounds: identification of the total bacterial community's structure[END_REF].

Nutrients were supplied by a nutritive solution dispersed on the surface of the biofilter. In order to follow the impact of nutrients supply on the degradation of a mixture of pollutants, the three columns were watered by different solutions with a fixed rate of 0.344 L day -1 . The first biofilter (BIO1) operated without nutrient supply; it was just sprinkled by deionized water. Different nutritive solutions were supplied to the second (BIO2) and third (BIO3)

biofilters. The composition of these solutions is shown in Table 3. Na2CO3 was added to supply for inorganic carbon since regarding sulfur compounds, even if both autotrophic [START_REF] Smet | Treatment of waste gases contaminated with odorous sulfur compounds[END_REF] and heterotrophic [START_REF] Kennes | Fungal biocatalysts in the biofiltration of VOC-polluted air[END_REF][START_REF] Rappert | Odor compounds in waste gas emissions from agricultural operations and food industries[END_REF] pathways can be considered, the autotrophic pathway is the most often reported in the available literature [START_REF] Smet | Treatment of waste gases contaminated with odorous sulfur compounds[END_REF][START_REF] Cha | Removal of organo-sulphur odour compounds by Thiobacillus novellus SRM, sulphur-oxidizing microorganisms[END_REF]; while composted wood mulch was expected to supply for organic carbon. Regarding K 2 HPO 4 and (NH 4 ) 2 SO 4 they were added to supply for nitrogen, phosphate and sulfate for microbial growth.

The removal efficiency RE of each bed was calculated as follows:

RE (%) = ( C ¿ -C out ) C ¿ ×100
Where RE: removal efficiency (%); C in : inlet concentration (mg m -3 ); C out : outlet concentration (mg m -3 ).

Results and Discussion

Influence of nutrients supply on treatment performance

During the first phase (until 14 days), it was observed that H 2 S was totally eliminated from the first day (Fig. 2a), showing that nutrients supply had no influence on the abatement of H 2 S in the considered conditions. Contrarily and even at low pollutants concentration (first phase), about 3 and 5 days were needed to achieve total EtSH and DMDS removal (Figs. 2c andb).

The negligible impact of nitrogen and phosphorus addition on H 2 S removal should be most likely related to their presence in the organic packing material used, composted wood mulch and pine bark. In addition, it was shown that in the presence of a mixture of sulfur compounds, H 2 S is preferentially metabolized [START_REF] Iranpour | Literature review of air pollution control biofilters and biotrickling filters for odor and volatile organic compound removal[END_REF][START_REF] Chung | Structure of the bacterial community in a biofilter during dimethyl sulfide (DMS) removal processes[END_REF].

An increase of the pollutant concentrations (7.6 mg m -3 -phase 2 -from day 14 to day 79)

clearly highlighted the influence of nutrients supply on the degradation of EtSH and DMDS.

Indeed, the removal efficiencies of these compounds in the columns sprinkled by nutrient solutions (BIO2 and BIO3) were higher than those in the first biofilter operating just with deionized water. The positive impact of nutrients appeared however limited in the case of DMDS, while it was especially pronounced in the case of EtSH. Indeed, final removal efficiencies were 95 and 82% for DMDS and EtSH in the absence of nutrients (BIO1) (Figs. 2b andc), while total DMDS removal (Fig. 2b) and 95% EtSH removal (Fig. 2c) were observed for a low supplementation (BIO2) and total removal of both pollutants for the high supplementation (BIO3) (Figs. 2b andc). Therefore and in most cases, the presence of nutrients is required to improve the biodegradation of recalcitrant sulfur compounds, in agreement with the available literature showing that nutrient deficiency may become a limiting factor for biofiltration efficiency [START_REF] Kumar | Biofiltration of Volatile Organic. Compounds (VOCs) -An Overview[END_REF]. These results clearly showed that microorganisms consumed preferentially the most easily biodegradable compound (H 2 S) in the biofilter to the detriment of the most recalcitrant pollutants (DMDS and EtSH). When treating a mixture of sulfur compounds, the presence of H 2 S inhibits the removal of some organic reduced compounds (EtSH and DMDS) which are more recalcitrant to biodegradation [START_REF] Ramirez | Biofiltration of reduced sulphur compounds and community analysis of sulphur-oxidizing bacteria[END_REF].

From these results, an order of degradation of these compounds can be established as follows:

H 2 S > DMDS > EtSH.
It should be observed that analysis of the impact of nutrients addition can only be qualitative but not quantitative owing on the one hand to the involvement of the packing material, and especially composted wood mulch, in the supply of nutrients, and on the other hand to the pathway involved in organic sulphur compounds assimilation, such as EtSH and DMDS, which can be both autotrophic [START_REF] Smet | Treatment of waste gases contaminated with odorous sulfur compounds[END_REF] and heterotrophic [START_REF] Kennes | Fungal biocatalysts in the biofiltration of VOC-polluted air[END_REF][START_REF] Rappert | Odor compounds in waste gas emissions from agricultural operations and food industries[END_REF].

To examine system stability and reliability, experiments were ceased for two weeks (Fig. 2). A clear impact on DMDS and EtSH removals was shown, while it appeared negligible on H 2 S elimination, since no decrease of its removal yield can be observed (Fig. 2a). Regarding DMDS, even if a decrease was noted after cessation of operation, the abatement rapidly increased again until recovering stable removal efficiencies, close to those observed before cessation of operation (Fig. 2b). Contrarily, a more significant impact on EtSH abatement can be noted, since an important lapse of time was needed (until approximately day 70) before recovering removal efficiencies observed before cessation of operation (Fig. 2c).

Influence of pH on the removal of H 2 S, DMDS and EtSH

The pH has a significant impact on treatment efficiency; microbial activity and colonization is disturbed by any modification of pH in the biofilter [START_REF] Anet | Packing material evolutions and odorous abatement of peat and heather biofilters operating in rendering industry[END_REF][START_REF] Hernandez | Simultaneous removal of H 2 S, NH 3 , and Ethyl mercaptan in biotrickling filters packed with poplar wood and polyurethane foam: Impact of pH during startup and crossed effects evaluation[END_REF]. The pH was therefore monitored in the middle of each stratum, namely at 25 and 75 cm height, and showed a decrease in the two stratums (Figure 3). At the end of experiment, in all biofilters the pH was lower than 6.5

in the stratum of composted wood mulch and less than 3 in the stratum of pine bark. From this, the nutrient solution or water added to biofilter can maintain pH on the top of the bed, while the decrease was more pronounced at the bottom, from about 7 to less than 3, versus from about 8.5 to 6.5 at the top.

It should however be observed the similar pH time-courses observed for the three biofilters, close pH in all biofilters can be noted for a given time (Fig. 3). The impact of the nutrients solution on the pH appeared therefore negligible.

In addition, at the beginning of the experiments, the pH of the aqueous phase collected at the bottom of each column was neutral. Over time these leachates became acidic (Fig. 4); from the 10 th day, the pH was less than 6 and decreased until final values close to 2.2 (Fig. 4).

Both autotrophic [START_REF] Smet | Treatment of waste gases contaminated with odorous sulfur compounds[END_REF] and heterotrophic [START_REF] Kennes | Fungal biocatalysts in the biofiltration of VOC-polluted air[END_REF][START_REF] Rappert | Odor compounds in waste gas emissions from agricultural operations and food industries[END_REF] pathways can be considered for the assimilation of organic sulphur compounds, such as EtSH and DMDS. However, the autotrophic pathway, as it was also the case for H 2 S assimilation, is the most often reported in the available literature [START_REF] Smet | Treatment of waste gases contaminated with odorous sulfur compounds[END_REF][START_REF] Cha | Removal of organo-sulphur odour compounds by Thiobacillus novellus SRM, sulphur-oxidizing microorganisms[END_REF]. Sulfur compounds assimilation as energy sources leads to the production of H 2 SO 4 inducing medium acidification [START_REF] Smith | Mechanism of oxidation of dimethyl disulphide by Thiobacillus thioparus strain E6[END_REF], as experimentally confirmed since the amount of sulfates accumulated in the leachate at the end of the first running period was close to 1 g. H 2 S removal was not affected by pH; it was completely removed in the biofilter.

Regarding EtSH and DMDS and even if their removal was reported to be optimal when the pH of the bed remains between 6 and 8 [START_REF] Premkumar | Effect of secondary parameters on biofilter treating industrial effluent[END_REF], biofilter acidification with time did not seem to have an impact on their assimilation, since and even after a disturbance phenomenon (cessation of operation) optimal removal efficiencies were recovered.

Removal of H 2 S, DMDS and EtSH along the biofilter

The elimination of H 2 S, EtSH and DMDS was followed all along the biofilters height.

Removal efficiencies values at days 3 and 72, namely for the low and high pollutant amounts are given in Figs. 5 and6 for instance. In addition, the mean pH values in the stratum of pine bark and that of composted wood mulch are given in Table 4.

As indicated above (see 3.2), the degradation of H 2 S was not affected by the pH of the packing materials; microorganisms involved in H 2 S degradation can tolerate acidic conditions [START_REF] Li | Field measurements and modeling of two-stage biofilter that treats odorous sulfur air emissions[END_REF]. For the low pollutant concentrations, the influence of nutrients supply appeared also almost negligible since no real difference between the three biofilters can be underlined (Fig. 5a). The total bed height was not needed for a complete H 2 S removal, since it was totally eliminated on 87.5 cm height irrespective of the considered biofilter.

About half of the H 2 S was removed in the pine bark (first 50 cm), namely in the most acidic part of the bed (mean pH of 6.63 -Table 4) and the rest was removed in the composted woold much, namely at alkaline pH (mean pH of 8.54 -Table 4).

Contrarily to H 2 S, EtSH and DMDS removals were clearly affected by nutrients supply in the aqueous solution, and this impact was especially pronounced in the case of a total absence of nutrients in the spraying water (Figs. 5b andc). Only about 30% of EtSH and DMDS were removed in the first part of the bed, pine bark, and only 63 and 72% of DMDS and EtSH on the total bed height in the absence of nutrients supply, respectively; while in the presence of nutrients in the watering solution removal yields were in the range 90-97% (Figs. 5b andc).

During the second phase, namely at high pollutant concentrations, possible pH impact on the removal efficiencies can be more clearly highlighted, since medium acidification continued and from the 65 th day pH below 3 were measured in the three columns in the first stratum of the bed, pine bark. As observed, H 2 S removal was not affected by the pH, since even in the absence of nutrients supply only 62.5 cm of the bed were needed for its complete elimination and about 82% was removed in the pine bark (Fig. 6a), where the mean pH was 2.9 (Table 4).

Microbial populations involved in H 2 S degradation appeared therefore active in a large pH range, from less than 3 in the pine bark at the end of the second phase to close to 9 in the composted wood much in the beginning of the first phase (Table 4).

The impact of nutrients supply on the one hand and that of pH on the other hand were clearly andc). Regarding biofiltration at high pollutant amounts (BIO3), it should be observed that less than 20% pollutants degradation were obtained on the first 40 cm of pine bark (Figs. 6b andc). The main part of DMDS and EtSH was therefore removed in the second stratum, composted wood mulch, namely at neutral pH (mean pH close to 6 -Table 4), since an almost total DMDS removal was observed in all biofilters (Fig. 6b) and an almost total EtSH removal in BIO2 and BIO3 (Fig. 6c). The pH impact was therefore clearly shown; DMDS and EtSH removal occurred mainly at neutral pH in agreement with the relevant literature [START_REF] Silva | Modelling of the biofiltration of reduced sulphur compounds through biotrickling filters connected in series: Effect of H2S[END_REF].

The impact of nutrients supply was also clearly underlined, and appeared significant for the three pollutants. Indeed, pine bark was almost enough for a total H 2 S degradation for the high supplementation (Fig. 6a), while about 62 cm of bed was needed for its total removal during the other experiments (BIO1 and BIO2 -Fig. 6a). An improvement of the DMDS and EtSH degradation rates were also shown at high supplementation (Figs. 6b andc); while a total absence of nutrients had also an impact on the rate of EtSH degradation, as well as on its final removal yield, less than 80% versus about 95% for the other experiments, BIO2 and BIO3 (Fig. 6c).

Conclusions

The obtained results showed a significant impact of nutrients supply and pH on sulfur compounds biodegradation. An improvement of the DMDS and EtSH degradation rates were shown at high nutrients supply; while a total absence of nutrients had an impact on both the rate of EtSH degradation, and its final removal yield, less than 80% versus an almost total removal for the other experiments. Nitrogen and/or phosphorus limitations should be most likely considered to account for this nutrients effect; this assumption should be subsequently validated. An effect of the nutrients addition was also shown for H 2 S at high supplementation, even if it was always totally removed.

Microbial populations involved in H 2 S degradation appeared to be active in a large pH range, from less than 3 to close to 9; while and in agreement with the relevant literature, DMDS and EtSH removal occurred mainly at neutral pH. As expected, the most easily biodegradable compound (H 2 S) was preferentially consumed to the detriment of the most recalcitrant pollutants (DMDS and EtSH).

  shown regarding the degradation of the other sulfur compounds, DMDS and EtSH. Indeed, their degradation in the pine bark remained limited, since DMDS and EtSH removal efficiencies in the stratum of pine bark were 10 and 15% in the absence of nutrients supply, 10 and 15% for the low supplementation, and 40 and 25% for the high supplementation(Figs.6b

aTable 2 .Table 3 .Table 4 .Figure 1 .

 2341 Figure legendsFigure 1. Schematic representation of the experimental setup: 1 air inlet, 2 humidifier, 3 reserve of liquid pollutants (DMDS+EtSH), 4 volatilization chamber, 5 H 2 S gas cylinder, 6 flow meter, 7 water or nutrient solution reservoir, 8 pump, 9-11 pH measuring points, 12-18 points of gas sampling, 19 leachate outlet.

Figure 2 .

 2 Figure 2. Elimination efficiency (%) of H 2 S (a), DMDS (b) and EtSH (c) in the three biofilters.

Figure 3 .

 3 Figure 3. pH time-courses in the three biofilters in the middle of the first and second stratum, 25 and 75 cm bed height respectively.

Figure 4 .

 4 Figure 4. pH time-courses in the leachates of the three biofilters.

Figure 5 .

 5 Figure 5. Removal of H 2 S (a), DMDS (b) and EtSH (c) along the biofilter height (on the 3 rd day).

Figure 6 .

 6 Figure 6. Removal of H 2 S (a), DMDS (b) and EtSH (c) along the biofilter height (on the 72 nd day).

Table 1 .

 1 Characteristics of the packing materials used.

	Properties	Composted wood mulch	Pine bark
	pH	8.62	4.75
	%C	36.6%	48.1%
	%O	25.5%	45.4%
	%H	3.5%	5.9%
	%N ε (-) a	1.1% 45.5%	0.1% 56.3%
	ρ (kg m -


Florian Lemarié

Patrick Marchesiello

Laurent Debreu

Eric Blayo

Sensitivity of ocean-atmosphere coupled models to the coupling method : example of tropical cyclone Erica

Keywords: ocean-atmosphere coupling, domain decomposition methods, air-sea interface, Schwarz algorithms Couplage océan-atmosphère, méthodes de décomposition de domaine, interface air-mer, Algorithmes de Schwarz

In this paper, the sensitivity of Atmospheric and Oceanic Coupled Models (AOCMs) to the coupling method is investigated. We propose the adaptation of a Schwarz-like domain decomposition method to AOCMs. We show that the iterative process of the method ensures consistency of the coupled solution across the air-sea interface, contrarily to usual ad-hoc algorithmic approaches. The latter are equivalent to only one iteration of a Schwarz-like iterative method, which does not provide a converged state. It is generally assumed that this lack of consistency does not affect significantly the physical properties of the solution. The relevancy of this statement is first assessed in a simplified problem, then in the realistic application of a mesoscale atmospheric model (WRF) coupled with a regional oceanic model (ROMS) to simulate the genesis and propagation of tropical cyclone Erica. Sensitivity tests to the coupling method are carried out in an ensemble approach. We show that with a mathematically consistent coupling the spread of the ensemble is reduced, suggesting that there is room for further improvements in the formulation of AOCMs at a mathematical and numerical level.

1 Introduction

Context

Many applications in global/regional oceanography and meteorology require Atmospheric and Oceanic Coupled Models (AOCMs) which account for important air-sea feedbacks. Separate integrations of the oceanic and atmospheric numerical models in forced mode (i.e. without feedback from one component to the other) may be satisfactory for numerous applications and process studies (e.g. [START_REF] Marchesiello | Equilibrium structure and dynamics of the California Current System[END_REF][START_REF] Colas | Heat balance and eddies in the Peru-Chile Current System[END_REF][START_REF] Lemarié | Are There Inescapable Issues Prohibiting the use of Terrain-Following Coordinates in Climate Models[END_REF]. However, two-way coupling is essential for analyzing energetic and complex phenomena like tropical cyclones (e.g. [START_REF] Bao | Numerical simulations of air-sea interaction under high wind conditions using a coupled model : A study of hurricane development[END_REF][START_REF] Jullien | Ocean feedback to tropical cyclones: climatology and processes[END_REF], eastern boundary upwellings (e.g. [START_REF] Perlin | Numerical simulation of air-sea coupling during coastal upwelling[END_REF][START_REF] Capet | Eddies in easternboundary subtropical upwelling systems[END_REF], and more generally climate trends (e.g. [START_REF] Large | Attribution and impacts of upper-ocean biases in CCSM3[END_REF][START_REF] Terray | The role of the frequency of SST coupling in the Indian Monsoon variability and monsoon-ENSO-IOD relationships in a global coupled model[END_REF]. In that case, connecting the two model solutions at the air-sea interface is an arduous task which jointly raises mathematical, physical and computational issues.

Besides the oceanic and atmospheric models, the formulation of AOCMs is based on several components: a parameterization of the turbulent air-sea fluxes, a coupling infrastructure, and a coupling algorithm. Several parameterizations for surface atmospheric flow dynamics under oceanic influence have been designed. Those parameterizations are derived at a semi-empirical level and are based on field and laboratory experiments designed to carefully tune the parameter values [START_REF] Fairall | Bulk parameterization of air-sea fluxes: updates and verification for the COARE algorithm[END_REF][START_REF] Large | Surface fluxes for practitioners of global ocean data assimilation[END_REF]. As for computational issues, numerous coupling softwares were developed in the last decade [START_REF] Hill | The Architecture of the Earth System Modeling Framework[END_REF][START_REF] Joppich | MpCCI -a tool for the simulation of coupled applications[END_REF][START_REF] Redler | OASIS4 -a coupling software for next generation earth system modelling[END_REF]. Those tools are needed to handle message passing, synchronization in time and regridding procedures (i.e., interpolation/extrapolation) between the computational grids. The last ingredient, of numerical nature, in the design of an AOCM is a consistent coupling algorithm; this is the main subject of this paper. The notion of consistency associated with the ocean atmosphere coupling problem will be clarified below in Sec. 1.2. It is known that coupled model solutions exhibit a strong sensitivity to model parameters [START_REF] Bengtsson | From short-range barotropic modelling to extended-range global weather prediction: a 40-year perspective[END_REF][START_REF] Mcwilliams | Irreducible imprecision in atmospheric and oceanic simulations[END_REF]. More specifically, sensitivity to perturbations in the initial conditions [START_REF] Ploshay | Large sensitivity to initial conditions in seasonal predictions with a coupled ocean-atmosphere general circulation model[END_REF], coupling frequency (Lebeaupin Brossier et al., 2009;[START_REF] Terray | The role of the frequency of SST coupling in the Indian Monsoon variability and monsoon-ENSO-IOD relationships in a global coupled model[END_REF][START_REF] Masson | Impact of intra-daily SST variability on ENSO characteristics in a coupled model[END_REF], and to air-sea flux formulation (Lebeaupin Brossier et al., 2008) have been reported, in addition to the sensitivity that is inherent to each component of the coupled system [START_REF] Tribbia | The reliability of improvements in deterministic shortrange forecasts in the presence of initial state and modeling deficiencies[END_REF]. This is arguably a source of concern when it comes to assess the AOCMs solution. Uncertainties in the specification of air-sea fluxes is strongly related to the empirical nature of atmospheric surface layer parameterizations resulting from the extreme complexity of the physical processes that we wish to represent. In general, a given parameterization is designed and tuned for a particular geographical region or a particular range of static stability. A lot of efforts are being directed toward improving physical parameterizations because these are usually considered as the major source of errors in AOCMs. However, it is arguably crucial to keep working on other model developments and identify other possible sources of error/deficiency. Preliminary studies of the fifth phase of the Coupled Model Comparison Project (CMIP5) do not report significant improvements in biases in present-day climate compared to CMIP3 (e.g. [START_REF] Roehrig | The present and future of the west african monsoon: a process-oriented assessment of cmip5 simulations along the amma transect[END_REF][START_REF] Gupta | Climate drifts in the cmip5 models[END_REF]. This suggests that finer horizontal resolution is not sufficient to cure the aforementioned problems. There is clearly a need for a more complete understanding of what goes on in numerical coupled models, notably through a finer consideration of numerical methods and modeling assumptions. To our knowledge, no systematic sensitivity study of the algorithmic aspects of air-sea coupling has been reported. We investigate this sensitivity in the present paper, starting with the definition of the coupling problem.

Ω oce Ω atm L oce U o = f oce L atm U a = f atm B atm U a = g atm B oce U o = g oce Γ @Ω ext atm @Ω ext oce surface layer ABL OBL habl hobl K o m (z) K a m (z) ∆U = U a (z a ) -U o (z o ) z z x, y

Coupling Problem

We symbolically describe the oceanic and atmospheric circulation models by partial differential operators L oce and L atm corresponding to the systems of equations solved by numerical models.

Traditionally, we consider the primitive equations in the ocean and the fully-compressible Euler equations in the atmosphere. On the computational domain Ω = Ω atm ∪ Ω oce (with external boundaries ∂Ω ext atm and ∂Ω ext oce , see Fig. 1), the integration over a time period

[0, T ] corresponds to    L atm U a = f atm , in Ω atm × [0, T ], B atm U a = g atm , in ∂Ω ext atm × [0, T ], F atm U a = F oa (U o , U a , R), on Γ × [0, T ], (1)
   L oce U o = f oce , in Ω oce × [0, T ], B oce U o = g oce , in ∂Ω ext oce × [0, T ], F oce U o = F oa (U o , U a , R), on Γ × [0, T ], (2)
with appropriate boundary conditions provided through the boundary operators B oce and B atm (the initialization of coupled models is an open problem that is beyond the scope here). In ( 1) and

(2), U a = (u a h , T a ) t and U o = (u o h , T o ) t are the state variables with u h the horizontal velocity and T the (potential) temperature, f atm and f oce are forcing terms. For the sake of simplicity, we focus on temperature and momentum, we do not include salinity and humidity in the continuous formulation of the problem (it is however included in the practical application discussed in Sec. 5). The oceanic domain Ω oce and the atmospheric domain Ω atm have a common interface Γ (Fig. 1). F oa is a function allowing the computation of air-sea fluxes. This function, generally based on the atmospheric surface layer similarity theory (e.g. [START_REF] Large | Surface fluxes for practitioners of global ocean data assimilation[END_REF], depends on U a and U o in the vicinity of the air-sea interface (as shown in Fig. 1), and on a set of non-turbulent radiative fluxes R.

In ( 1)-( 2), the interface operators are defined as

F atm • = ρ a K a ∂ z •, F oce • = ρ o K o ∂ z •,
where z is positive upward, ρ a and ρ o are the densities of the fluids, and

K a =   K a m K a m c a p K a t   , K o =   K o m K o m c o p K o t   .
We note K a m and K o m the eddy viscosities, and K a t and K o t the eddy diffusivities. The c p terms correspond to the specific heat of the fluid.

In forced mode, U o (resp. U a and R) in (1) (resp. ( 2)) is provided offline by existing satellitebased or reanalysis products. In coupled mode, both models are run either simultaneously or successively on the same space-time interval. In this case, the consistency required at the air-sea interface is the continuity of vertical fluxes1

ρ a K a m ∂ z u a h = ρ o K o m ∂ z u o h = τ on Γ × [0, T ] ρ a c a p K a t ∂ z T a = ρ o c o p K o t ∂ z T o = Q net on Γ × [0, T ] Q net = R + Q S
where the surface wind stress τ and the sensible heat flux Q S (consequently the net heat flux Q net ) are computed using the function F oa = (τ , Q net ) t previously introduced. The turbulent air-sea fluxes are given by a parameterization of the atmospheric surface layer. They usually take the form

τ = ρ a C D ∆U ∆U, Q S = ρ a c a p C H ∆U ∆T, (3)
where C D and C H are exchange coefficients that depend on surface roughness and local stability. ∆U (resp. ∆T) correspond to the velocity (resp. temperature) jump across the air-sea interface which is defined, in a bulk way, as the region between the lowest vertical level in the atmospheric model and the shallowest vertical level in the oceanic model. At this point we have defined all the necessary notations to formulate the coupled problem under study:

Find U a and U o that satisfy

               L atm U a = f atm in Ω atm × [0, T ] B atm U a = g atm in ∂Ω ext atm × [0, T ] L oce U o = f oce in Ω oce × [0, T ] B oce U o = g oce in ∂Ω ext oce × [0, T ] F atm U a = F oce U o on Γ × [0, T ] = F oa (U o , U a , R) (4)
for given initial and boundary conditions. It is worth mentioning that [START_REF] Lions | Mathematical theory for the coupled atmosphere-ocean models[END_REF] proved the existence of a global weak solution to problem (4) with interface conditions given by (3), considering that L atm and L oce are the primitive equations of the atmosphere and the ocean.

This paper is organized as follows. The existing coupling methods currently in use to solve (4) are presented in Sec. 2. In Sec. 3, we provide critical comments about these methods and we introduce a new theoretical framework to solve the coupling problem. In Sec. 4 we present the formulation of our AOCM based on the Weather Research and Forecasting (WRF, Skamarock and Klemp, 2008) model and the Regional Oceanic Modeling System (ROMS, [START_REF] Shchepetkin | The Regional Oceanic Modeling System (ROMS): A split-explicit, free-surface, topography-following-coordinate ocean model[END_REF]. In Sec. 5 we assess the sensitivity of this model to the algorithmic approach for the simulation of a tropical cyclone. We show that a mathematically consistent algorithm leads to reduced stochasticity, i.e., less sensitivity of the coupled solution to perturbations in the initial condition and coupling frequency. We summarize and discuss our findings in Sec. 6.

Inria

Ω oce t i-1 t i-2 t i Ω atm t i+1 F oa (hU oce i i-1 , U atm | i ) F oa (hU o i i-2 , U a | i-1 ) i-1 F oa (hU o i i-3 , U a | i-2 ) i-2 F oa (hU o i i-1 , U a | i ) i F oa (hU o i i-3 , U a | i-2 ) F oa (hU o i i-2 , U a | i-1 ) . i . i-2 . i-1 i i -1 i -2
Figure 2: Schematic view of the coupling by time windows. Three time windows, numbered from i -2 to i, are considered here. The arrows correspond to an exchange of information between the two models. The operator . i denotes a temporal average over the time window i.

A Classification of Existing Methods

Coupling methods can be categorized in two different groups : 1-way algorithms, when one model receive informations but does not send anything back to the other (e.g. [START_REF] Muller | High-resolution atmospheric forcing for regional oceanic model: the iroise sea[END_REF], and 2-way algorithms, when a feedback is introduced. We focus here on the latter. The usual coupling strategies can be found in [START_REF] Bryan | The ncar csm flux coupler[END_REF] mostly for global configurations and in [START_REF] Bao | Numerical simulations of air-sea interaction under high wind conditions using a coupled model : A study of hurricane development[END_REF] and [START_REF] Perlin | Numerical simulation of air-sea coupling during coastal upwelling[END_REF] for regional high-resolution studies. A first algorithmic approach used for global problems is based on the exchange of averaged-in-time fluxes between the models (usually referred to as asynchronous coupling) whereas a second one deals with instantaneous fluxes (usually referred to as synchronous coupling). In this section, we expose the key differences between the two strategies.

Asynchronous Coupling by Time Windows (Based on Mean Fluxes)

Asynchronous coupling is the strategy used in most climate models of the IPCC 2 . For this method, the total simulation time [0, T ] is split into M smallest time windows

[t i , t i+1 ], i.e. [0, T ] = ∪ M i=1 [t i , t i+1 ].
The length of those time windows are typically between 1 hour and 1 day depending on target applications and the need to resolve or not the diurnal cycle. On a given time window, atmospheric and oceanic models only exchange time-averaged quantities. It has the advantage of requiring few communications between models. Noting . i a temporal average on the time window [t i , t i+1 ], the coupling algorithm is

L atm U a = f atm in Ω atm × [t i , t i+1 ] F atm U a = F oa ( U o i-1 , U a , R) on Γ × [t i , t i+1 ] then L oce U o = f oce in Ω oce × [t i , t i+1 ] F oce U o = F atm U a i on Γ × [t i , t i+1 ]
(5) This algorithm is schematically described in Fig. 2. First the atmospheric model is advanced from t i to t i+1 using the averaged ocean state computed on the previous time window. Then, 2 Intergovernmental Panel on Climate Change RR n°8651 the fluxes used to force the atmospheric model are averaged and applied to the oceanic model on the same time interval (the fluxes are generally piecewise constant on each time window). This methodology ensures that over a time interval [t i , t i+1 ] both models are forced by the exact same mean fluxes, implying a strict conservation of the quantities. Indeed, in both models the time integral of surface fluxes over the simulation equals F atm U a [0,T ] . However, the solution of algorithm ( 5) is not rigorously solution of the original problem (4) because there is a synchronicity issue. The modification of the oceanic state U o on [t i , t i+1 ] is not provided to the atmospheric component on the proper time interval

[t i , t i+1 ] but only on [t i+1 , t i+2 ].
It is usually assumed that the loss of synchronization does not significantly impair the coupled solution, but we will show on the contrary that algorithm (5) leads to significant numerical errors (Sec. 33.2).

Synchronous Coupling at the Time Step Level (Based on Instantaneous Fluxes)

In nature, the ocean and the atmosphere continuously exchange fluxes on scales ranging from global to micro scales. Therefore, proper coupling frequency between numerical models should be as small as possible, typically the largest time step between the oceanic and the atmospheric model. In this regard, a relevant algorithm would consist in exchanging instantaneous fluxes.

Let us consider that the oceanic time step ∆t o is such that ∆t o = N ∆t a . The corresponding algorithm reads

L atm U a = f atm in Ω atm × [t i , t i + N ∆t a ] F atm U a = F oa (U o (t i ), U a (t), R(t)) on Γ × [t i , t i + N ∆t a ] L oce U o = f oce in Ω oce × [t i , t i + ∆t o ] F oce U o = F oa (U o (t i ), U a (t i ), R(t i )) on Γ × [t i , t i + ∆t o ] (6)
The oceanic and atmospheric components are integrated forward for a time period corresponding to ∆t o (or equivalently N ∆t a ). Data exchange of instantaneous values is then performed and model integration continues for another ∆t o period of time. This process is repeated until the final forecast time (Fig. 3). In (6) the oceanic component receives instantaneous values from the atmosphere but integrated values can also be considered between t i and t i + N ∆t a to avoid aliasing errors (Fig. 3,b). Both choices raise either conservation, aliasing or synchronization problems. In addition, algorithm ( 6) is difficult to implement efficiently from computational and numerical view points. Code communications are extremely frequent and time integration schemes must be carefully considered for consistent interfacial conditions.

At first glance, (6) may appear as a solution of the full problem (4). This is, however, formally true only in the limit ∆t o , ∆t a → 0. Indeed, in most numerical models vertical diffusion is treated implicitly-in-time; air-sea fluxes are thus provided at time t i + ∆t rather than t i as in (6). As for algorithm (5), the loss of mathematical consistency due to numerical implementation may lead to significant errors (Sec. 33.2).

Inria

Ω oce Ω atm t t t + ∆t o t + N ∆t a t +2N ∆t a t +2∆t o F oa F oa F oa Ω oce Ω atm t t t + ∆t o t + N ∆t a t +2N ∆t a t +2∆t o F oa F oa a) b)

Numerical and Physical Considerations

In the previous section, we showed that the numerical methods generally used to solve the coupled problem (4) are not fully satisfying from a mathematical point of view. In this section, we specificy further the ocean-atmosphere coupling problem based on physical considerations, then we propose an alternative method to solve it.

Physical Constraints

Coupling methods used in the context of AOCMs differ from the methods usually used for multiphysics (fluid-fluid) coupling. However, the latter can not be easily applied to AOCMs because they often assume a full representation of all the spatio-temporal scales involved in the problem.

The oceanic and atmospheric variability is spectrally broad band across all scales from micro O(10 -3 m, 10 -3 s) to global O(10 7 m, 10 10 s) scales. AOCMs contain essential parameterization to account for unresolved processes. The fluxes exchanged by the models are not the result of a discrete derivative in the vicinity of the air-sea interface but are given by atmospheric surface layer parameterizations based on the so-called bulk aerodynamic formulae. Bulk formulations are symbolically represented by the function F oa in (4). They are defined and calibrated semiempirically using measurements averaged in time over about an hour or more, and for a restricted range of stability values [START_REF] Large | Surface fluxes for practitioners of global ocean data assimilation[END_REF]. There is little knowledge and observations of air-sea fluxes at high temporal scales (see discussion in [START_REF] Danabasoglu | Diurnal coupling in the tropical oceans of ccsm3[END_REF], Sec. 2) and the sign of airsea fluxes is uncertain on time scales less than 10 minutes. In addition, high frequency physical processes are associated with the wavy boundary layer that requires specific time-dependent equations to be solved [START_REF] Bao | Numerical simulations of air-sea interaction under high wind conditions using a coupled model : A study of hurricane development[END_REF]. The wavy boundary layer tends to act as a low-pass filter on air-sea exchanges. Therefore and for all reasons mentioned above, we consider mean hourly fluxes preferable when using bulk formulations (see [START_REF] Large | Surface fluxes for practitioners of global ocean data assimilation[END_REF], for a discussion). A timeaveraging procedure, acting as a smoother of smallest time scales, is thus adopted in algorithm (5). We consider that algorithm ( 6) is relevant only if additional physical processes predominant on short time-scales are explicitly addressed in the flux computation.

Numerical Considerations

We showed in Sec. 2 that usual coupling methods used in AOCMs are not entirely satisfactory with respect to consistency, conservation or synchronization. An illustration of numerical errors is given here based on a simplified one-dimensional coupling problem with an interface at z = 0

:    ∂ t q 2 -∂ z (ν 2 ∂ z q 2 ) = f 2 , in ]0, L 2 [×[0, T ], q 2 (L 2 , t) = 0, t ∈ [0, T ], ν 2 ∂ z q 2 (0, t) = ν 1 ∂ z q 1 (0, t), t ∈ [0, T ], (7)
   ∂ t q 1 -∂ z (ν 1 ∂ z q 1 ) = f 1 , in ] -L 1 , 0[×[0, T ], q 1 (-L 1 , t) = 0, t ∈ [0, T ], q 1 (0, t) = q 2 (0, t), t ∈ [0, T ]. (8)
The complete setup for this simplified test-case is described in App. A, and a thorough mathematical and numerical study of this problem can be found in Lemarié et al. (2013b,c). We impose continuity of tracer q and its vertical flux at the interface. Fig. 4 (left panel) shows the ℓ 2 -norm of errors associated with synchronous and asynchronous methods. Note that we chose a testcase forced with right hand sides f 1 and f 2 which are identical whatever the coupling method. That is the reason why the coupling error does not grow significantly with time. It must be clear that setting f 1 = f 2 = 0 in (7-8) would lead to much larger errors constantly growing with time. Note that, whatever the method used, the error ℓ 2 -norm is initially zero because the initial condition satisfies the original equation.

Based on purely numerical arguments, the synchronous method is preferable over the asynchronous one as it leads to much smaller error. However, because the synchronous method does not satisfy the aforementioned physical constraints (Sec. 33.1) for coupling and also because of computational burden, we strive to propose an alternative approach which would theoretically prove more consistent.

Global in Time Schwarz Method

The Schwarz-like domain decomposition methods (see [START_REF] Gander | Schwarz methods over the course of time[END_REF], for a review) are widely used for coupling problems with different physics and/or different numerical treatment. Originally introduced for stationary problems, those methods have been recently extended to time-dependent problems to provide a global-in-time Schwarz method, a.k.a. Schwarz waveform relaxation (e.g. [START_REF] Gander | Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF]. The idea is to separate the original problem on Ω = Ω atm ∪ Ω oce into subproblems on Ω atm and Ω oce , which can be solved separately. An iterative process is then applied to achieve convergence to the solution of the original problem. The main drawback of this approach is the iterative procedure which increases the computational cost of coupling, especially when convergence is slow (note that there is currently an active research aiming at optimizing the convergence speed of Schwarz-like methods; see discussion in Sec. 6). These methods have already been applied to oceanic models for improving Poisson and Helmholtz solvers, open boundary specification and nesting techniques [START_REF] Debreu | On the schwarz alternating method for oceanic models on parallel computers[END_REF][START_REF] Blayo | Nesting ocean models[END_REF][START_REF] Cailleau | Comparison of different numerical methods used to handle the open boundary of a regional ocean circulation model of the bay of biscay[END_REF]. They also seem well suited to our ocean-atmosphere coupling problem (4). Using the notations introduced previously, the iterative algorithm on a time window [t i , t i+1 ] can be written as follows (for a given initial condition at t = t i ) :

Loop over k until convergence L atm U a k = f atm , in Ω atm × [t i , t i+1 ] F atm U a k = F oa (U o k-1 , U a k , R k ), on Γ × [t i , t i+1 ] L oce U o k = f oce , in Ω oce × [t i , t i+1 ] F oce U o k = F atm U a k , on Γ × [t i , t i+1 ] (9)
where the subscripts k denotes the iteration number. The first guess U o k=0 on Γ × [t i , t i+1 ] is generally taken from the converged solution on the previous time window [t i-1 , t i ]. The two models, at each iteration, are run successively: this is the so called multiplicative form of the algorithm. If the condition

F oce U o k = F atm U a k is replaced by F oce U o k = F atm U a
k-1 both models can be run in parallel over the whole time window [t i , t i+1 ]: this is the parallel form of the algorithm. When convergence is reached, this algorithm gives the exact solution to (4). Note that with algorithm (9) the solution is independent of the size of time window [t i , t i+1 ], as opposed to the asynchronous coupling method (akin to performing only one iteration of the Schwarz algorithm). However, for physical constraints on high-frequency treatment mentioned in Sec. 33.1, algorithm (9) should be modified to include time-averaging of the quantities near the air-sea interface as in (5):

Loop over k until convergence L atm U a k = f atm , in Ω atm × [t i , t i+1 ] F atm U a k = F oa U o k-1 i , U a k , R k , on Γ × [t i , t i+1 ] L oce U o k = f oce , in Ω oce × [t i , t i+1 ] F oce U o k = F atm U a k i , on Γ × [t i , t i+1 ] (10)
Fig. 4 (right panel) shows that numerical errors using algorithm (10) are significantly reduced compared to the asynchronous method (5). Moreover, we numerically checked on a simple test case that algorithm (10) converges.

Partial Conclusion

We have so far discussed numerous subtleties of numerical and physical nature involved in the design of AOCMs. Before proceeding to a real-case study, we draw a few remarks based on our survey of the algorithmic aspects of ocean-atmosphere coupling. Looking at (5) and ( 10), it appears that the asynchronous coupling method currently in use in global climate models corresponds to one iteration of the multiplicative form of a global-in-time Schwarz algorithm [START_REF] Lemarié | Algorithmes de Schwarz et couplage océan-atmosphère[END_REF]. In this respect, asynchronous coupling is mathematically inconsistent because it does not give the solution of the coupling problem (4) but an approximation. It can also easily be shown that the synchronous coupling ( 6) is equivalent to one iteration of a local-intime Schwarz algorithm [START_REF] Cai | Local multiplicative schwarz algorithms for steady and unsteady convection-diffusion equations[END_REF]. As described above, it is possible to satisfy the required consistency at the expense of an iterative process.

To avoid the burden of the iteration process associated with Schwarz methods, a monolithic scheme (i.e. a single model including both the oceanic and atmospheric physics) may be used, but that would impose the ocean and atmosphere to be advanced on the same horizontal grid with the same time-step. This is unnatural in the case of space-time multi-physics problems. It must be clear that the Schwarz method does not interfere with the time-scale of physical processes. Convergence of the Schwarz method does not imply any instantaneous physical adjustment, e.g. of the oceanic boundary layer to the overlying atmospheric conditions. Its mode of action is on the mathematical consistency of the solution at the air-sea interface. It should not be Inria confused with other iterative methods sometimes used in convective adjustment schemes or bulk formulations.

When using an iterative method, there are intertwined concerns : the computational cost and the convergence speed. The present paper is a preliminary study where attention is given to the relevancy of consistent methods, rather than their computational cost. However, one important question that first needs addressing is whether the Schwarz method actually converges in a practical application. We address this question in the next section by applying the multiplicative Schwarz algorithm to model the genesis and movement of tropical cyclone (TC) Erica in the southwest Pacific ocean during March 2003.

Model description and Experimental Setup

In this section we briefly present the numerical models composing our AOCM and our coupling strategy based on algorithm (10) and applied to tropical cyclone Erica.

Oceanic Model

The numerical oceanic model is ROMS [START_REF] Shchepetkin | The Regional Oceanic Modeling System (ROMS): A split-explicit, free-surface, topography-following-coordinate ocean model[END_REF] in its AGRIF-IRD3 version, see [START_REF] Shchepetkin | Correction and commentary for "Ocean forecasting in terrain-following coordinates: Formulation and skill assessment of the Regional Ocean Modelling System[END_REF] for a description of the various ROMS kernels. ROMS is a split-explicit time stepping, hydrostatic, Boussinesq, free-surface primitive equation model specifically designed for regional applications. ROMS equations are formulated using a generalized terrain-following σ-coordinate, that can be configured to enhance resolution in the surface boundary layer at the air-sea interface. Our experiment is setup in a configuration with four open boundaries, 50 σ-coordinate levels (grid parameters are chosen to ensure that the surface grid-box depth is at most 1 m), a horizontal grid with ∆x = 1/3 o , and ∆t = 1800s. The domain roughly extends from 6 • S to 25 • S in latitude and from 142.5 • E to 172.5 • E in longitude (Fig. 5). Note that special care must be given to the numerical schemes for tracer transport to properly simulate this area of complex bathymetry [START_REF] Marchesiello | Spurious diapycnal mixing in terrain-following coordinate models: The problem and a solution[END_REF]. Vertical mixing of tracers and momentum to predict K o m (z) and K o t (z) is done with the KPP parameterization [START_REF] Large | Oceanic vertical mixing: A review and a model with a nonlocal boundary layer parameterization[END_REF]. The coupled simulations are conducted without any flux correction scheme nor sea-surface temperature or salinity restoring. Model initialization results from a ten-year spin-up simulation forced by climatological atmospheric fluxes provided by the Comprehensive Ocean-Atmosphere Dataset (COADS) fields and by the QuikSCAT Climatology of Ocean Winds (QuikCOW). The oceanic boundary conditions are interpolated from a climatology based on the SODA reanalysis (the forcing fields are constructed using the Romstools utilities, [START_REF] Penven | Software tools for pre-and postprocessing of oceanic regional simulations[END_REF]). Note that the spin-up simulation develops large intrinsic variability at the mesoscale that is uncorrelated at this scale with the actual motion (no data assimilation techniques were used). Our goal here is to document the solution sensitivity to the coupling methodology rather than provide the best possible hindcast.

Atmospheric Model

The atmospheric model used in our experiment is the WRF-ARW4 solver (Skamarock and [START_REF] Inria Skamarock | A time-split nonhydrostatic atmospheric model for weather research and forecasting applications[END_REF]. WRF integrates the fully compressible nonhydrostatic Euler equations formulated using a terrain-following mass vertical coordinate. The model grid has a horizontal resolution of 1/3 • with 31 vertical levels and the time step is 180s. The meteorological data used for model initialization and boundary conditions are the NCEP2 reanalysis5 . Note that there was no need for bogus injection to initialize the generation of TC Erica, as is sometimes done for cyclone studies, because the initial perturbation was captured by the NCEP2 reanalysis. The physical options used for the present study are the WSM3 (WRF Single-Moment 3-class) scheme for microphysics, the Rapid Radiative Transfer Model for longwave radiation, the Dudhia shortwave radiation scheme, the 5-layer thermal-diffusion land-surface model, and the Betts-Miller-Janjic cumulus parameterization. The Planetary Boundary Layer scheme used to compute K a m (z) and K a t (z) is a non-local K-profile scheme (the Yonsei University (YSU) scheme, [START_REF] Hong | A new vertical diffusion package with an explicit treatment of entrainment processes[END_REF]). Similar WRF and ROMS configurations were used in an uncoupled mode in [START_REF] Jourdain | Mesoscale simulation of tropical cyclones in the south pacific : climatology and interannual variability[END_REF] and [START_REF] Jullien | Impact of tropical cyclones on the heat budget of the South Pacific Ocean[END_REF] respectively to study the statistics of cyclonic activity over the South-Pacific, and the oceanic response to tropical cyclones. A following ROMS-WRF coupled study using our coupling method has recently been addressed [START_REF] Jullien | Ocean feedback to tropical cyclones: climatology and processes[END_REF]. The present study is focused on coupling rather than TC science.

i-1 t i-2 t i t i+1 i i -1 i -2 time F oa k=1 k=3 k=2 k=1 k=2 k=3 k=2 k=1 k=3

Model Coupling Strategy

We implement a global-in-time Schwarz method to couple WRF and ROMS. In practice, the different steps corresponding to algorithm (10) on a given time window [t i , t i+1 ] and for iteration k are :

1. Compute the atmospheric solution from t i to t i+1 using

U o k-1 (or U o [ti-1,ti] for k = 1)
2. Send averaged air-sea fluxes F atm U a k i on [t i , t i+1 ] to the oceanic model 3. Integrate the oceanic model on the same time period (using F atm U a k (t i ) and F atm U a k i , fluxes from t i to t i+1 are linearly reconstructed; see Fig. 6)

Send the newly computed U o

k to the atmospheric model.

Steps 1 to 4 are iteratively applied until convergence (or until a fixed number of iterations is attained). On each time-window [t i , t i+1 ] the initial condition at t = t i corresponds to the converged solution from the previous time-window

[t i-1 , t i ].
As mentioned earlier, one sequence of steps 1 to 4, without iterating, corresponds to (5).

In our practical implementation, the turbulent components of air-sea fluxes are computed in the WRF surface layer scheme based on a classical similarity theory which uses stability functions from [START_REF] Paulson | The mathematical representation of wind speed and temperature profiles in the unstable atmospheric surface layer[END_REF], [START_REF] Dyer | Flux-gradient relationships in the constant flux layer[END_REF][START_REF] Dyer | Flux-gradient relationships in the constant flux layer[END_REF][START_REF] Webb | Profile relationships: The log-linear range, and extension to strong stability[END_REF] to compute the surface transfer coefficients for heat, moisture, and momentum. For our experiments we consider the same horizontal grids in WRF and ROMS. This is not a limitation of our method but a choice to ensure that the focus is on coupling rather than interpolation/extrapolation errors.

Inria

Experimental Setup

Erica, a category 4 cyclone (on the Australia and Fiji scale), was generated in March 2003 off Australia and reached New Caledonia a few days later causing human and material damage on the island. This event was poorly forecast as various attempts gave a totally different track than observed. As a result the local alert system failed. Cyclone track prediction has improved in the last decade but large uncertainties remain, in part due to the neglect of ocean and atmosphere coupling [START_REF] Bender | Numerical simulations of tropical cyclone-ocean interaction with a high-resolution coupled model[END_REF]. Latent heat release to the atmosphere is the source of cyclone intensity but TC-induced cooling of the surface ocean is a powerful negative feedback affecting both intensity and trajectory [START_REF] Jullien | Impact of tropical cyclones on the heat budget of the South Pacific Ocean[END_REF][START_REF] Jullien | Ocean feedback to tropical cyclones: climatology and processes[END_REF]. Representing the correct two-way interaction between sea surface temperature (SST) and tropical cyclones is thus of primary importance to cyclogenesis.

The present study is designed to check the viability of the Schwarz method in a strong coupling case and to address two important questions, namely the convergence properties of the method in a fully realistic framework and the impact of the coupling method on the coupled solution. Our analysis focuses on the inter-comparison between coupling strategies rather than the direct comparison with observations. The coupled model is not optimally tuned for that. The default bulk coefficients in WRF are not set for extreme events and wave age effects are neglected in the absence of a coupled wave model. In addition, TC intensity is underestimated here as horizontal grid spacing of more than 2-3 km cannot explicitly solve eyewall dynamics (including vortex Rossby waves and mesovortices) with consequences on the intensity of ascent surrounding the eye, warming of the core and associated cyclonic intensity [START_REF] Gentry | Sensitivity of simulated tropical cyclone structure and intensity to horizontal resolution[END_REF][START_REF] Hill | Analysis of idealized tropical cyclone simulations using the weather research and forecasting model: Sensitivity to turbulence parameterization and grid spacing[END_REF]. Nevertheless, our simulations of Erica are realistic enough for the task at hand as it reproduces the observed intensity and track with some accuracy (see Fig. 11).

Ensemble Design and Simulation Strategy

Let us now describe the ensemble design. We restrict our model comparison to that between the asynchronous (5) and the Schwarz (10) methods. The synchronous coupling method (at the time-step level) is cumbersome and tedious to implement. In addition, as explained earlier, current parameterizations are inadequate for this type of coupling.

AOCMs tends to exhibit strong sensitivity to model parameters like the initial condition and/or the coupling frequency (Sec. ??). To check whether part of this sensitivity is related to coupling consistency errors (rather than stochasticity), we design two ensembles of 18 members, one for each coupling method. The 18 members are generated through perturbations of initial conditions and coupling frequency. We consider time windows of 3 hours and 6 hours, and three different initial conditions are chosen for the atmospheric and the oceanic models. The atmospheric model is started either from Feb. 28, March 1 or March 2, corresponding to distinct synoptic setting (Fig. 7). As for the oceanic model, three initial conditions for the month of March are selected from the 10 year climatological spin-up. All coupled simulations cover the duration of cyclonic event, ending in March 16, 2003. For the ensemble integrated using the Schwarz method we systematically proceed to M = 9 iterations. This is done to avoid choosing an arbitrary stopping criterion, and to keep the computational cost at a reasonable level.

Figure 7: Initial conditions for the surface winds on 28 Feb. 2003(a), 01 Mar. 2003 (vectors, b), and 02 Mar. 2003 (vectors, c). Difference in wind speed between Mar. 01 and Feb. 28 (shaded,b), and between Mar. 02 and Feb. 28 (shaded, c) Inria 5 Numerical results

Convergence Properties

First, it is natural to check the behavior of the coupled solution with respect to the iterations when the Schwarz method is used. To do so, we introduce a convergence rate R k using the sea surface temperature T o (z = 0) averaged in time on a given time window. R k is defined as the ratio between the ℓ 2 -norm of the error at two successive iterations, considering that the solution at iterate k = M is the true solution of the problem : andn x (resp. n y ) the number of grid points in the zonal (resp. meridional) direction. For the algorithm to converge R k must remain smaller than unity, and the smaller R k the faster the convergence. In this study, it turned out that four iterations were usually sufficient for convergence (Fig. 8). However, the convergence speed of the method should depend on the model configuration and grid resolution. Finer resolution would lead to more energetic nonlinear effects which may impair the convergence speed of the method. These complications are left for a future study.

R k = e k 2 e k-1 2 = T o k (z = 0) -T o M (z = 0) 2 T o k-1 (z = 0) -T o M (z = 0) 2 , e k 2 = nx×ny i=1 T o k i (z = 0) -T o M i (z = 0) 2 with • 2 the ℓ 2 -norm,
From a theoretical point of view, it is straightforward to show that Schwarz methods, and more generally iterative methods, are slower to converge in presence of low frequencies in the error as long as the problem under consideration is only weakly nonlinear. We can thus anticipate that the main differences in the coupled solutions integrated with Schwarz and asynchronous methods occur at low frequencies. Fig. 9 shows the temporal spectrum

Sp ω { u a h (z = 10 m) k=M -u a h (z = 10 m) asyn }
of the difference between the coupling methods ( u a h asyn being the solution obtained with the asynchronous method, i.e. one iteration only), with ω the frequency. As expected, the effect of iterations is primarily visible for low time-frequencies, which suggests that the Schwarz method would have a significant impact on large-scale and climate processes. Note that it is theoretically possible to design a Schwarz method which could converge equally fast at all frequencies (this is the so-called Optimized Schwarz Method). However, this method is tractable only for relatively simple academic problems and is not straightforward to generalize to the ocean-atmosphere coupling problem (see discussion in Sec. 6).

Ensemble Spread

We checked that the Schwarz method converges at least in one realistic application of our coupled model system. Let us now illustrate the impact of the method on the robustness of the coupled solutions. Our analysis will focus on the cyclone track and intensity of each member of the two ensembles corresponding to the two coupling algorithms. We identify the cyclone track using minimum pressure at the first vertical level of the atmospheric model. The track envelop shown in Fig. 10 approximately follows the observed track with motion speed and underlying sea surface temperature in the same range (Fig. 11). It is striking that the dispersion of trajectories is significantly smaller for the ensemble based on the Schwarz method. In this case, the maximum standard deviation is 164 km and the standard deviation is around 92 km compared to respectively 212 km and 125 km with the asynchronous method. Note that the track spread in both ensembles tend to reduce when approaching the south-east corner of the domain, presumably because all simulations share identical boundary conditions.

The same behavior as for cyclone track is manifest in time histories of the maximum surface wind speed (not shown), with a maximum deviation of 6.2 m s -1 for the asynchronous method compared to 3.7 m s -1 for the Schwarz method (3.7 m s -1 versus 2.2 m s -1 for the standard deviation). Therefore, two very important features of tropical cyclone events, namely their trajectory and intensity, are more robustly simulated when improving the consistency of coupling method.

Discussion and Conclusion

We have emphasized in this paper the key role of the coupling algorithm in the design of an atmospheric and oceanic coupled model. Very popular coupling methods used in regional and global climate models do not provide the exact solution to the ocean-atmosphere coupling problem (4), but an approximation. We introduced a natural and non-intrusive method to solve this problem and we showed its relevance. This method called Global-in-Time Schwarz Method is based on an iterative process. It can easily be shown that the usual synchronous (6) and asynchronous (5) methods prominent in climate models correspond to only one iteration of a Schwarz algorithm. However, the iterative procedure can substantially increase the computational cost of the AOCM (by a factor corresponding to the number of iterations to reach convergence).

A realistic application is produced by an ensemble of simulations through perturbations of the initial conditions and coupling frequency of a regional AOCM. The results suggest that part of the model sensitivity to perturbed parameters can be attributed to inaccuracies in the coupling method. Specifically, coupling inconsistencies can spuriously increase the physical stochasticity of atmospheric and oceanic events (materialized here by the ensemble spread). For our particular case, three iterations of the Schwarz method are sufficient to improve the coupled solutions with respect to the ensemble spread. We show that the iterative process should significantly modify the low-frequency component of the solution, whereas high-frequencies converge very rapidly. Therefore, we anticipate that the iterative procedure would be particularly relevant to long-term climate studies. However, it remains to be seen how robust is the convergence regardless of the model formulation, particularly the boundary layer parameterizations at the air-sea interface. In this regard, the Schwarz algorithm can also be used as diagnostic tool to assess consistency between atmospheric and oceanic boundary layer parameterizations: two schemes could be recognized as consistent if they lead to a converging Schwarz algorithm.

For some applications, the computational cost of the Schwarz algorithm may be unaffordable and one would thus proceed to only one iteration of the method, thereby retrieving the asynchronous method. Based on our study and on our earlier theoretical work (Lemarié et al., 2013a,b,c), we can formulate some recommendations when the asynchronous coupling method is used :

• The coupling frequency (i.e. the size of time window between two exchanges) can be adequately set to minimize coupling errors and that of surface flux parameterizations.

Because the Schwarz algorithm converges slowly at low-frequency, the shorter the time window the smaller the coupling error. However, because of large uncertainties in the specification of surface fluxes, the larger the time window the smaller the parameterization error. With those constraints and the requirement of a proper representation of the diurnal cycle, a coupling frequency between 1 and 3 hours would provide a good compromise. In addition, it is generally observed that the wave field acts as a low-pass filter on air-sea exchanges, with a cutoff around 2 hours, which suggests that a time window of 2 hours may be considered as a relevant choice.

• Lemarié et al. (2013b,c) show that the mathematical type of interface condition has great impact on the convergence speed of the Schwarz algorithm. In the case of a diffusion problem, replacing Dirichlet-Neumann 6 conditions by a Robin (a.k.a. Fourier) interface condition (i.e., a linear combination of Dirichlet and a Neumann conditions) can dramatically improve the convergence speed if the weights in the linear combination are properly chosen. Even in the case of one single iteration of the method (i.e., the asynchronous method) coupling errors can be reduced by simply working on the type of interface condition. However, the theoretical problem is not straightforward for AOCMs because of bulk formulations. It is a work in progress.

A One-Dimensional Diffusion Test Problem

In this section, we compare different coupling methods using a simple test case, which is assumed representative of an ocean-atmosphere type problem. To do so, we define two subdomains

Ω 1 = ]-L 1 , 0[ and Ω 2 =]0, L 2 [ with L 1 = L 2 = 250 m.
In this example, we consider the one dimensional diffusion equation of a scalar quantity q

∂ t q -∂ z (ν∂ z q) = f, (11)
with ν a diffusion coefficient such that ν = ν 1 in Ω 1 and ν = ν 2 in Ω 2 . We choose ν 1 = ν 2 to model the heterogeneous physical properties between the two subdomains. Uniqueness of solutions for the coupling problem is obtained with Dirichlet-Neumann conditions at z = 0 (i.e. we require the equality of the subproblems solutions and of their fluxes). The corresponding coupling problem reads   

∂ t q 2 -∂ z (ν 2 ∂ z q 2 ) = f 2 , in ]0, L 2 [×[0, T ], q 2 (L 2 , t) = q 2 (t), t ∈ [0, T ], ν 2 ∂ z q 2 (0, t) = ν 1 ∂ z q 1 (0, t), t ∈ [0, T ], (12)
   ∂ t q 1 -∂ z (ν 1 ∂ z q 1 ) = f 1 , in ] -L 1 , 0[×[0, T ], q 1 (-L 1 , t) = q 1 (t),
t ∈ [0, T ], q 1 (0, t) = q 2 (0, t), t ∈ [0, T ],

(13) for a given initial condition. We pose the problem to get analytical solutions of the form q 2 (z, t) = q 0 8 3 -exp -z α 2 1 + cos 2 πt T q 1 (z, t) = q 0 8 1 + exp z α 1 1 + cos 2 πt T where the condition ν 2 α 1 = ν 1 α 2 is sufficient to ensure the proper regularity of the solution across the interface at z = 0. To obtain such analytical solutions, we set the right hand sides f 1 and f 2

f 1 (z, t) = - q 0 8 π T sin 2πt T 1 + exp z α 1 + ν 1 α 2 1 exp z α 1 1 + cos 2 πt T , f 2 (z, t) = - q 0 8 π T sin 2πt T 3 -exp - z α 2 - ν 2 α 2 2 exp - z α 2 1 + cos 2 πt T .
Equation ( 11) is discretized using a backward Euler scheme in time and a second-order scheme on a staggered grid in space. The parameter values are α 1 = 50, α 2 = 10, q 0 = 15, T 0 = 24 × 3600 s, ν 1 = 1 m s -1 , ν 2 = 0.2 m s -1 , ∆z = 1 m, ∆t = 900 s, and the total simulation time is T = 2000 days. To numerically solve the coupling problem (12-13), we use the six different methods described in the paper : the asynchronous method with instantaneous and averaged fluxes (the fluxes are reconstructed using a linear conservative scheme), the global-in-time Schwarz method with instantaneous and averaged fluxes, the synchronous method and the local-in-time Schwarz method. The global-in-time Schwarz method and the asynchronous method are integrated with time windows of 6 hours. As expected, the Schwarz methods (local or global in time) with instantaneous fluxes provide strictly the same solution q ⋆ (z, t)(z ∈ Ω 1 Ω 2 , t ∈ [0, T ]) which is used as a reference to compute the ℓ 2 -norm of the error for other methods

ε(t i ) 2 = k=N k=1 |q(z k , t i ) -q ⋆ (z k , t i )| 2 ,

Inria

with N the number of grid points between -L 1 and L 2 . The time evolution of ε 2 is represented in Figure 4, and discussed in Sec. 33.2 and 33.3.

Figure 1 :

 1 Figure 1: Schematic view of the two non-overlapping subdomains Ω atm and Ω oce . The mathematical notations correspond to those introduced in Sec. 1.2. The boundary layers encompassing the air-sea interface Γ and the atmospheric surface layer (gray shaded) are of particular importance for the coupling problem. The quantity ∆U = U a (z a ) -U o (z o ) is directly involved in the computation of air-sea fluxes, with z a (resp. z o ) the location of the lowest (resp. shallowest) vertical level in the atmospheric (resp. oceanic) model.

Figure 3 :

 3 Figure 3: Example of two coupling strategies at the time step level. ∆t o and ∆t a denote the (baroclinic) time steps respectively of the oceanic and the atmospheric models, with ∆t o = N ∆t a (N = 6 here). The arrows represent an exchange of information with the surface layer parameterization function F oa . For the atmospheric component this exchange is based on instantaneous values for algorithm a) and on time integrated values for algorithm b).

Figure 4 :

 4 Figure4: Time history of the error ℓ 2 -norm for different coupling methods : synchronous and asynchronous methods (left), asynchronous and global-in-time Schwarz methods with flux averaging (right). Note that, whatever the method used, the error ℓ 2 -norm is initially zero because the initial condition satisfies the original equation.

Figure 5 :

 5 Figure 5: Snapshots (March 12, 2003 at 8 p.m. GMT) of (a) ROMS sea surface temperature (b) WRF 10 meter winds during a coupled simulation.

  t

Figure 6 :

 6 Figure 6: Piecewise linear reconstruction of the air-sea fluxes for the oceanic model in the case of N = 3 iterations. The initial value (grey circles) and the average on a given time-window are sufficient to reconstruct the linear function.

Figure 8 :Figure 9 :Figure 10 :Figure 11 :

 891011 Figure 8: Ensemble mean of the evolution of the error ℓ 2 -norm (left) and of the convergence rate R k with respect to the iterates k in the case of time windows of 3 hours and 6 hours. The computation of R k and e k 2 is done using sea surface temperature.

3

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Coupling Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2 A Classification of Existing Methods 7 2.1 Asynchronous Coupling by Time Windows (Based on Mean Fluxes) . . . . . . . 7 2.2 Synchronous Coupling at the Time Step Level (Based on Instantaneous Fluxes) . 8 Numerical and Physical Considerations 10 3.1 Physical Constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.2 Numerical Considerations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.3 Global in Time Schwarz Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.4 Partial Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 4 Model description and Experimental Setup 14 4.1 Oceanic Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 4.2 Atmospheric Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 4.3 Model Coupling Strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 4.4 Experimental Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 4.5 Ensemble Design and Simulation Strategy . . . . . . . . . . . . . . . . . . . . . . 17 5 Numerical results 19 5.1 Convergence Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 5.2 Ensemble Spread . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

RR n°8651

Here, we assume wind-wave equilibrium, i.e., the atmospheric momentum flux to the wave field is immediately transferred to the ocean through wave breaking. The addition of a wave component to the coupling system would act as a low-pass filter to air-sea exchanges (see Sec.33.1 and Sec.6) RR n°8651

Adaptive Grid Refinement in Fortran-Institut de Recherche pour le Développement http://www.romsagrif. org/

Advanced Research WRF, http://www.mmm.ucar.edu/wrf/users/

National Center for Environmental Prediction, http://nomad3.ncep.noaa.gov/ncep_data/index.html Inria

A Dirichlet condition amounts to specifying a value in the interface, whereas a Neumann condition amounts to specifying a flux. Inria