Xiaoxuan Liu

Yang Wang

Chao Chen

Aura Tintaru

Yu Cao

Juan Liu

Fabio Ziarelli

Jingjie Tang

Hongbo Guo

Roseline Rosas

Suzanne Giorgio

Laurence Charles

Palma Rocchi

Ph.D Ling Peng
email: ling.peng@univ-amu.fr

A Fluorinated Bola-Amphiphilic Dendrimer for On-Demand Delivery of siRNA, via Specific Response to Reactive Oxygen Species

Keywords: bola-amphiphile, stimuli-responsive, siRNA delivery, gene therapy, 19 F-NMR

published or not. The documents may come

INTRODUCTION

Molecular science has revolutionized the central paradigm of drug delivery, especially for establishing smart or intelligent materials to deliver therapeutic agents on-demand. [1][2][3] In particular, many of the advances in the design of delivery systems for anticancer therapeutics has been inspired by a growing understanding of tumor microenvironments and exploitation of specific characteristics of cancer cells and subtle difference in tumor lesions. 1,2,4 Tumor microenvironments are characterized by a variety of atypical features such as abnormal tumor vasculature, absence of lymphatic drainage, altered redox environment, hypoxia and lower pH gradient etc. 5 It is of note that an intrinsic factor specifically linked to cancer cells is the high level of reactive oxygen species (ROS), as ROS are constantly generated within cancer cells because of the highly stressful environment caused by rapid and uncontrollable cancer cell proliferation. 6,7 Despite the obvious interest in harnessing this factor for specific targeting in cancer treatment, few strategies have been explored to develop delivery systems that are capable of ROS-triggered controlled release, in particular for small interfering RNA (siRNA) therapeutics. 8,9 Therapeutics based on siRNA provide an enormous opportunity for cancer treatment by virtue of the ability of siRNA to specifically and efficiently turn off the expression of genes associated with cancer development and drug resistance. [10][11][12] However, the main challenge facing siRNA therapeutics is their safe and efficient delivery. 13,14 This is because siRNA is too negatively charged to spontaneously cross biomembranes; at the same time, siRNA is vulnerable to nuclease attack. Consequently, safe and efficient on-demand delivery vectors, which can prevent siRNA degradation and convey functional siRNA into cancer cells, are in high demand.

During the past years, various natural and synthetic delivery systems have been developed for siRNA delivery. Among these, the most effective are viral vectors. However, because of concerns about the immunogenicity and safety of viral delivery systems, there is an urgent demand for developing the alternative nonviral vectors. 13,14 The two main classes of nonviral vectors are based on lipids and polymers. [14][15][16][17] Dendrimers, a special type of synthetic polymer, have recently emerged as a promising delivery platform for siRNA therapeutics by virtue of their well-defined and precisely-controlled molecular structures as well as the unique multivalent cooperativity confined within a nanoscale volume. 18,19 In particular, amphiphilic dendrimers with judiciously tailored hydrophilic and hydrophobic components are able to benefit from the combined advantages of lipid and dendrimer vectors for effective siRNA delivery. [20][21][22] Here, we report an innovative bola-amphiphilic dendrimer bola4A for on-demand delivery of siRNA in response to the high level of ROS in cancer cells (Figure 1). This dendrimer features a bola-lipid/dendrimer hybrid harboring a ROS-sensitive thioacetal group at the central hydrophobic core and two poly(amidoamine) (PAMAM) dendrons as the polar terminals (Figure 1A). The PAMAM dendrons, which bear amine terminals and are positively charged at physiological pH, have been devised to interact with the negatively charged siRNA via electrostatic interactions, 23 while the hydrophobic "bola-lipid" core scaffold is to mimic the robust and strong assembly properties of bola-amphiphiles observed in extremophile archae bacteria. 24 The thioacetal entity at the bola-lipid core is supposed to decompose upon exposure to the high level of ROS in cancer cells, 8,9 thereby promoting siRNA unpacking for potent gene silencing within cancer cells (Figure 1B). Also of note is the presence of fluorine tags within this vector, which allow tracking of the ROS-responsive delivery process by 19 F-NMR. 25 Last but not least, this distinctive and ingenious bola-amphiphilic vector, by combining the advantages of both lipid and dendrimer vectors, will provide a new perspective on the design of ROS-responsive vectors for targeted siRNA delivery.

RESULTS AND DISCUSSION

Dendrimer bola4A is readily synthesized via "click" chemistry. We synthesized bola4A according to the plan illustrated in Scheme 1 (and in Scheme S1 in supplementary materials).

The bola-core 1 was prepared via the condensation of 3,5-difluorobenzaldehyde with the freshly prepared corresponding azido-bearing thiol in the presence of the Lewis acid catalyst trifluoroboronetherate. 26,27 Using the Cu-catalyzed Huisgen "Click" reaction, 20,21,28 we successfully conjugated the bola-core 1 with the alkynyl-terminated dendron 2 to yield 3. The bola-dendrimer 3 was subsequently subjected to amidation with ethylenediamine to provide the desired bola4A. After purification via dialysis, bola4A was obtained in excellent purity with yield exceeding 90%.

Bola4A responds to ROS. With the bola-dendrimer bola4A in hand, we investigated its response to ROS under simulated conditions, that is, in the presence of H2O2 (Figure 2). We incubated bola4A with H2O2 and found that the sharp 19 F-NMR signal originating from the intrinsic fluorine atoms of bola4A ( = -110.5 ppm) disappeared rapidly, while new 19 F-NMR peaks associated with the degradation products of bola4A progressively appeared (Figure 2A and Figure S1). This finding implies that bola4A was decomposed following exposure to H2O2.

Subjecting the same samples to electrospray ionization mass spectrometry (ESI-MS) further confirmed the degradation of bola4A upon treatment with H2O2 (Figure 2B). On the one hand, signals assigned to bola4A in multiple protonation states (from 3+ to 6+) were observed to decrease and disappear as a function of the H2O2 exposure time (Figure 2B). On the other hand, new signals detected during the H2O2 treatment indicated the decomposition of bola4A and the formation of degradation products (Figure 2B). Based on the elemental composition derived from the accurate mass measurements (Table S1), structures could be proposed for these degradation products (Figure 2C), highlighting the effective disintegration of the thioacetal function after exposure to H2O2. Collectively, these results demonstrate that bola4A is readily decomposed under ROS conditions, and hence possesses favorable properties for potential controlled release in response to ROS.

Bola4A forms nanoparticles with siRNA, protects siRNA and promote cellular uptake.

For siRNA delivery, it is important that the delivery vector is able to bind and condense the siRNA into nanosized complexes and protect it from degradation before promoting its cellular uptake. With this in mind, we examined the formation of siRNA/bola4A complexes using gel shift analysis. As shown in Figure 3A, bola4A was able to form stable complexes with siRNA and completely retard the migration of siRNA on the gel at an N/P ratio ≥ 2.5. Further results from transmission electron microscopy (TEM) and scanning electron microscopy (SEM) revealed that the resulting siRNA/bola4A complexes formed uniform, compact and spherical nanoparticles (Figures 3B and3C). Additional dynamic light scattering (DLS) analysis confirmed that the siRNA/bola4A nanoparticles were well dispersed with a size average of 50 nm (Figure 3D). The surface charge of these nanoparticles was characterized by a -potential of +28 mV, implying stable colloidal nanoparticles. Indeed, the so-formed siRNA/bola4A nanoparticles effectively protected the siRNA from degradation by the enzyme RNase (Figure 3E), further indicating the formation of stable siRNA/bola4A complexes. Finally, these siRNA/bola4A nanoparticles were rapidly and efficiently internalized by cells (Figure 3F), a benefit and advantageous prerequisite for effective siRNA delivery.

Bola4A/siRNA complexes are responsive to ROS. Importantly, the siRNA/bola4A nanoparticles were readily disassembled in response to high ROS levels in cancer cells. We first studied the siRNA/bola4A nanoparticles in response to H2O2, which simulates ROS conditions. As shown by TEM imaging, the siRNA/bola4A nanoparticles collapsed upon exposure to H2O2, accompanied by a significant change in morphology, suggesting a ROStriggered disassembly of the siRNA/bola4A complexes (Figure 4A and Figure S2). In line with the TEM imaging, further results from DLS analysis (Figure S3) also indicated the destruction of the siRNA/bola4A nanoparticles upon treatment with H2O2.

Moreover, we examined the ROS-triggered decomposition of bola4A once the siRNA/bola4A complexes were internalized into ROS-abundant human prostate cancer PC-3 cells. To do this, we used 19 F high-resolution magic angle spinning (HRMAS) NMR, a nondestructive method for in situ analysis of 19 F-containing compounds within intact cells. 19 F HRMAS NMR is particularly powerful and sensitive for investigating fluorinated bola4A in cells since there are no endogenous fluorinated compounds in cells, and hence no 19 F-NMR signal in the control background. 25,29,30 As shown in Figure 4B, only weak 19 F-NMR signals were observed in ROS-abundant prostate cancer PC-3 cells, although FACS flow cytometry demonstrated that the siRNA/bola4A complexes were rapidly and effectively taken up by PC-3 cells (Figure 3F). In contrast, a sharp, strong and singular 19 F-NMR signal was detected in ROS-depleted PC-3 cells which were pretreated with the antioxidant N-acetyl cysteine (NAC) (Figure 4C). It is of note that treating PC-3 with the antioxidant NAC effectively downregulated the ROS level significantly (Figure 5A). In addition, the chemical shift of the 19 F-NMR signal detected in the NAC-treated PC-3 cells was a perfect match with that of bola4A, indicating that bola4A maintained intact in ROS-poor PC-3 cells pretreated with the antioxidant NAC. Similar results were also observed with ROS-poor Chinese hamster ovary (CHO) cells (Figure 4D). Collectively, these results indicate that bola4A was indeed decomposed or metabolized in ROS-rich PC-3 cells, but not in ROS-poor CHO or antioxidanttreated low-ROS PC-3 cells. Thus, bola4A demonstrates effective ROS-responsiveness for cancer cells.

Bola4A-mediated specific ROS-responsive delivery of siRNA and gene silencing.

Encouraged by the ROS-responsive properties of our bola4A dendrimer, we evaluated its ability to deliver siRNA and inhibit gene expression in two high-ROS cell lines, human prostate cancer PC-3 cells and breast cancer MCF-7 cells, and three low-ROS cell lines, human embryonic kidney (HEK) cells, Chinese hamster ovary (CHO) cells and antioxidant-pretreated PC-3 cells (Figure 5). The siRNA molecules used in this study were devised to target either heat shock protein 27 (Hsp27) 31,32 or translationally controlled tumor protein (TCTP) 33,34 , both of which are actively involved in cancer development and drug resistance. Following bola4Amediated siRNA delivery, expression of Hsp27 and TCTP was considerably suppressed in the two high-ROS content cancer cell lines PC-3 (Figure 5B) and MCF-7 (Figure 5C), whereas no noticeable gene silencing was observed in the low-ROS content HEK (Figure 5D) and CHO cells (Figure 5E). This can be reasonably ascribed to the inherently higher level of ROS in cancer cells, which leads to the decomposition of bola4A and consequently the disassembly of the siRNA/bola4A complexes, thereby enhancing siRNA release and gene silencing.

We further demonstrated that down-regulation of the ROS level in PC-3 cells by pretreatment with the antioxidant NAC (Figure 5A) led to a dramatic decrease in gene silencing (Figure 5F). This provides additional evidence that gene silencing is specifically triggered by ROS. Together, our results indicate that bola4A is able to mediate specific and efficient siRNA delivery and gene silencing in response to a ROS-rich environment, in perfect agreement with our design concept for bola4A as a ROS-responsive vector.

Bola4A benefits from the integrated delivery advantages of lipid and dendrimer vectors.

In addition to the ROS-responsive feature of bola4A, we wanted to forge a strong and stable vector based on our bola-amphiphile by combining the benefits of both lipid and PAMAM dendrimer vectors for siRNA delivery. 20,21 Our bola4A is a lipid/dendrimer hybrid bearing a hydrophobic chain entity and two hydrophilic PAMAM dendrons (Figure 1A). When we compared the gene silencing activity of bola4A with mono4A (the amphiphilic dendrimer without the bola-lipid core), dendron4A (the dendron entity alone) or the bola-core 1, only bola4A was effective (Figure 6A). Thus, the unique bola-amphiphilic architecture did indeed endow bola4A with the ability to deliver siRNA effectively. It is to mention that the hydrophobic bola-core in bola4A is shorter than the native phospholipid bilayer. We designed bola4A with this unique molecular architecture to avoid insertion into the cell membrane, thus obviating deleterious effects on cell membrane integrity and cell viability. The biocompatibility of bola4A was confirmed by the absence of serum hemolysis (Figure S4) and cytotoxicity using both lactate dehydrogenase (LDH) and MTT tests (Figure S5). We further examined bola4A-mediated siRNA delivery and gene silencing in the presence of dioleoylphosphatidylethanolamine (DOPE). DOPE is a fusogenic lipid, which promotes fusion and is frequently used to enhance the delivery efficacy of lipid vectors. Our results showed that gene silencing was significantly increased in the presence of DOPE (Figure 6B), confirming that bola4A is endowed with the delivery characteristics of lipid vectors.

We also surmised that bola4A might profit from the "proton sponge effect" 35 of PAMAM dendron entities for effective nucleic acid delivery. The proton sponge phenomenon occurs for PAMAM dendrimers in acidic environments such as endosomes, and is thought to be important for release of the cargo from endosomes into the cytoplasm. The tertiary amine groups in the interior of the PAMAM dendron are ready to mop up protons within endosomes (Figure S6), leading to an ionic imbalance which results in endosomal lysis and cargo release. We used the proton pump inhibitor bafilomycin A1 to impede endosome acidification, in order to test whether endosomal acidification affects bola4A-mediate siRNA delivery and gene silencing. Indeed, treatment with bafilomycin A1 significantly lowered the level of Hsp27 expression (Figure 6C), suggesting that acidic endosomes are necessary for effective siRNA release and delivery by bola4A. This result is consistent with the hypothesis that the PAMAM dendrons act as a proton sponge, and that this property is important for siRNA delivery.

CONCLUSION

In conclusion, we have established a ROS-sensitive bola-amphiphilic dendrimer bola4A as an innovative on-demand vector for stimulus-responsive siRNA delivery and gene silencing.

The so-devised bola4A is also able to benefit from the combined advantages of lipid and dendrimer vectors. The distinctive ROS-sensitive thioacetal motif within bola4A allows efficient disassembly of the siRNA/bola4A complexes under ROS-rich conditions for effective siRNA delivery and potent gene silencing in cancer cells. In addition, the presence of fluorine atoms within this vector makes it possible to study the ROS-responsive delivery process by 19 F-NMR. Collectively, our results show that this ROS-sensitive bola-amphiphilic dendrimer offers a unique opportunity to achieve controllable release of siRNA for effective gene silencing in cancer cells by capitalizing on the high level of intracellular ROS. Our study provides a new perspective on the design of tailor-made stimulus-responsive materials for ondemand drug delivery.

ASSOCIATED CONTENT

Supporting information: supplementary figures, materials and methods, dendrimer synthesis and characterization as well as all experimental protocols for NMR, MS, TEM, DLS, cell uptake, siRNA delivery, gene silencing etc. This information is available free of charge via Internet.

Figure legends:

Scheme 1: Synthesis of the bola-amphiphilic dendrimer bola4A. S1), of degradation products of bola4A formed upon treatment with H2O2. S1), of degradation products of bola4A formed upon treatment with H2O2.

Figure 1 :

 1 Figure 1: The bola-amphiphilic dendrimer bola4A, designed for ROS-responsive siRNA delivery. (A) Chemical structure of bola4A. (B) Schematic representation of bola4A for the ROS-triggered delivery of siRNA and consequential gene silencing. Bola4A dendrimers form complexes with siRNA molecules, which can be internalized by the cancer cell before releasing siRNA in response to ROS, leading to effective gene silencing.

Figure 2 :

 2 Figure 2: Study of bola4A decomposition upon treatment with 200 mM H2O2 (simulated ROS conditions) for 0, 2 and 24 h, using (A) 19 F-NMR and (B) ESI-MS analysis. (C) Proposed structures, based on accurate mass measurements (TableS1), of degradation products of

Figure 3 :

 3 Figure 3: Bola4A is able to form nanocomplexes with siRNA, protect siRNA from degradation and promote cellular uptake. (A) Agarose gel migration of siRNA (200 ng per well) in the presence of bola4A dendrimer at N/P charge ratios of 1/5 -10/1. (B) TEM image of the siRNA/bola4A complexes using 5 ng/µL siRNA and bola4A at an N/P ratio of 10. (C) SEM image of the siRNA/bola4A complexes using 5 ng/µL siRNA and bola4A at an N/P ratio of 10. (D) Size distribution of the siRNA/bola4A complexes (at an N/P ratio of 10 with 1 M siRNA) determined using DLS. (E) Protection of siRNA by bola4A against enzymatic degradation. Compared to the naked siRNA (200 ng per well) which was degraded within 5 min in the presence of RNase, siRNA complexed with bola4A at an N/P ratio of 10 was resistant to RNase and remained stable even after 1 h incubation. (F) Uptake of Dy647-labeled Hsp27 siRNA/bola4A complexes (at an N/P ratio of 10 with 50 nM siRNA) by human prostate cancer PC-3 cells evaluated using flow cytometry. Experiments were carried out in triplicate.

Figure 4 :

 4 Figure 4: (A) TEM images of siRNA/bola4A complexes (N/P=10) before and after incubation with H2O2 (1.06 M) for 2, 4 and 24 h at 37 ºC. Scale bars indicate 200 nm. 19 F HRMAS NMR recording of bola4A at different time points (0, 1, 2 and 4 h) in (B) ROS-rich normal prostate cancer PC-3 cells, (C) ROS-poor PC-3 cells pretreated with the antioxidant N-acetyl cysteine (NAC) and (D) ROS-poor Chinese hamster ovary (CHO) cells after treatment with the

Figure 5 :

 5 Figure 5: Bola4A-mediated specific and efficient ROS-responsive siRNA delivery and gene silencing. (A) ROS levels in Chinese hamster ovary (CHO) cells, human embryonic kidney (HEK) cells, breast cancer MCF-7 cells, prostate cancer PC-3 cells, and PC-3 cells pretreated with the antioxidant N-acetyl-cysteine (NAC) (10 mM) quantified using CellROX ® orange reagent by flow cytometry. Bola4A-mediated siRNA delivery and gene silencing in ROSabundant (B) PC-3 cells and (C) MCF-7 cells, as well as in ROS-poor (D) HEK cells, (E) CHO cells and (F) PC-3 cells pretreated with NAC (50 nM siRNA at N/P=10). SiRNAs targeting heat shock protein 27 (Hsp27) and translationally controlled tumor protein (TCTP) were used.

Figure 6 :BFigure 1 :

 61 Figure 6: Bola4A-mediated siRNA delivery benefits from both the distinctive bolaamphiphilic structure and the delivery features of lipid and dendrimer vectors. (A) Compared to bola4A, neither the amphiphilic dendrimer mono4A without the bola-lipid core, nor the dendron entity dendron4A, nor the bola-core 1 led to any gene silencing (50 nM siRNA at N/P=10). (B) Dioleoylphosphatidylethanolamine (DOPE) enhanced the bola4A-mediated siRNA delivery and gene silencing (20 nM siRNA at N/P=10). (C) Bafilomycin A1 decreased the bola4A-mediated gene silencing (50 nM siRNA at N/P=10). PC-3 cells and Hsp27 siRNA were used.

Figure 2 :

 2 Figure 2: Study of bola4A decomposition upon treatment with 200 mM H2O2 (simulated ROS conditions) for 0, 2 and 24 h, using (A) 19 F-NMR and (B) ESI-MS analysis. (C) Proposed structures, based on accurate mass measurements (TableS1), of degradation products of bola4A formed upon treatment with H2O2.

Figure 3 :

 3 Figure 3: Bola4A is able to form nanocomplexes with siRNA, protect siRNA from degradation and promote cellular uptake. (A) Agarose gel migration of siRNA (200 ng per well) in the presence of bola4A dendrimer at N/P charge ratios of 1/5 -10/1. (B) TEM image of the siRNA/bola4A complexes using 5 ng/µL siRNA and bola4A at an N/P ratio of 10. (C) SEM image of the siRNA/bola4A complexes using 5 ng/µL siRNA and bola4A at an N/P ratio of 10. (D) Size distribution of the siRNA/bola4A complexes (at an N/P ratio of 10 with 1 M siRNA) determined using DLS. (E) Protection of siRNA by bola4A against enzymatic degradation. Compared to the naked siRNA (200 ng per well) which was degraded within 5 min in the presence of RNase, siRNA complexed with bola4A at an N/P ratio of 10 was resistant to RNase and remained stable even after 1 h incubation. (F) Uptake of Dy647-labeled Hsp27 siRNA/bola4A complexes (at an N/P ratio of 10 with 50 nM siRNA) by human prostate cancer PC-3 cells evaluated using flow cytometry. Experiments were carried out in triplicate.

Figure 4 :

 4 Figure 4: (A) TEM images of siRNA/bola4A complexes (N/P=10) before and after incubation with H2O2 (1.06 M) for 2, 4 and 24 h at 37 ºC. Scale bars indicate 200 nm. 19 F HRMAS NMR recording of bola4A at different time points (0, 1, 2 and 4 h) in (B) ROS-rich normal prostate cancer PC-3 cells, (C) ROS-poor PC-3 cells pretreated with the antioxidant N-acetyl cysteine (NAC) and (D) ROS-poor Chinese hamster ovary (CHO) cells after treatment with the siRNA/bola4A complexes.

Figure 5 :

 5 Figure 5: Bola4A-mediated specific and efficient ROS-responsive siRNA delivery and gene silencing. (A) ROS levels in Chinese hamster ovary (CHO) cells, human embryonic kidney (HEK) cells, breast cancer MCF-7 cells, prostate cancer PC-3 cells, and PC-3 cells pretreated with the antioxidant N-acetylcysteine (NAC) (10 mM) quantified using CellROX ® orange reagent by flow cytometry. Bola4Amediated siRNA delivery and gene silencing in ROS-abundant (B) PC-3 cells and (C) MCF-7 cells, as well as in ROS-poor (D) HEK cells, (E) CHO cells and (F) PC-3 cells pretreated with NAC (50 nM siRNA at N/P=10). SiRNAs targeting heat shock protein 27 (Hsp27) and translationally controlled tumor protein (TCTP) were used.

Figure 6 :

 6 Figure 6: Bola4A-mediated siRNA delivery benefits from both the distinctive bola-amphiphilic structure and the delivery features of lipid and dendrimer vectors. (A) Compared to bola4A, neither the amphiphilic dendrimer mono4A without the bola-lipid core, nor the dendron entity dendron4A, nor the bola-core 1 led to any gene silencing (50 nM siRNA at N/P=10). (B) Dioleoylphosphatidylethanolamine (DOPE) enhanced the bola4A-mediated siRNA delivery and gene silencing (20 nM siRNA at N/P=10). (C) Bafilomycin A1 decreased the bola4A-mediated gene silencing (50 nM siRNA at N/P=10). PC-3 cells and Hsp27 siRNA were used.

ACKNOWLEDGEMENTS

Financial support from La Ligue Nationale Contre le Cancer (LP), Association pour la Recherche sur les Tumeurs de la Prostate (LP, XL), Association Franç aise contre les Myopathies (XL), Fondation pour la Recherche Mé dicale (YC), the international ERA-Net EURONANOMED European Research project "Target4Cancer" (LP), China Scholarship Council (YW, CC, JL, JT), PACA Canceropôle, INCa, CNRS and INSERM is gratefully acknowledged. We thank Serge Netische and Damien Chaudanson for SEM experiments,

AUTHOR INFORMATION

LP conceived the project. XL, YW, LP designed experiments, YW, CC, YC and JT synthesized the dendrimer, AT, FZ, RR and LC executed MS and NMR analysis, CC, JT, SG, HG performed TEM and DLS experiments, XL carried out all biological experiments with the help of PR and JL. XL and LP wrote the manuscript with comments from all the other authors.

REFERENCES:

Scheme 1: Synthesis of the bola-amphiphilic dendrimer bola4A.


Guillaume Bertrand

Thomas Gutierrez

Bruno Tourlière

Jérémie Melleton

Eric Gloaguen

Florent Cheval-Garabédian

Field validation of mineral prospectivity approaches, a first test of the CBA method

Field validation of mineral prospectivity approaches, a first test of the CBA method

Bertrand Guillaume 1,2 , Gutierrez Thomas 1 , Tourlière Bruno 1 , Melleton Jérémie 1 , Gloaguen Eric 1,2 and Florent Cheval-Garabédian 2 1 -BRGM, Orléans, France 2 -ISTO, UMR7327, Orléans, France Numerous methods of mineral prospectivity mapping have been developed during the last decades and it is often difficult to evaluate their reliability and adequacy. Some objective a posteriori approaches allow calculating their performance such as, for instance, AUC values on ROC curves. However, we believe that field control is a necessary step, especially in the context of mineral targeting in greenfield exploration.

Mineral prospectivity mapping generally relies on unequivocal associations between discrete data (e.g., known deposits) and polygons (e.g., lithology, geophysics, geochemistry, etc.). Consequently, the quality of results obtained from 'classical' methods (e.g., Weight of Evidence, Fuzzy Logic, Logistic Regression, Neural Network, etc.) strongly depends on the accuracy of input data (location of points and contours of polygons) which, in many cases, may be questionable. To address this issue of geographic inaccuracy of map elements, BRGM has developed the CBA (Cell Based Association) method. Its base principle is not to rank polygons in which known deposits are located, but to identify favorable association neighboring these deposits (i.e. in cells of a regular grid containing them).

In order to assess the adequacy of CBA for mineral targeting in early exploration phases, we have done a first field validation campaign for Sb in the Vendée region (western France). 107 soil samples have been collected (auger drill) in both favorable and unfavorable cells (according to CBA results) and across know Sb-bearing veins. They have been analyzed at ultra-trace level using ICP-MS/ICP-AES (detection range of 0.05 to 10 000 ppm for Sb). In this contribution, we present an overview of the CBA method, its resulting prospectivity map for Sb in Vendée, the sampling and analyzing procedure, the preliminary results we have obtained and the first conclusion they allow.


Hortense Laeuffer
email: hortense.laeuffer@u-bordeaux.fr

Brice Guiot

Jean-Christophe Wahl

Nicolas Perry
email: nicolas.perry@ensam.eu

Florian Lavelle

Christophe Bois
email: christophe.bois@u-bordeaux.fr

RELATIONSHIP BETWEEN DAMAGE AND PERMEABILITY IN LINER-LESS COMPOSITE TANKS: EXPERIMENTAL STUDY AT THE LAMINA LEVEL

The aim of this study is to provide a relevant description of damage growth and the resultant crack network to predict leaks in liner-less composite vessels. Tensile tests were carried out on three different laminates: [0 2 /90 n /0 2 ], [+45/-45] 2s and [0/+67.5/-67.5] s . Number n varies from 1 to 3 in order to study the effect of ply thickness. Transverse crack and delamination at crack tips were identified with an optical microscope during tensile loading. A length of 100 mm was observed for several loading levels to evaluate statistical effects. Results highlight a preliminary step in the damage scenario with small crack densities before a second step where the crack growth speeds up. In bulk, cross-section examinations showed that no delamination occurred at crack tip in the material of the study (M21 T700). Cross-section examinations were also performed on [+45/-45] 2s and [0/+67.5/-67.5] s layups in order to bypass the issue of free edge effects. Damage state in those layups was shown to be significantly different in the bulk than at the surface. Observations of the damage state in bulk for those layups demonstrated that there is no transverse crack in [+45/-45] 2s specimens subjected to shear strains up to 4%, and that interactions between damage of consecutive plies strongly impact both the damage kinetics and the arrangement of cracks. These elements are fundamental for the assessment of permeability performance, and will be introduced in the predictive model.

INTRODUCTION

Designing liner-less composite vessels for launch vehicles enables to save both cost and weight. Therefore, this is a core issue for aerospace industry. In usual composite vessels, the liner is a metallic or polymer part in charge of the gas barrier function. One challenge when designing a liner-less vessel is to reach the permeability requirement with the composite wall itself. Pristine composite laminates meet the permeability requirement, but as those materials are heterogeneous, damage growth may occur in cases of low thermo-mechanical loads. Transverse cracks and micro-delamination in adjacent plies grow and may connect together (Fig. 1), resulting in a leakage network through the composite wall.

On the ply scale, transverse cracks are generated by transverse stress (mode I) and shear stress (mode II). Developing a model that predicts damage densities and arrangement for any laminate requires reliable experimental data for both modes. The aim of this study is to provide a relevant description of damage growth and the resultant network to predict leaks.

The first part of this experimental study describes the method for damage characterisation. The second part presents the results obtained on the evolution of damage densities. Finally, the third part is devoted to the morphology of the damage network.

METHOD FOR DAMAGE CHARACTERISATION

Damage description

The meso-damage state of each ply of a laminate can be described by two damage densities, as illustrated in Fig. 2: the crack density ρ, which is the average number of transverse cracks over the observed length L, and the micro-delamination length µ, which is the average length of microdelamination at each crack tip. The corresponding dimensionless variables are defined by: the crack rate :

ρ " ρh " N L h (1)
the delamination rate :

μ " µρ " µ N L (2)
where h denotes the ply thickness and N the number of cracks observed on the length L. Crack and micro-delamination rates are average variables, so that the damage stage is considered homogeneous in each ply. When used as damage variables in a model, the damage densities define the crack pattern which is commonly considered as periodic [START_REF] Ladevèze | Relationships between 'micro' and 'meso' mechanics of laminated composites[END_REF]. Experimental values and evolution of damage densities can be obtained by observation by optical microscopy [START_REF] Huchette | Sur la complémentarité des approches expérimentales et numériques pour la modélisation des mécanismes d'endommagement des composites stratifiés[END_REF][START_REF] Malenfant | Étude de l'influence de l'endommagement sur la perméabilité des matériaux composites, application à la réalisation d'un réservoir cryogénique sans liner[END_REF][START_REF] Bois | A multiscale damage and crack opening model for the prediction of flow path in laminated composite[END_REF] or X-ray tomography as well.

Instrumentation and method

Transverse cracks are generated by transverse stress (mode I) and shear stress (mode II). This requires testing different laminates in order to characterise damage growth in the ply in both modes or in mixed mode I + II. In this study, three different laminates were tested: [0 2 /90 n /0 2 ], [+45/-45] 2s and [0/+67.5/-67.5] s . The effect of ply thickness is also studied by making number n vary from 1 to 3. In the present work, the assessment of the damage densities was carried out by optical microscopy for several loading levels. Tensile test were performed on specimens polished on one edge. Transverse crack and delamination at crack tips were identified with a travelling optical microscope (Fig. 3). The observation area was chosen quite large, i.e. about 100 mm, in order to evaluate statistical effects.

The material is carbon fibre and thermoset matrix M21 T700. Plate samples were manufactured by Automated Fibre Placement (AFP). The thickness of the elementary layer is 0.26 mm. Several layers of the same orientation can be stacked together in order to obtain a thicker ply, e.g. in a [0 2 /90 3 /0 2 ] laminate the thickness of the 90 0 ply is 0.78 mm. A specific protocol was applied in order to bypass the issue of free edge effects. Observation of the surface under loading was combined to cross-section examination : polishing of the edge was performed after loading in order to remove from 15 µm to 2 mm of the edge surface and observe the damage state in bulk. After polishing, the sample was loaded to a lower level so that damage did not propagate but cracks opened and became more easy to distinguish. Those steps were repeated for increasing maximal loading to obtain the evolution laws of damage in bulk. Beside crack growth characterisation, performing several cross-section examinations through the width of a specimen also allows to study the arrangement of cracks in three dimensions.

EVOLUTION OF DAMAGE DENSITIES

Edge effects

Cross section examinations were carried out on damaged tensile test specimens of [0 2 /90 n /0 2 ] and [+45/-45] 2s layups. Measurements were performed on representative areas of about 50 mm long. Fig. 4 shows a transverse crack and the associated micro-delamination at crack tips on the edge of the sample before and after removing 15 microns by a polishing procedure. 15 microns under the edge surface, no visible delamination remains. Fig. 5(a) presents the evolution of micro-delamination rates before and after polishing for three ply thicknesses. Independently of the ply thickness, the micro-delamination rate tends quickly to zero. Cross section examinations were continued to verify the length of cracks. The crack rate with respect to the polishing depth is reported in Fig. 5(b). Cracks were shown to be continuous for double and triple 90 plies, while a few cracks disappeared in the [0 2 /90 1 /0 2 ] lay-up. This is consistent with edge effects in this layup. This variation of crack rate for the simple ply is about 15% of the crack rate at the surface (see Fig. 6). The [+45/-45] 2s layup was subjected to transverse strains up to ε 12 " 4%, polished and then observed. For the highest strain and after removing 1.5 mm, only a few transverse cracks remained in the central (and double) ply and in external plies. Deeper, after removing 3 mm of the surface, no crack remained at all. To exclude the eventuality of cracks being too closed to be observable, the sample was observed under loading. The results demonstrates that no transverse cracks nor micro-delamination occur in the bulk under shear load. This phenomenon can be due to the use of a toughened matrix (M21) with thermoplastic nodules. For another material, e.g. with a more fragile matrix, delamination would likely be lower in bulk than at the surface but may persist in bulk. Although no meso-damage is observable, this layup nevertheless undergoes irreversible strains and stiffness loss due to diffuse damage at the micro (fibre) scale. However, if shear loading alone does not lead to transverse cracking in this material, the shear component associated to a tensile loading will likely contribute to crack growth. This contribution may require additional testing on other layups.

y (µm) μ r0 2 90 1 0 2 s, ε " 1.5% r0 2 90 2 0 2 s, ε " 1.3% r0 2 90 2 0 2 s, ε " 1.5% r0 2 90 3 0 2 s, ε " 1.1% r0 2 90 3 0 2 s, ε " 1.5% (
y (mm) ρ r0 2 90 1 0 2 s, ε " 1.3% r0 2 90 1 0 2 s, ε " 1.5% r0 2 90 2 0 2 s, ε " 1.2% r0 2 90 3 0 2 s, ε " 1.3%
(b) Crack rate ρ vs polishing depth. Fig. 5. Damage rates with respect to the polishing depth y after the transverse strain ε was applied.

ε 22 (%) ρ r0 2 {90 1 {0 2 s edge r0 2 {90 1 {0 2 s bulk
Fig. 6. Crack rate ρ with respect to the transverse strain ε 22 at the surface of the edge and in bulk.

Evolution of transverse cracking

Evolution of transverse cracking was observed previously described in Section 2.2 for [+45/-45] 2s and [0 2 /90 n /0 2 ] layups with n varying from 1 to 3. As the first one was shown to be not subjected to transverse cracking, this section is devoted to the latter. The position of cracks measured for several loading levels on a [0 2 /90 1 /0 2 ] specimen is plotted in Fig. 7(a). The corresponding crack rate has been computed for the whole observed area and also for each separate third of the same area, see Fig. 7(b). The two figures highlight the benefit of observing a large enough and hence representative area. Depending on the piece of sample one chose to focus on, initiation and slope of crack growth may be very different, at least concerning the beginning of crack growth. It is likely due to weak areas (defect locations) that drive the position of the first cracks. The effect of defects on the damage threshold vanishes when damage increases: for higher strains, the distance between consecutive cracks becomes more homogeneous, and no major difference can be observed between the three small areas. The average crack rate computed on the whole length reveals a preliminary step with a progressive beginning of cracking.

The evolution of damage densities according to the applied strain for two ply-thicknesses are presented in Fig. 8. In all cases, the preliminary step in the damage scenario with small cracking rates before a second step where the crack growth speeds up was observed. The curves shows an increase in cracking threshold when the ply thickness decreases. This phenomenon is explained by the energy released by the fracture being lower for a thinner ply. This point is widely described in the literature [START_REF] Parvizi | Constrained cracking in glass fibre-reinforced epoxy cross-ply laminates[END_REF][START_REF] Gudmundson | Initiation and growth criteria for transverse matrix cracks in composite laminates[END_REF][START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF] and obviously makes thin ply very competitive for liner-less composite vessel. The modelling of the preliminary step is fundamental for the prediction of first leak path. This phenomenon due to the variability of the material properties can be introduced in the damage model through a probabilistic approach [START_REF] Nairn | Matrix microcracking in composites[END_REF]. Evolution laws for transverse cracks and micro-delamination were built based on Finite Fracture Mechanics and energy release rates [START_REF] Tual | Multiscale model based on a finite fracture approach for the prediction of damage in laminate composites[END_REF][START_REF] Laeuffer | A model for the prediction of transverse crack and delamination density based on a strength and fracture mechanics probabilistic approach[END_REF]. Probability density functions are defined for fracture toughness and strength threshold. A set of simulations are then performed and the average response of the model are computed according to the weight associated to the density functions. Results of the simulation for both ply thicknesses are presented in Fig. 8 by doted and dashed lines. Preliminary step in the damage scenario is well described even though the densities of probability were not accurately identified yet.

MORPHOLOGY OF THE DAMAGE NETWORK

To study the shape of a damage network, damage observation under loading and cross-section examination were applied to [0/+67.5/-67.5] s specimens. The interest of this lay-up is that transverse cracking can occur in three different and consecutive plies, leading to the creation of a network. The results presented here concern one damage network (without damage growth) observed at several polishing depth from 0.2mm to 5mm.

At the scale of the specimen, the results provide an overview of the network. The position of cracks is schematically drawn in Fig. 9 from the observations. The centre ply, ´67.5 0 , is also a double ply, and thus has a lower cracking threshold. In this ply, cracks occurred first and were continuous through the width. Conversely, cracks in the simple `67.5 0 plies were short and located around the cracks of the double ply. Hence, existence of cracks in an adjacent ply drives the position and the length of the new cracks. Moreover, the cracking threshold of the simple plies is also modified: cracks occur almost simultaneously in the three plies despite their different thicknesses. Tests on [0/+67.5/0/67.5] s specimens in which cracked plies are isolated are scheduled. Comparing the results obtained with isolated and not-isolated plies will make it possible to quantify the effect of the interaction between cracks in adjacent plies. Fig. 10 focuses on the intersection of three cracks. At the interface between two plies, delamination connects the cracks. This makes the connection area larger than the Crack Opening Displacement (COD) at crack tip and may increase the leakage rate induced by a leak path.

CONCLUSION

Transverse crack and delamination at crack tips were identified with an optical microscope during tensile loading. A length of 100 mm was observed for several loading levels. This allowed to highlight a preliminary step in the damage scenario with small crack densities and progressive growth before a second step with a steeper growth. In bulk, cross-section examinations showed that no delamination occurred at crack tip in the material of the study. Cross-section examinations were also performed for the observation of [+45/-45] 2s and [0/+67.5/-67.5] s layups in order to bypass the issue of free edge effects. Damage state in those layups was shown to be significantly different through the width of the specimens than at the surface of the edges. Particularly, there is no transverse crack in [+45/-45] 2s specimens subjected to shear strains up to 4%. It was also observed that crack growth and crack length are modified by the damage state of adjacent plies. These elements are fundamental for the assessment of permeability performance, and thus will be introduced in the model.

Predicting the percolation of the network also requires to describe the network in terms of number of connections between two adjacent plies. It could be achieved by using only crack density and ply angle, but this is not trivial any more since crack growth and crack length are modified by the damage state of neighbouring plies. Cross-section examinations give an insight into the network pattern, nevertheless this method is restrictive because it is destructive, not very accurate and time-consuming. Additional experiments involving X-ray tomography are required to accurately characterise the crack network pattern. This kind of experiment remains a challenge because of the mismatch between the size of the crack pattern (3mm in the case of the [0/+67.5/-67.5] s specimen) and the size of its elements (COD « 2µm when the specimen is unloaded). Those experiments will also allow to assess the effect of the interactions between damage of consecutive plies.

Fig. 1 .

 1 Fig. 1. Transverse crack and delamination : crack network in two damaged plies and micrograph of one transverse crack with delamination at crack tip.

3 Fig. 2 .

 32 Fig. 2. Measurement of the crack density ρ and average delamination at crack tip µ in a [0 2 /90/0 2 ] lay-up.

Fig. 3 .

 3 Fig. 3. Damage observation under tensile loading test.

Fig. 4 .

 4 Fig. 4. Cross-section examinations on [0 2 /90 2 /0 2 ]: diagram of the cross-sections and micrographs at the surface (y " 0µm) and after removing 15 microns (y " 15µm).

  a) Delamination rate μ vs polishing depth.

L

  55 σ u ε " 0.9% σ " 0.63 σ u ε " 1.0% σ " 0.69 σ u ε " 1.1% σ " 0.82 σ u ε " 1.3% σ " 0.99 σ u ε " Position of cracks observed on the edge for 5 loading levels. = 60 -90 mm L = 30 -60 mm L = 0 -30 mm L = 0 -90 mm (b) Crack density ρ at the surface.

Fig. 7 .

 7 Fig. 7. Damage measurements on a [0 2 /90 1 /0 2 ] specimen, ply thickness h 90 " 0.26 mm.

r0 2 Fig. 8 .

 28 Fig. 8. Measured and predicted crack rates.

Fig. 9 .

 9 Fig. 9. Crack network in a [0/+67.5/-67.5] s : diagram of the network and micrograph of the three damaged plies.

Fig. 10 .

 10 Fig. 10. Intersection between the cracks of three plies in a [0/+67.5/-67.5] s for several polishing depth.

ACKNOWLEDGEMENTS

The authors acknowledge the council of Region Aquitaine and the French space agency CNES for their support.


Masoud Rostami

Vladimir Zeitlin
email: zeitlin@lmd.ens.fr

Improved moist-convective rotating shallow water model 1 "Improved moist-convective rotating shallow water model and its application to instabilities of hurricane-like vortices"

Keywords: Moist Convection, Rotating Shallow Water, Tropical Cyclones, Baroclinic Instability

We show how the two-layer moist-convective rotating shallow water model (mcRSW), which proved to be a simple and robust tool for studying effects of moist convection on large-scale atmospheric motions, can be improved by including, in addition to the water vapour, precipitable water, and the effects of vaporisation, entrainment, and precipitation. Thus improved mcRSW becomes cloud-resolving. It is applied, as an illustration, to model the development of instabilities of tropical cyclone-like vortices.

Introduction

Massive efforts have been undertaken in recent years in order to improve the quality of weather and climate modelling, and significant progress was achieved. Nevertheless, water vapour condensation and precipitations remain a weak point of weather forecasts, especially long-term ones. Thus, predictions of climate models significantly diverge in what concerns humidity and precipitations [START_REF] Stevens | What climate models miss?[END_REF] . The complexity of thermodynamics of the moist air, which includes phase transitions and microphysics, is prohibitive.

That is why the related processes are usually represented through simplified parameterisations in the general circulation models.

However, the essentially non-linear, switch character of phase transitions poses specific problems in modelling the water cycle.

Parametrisations of numerous physical processes in general circulation models often obscure the role of the water vapour cycle upon the large-scale atmospheric dynamics. The moist-convective rotating shallow water (mcRSW) model was proposed recently, precisely, in order to understand this role in rough but robust terms. The model is based on vertically averaged primitive equations with pseudo-height as vertical coordinate. Instead of proceeding by a direct averaging of the complete system of equations with full thermodynamics and microphysics, which necessitates a series of specific ad hoc hypotheses, a hybrid approach is used, consisting in combination of vertical averaging between pairs of isobaric surfaces and Lagrangian conservation of the moist enthalpy [START_REF] Bouchut | Fronts and nonlinear waves in a simplified shallow-water model of the atmosphere with moisture and convection[END_REF][START_REF] Lambaerts | Simplified two-layer models of precipitating atmosphere and their properties[END_REF]. Technically, convective fluxes, i.e. an extra vertical velocity across the material surfaces delimiting the shallow-water layers, are added to the standard RSW model, and are linked to condensation. For the latter a relaxation parametrisation in terms of the bulk moisture of the layer, of the type applied in general circulation models, is used. Thus obtained mcRSW model combines simplicity and fidelity of reproduction of the moist phenomena at large scales, and allows to use efficient numerical tools available for rotating shallow water equations. They also proved to be useful in understanding moist instabilities of atmospheric jets and vortices [START_REF] Lambaerts | Moist versus dry baroclinic instability in a simplified two-layer atmospheric model with condensation and latent heat release[END_REF][START_REF] Lahaye | Understanding instabilities of tropical cyclones and their evolution with a moist-convective rotating shallow-water model[END_REF]Rostami and Zeitlin 2017;Rostami et al. 2017).

The mcRSW model, however, gives only the crudest representation of the moist convection. The water vapour can condense, but after that the liquid water is dropped off, so there are no co-existing phases and no inverse vaporisation phase transition in the model. Yet, it is rather simple to introduce precipitable water in the model, and link it to the water vapour through bulk condensation and vaporisation. At the same time, the convective fluxes present in mcRSW can be associated with entrainment of precipitable water, and its exchanges between the layers, adding more realism in representing the moist convection. Below, we will make these additions to the mcRSW model, and thus obtain an "improved" mcRSW, which we call imcRSW. We will illustrate the capabilities of the new model on the example of moist instabilities of hurricane-like vortices. Multi-layer modelling of tropical cyclones goes back to the pioneering paper [START_REF] Ooyama | Numerical simulation of the life cycle of tropical cyclones[END_REF], which had, however, a limited range due to the constraint of axisymmetry. Strictly barotropic models were also used, e.g. [START_REF] Guinn | Hurricane spiral bands[END_REF], as well as shallow water models with ad hoc parametrisations of latent heat release, e.g. [START_REF] Hendricks | Hurricane eyewall evolution in a forced shallow-water model[END_REF]. The imcRSW model is a logical development of such approach.

Derivation of the improved mcRSW

Reminder on mcRSW and its derivation

Let us recall the main ideas and the key points of derivation of the 2-layer mcRSW model. The starting point is the system of "dry" primitive equations with pseudo-height as vertical coordinate [START_REF] Hoskins | Atmospheric frontogenesis models: Mathematical formulation and solution[END_REF]. We recall that pseudo-height is the geopotential height for an atmosphere with an adiabatic lapse rate: z = z 0 1 -(p/p 0 ) R/cp , where z 0 = cpθ 0 /g, and the subscript 0 indicates reference sea-level) values. Horizontal momentum and continuity equations are vertically averaged between two pairs of material surfaces z 0 , z 1 , and z 1 , z 2 , where z 0 is at the ground, and z 2 is at the top. The pseudo-height z being directly related to pressure, the lower boundary is a "free surface" and the upper boundary is considered to be at a fixed pressure ("rigid lid"). The meanfield approximation is then applied, consisting, technically, in replacing averages of the products of dynamical variables by products of averages, which expresses the hypothesis of columnar motion. In the derivation of the "ordinary" RSW the fact that material surfaces z i , i = 0, 1, 2 are moving, by definition, with corresponding local vertical velocities w i allows to eliminate these latter. The main assumption of the mcRSW model is that there exist additional convective fluxes across z i , such that

w 0 = dz 0 dt , w 1 = dz 1 dt + W 1 , w 2 = dz 2 dt + W 2 , (1)
where W 1,2 are contributions from the extra fluxes, whatever their origin, cf. Figure 1. The resulting continuity equations for the thicknesses of the layers h 2 = z 2 -z 1 , h 1 = z 1 -z 0 are modified in a physically transparent way, acquiring additional source and sink terms:

∂ t h 1 + ∇ • (h 1 v 1 ) = -W 1 , ∂ t h 2 + ∇ • (h 2 v 2 ) = +W 1 -W 2 .
(2)

The modified momentum equations contain the terms of the form W i v at the boundaries z i of the layers. An additional assumption is, hence, necessary, in order to fix the value of the horizontal velocity at the interface. In the layered models the overall horizontal velocity, by construction, has the form

v(z) = N i=1 v i H(z i -z)H(z -z i-1 )
, where H(z) is Heaviside (step-) function. Assigning a value to velocity at z i means assigning a value to the Heaviside function at zero, where it is not defined. This a well-known modelling problem, and any value between zero and one can be chosen, depending on the physics of the underlying system. In the present case this choice would reflect the processes in an intermediate buffer layer interpolating between the main layers, and replacing the sharp interface, if a vertically refined model is used. The "asymmetric" (non-centred) assignment H(0) = 1 was adopted in previous works. The "symmetric" (centred) assignment H(0) = 1/2 will be adopted below. This choice does not affect qualitatively the previous results obtained with mcRSW, however it does affect the forcing terms in conservation laws. It corresponds to a choice of efficiency of momentum transport between the layers. In this way, the vertically averaged momentum equations become:

∂ t v 1 + (v 1 • ∇)v 1 + f k × v 1 = -∇φ(z 1 ) + g θ1 θ0 ∇z 1 + v1-v2 2h1 W 1 , ∂ t v 2 + (v 2 • ∇)v 2 + f k × v 2 = -∇φ(z 2 ) + g θ2 θ0 ∇z 2 + v1-v2 2h2 W 1 + v1 2h2 W 2 , (3)
Note that, whatever the assignment for Heaviside function, the total momentum of the two-layer system (z

1 -z 0 )v 1 + (z 2 -z 1 )v 2 is
locally conserved (modulo the Coriolis force terms). In what follows, will be assuming that W 2 = 0.

The system is closed with the help of hydrostatic relations between geopotential and potential temperature, which are used to express the geopotential at the upper levels in terms of the lower-level one:

φ(z) = φ(z 0 ) + g θ1 θ0 (z -z 0 ) if z 0 ≤ z ≤ z 1 , φ(z 0 ) + g θ1 θ0 (z 1 -z 0 ) + g θ2 θ0 (z -z 1 ) if z 1 ≤ z ≤ z 2 , (4)
The vertically integrated (bulk) humidity in each layer

Q i = zi zi-1 qdz, i = 1, 2
, where q(x, y, z, t) is specific humidity, measures the total water vapour content of the air column, which is locally conserved in the absence of phase transitions. Condensation introduces c 2017 Royal Meteorological Society Prepared using qjrms4.cls a humidity sink:

∂ t Q i + ∇ • (Q i v i ) = -C i , i = 1, 2.
(5)

In the regions of condensation (C i > 0) specific moisture is saturated q(z i ) = q s (z i ) and the potential temperature θ(z i ) + (L/cp)q s (z i )

of an elementary air mass W i dt dx dy, which is rising due to the latent heat release, is equal to the potential temperature of the upper layer θ i+1 :

θ i+1 = θ(z i ) + L cp q(z i ) ≈ θ i + L cp q(z i ), (6)
If the background stratification, at constant θ(z i ) and constant q(z i ), is stable θ i+1 > θ i , by integrating the three-dimensional equation of moist-adiabatic processes

d dt θ + L cp q = 0. (7)
we get

W i = β i C i , β i = L cp(θ i+1 -θ i ) ≈ 1 q(z i ) > 0. (8)
In this way the extra vertical fluxes in (3), (2) are linked to condensation. For the system to be closed, condensation should be connected to moisture. This is done via the relaxation parametrisation, where the moisture relaxes with a characteristic time τc towards the saturation value Q s , if this threshold is crossed:

C i = Q i -Q s i τc H(Q i -Q s i ). (9)
Essentially nonlinear, switch character of the condensation process is reflected in this parameterisation, which poses no problem in finite-volume numerical scheme we are using below. For alternative, e.g. finite-difference schemes smoothing of the Heviside function could be used. In what follows we consider the two-layer model assuming that the upper layer is dry, and even with entrainment of water from the lower moist layer, water vapour in this layer is far from saturation, so the convective flux W 2 is negligible. In this way we get the mcRSW equations for such configuration:

             ∂ t v 1 + (v 1 • ∇)v 1 + f k × v 1 = -g∇(h 1 + h 2 ) + v1-v2 2h1 βC, ∂ t v 2 + (v 2 • ∇)v 2 + f k × v 2 = -g∇(h 1 + sh 2 ) + v1-v2 2h2 βC, ∂ t h 1 + ∇ • (h 1 v 1 ) = -βC, ∂ t h 2 + ∇ • (h 2 v 2 ) = +βC, ∂ t Q + ∇ • (Qv 1 ) = -C, C = Q-Q s τc H(Q -Q s ) , (10)
where s = θ 2 /θ 1 > 1 is the stratification parameter, v 1 = (u 1 , v 1 ) and v 2 = (u 2 , v 2 ) are the horizontal velocity fields in the lower and upper layer (counted from the bottom), with u i zonal and v i meridional components, and h 1 , h 2 are the thicknesses of the layers, and we will be considering the Coriolis parameter f to be constant.

As in the previous studies with mcRSW, we will not develop sophisticated parameterisations of the boundary layer and of fluxes across the lower boundary of the model. Such parameterisations exist in the literature [START_REF] Schecter | Hurricane formation in diabatic ekman turbulence[END_REF], and may be borrowed, if necessary. We will limit ourselves by the simplest version of the exchanges with the boundary layer, with a a source of bulk moisture in the lower layer due to surface evaporation E. The moisture budget thus becomes:

∂ t Q + ∇ • (Qv 1 ) = E -C (11)
The simplest parametrisations being used in the literature are the relaxational one

E = Q -Q τ E H( Q -Q), (12)
and the one where surface evaporation is proportional to the wind, which is plausible for the atmosphere over the oceanic surface:

E ∝ |v|; (13)
The two can be combined, in order to prevent the evaporation due to the wind to continue beyond the saturation:

Es = Q -Q τ E |v|H( Q -Q). (14)
The typical evaporation relaxation time τ E is about one day in the atmosphere, to be compared with τc, which is about an hour. Thus

τ E
τc. Q can be taken equal, or close to Qs, as we are doing, but not necessarily, as it represents complex processes in the boundary layer, and can be, in turn, parametrised.

c 2017 Royal Meteorological Society Prepared using qjrms4.cls

Improving the mcRSW model

An obvious shortcoming of the mcRSW model presented above is that, although it includes condensation and related convective fluxes, the condensed water vapour disappears from the model. In this sense, condensation is equivalent to precipitation in the model. Yet, as is well-known, condensed water remains in the atmosphere in the form of clouds, and precipitation is switched only when water droplets reach a critical size. It is easy to include precipitable water in the model in the form of another advected quantity with a source due to condensation, and a sink due to vaporisation, the latter process having been neglected in the simplest version of mcRSW. We thus introduce a bulk amount of precipitable water, W (x, y, t), in the air column of a given layer. It obeys the following equation in each layer:

∂ t W + ∇.(W v) = +C -V, (15)
where V denotes vaporisation. Vaporisation can be parametrised similarly to condensation:

P = Qs -Q τv H(Qs -Q). (16)
Opposite to the condensation, vaporisation engenders cooling, and hence a downward convective flux, which can be related to the background stratification along the same lines as upward flux due to condensation:

Wv = -β * V, β * = L * Cv(θ 2 -θ 1 ) , ( 17
)
where L * is the latent heat absorption coefficient, Cv is specific heat of vaporisation. β * is an order of magnitude smaller than β. There is still no precipitation sink in (15). Such sink can be introduced, again as a relaxation with a relaxation time τp, and conditioned by some critical bulk amount of precipitable water in the column:

P = W -Wcr τp H(W -Wcr). (18)
The extra fluxes (17) due to cooling give rise to extra terms in mass and momentum equations of the model in each layer. Another important phenomenon, which is absent in the simplest version of mcRSW is the entrainment of liquid water by updrafts. This process can be modelled in a simple way as a sink in the lower-layer precipitable water equation, which is proportional, with some coefficient γ, to the updraft flux, and hence, to condensation, and provides a corresponding source of precipitable water in the upper-layer.

Including the above-described modifications in the mcRSW models, and neglecting for simplicity 1) condensation and precipitations in the upper layer, by supposing that it remains far from saturation, 2) vaporisation in the lower layer, which is supposed to be close to saturation, we get the following system of equations:

                               d 1 v 1 dt + f ẑ × v 1 = -g∇(h 1 + h 2 ) + ( βC -β * V h 1 )( v 1 -v 2 2 ), d 2 v 2 dt + f ẑ × v 2 = -g∇(h 1 + sh 2 ) + ( βC -β * V h 2 )( v 1 -v 2 2 ), ∂ t h 1 + ∇.(h 1 v 1 ) = -βC + β * V, ∂ t h 2 + ∇.(h 2 v 2 ) = +βC -β * V, ∂ t W 1 + ∇.(W 1 v 1 ) = +(1 -γ) C -P, ∂ t W 2 + ∇.(W 2 v 2 ) = +γ C -V, ∂ t Q 1 + ∇.(Q 1 v 1 ) = -C + E, ∂ t Q 2 + ∇.(Q 2 v 2 ) = V, (19)
where

d i .../dt = ∂ t ... + (v i •∇)..., i = 1, 2.
Here C is condensation in the lower layer considered to be close to saturation, W i is the bulk amount of precipitable water and Q i bulk humidity in each layer, γ is the entrainment coefficient, V is vaporisation in the upper layer, considered as mostly dry. C, V , and P obey ( 9), ( 16), ( 18), respectively. Note that if the above-formulated hypotheses of mostly dry upper layer, and almost saturated lower layer are relaxed (or get inconsistent during simulations), the missing condensation, precipitation, and vaporisation in the corresponding layers can be easily restituted according to the same rules.

Conservation laws in the improved mcRSW model

As was already said, the total momentum of the system is locally conserved in the absence of the Coriolis force (f → 0), as can be seen by adding the equations for the momentum density in the layers:

(∂ t + v 1 .∇)(h 1 v 1 ) + h 1 v 1 ∇.v 1 + f ẑ × (h 1 v 1 ) = -g∇ h 2 1 2 -gh 1 ∇h 2 -( v 1 + v 2 2 )(βC -β * V ) (20a) (∂ t + v 2 .∇)(h 2 v 2 ) + h 2 v 2 ∇.v 2 + f ẑ × (h 2 v 2 ) = -gs∇ h 2 2 2 -gh 2 ∇h 1 + ( v 1 + v 2 2 )(βC -β * V ) (20b)
The last term in each equation corresponds to a Rayleigh drag produced by vertical momentum exchanges due to convective fluxes.

The total mass (thickness) h = h 1 + h 2 is also conserved, while the mass in each layer h 1,2 is not. However, we can construct a moist enthalpy in the lower layer

m 1 = h 1 -βQ 1 -β * W 2 , (21)
∂ t m 1 + ∇.(m 1 v 1 ) = 0, i = 1, 2. ( 22
)
The inclusion of precipitable water in the upper layer in ( 21) is necessary to compensate the downward mass flux due to vaporisation.

The dry energy of the system E = dxdy(e 1 + e 2 ) is conserved in the absence of diabatic effects, where the energy densities of the layers are:

       e 1 = h 1 v 2 1 2 + g h 2 1 2 , e 2 = h 2 v 2 2 2 + gh 1 h 2 + sg h 2 2 2 .
In the presence of condensation and vaporisation, the energy budget changes and the total energy density e = e 1 + e 2 is not locally conserved, acquiring a sink/source term:

∂ t e = -∇ • fe -(βC -β * V )g(1 -s)h 2 , ( 23
)
where fe is the standard energy density flux in the two-layer model. For the total energy E = dxdy e of the closed system we thus get

∂ t E = (βC -β * V )g(s -1) dxdy h 2 . ( 24
)
For stable stratifications s > 1, the r.h.s. of this equation represents an increase (decrease) of potential energy due to upward (downward) convective fluxes due to condensation heating (vaporisation cooling). Note that with "asymmetric" assignment of Heaviside function at zero, an extra term corresponding to kinetic energy loss due to Rayleigh drag would appear in the energy budget, cf [START_REF] Lambaerts | Simplified two-layer models of precipitating atmosphere and their properties[END_REF].

Potential vorticity (PV) is an important characteristics of the flow. In the presence of diabatic effects it ceases to be a Lagrangian invariant, and evolves in each layer as follows:

d 1 dt ( ζ 1 + f h 1 ) = ( ζ 1 + f h 1 ). (βC -β * V ) h 1 + ẑ h 1 • ∇ × v 1 -v 2 2 . βC -β * V h 1 , (25a)
d 2 dt ( ζ 2 + f h 2 ) = -( ζ 2 + f h 2 ). (βC -β * V ) h 2 + ẑ h 2 • ∇ × v 1 -v 2 2 . βC -β * V h 2 , (25b)
where

ζ i = ẑ•(∇ × v i ) = ∂xv i -∂yu i (i = 1, 2
) is relative vorticity, and q i = (ζ i + f )/h i is potential vorticity in each layer. One can construct a moist counterpart of potential vorticity in the lower layer with the help of the moist enthalpy (21), cf. Lambaerts et al.

(2011):

q 1m = ζ 1 + f m 1 . ( 26
)
The moist PV is conserved in the lower layer, modulo the Rayleigh drag effects:

d 1 dt ( ζ 1 + f m 1 ) = + ẑ • ∇ × v 1 -v 2 2 . βC -β * V m 2 1 . (27)
Note that the "asymmetric" assignment of the value of the step-function, which was discussed above, renders the moist PV in the lower layer conserved.

3. Illustration: application of improved mcRSW model to moist instabilities of hurricane-like vortices 3.1. Motivations

We will illustrate the capabilities of the improved mcRSW, the imcRSW, on the example of moist instabilities of hurricane-like vortices.

The mcRSW model, in its simplest one-layer version, captures well the salient properties of moist instabilities of such vortices, and clearly displays an important role of moisture in their development [START_REF] Lahaye | Understanding instabilities of tropical cyclones and their evolution with a moist-convective rotating shallow-water model[END_REF]). Below we extend the analysis of [START_REF] Lahaye | Understanding instabilities of tropical cyclones and their evolution with a moist-convective rotating shallow-water model[END_REF] to baroclinic tropical cyclones (TC), and use the imcRSW to check the role of new phenomena included in the model. Some questions which remained unanswered will be addressed, as well as new ones, possible to answer with the improved version of the model. In particular, we will investigate the influence of the size of TC (the radius of maximum wind) upon the structure of the most unstable mode, the role of vertical shear, and the evolution of inner and outer cloud bands at nonlinear stage of instability.

Fitting velocity and vorticity distribution of the hurricanes

We begin with building velocity and vorticity profiles of a typical TC within the two-layer model. An analytic form of the velocity profile is convenient both for the linear stability analysis, and for initialisations of the numerical simulations, so we construct a simple analytic fit with a minimal number of parameters.

V i (r) =              i (r -r 0 ) αi e -mi(r-r0) β i max (r -r 0 ) αi e -mi(r-r0) β i
, r ≥ r 0 , m 0 r, r r 0 .

(28)

Here i = 1, 2 indicate the lower and upper layer, respectively, r is the non-dimensional distance from the center, i measures the intensity of the velocity field, r 0 sets the non-dimensional distance of maximum wind from the centre, and other parameters allow to fit the shape of the distribution. A cubic Hermite interpolation across r = r 0 is made to prevent discontinuity in vorticity. Here and below we use a simple scaling where the distances are measured in units of barotropic deformation radius R d = √ gH/f , and velocities are measured in units of √ gH, where H is the total thickness of the atmospheric column at rest. (Hence, the parameter acquires a meaning of Froude number). Under this scaling the Rossby number of the vortex is proportional to the inverse of the non-dimensional radius of maximum wind (RM W ). A useful property of this parametrisation is a possibility to tune the ascending or descending trends of the wind near and far from the velocity peak. Velocity is normalised in a way that the maximum velocity is equal to . We suppose that velocity profile (28) corresponds to a stationary solution of "dry" equations of the model. Such solutions obey the cyclo-geostrophic balance in each layer:

V 1 r + f V 1 = g ∂ ∂r (H 1 + H 2 ) , (29a)
V 2 r + f V 2 = g ∂ ∂r (H 1 + αH 2 ) , (29b)
so the related H i (r) are obtained by integrating these equations using (28). The radial distribution of the relative vorticity in the vortex is given by (1/r)d [rV (r)] /dr. It should be emphasised that the radial gradient of the PV corresponding to the profile (28) has sign reversal, and hence the instability of the vortex is expected. Typical velocity and vorticity fields of an intense (category 3) vortex are presented in Figure 2. 1.

In what follows, we will be studying instabilities of thus constructed vortices, and their nonlinear saturation. The strategy will be the same as in [START_REF] Lahaye | Understanding instabilities of tropical cyclones and their evolution with a moist-convective rotating shallow-water model[END_REF]: namely, we identify the unstable modes of the vortex by performing detailed linear stability analysis of the "dry" adiabatic system, with switched-off condensation and vaporisation, and then use the unstable modes to initialise numerical simulations of nonlinear saturation of the instability, by superimposing them, with small amplitude, onto the background vortex. We will give below the results of numerical simulations of developing instabilities for three typical configurations which are presented in Table 1: weak barotropic (BTW) and baroclinic (BCW), and strong baroclinic (BCS) vortices.

Results of the linear stability analysis: the most unstable mode and its dependence on the radius of maximal wind

By applying the standard linearisation procedure and considering small perturbations to the axisymmetric background flow, we determine eigenfrequencies and eigenmodes of the "dry" system linearised about the background vortex. The technicalities of such c 2017 Royal Meteorological Society Prepared using qjrms4.cls analysis are the same as in [START_REF] Lahaye | Centrifugal, barotropic and baroclinic instabilities of isolated ageostrophic anticyclones in the two-layer rotating shallow water model and their nonlinear saturation[END_REF] extended to the two-layer configuration as in Rostami and Zeitlin (2017), and we pass directly to the results, which we will not give in detail either, limiting ourselves by what is necessary for numerical simulations in the next section.

The most unstable mode with azimuthal wavenumber l = 3 of the BCS vortex is presented in Fig. 3. The unstable mode of Figure 3 Pressure is clearly of mixed Rossby -inertia gravity wave type. Such unstable modes of hurricane-like vortices are well documented in literature, both in shallow-water models [START_REF] Zhong | An eigenfrequency analysis of mixed rossby-gravity waves on barotropic vortices[END_REF] and full 3D models [START_REF] Menelau | On the relative contribution of inertia-gravity wave radiation to asymmetric intsabilities in tropical cyclone-like vortices[END_REF]). It should be stressed that at Rossby numbers which are about 40 and small Burger number, as can be inferred from the right panel of Fig. 2, the vortex Rossby wave part of the unstable mode, which is known since the work [START_REF] Montgomery | A theory for vortex Rossby-waves and its application to spiral bands and intensity changes in hurricanes[END_REF] and is clearly seen in the upper panels of Fig. 3, is inevitably coupled to inertia-gravity wave field through Lighthill radiation mechanism, cf. [START_REF] Zeitlin | Decoupling of balanced and unbalanced motions and inertia -gravity wave emission: Small versus large rossby numbers[END_REF].

The most unstable modes of BCW and BTW vortices have wavenumber l = 4. With our scaling the strength of the vortex is inversely proportional to its non-dimensional RM W , and thus the structure of the most unstable mode depends on RM W . Yet, as follows from

Figure 4, the mode l = 4 is dominant through the wide range of RM W . In general, higher values of RM W correspond to higher azimuthal wavenumbers and lower growth rates.

Nonlinear evolution of the instability

We now use the unstable modes identified by the linear stability analysis to initialise numerical simulation of nonlinear evolution of the instability. We superimpose the unstable modes with weak amplitude (several per cent with respect to the background values) onto the vortex and trace the evolution of the system, as follows from numerical simulations with finite-volume well-balanced scheme developed for moist-convective RSW model [START_REF] Bouchut | Fronts and nonlinear waves in a simplified shallow-water model of the atmosphere with moisture and convection[END_REF]). Numerical simulations with each of the vortex configurations of Table 1 were performed both in "dry" (M), with diabatic effects switched off, and moist-convective (MCEV) environments.

The values of parameters controlling condensation, evaporation, vaporisation, and precipitation in MCEV environment are given in Table 2. The values of parameters controlling condensation: τc, Q s , stratification: s = θ 2 /θ 1 , evaporation: τ E , vaporisation: τv, and precipitation: Wcr, τp. ∆t is the time-step of the code.

τc Q s Q s τ E τv Wcr τp γ 5∆t 0.75 ≈ Q s 1.5 200∆t ≈1 day 10τc 0.01 3∆t 0.3 content.
We present below some outputs of the simulations, illustrating different aspects of moist vs "dry" evolution, and difference in behaviour of baroclinic and barotropic vortices.

Evolution of potential vorticity

We start by evolution of the PV field of the weak cyclone, as it is understandably slower than the evolution of the strong one, and different stages can be clearly distinguished. The evolution of potential vorticity in both layers during nonlinear saturation of the instability of the BCW vortex in MCEV environment is presented in Fig. 5,6. The simulations show formation of a transient polygonal pattern inside the RM W at initial stages, with a symmetry of the most unstable linear mode. The patterns of this kind are frequently observed [START_REF] Muramatsu | The structure of polygonal eye of a typhoon[END_REF][START_REF] Lewis | Polygonal eye walls and rainbands in hurricanes[END_REF][START_REF] Muramatsu | The structure of polygonal eye of a typhoon[END_REF][START_REF] Kuo | A possible mechanism for the eye rotation of typhoon herb[END_REF]). The polygon is further transformed into an annulus of high PV. Such annuli of elevated vorticity (the so-called hollow PV towers (Hendricks and Schubert 2010) ) are found in both moist-convective and dry cases. It is worth mentioning that the growth of the primary unstable mode is accompanied by enhancement of outer gravity-wave field, as follows from the divergence field presented in Fig. 7.

As follows from Figs. 5, 6 the polygon loses its shape at t ≈ 17. At this time the modes with azimuthal wavenumbers l = 1, 2 are being produced by nonlinear interactions, and start to grow and interact with the polygonal eye-wall, which leads to symmetry loss by the core. A secondary, dipolar instability of the core thus develops, and gives rise to formation of an elliptical central vortex, corresponding to azimuthal mode l = 1, and of a pair of satellite vortices indicating the presence of l = 2 mode. The interaction of initial l = 4 mode with emerging l = 1 and l = 2 modes is accompanied by inertia-gravity wave (IGW) emission, and enhancement of water vapour condensation that will be discussed below. It should be emphasised that interaction between l = 2 mode and elliptical eye, of the kind we observe in simulations, was described in TC literature, e.g. [START_REF] Kuo | A possible mechanism for the eye rotation of typhoon herb[END_REF], where reflectivity data from a Doppler radar were used to hypothesise that it was due to azimuthal propagation of l = 2 vortex Rossby waves around the eye-wall.

Further nonlinear evolution consists in breakdown of the central ellipse with subsequent axisymmetrisation of the PV field, and its intensification at the center. This process characterises the evolution of both BTW (not shown) and BCW vortices, but is more efficient in the baroclinic case, as follows from Fig. 8. As seen in the Figure, the azimuthal velocity in the core region with r < 0.5RM W is subject to strong intensification. The exchanges of PV between the eye-wall and the eye, and intensification are known from the barotropic simulations [START_REF] Montgomery | A theory for vortex Rossby-waves and its application to spiral bands and intensity changes in hurricanes[END_REF][START_REF] Schubert | Polygonal eyewalls, asymmetric eye contraction, and potential vorticity mixing in hurricanes[END_REF][START_REF] Lahaye | Understanding instabilities of tropical cyclones and their evolution with a moist-convective rotating shallow-water model[END_REF]). As we see in Fig. 8, the intensification is enhanced by baroclinicity of the background vortex. This is confirmed by Fig. 9, and by Fig. 10, which illustrate the enhancement effect of both moist convection and baroclinicity upon palinstrophy, which is defined as

P(t) = 1 2 ∇ζ.∇ζdxdy, (30)
in each layer, and which diagnoses the overall intensity of vorticity gradients.

It is worth emphasising that, because of higher vorticity and smaller RMW, the axisymmetric steady state is achieved in the lower layer more rapidly than in the upper one in the case of baroclinic vortices.

In the case of intense vortex, nonlinear evolution of the instability follows similar scenario, but is considerably accelerated, as follows from Fig. 11.

Spiral cloud bands

Tropical cyclones exhibit specific cloud patterns. The new version of the model gives a possibility to follow clouds and precipitation, and it is interesting to test its capability to produce realistic cloud patterns. There are two types of cloud and rain bands associated with tropical cyclones, as reported in literature: the inner bands, that are situated close to the vortex core, within ≈ 2 RM W , and the outer spiral bands located farther from the centre and having larger horizontal scales [START_REF] Guinn | Hurricane spiral bands[END_REF][START_REF] Wang | How do outer spiral rainbands affect tropical cyclone structure and intensity[END_REF]. Fig. 12 shows formation of inner and outer cloud bands, the latter having characteristic spiral form, during nonlinear evolution of the instability. Spiral cloud bands are related to inertia-gravity "tail" of the developing unstable mode. The link of spiral bands to inertia-gravity waves in "dry" RSW model of hurricanes was discussed in literature [START_REF] Zhong | A theory for mixed rossby-gravity waves in tropical cyclones[END_REF]). Here we see it in "cloud-resolving" imcRSW. It is to be stressed that amount of clouds strongly depends on the initial water vapour content. If it is closer to the saturation value, the amount of clouds and precipitation obviously increases and eventually covers the whole vortex.

Conclusions and discussion

Thus, we have shown that the moist-convective rotating shallow water model "augmented" by adding precipitable water, and relaxational parametrisations of related processes of vaporisation, precipitation, together with entrainment, is capable to capture some salient features of the evolution of instabilities of hurricane-like vortices in moist-convective environment, and allows to analyse the importance of moist processes on the life-cycle of these instabilities. There exist extended literature on the dynamics of the hurricanes eyewall, with tentative explanations in terms of transient internal gravity waves, which form spiral bands, cf. [START_REF] Lewis | Polygonal eye walls and rainbands in hurricanes[END_REF]Hawkins (1982) Willoughby (1978), [START_REF] Kurihara | On the development of spiral bands in a tropical cyclone[END_REF], or alternative explanations [START_REF] Guinn | Hurricane spiral bands[END_REF] in terms of PV dynamics and vortex Rossby waves. Thus [START_REF] Schubert | Polygonal eyewalls, asymmetric eye contraction, and potential vorticity mixing in hurricanes[END_REF] obtained formation of polygonal eyewalls as a result of barotropic instability near the radius of maximum wind in a purely barotropic model, without gravity waves. A detailed analysis of instabilities of tropical cyclones was undertaken with a cloud-resolving model in [START_REF] Naylor | Evaluation of the impact of moist convectionon the developmentof asymmetric inner core instabilities in simulated tropical cyclones[END_REF], and showed that the results of [START_REF] Schubert | Polygonal eyewalls, asymmetric eye contraction, and potential vorticity mixing in hurricanes[END_REF] gave a useful first approximation for the eyewall instabilities. As was already mentioned in section 3.3, at high Rossby numbers the vortex Rossby wave motions are inevitably coupled to inertia-gravity waves, and our linear stability analysis confirms this fact. The mixed character of the wave perturbations of axisymmetric hurricane-like vortices was abundantly discussed in literature, e.g. [START_REF] Zhong | A theory for mixed rossby-gravity waves in tropical cyclones[END_REF]. A detailed analysis of instabilities of hurricane-like vortices in continuously stratified fluid was given recently by [START_REF] Menelau | On the relative contribution of inertia-gravity wave radiation to asymmetric intsabilities in tropical cyclone-like vortices[END_REF], where it was shown that the inertia-gravity part of the unstable modes intensifies with increasing Froude number. The vortex profiles used above in section 3.2 have moderate Froude numbers, and the corresponding unstable modes have weak inertia-gravity tails. They are, however, sufficient to generate spiral cloud patterns, as we showed. The development of the instability of the eyewall proper at early stages (up to ≈ 40f -1 ) is only weakly influenced by moist convection, in accordance with findings of [START_REF] Naylor | Evaluation of the impact of moist convectionon the developmentof asymmetric inner core instabilities in simulated tropical cyclones[END_REF]. This can be seen from comparison of the right and left panels of Fig. 9 and from Fig. 10.

An advantage of our model, as compared to simple barotropic models, is its ability to capture both vorticity and inertia gravity waves Although we limited ourselves above by an application to tropical cyclones, the model can be used for analysis of various phenomena in mid-latitudes and tropics. The passage to the equatorial beta-plane is straightforward in the model, and it can be easily extended to the whole sphere. An important advantage of the model is that it allows for self-consistent inclusion of topography in the numerical scheme, giving a possibility to study a combination of moist and orographic effects. As was already mentioned, more realistic parametrisations of the boundary layer are available, and generalisations to three-layer versions are straightforward.

c 2017 Royal Meteorological Society Prepared using qjrms4.cls

Figure 1 .

 1 Figure1. Notations for the simplified two-layer scheme with mass flux across material surfaces. From[START_REF] Lambaerts | Simplified two-layer models of precipitating atmosphere and their properties[END_REF], with permission of AIP.

Figure 2 .

 2 Figure2. Normalised radial structure of azimuthal tangential wind with a fixed slope close to the centre (left panel) and relative vorticity (right panel) in both layers corresponding to the BCS vortex in Table1.

Figure 3 .

 3 Figure 3. Upper row: Pressure and velocity fields in the x -y plane in the lower (left panel) and upper (right panel) layers corresponding to the most unstable mode with azimuthal wavenumber l = 3 of the BCS vortex of Table1. Lower row: left panel-corresponding divergence field of the most unstable mode, right panel -radial structure of three components of the most unstable mode: pressure anomaly η, and radial (v) and tangential (u) components of velocity ; dashed (solid) lines: imaginary (real) part, thick (thin) lines correspond to upper (lower) layer. Note that the domain in the lower left panel is ≈ ten times larger than that of the upper panels.

Figure 4 .

 4 Figure4. Dependence of the linear growth rates (in units of f -1 ) of the unstable modes with different azimuthal wavenumbers on the radius of maximum wind (RMW ).

Figure 5 .

 5 Figure5. Nonlinear evolution of the most unstable l = 4 mode superimposed onto the background BCW vortex in MCEV environment, as seen in potential vorticity field in the lower layer. Formation of meso-vortices (zones of enhanced PV in the vorticity ring) is clearly seen, giving way to axisymmetrisation and monotonisation of the PV profile. Time is measured in units of f -1 .

Figure 6 .Figure 7 .Figure 8 .

 678 Figure6. Same as in Fig.5, but for the upper layer.

Figure 9 .Figure 10 .Figure 11 .Figure 12 .

 9101112 Figure9. Effect of vertical shear on intensification of vorticity in the vortex core in environments without (M ) and with (MCEV ) moist convection and surface evaporation. The vorticity is normalised by its initial value.

Table 1 .

 1 Parameters of the background vortices. BCW(S): weak(strong) baroclinic, BTW: weak "barotropic", without vertical shear, l: the most unstable azimuthal mode

	conf ig. l	1	2	α 1	α 2	β 1	β 2	m 1	m 2	r 01	r 02
	BCS	3 0.41 0.49	4.5	4.5	0.180 0.178 48 47.5 0.01 0.0101
	BTW	4	0.4	0.40 2.25 2.25	0.25	0.25	14	14	0.1	0.1
	BCW	4	0.4	0.36 2.25 2.25	0.25	0.237 14 12.6	0.1	0.115
	observations were collected from Atlantic and eastern Pacific storms during 1977 -2001:

It has the form which is consistent with Mallen et al. (2005) where flight-level c 2017 Royal Meteorological Society Prepared using qjrms4.cls

Table 2

 2

	Prepared using qjrms4.cls

. It must be stressed that amount of precipitable water in each layer is highly sensitive to the values of parameters, especially to the intensity of surface evaporation. Condensation and precipitation time scales are chosen to be short, just few time steps ∆t, while vaporisation and surface evaporation time-scales are much larger which is consistent with physical reality. Changing these parameters within the same order of magnitude does not lead to qualitative changes in the results. Wcr is an adjustable parameter that controls precipitation and γ controls entrainment of condensed water. The MCEV simulations were initialised with spatially uniform moisture c 2017 Royal Meteorological Society