MS Audrey Bouchet

MS †, † †, §,• Benjamin Lemasson

PhD Géraldine Le Duc

MS †, † † Cécile Maisin

MS Elke Bräuer-Krisch

PhD || Erik Albert Siegbahn

Luc Renaud

PhD Enam Khalil

PhD †, † † Chantal Rémy

MS †, † † Cathy Poillot

Renaud ¶ Luc

PhD Alberto Bravin

MD $ Jean A Laissue

PhD Emmanuel L Barbier

PhD Raphaël Serduc
email: raph.serduc@gmail.com

Preferential effect of synchrotron microbeam radiation therapy on intracerebral 9L gliosarcoma vascular networks

Keywords: Synchrotron Microbeam radiation therapy, intracerebral 9L gliosarcoma, tumor vasculature

Purpose

Tumor vasculature is deemed to be a promising target for the treatment of solid cancerous lesions. Different therapeutic approaches, alone or combined, such as antiangiogenic therapy, chemo-, radio-and immunotherapy are currently refined in order to specifically interfere with blood vessel growth in the tumor [START_REF] Jain | Angiogenesis in brain tumours[END_REF]. This objective might be also achieved by a novel form of radiosurgery, known as Microbeam Radiation Therapy (MRT) [START_REF] Slatkin | Subacute neuropathological effects of microplanar beams of x-rays from a synchrotron wiggler[END_REF], which uses a particular irradiation geometry and X-rays generated in a synchrotron light source; the beam is spatially fractionated into an array of quasi-parallel planar microbeams which are a few tens microns wide and separated by a few hundred microns. This unique irradiation geometry allows very high dose (hectograys) deposition efficient for the treatment of brain tumors in small animals while ensuring an unusually high normal tissue tolerance, even in developing larger animals [START_REF] Brauer-Krisch | Effects of pulsed, spatially fractionated, microscopic synchrotron X-ray beams on normal and tumoral brain tissue[END_REF][START_REF] Laissue | Prospects for microbeam radiation therapy of brain tumours in children to reduce neurological sequelae[END_REF].

The preferential effect on malignant tissues has been mainly attributed to selective effects of microbeams on immature tumor vessels versus lack of microbeam effects on the differentiated normal vasculature [START_REF] Dilmanian | Response of rat intracranial 9L gliosarcoma to microbeam radiation therapy[END_REF], although other factors might be operative. Lesions induced by the microbeams in the normal capillary network can be rapidly repaired [START_REF] Serduc | Characterization and quantification of cerebral edema induced by synchrotron x-ray microbeam radiation therapy[END_REF][START_REF] Serduc | In vivo two-photon microscopy study of short-term effects of microbeam irradiation on normal mouse brain microvasculature[END_REF] and normal arteries are microbeam-resistant too [START_REF] Van Der Sanden | Tolerance of arteries to microplanar x-ray beams[END_REF]. Thus, there is no decrease in vessel density, blood volume, no hemorrhage and no ischemic necrosis after MRT in normal mammalian tissues. Conversely, no clear experimental proof of the effect of MRT on tumor vessels has been given in the literature. Our recent magnetic resonance imaging study failed to reveal important damage to tumor vessels after MRT [START_REF] Serduc | Brain tumor vessel response to synchrotron microbeam radiation therapy: a short-term in vivo study[END_REF]. Crossfired MRT slowed the growth of the heterologously implanted 9L gliosarcoma (9LGS) and the development of its vasculature. We now characterize the biological effects of 50 µm-wide microbeams, efficient in terms of tumor control and normal tissue sparing [START_REF] Serduc | Synchrotron microbeam radiation therapy for rat brain tumor palliation-influence of the microbeam width at constant valley dose[END_REF], on normal rat brain and intracerebral 9LGS, particularly on growth and morphometric parameters of blood vessels from day 2 to day after irradiation, and on the expression of the Vascular Endothelial Growth Factor.

Methods and Materials

All procedures related to animal care conformed to the Guidelines of the French Government with licenses 380325 and 380321. Rats were anesthetized with an intraperitoneal injection of xylazine/ketamine (64.5/5.4mg.kg -1 for irradiations) and 2.5% isoflurane (MRI follow-up) for maintenance.

Tumor inoculation

Transplantable 9LGS cells were implanted in 59, 10-week-old-Fisher rats weighing between 180 and 220g. Implantation procedure followed the protocol described in [START_REF] Serduc | Synchrotron microbeam radiation therapy for rat brain tumor palliation-influence of the microbeam width at constant valley dose[END_REF]. Briefly, 10 4 9LGS cells suspended in a 1μl DMEM medium were injected into the right caudate nucleus (3.5mm from the bregma) at a depth of 5.5mm from the dura.

Magnetic resonance imaging and group sorting

We used a 4.7T Bruker Avance console and volume/surface cross coil configuration.

Seven days after 9LGS implantation, rats underwent anatomical MRI (T 2 -weighted image) in order to sort them two groups with similar mean tumor size. Those 59 rats were divided into a control group (n=25) and a group of 34 rats and irradiated by MRT 10 days after 9LGS implantation. Rats (n=4 to 8) underwent MRI and 4 of them were culled at each delay for biological studies. The survival curves were established using 9 rats of each group. The number of days elapsed after implantation (n) and treatment (m) is designated as D n T m .

Radiation source and irradiation geometry

Irradiations were performed on the ID17 biomedical beamline at the European Synchrotron Radiation Facility using X-rays emitted tangentially from relativistic electron bunches circulating in a storage ring. The wiggler produces a wide spectrum of photons which extends from 50 over 350keV (median energy: 90keV). The quasi-laminar beam was shaped into an array (8mm wide) of rectangular, microscopically thin and quasi-parallel microbeams using a multislit collimator [START_REF] Bräuer-Krisch | New technology enables precision multi slits collimators for MRT (Microbeam Radiation Therapy)[END_REF]. Rats in the treated group were irradiated (D 10 T 0 ) after tumor inoculation using two, 10mm high, arrays of 40 vertical, planar microbeams (width 50μm, 200μm on-center spacing), one anteroposterior, the other lateral, both focused for orthogonal intersection (crossfiring) in the anterior part of the right hemisphere, where the 9LGS was located. Before reaching and after leaving the volume of intersection enclosing the tumor, each array alone covered a volume of the rat's head in which it was unidirectional. Animals were aligned into the beam using 3 live cameras and the coordinates of the center of the microbeam arrays were 3.5mm left from the medio-lateral line, -6mm in the dorso-ventral direction from the skin surface and at the bregma point in the anterioposterior direction. The in-microbeam entrance dose was 400Gy per array, the mean dose rate ~16000Gy.s -1 .

Monte Carlo simulation and doses calculations

The doses were calculated by means of the Monte Carlo method normally used for MRT dose planning and quality assurance at the ESRF [START_REF] Siegbahn | Determination of dosimetrical quantities used in microbeam radiation therapy (MRT) with Monte Carlo simulations[END_REF]. Further details of our MRT dosimetry protocols are provided as Supplementary Material.

Rat follow-up by MRI

Eight rats were imaged one day before (D 9 T -1 ) and animals from all groups at 2, 5, 8 15 and 45 days after MRT treatment. For each group, 4 rats were followed up and imaged at every time point. Four rats per groups were also imaged at every delay after irradiation and culled for immunohistological studies and western blot analysis.

Each MRI session (duration 1h15), detailed in the Supplementary material, was organized as follow:

-Anatomical T 2 -weighted images -BVf (blood volume fraction) and VSI (vessel size index) measurements -Vascular permeability of P846 (Gd-based contrast agent) measurements

MRI data analysis

BVf, VSI and vessel permeability maps were computed pixel by pixel using an inhouse program developed under Matlab (MathWorks, NA). BVf and VSI were calculated as described previously [START_REF] Tropres | Vessel size imaging[END_REF][START_REF] Tropres | In vivo assessment of tumoral angiogenesis[END_REF]. Null pixels were excluded in the mean calculation in order to give information only on perfused brain regions. Vessel permeability was estimated by subtracting the peak MR signal obtained after the injection of Gd-DOTA from the peak MR signal measured before injection and expressed as a percentage of contrast enhancement.

Three regions of interest (ROIs) were manually delineated on the T 2 -weighted image on which the tumor appeared the largest (see fig. 1) i) the whole tumor, ii) the ipsilateral (right) hemisphere without tumor and iii) the contralateral hemisphere, and reported on the other maps. Tumor and surrounding parts of the normal brain in the right hemisphere were exposed to the two intersecting arrays, normal brain parts of the left hemisphere unidirectionally to the lateral array alone. Tumor volumes were determined by multiplying the pixel size by the number of "tumoral pixels" on all MR slices on which the tumor was visible. Two-way ANOVA followed by Bonferroni post-hoc tests were used for statistics. A Mann-Whitney test was used to compare 9LGS volumes. The statistical analyses were performed with the Prism program (GraphPad Software, USA). Values were considered significantly different when p<0.05.

Immuno / histological analyses (IHC)

Type-IV collagen (C IV ), a protein constituent of the vessel basal lamina, and RECA-1 (Rat Endothelial Cell Antigen), an antigen of endothelial cells, were labeled as previously described [START_REF] Serduc | In vivo two-photon microscopy study of short-term effects of microbeam irradiation on normal mouse brain microvasculature[END_REF][START_REF] Serduc | First trial of spatial and temporal fractionations of the delivered dose using synchrotron microbeam radiation therapy[END_REF]. Cell proliferation activity was assessed by means of Ki67 labeling. Further details of our IHC protocols and quantification are provided as Supplementary Material.

Western blot analysis of VEGF

The western blot protocol of VEGF is described in the Supplementary Material section.

Results

The main results of this study are summarized in table 1. Radiation doses were calculated in a rat head phantom using Monte Carlo simulations. The peak dose deposited at 1cm depth in rat brain was 350/700Gy, generating a valley dose of 12.5/25Gy in the unidirectionally/bidirectionally irradiated brain regions respectively.

Survival curves

Figure 1C shows the Kaplan-Meier representation of the survival of untreated and MRT-treated 9LGS bearing rats. MRT significantly increased the MST (mean survival time) of treated animals compared with the control group of non-irradiated rats (65 days versus 20, respectively, log rank test: p=0.003). Thereafter, the evaluation of tumor volume using on T 2 -weighted MR images became more problematic since 2 rats did not exhibit tumor anymore (the hyposignal as shown in figure 1Ah was taken into account for lesion size determination); 2 rats showed small tumors and 1 rat a recurrent tumor comparable to those measured in control rats at D 18 T 8 .

MRT effects on 9LGS growth

MRT effects on 9LGS proliferative activity

The analyses Ki67 immunolabeling and HE stained sections revealed that MRTinduced damages were detectable from D 15 T 5 (fig. 2): cell density decreased on H&E stained brain sections in irradiated tumors; the proliferative activity (proportion of Ki67 positive cells) of the tumor was distinctly reduced. Giant tumoral cells were observed. At D 25 T 15 , the 9LGS tissue was strongly disorganized, exhibiting cellular polymorphism and large lacunae and/or necrotic regions in the core. Only few dividing cells were found at the tumor periphery.

For two thirds of the animals observed at D 55 T 45 by immunohistology, it was not possible to delimit a tumor. Conversely, one third of the animals presented a tumor with numerous Ki67 positive cells which had invaded the whole ipsilateral hemisphere (fig. 2). Normal tissue exhibited astrocyte and neuron loss in the microbeam paths starting at day 5 after irradiation.

MRT effects on the evolution of normal and tumoral blood vessels of the brain

Blood volume fraction

Cerebral BVf (fig. 3) were significantly higher in the 9LGS than in the ipsi-and contralateral normal brain for all observation times. A decrease in tumoral blood volume fraction (TBVf) was detected between D 9 T -1 and D 18 T 8 in irradiated or non-irradiated tumors (p<0.001), whereas the CBVfs of the ipsi-and contralateral normal hemisphere did not change. Thereafter, in irradiated animals, TBVf decreased between D 18 T 8 and D 25 T 15 (<0.001).

Vessel size index

VSI values increased similarly in both irradiated and non irradiated tumors until D18, where they were significantly higher in tumors than in normal tissues (p<0.001). Thereafter, VSI values decreased in irradiated tumors between D 18 T 8 and D 25 T 15 (-23.8%, p<0.001). In normal tissues surrounding the tumor within the volume targeted by the intersecting, crossfired arrays (ipsilateral+MRT), VSI increased between D 18 T 8 and D 55 T 45 (8.5 to 13.7µm, p<0.001).

Contrast enhancement

Vessel permeability (indicated as contrast MR enhancement, figure 3) was significantly higher in 9LGS than in the normal ipsi-or contralateral hemisphere, and it was not modified in normal brain tissue after unidirectional or crossfired irradiations. Conversely, MRT led to an increase in tumor vessel permeability during the 5 days following the irradiation: The MR signal enhancement reached 26.6% at D15 while it was only 20.3% in nonirradiated tumors (p<0.001). Then, the blood vessel permeability of irradiated tumors decreased by 74% between D 15 T 5 and D 25 T 15 (p<0.001).

Blood vessels; immunohistological observations

Vascular networks were observed on brain sections labeled by C IV (basal membrane) and RECA-1 (endothelial cells, fig. 45). The co-expression ratio of these proteins in vessels is reported in on each image on Figure 4. During the entire experiment, no damage was detected in the unidirectionally irradiated volume of the contralateral hemisphere (fig. 4). All brain capillaries were double stained by C IV / RECA-1 antisera. Conversely, the morphology of blood vessels located in the normal brain tissue surrounding the tumor within the volume of intersecting arrays changed around D 55 T 45 : the normal vascular network showed an important disorganization and increased tortuosity. C IV immunolabeling revealed the presence of narrow (diameter 1-2 µm) entities without endothelial cells, as shown by a quasi-complete lack of reactivity for RECA-1 labeling. Few capillaries of normal size lacked homogenous RECA-1 labeling.

Tumor vessel density was lower than that observed in normal brain tissue. Spatial vessel distribution was irregular and vessel diameters were distinctly larger compared with normal brain blood vessels. The number of blood vessels per surface area in non-irradiated as well as in irradiated tumors decreased with time after implantation. However, beginning at D 12 T 2 , RECA-1 labeling intensity gradually faded in 9LGS vessels in the beam-crossing volume. In tumors observed from D 18 T 8 to D 55 T 45 , most of 9LGS vessels were not immunoreactive for RECA-1, as shown in figure 5 for a rat culled at D 18 T 8 .

MRT effects on VEGF expression

No VEGF expression was detected in contralateral (left) hemispheres of control rats during the whole experiment. Conversely, MRT induced an expression of VEGF increasing with time in normal brain tissues, i.e., from D 12 T 2 or D 15 T 5, in both hemispheres, whether irradiated by one or by both intersecting arrays (Fig. 6). At D 55 T 45 , the amount of VEGF detected in normal brain parts of the ipsilateral hemisphere situated in the beam-crossing region was >2.5-times higher (p<0.001) than the one measured in the left hemisphere irradiated unidirectionally by the lateral array (Contra+MRT).

In 9LGS, the relatively low expression of VEGF did not markedly change during their growth (fig. 6). Conversely, the irradiation led to a progressive and significant increase in VEGF expression in tumors. Indeed, the ratio VEGF/GADPH changed from 0.4 to 1.1 between D 12 T 2 and D 18 T 8 while it fluctuated from 0.5 to 0.8 in untreated tumors during the same period (NS for untreated tumors, p<0.05 in irradiated tumors). At D 25 T 15 and D 55 T 45, the relative expression of VEGF increased significantly in irradiated tumors (p<0.001), to values 1.35 and 1.8 times higher than the one measured at D 18 T 8 .

Discussion

The crossfired irradiation increased significantly the MST of rats bearing intracerebral 9LGS compared with unirradiated controls (p=0.0003). MRT inhibited tumor growth during few weeks after treatment as shown by MRI. Histopathologically, MRT induced a decrease in tumor cell density and reduced the proliferative activity of the lesion. The salient findings in this study were 1) the preferential loss of endothelial cells in tumor vessels exposed to crossfired arrays, from day 5 after treatment and 2) the significant decrease in tumoral blood volume fraction during the second week after irradiation; conversely, the microvasculature and the perfusion in the normal brain volumes irradiated unidirectionally in the contralateral hemisphere remained undamaged. Normal blood vessels in peritumoral normal brain tissue were damaged by both intersecting arrays.

The hypothesis of a preferential effect of MRT on tumoral blood vessels rather than on normal vasculature is well established [START_REF] Dilmanian | Response of rat intracranial 9L gliosarcoma to microbeam radiation therapy[END_REF][START_REF] Laissue | Neuropathology of ablation of rat gliosarcomas and contiguous brain tissues using a microplanar beam of synchrotronwiggler-generated X rays[END_REF]. However, Serduc et al. described a limited influence of the vascular component during 9LGS palliation in mice by MRT [START_REF] Serduc | Brain tumor vessel response to synchrotron microbeam radiation therapy: a short-term in vivo study[END_REF]. A significant increase in tumor vessel permeability on D 12 T 2 and D 15 T 5 (+44% and +23%) preceded tumor necrosis in the present experiment. No change in blood-brain barrier permeability was detected in ipsilateral (exposed to both arrays) and contralateral (exposed to one array only) normal brain tissues. This transient tumor-blood barrier might be used to deliver intravenously injected drugs to the tumor core. Adjuvant chemotherapy would be justified by the fact that in previous studies MRT did not sterilize intracerebral 9LGS, a highly resistant gliosarcoma model, even when adequate irradiation configurations and valley doses were used (45Gy (15), 36Gy [START_REF] Serduc | Microbeam radiation therapy parameters optimization for rat brain tumors palliation. Influence of the microbeam width at constant valley dose[END_REF]). In the present experiment, tumor growth was inhibited during several weeks after MRT in association with reduced proportions of proliferative Ki67-positive tumor cells, starting on D 15 T 5 . Two weeks after MRT, large necrotic regions were detected in the tumor core. The presence of proliferative cells at the tumor periphery may indicate tumor recurrence.

As several rats survived for a very long period after MRT for 9LGS in other experiments [START_REF] Dilmanian | Response of rat intracranial 9L gliosarcoma to microbeam radiation therapy[END_REF][START_REF] Laissue | Neuropathology of ablation of rat gliosarcomas and contiguous brain tissues using a microplanar beam of synchrotronwiggler-generated X rays[END_REF], despite the presence of ~2% surviving clonogenic 9LGS cells in the tumor [START_REF] Dilmanian | Response of rat intracranial 9L gliosarcoma to microbeam radiation therapy[END_REF], it was suggested that tumor ablation was due to the preferential effect of MRT on tumor blood vessels resulting in ischemic necrosis of the tumor. However, in the present study, morphometric parameters of the vascular network were not affected shortly after MRT, i.e. between D 12 T 2 and D 18 T 8 . Indeed, we observed no difference in BVf and VSI values between irradiated and non irradiated tumors. Conversely, during this period, RECA-1 labeled endothelia were not detected in all Type-IV collagen positive vascular structures of irradiated tumors, while their vascular counterparts in the time-matched non-irradiated tumors exhibited obvious double staining. Loss of RECA-1 immunoreactivity indicates a loss of endothelial cells. Alternatively, some vessels might have been occluded by thrombosis or following endothelial denudation. However, BVf and VSI values measured by MRI after Sinerem iv. injection did not indicate a decreased tumor perfusion after irradiation. Thus, the loss of endothelial cells might not have measurably impaired tumor blood supply after MRT. Loops of extracellular matrix, devoid of endothelia, seem able to temporarily irrigate tumors with blood [START_REF] Maniotis | Vascular channel formation by human melanoma cells in vivo and in vitro: Vasculogenic mimicry[END_REF]. Such C IV positive vessel-like structures lacking immunoreactive endotheliocytes, might thus permit a transitory perfusion of 9GLS; the parameters showing the evolving perfusion after irradiation in the same animal support this interpretation (Fig. 1A): after D 18 T 8, tumor size reduction until D 25 T 15 correlates with the significant decrease in TBVf and in VSI (respectively -24.6%, p<0.05 and -23.8%, p<0.01). Conversely, BVf and VSI values, and morphology of brain vessels in the normal brain parts exposed to crossfired arrays remained unchanged during the same period. These data persuasively support the existence of a preferential damaging effect of MRT on tumor vessels. During the second week after treatment, tumor necrosis can be associated with the delayed decrease in TBVf and VSI.

The use of 50 µm-wide microbeams significantly improved the MSTs for rats with 9LGS in our laboratory [START_REF] Serduc | First trial of spatial and temporal fractionations of the delivered dose using synchrotron microbeam radiation therapy[END_REF][START_REF] Serduc | Microbeam radiation therapy parameters optimization for rat brain tumors palliation. Influence of the microbeam width at constant valley dose[END_REF], notably in this study. This particularly radioresistant tumor requires very high radiation dose deposition for palliation or even ablation; a 40Gy radiosurgical dose controlled only 50% of 9LGS [START_REF] Kim | Fractionated radiosurgery for 9L gliosarcoma in the rat brain[END_REF]. Tumor responses to MRT on D 55 T 45 in our study were markedly heterogeneous Indeed, 2/5 rats did not exhibit tumors; 2 animals showed well-delimited small tumors and 1 rat had a large tumor invading the whole hemisphere. The cause for these different responses is not understood but these 5 rats displayed a common change in MRI, a cerebral edema in parts of the ipsilateral hemisphere exposed to both intersecting arrays. These results were supported by demonstration of infiltrates of monocytes/macrophages (Supplementary materials, Fig. e1) and of a strong disorganization and remodeling of the normal vascular network at this site. Normal vessels appeared tortuous and larger, an observation consistent with the significant increase in VSI at D 55 T 45 (+9.8%). Thin matricial entities, probably not perfused, were commonly observed.

RECA-1 labeling was inhomogeneous.

A gradual increase in VEGF expression was detected in tumoral and in normal tissues.

However, the kinetics and the amount of the protein produced seemed to depend on the radiation dose (exposure to uni-or bidirectional irradiation) and the status of the tissue (normal or tumoral).

VEGF, secreted e.g. in ischemia, stroke and after irradiation [START_REF] Kim | Upregulation of VEGF and FGF2 in normal rat brain after experimental intraoperative radiation therapy[END_REF][START_REF] Verhoeff | Tumour control by whole brain irradiation of anti-VEGF-treated mice bearing intracerebral glioma[END_REF], is commonly induced in presence of hypoxia; changes in tumoral and normal brain oxygen saturation after MRT should thus be evaluated. VEGF in the ipsilateral hemisphere was detectable from the second day after MRT and reached levels comparable to those measured in non-irradiated tumors. Thus, VEGF is implied in cerebral response in the early phase after MRT but its role is not defined. The overexpression of VEGF after MRT could have increased blood vessel permeability and led to severe brain edema.

Since the size of the tumors did not change significantly between D 25 T 15 and D 55 T 45 , we assume that MRT-induced brain edema in the crossfired regions was the main cause of animal deaths in our study. These qualitative findings are consistent observations ( 16) that major histopathologic brain and vascular damage was confined to MRT beam-crossing regions, whereas only minor damage was apparent in zones of the brain irradiated unidirectionally [START_REF] Laissue | Neuropathology of ablation of rat gliosarcomas and contiguous brain tissues using a microplanar beam of synchrotronwiggler-generated X rays[END_REF]. A comprehensive and detailed quantitative brain damage index for the present study was not established. However, the histopathologic findings were comparable, despite the different irradiation modalities.

Cross-fired MRT , particularly efficient for tumor control, may require corticosteroid injections [START_REF] Piette | Treating gliomas with glucocorticoids: from bedside to bench[END_REF] to control ensuing brain edema. Further, recovery from radiation injury may be influenced by growth factors involved in cell renewal and differentiation of neural tissue, e.g. by VEGF. Thus, a brief therapy with VEGF can prevent or delay radiation myelopathy after spinal cord injury in rats [START_REF] Andratschke | Potential role of growth factors in diminishing radiation therapy neural tissue injury[END_REF]. Experimentally, inhibition of VEGF may increase local tumor control after irradiation, possibly by preventing tumor angiogenesis [START_REF] Verhoeff | Tumour control by whole brain irradiation of anti-VEGF-treated mice bearing intracerebral glioma[END_REF].

The preferential effect of MRT on tumor vasculature may lead to further improvements.

Most importantly, the tumor targeting for MRT used preclinically must be improved and adapted to clinical standards. Irradiation fields must be limited to the tumor; image guided, non coplanar irradiations, which limit normal tissue damage, should be adopted [START_REF] Serduc | First trial of spatial and temporal fractionations of the delivered dose using synchrotron microbeam radiation therapy[END_REF][START_REF] Serduc | High-precision radiosurgical dose delivery by interlaced microbeam arrays of high-flux low-energy synchrotron X-rays[END_REF]. Extratumoral anterior parts of the right hemisphere were irradiated by the intersecting lateral and anteroposterior arrays (CROSS). The 9LGS tumor in the right anterior hemisphere was also cross-irradiated. The number of days (n) elapsed after implantation or irradiation is designated as DnTn, respectively. RECA = Rat Endothelial Cell Antigen.

Parameter

Figure

  Figure1Ashows a T 2 -weighted MRI follow-up of irradiated and non irradiated 9LGS which

Figure 1 :

 1 Figure 1: MRI follow-up, survival curves and lesion size evolution. T 2 -weighted MR

Figure 2 :

 2 Figure 2: Temporal immunohistological follow-up of untreated and MRT treated

Figure 3 :

 3 Figure 3: MRI vessel parameter measurements. BVf (a), VSI (b) and vessel permeability

Figure 4 :

 4 Figure 4: Temporal immunohistological follow-up of unirradiated and MRT-irradiated

Figure 5 :

 5 Figure 5: Confocal images of unirradiated and MRT-irradiated 9LGS blood vessels.

Figure 6 :

 6 Figure 6: VEGF expression evolution. Western blot for VEGF in the contralateral (Contra)

Figure 1 :

 1 Figure 1: MRI follow-up, survival curves and lesion size evolution. T 2 -weighted MR

Figure Click here to download Figure: Figure 1 Figure 2 :

 12 Figure Click here to download Figure: Figure 1 MRIFup Lesion size Surv_rev.doc

Figure Click here to download Figure: Figure 2 Figure 3 :

 23 Figure Click here to download Figure: Figure 2 HEKI.doc

Figure 4 :

 4 Figure Click here to download Figure: Figure 3 BVVSIVP.doc

Figure 5 :

 5 Figure Click here to download Figure: Figure 4 CollRECA.doc

Figure Click here to download Figure: Figure 5 Figure 6 :

 56 Figure Click here to download Figure: Figure 5 Confocal.doc

Table 1 :

 1 MRT effects on blood vessels located either in a 9LGS tumor or in normal brain tissue. Parts of left hemisphere were irradiated unidirectionally by the lateral array (UNI).

		Blood vessels in	Irradiation	Observations
	Blood vessel morphology	Normal Brain	UNIDIRECTIONAL (UNI) CROSSFIRED (CROSS)	No changes. Endothelial cells and vascular basement membranes are immunolabeled Long term changes (D 55 T 45 ) Increased tortuosity. Narrow channels without endothelia. Reduced RECA reactivity.
				Short Term changes (D 12 T 2 ). Gradual
		Tumor	CROSS	decrease of RECA reactivity. Loss of tumor endothelia
				(denudation)
		Normal	UNI	No changes
	Blood Volume Fraction (BVf)	Brain	CROSS	No changes
		Tumor	CROSS	Significant decrease from D 18 T 8
		Normal	UNI	No changes
	Vessel Size Index (VSI)	Brain	CROSS	Significant increase at D 55 T 45
		Tumor	CROSS	Significant decrease from D 18 T 8
		Normal	UNI	No changes
	Vessel Permeability	Brain	CROSS	No changes
		Tumor	CROSS	Short term increase at D 12 T 2 and D 15 T 5


Pierre Lubin
email: lubin@enscbp.fr

Stéphane Glockner
email: glockner@enscbp.fr

Olivier Kimmoun
email: olivier.kimmoun@centrale-marseille.fr

Hubert Branger
email: branger@irphe.univ-mrs.fr

Numerical study of the hydrodynamics of regular waves breaking over a sloping beach

Keywords: Navier-Stokes, Numerical simulation, Large Eddy Simulation, Two-phase flow, Breaking waves, Splash-up, Air entrainment, Vortices

In the last three decades, great improvements have been brought to the knowledge of the hydrodynamics and general processes occurring in the surf zone, widely affected by the breaking of the waves. Nevertheless, the turbulent flow structure is still very complicated to investigate. The aim of this work is to present and discuss the results obtained by simulating twodimensional breaking waves by solving the Navier-Stokes equations, in air and water, coupled with a dynamic subgrid scale turbulence model (Large Eddy Simulation, LES). First, the ability of the numerical tool to capture the crucial features of this complicated turbulent two-phase flow is demonstrated.

Numerical results are compared with experimental observations provided by Kimmoun and Branger [1]. Spilling/plunging breaking regular waves are con-

Introduction

Simulating the air entrainment phenomenon generated by breaking waves remains a major challenge for modern CFD tools. Numerous problems motivated by fundamental research and applications, from environmental and coastal engineering sciences, require an accurate description of wave breaking. Highly complex hydrodynamic features are usually encountered in the surf zone: transition from irrotational flow motion to high frequency turbulence, interacting with large-and small-scale interface deformations, from overturning and breaking of the waves to complex fractioning and coalescence of bubbles and droplets. A broad range of relevant length and time scales is thus involved in this multiphase turbulent flow, making it extremely complicated to investigate both experimentally and numerically.

The general knowledge and understanding of turbulence generated by breaking waves have greatly been improved in the last three decades. Broken waves involve motions of different types and scales, including large-scale coherent vortical motions and small-scale turbulence [START_REF] Battjes | Surf-zone dynamics[END_REF]. Once waves break, a large amount of energy is released and turned into turbulence [START_REF] Svendsen | vances in Coastal and Ocean Eng[END_REF]. Nevertheless, a lot of work still has to be done, some aspects suffering from a lack of efficiency or unsolved limitations. It is widely accepted that such complicated high-Reynolds number turbulent multiphase flow is inaccessible for Direct Numerical Simulations (DNS), Reynolds Averaged Navier-Stokes Equations (RANS) or Large Eddy Simulations (LES) therefore being effective alternatives.

Zhao and Tanimoto [START_REF] Zhao | Numerical simulation of breaking waves by large eddy simulation and vof method[END_REF] first applied the LES method to breaking waves and showed very promising results compared with experimental measurements, considering a two-dimensional configuration. Then some numerical works proved the ability of the LES method [START_REF] Christensen | Large eddy simulation of breaking waves[END_REF][START_REF] Hieu | Numerical simulation of breaking waves using a two-phase flow model[END_REF][START_REF] Zhao | Numerical simulation of breaking waves by multi-scale turbulence model[END_REF][START_REF] Watanabe | Three-dimensional vortex structures under breaking waves[END_REF][START_REF] Christensen | Large eddy simulation of spilling and plunging breakers[END_REF] to deal with spilling, plunging, strongly plunging and spilling/plunging breakers. The authors provided a great improvement in the numerical methods allowing a description with a very promising accuracy for both the free surface and the general behavior of the turbulent flow structures. Major contributions gave interesting insights into the hydrodynamics under broken waves. Nevertheless, as all the cited authors stated, air entrainment was not taken into account, although it was widely identified as an important feature to be considered. Indeed, most numerical surf zone studies are based on single phase flow assumption, ignoring the air phase for the convenience of computation. Moreover, comparisons led to the conclusion that the poor description of air/water mixing was responsible for the discrepancies observed due to coarse mesh grid resolutions.

Strong surface distortions, surface tension and air-water interactions were usually omitted. A generally good agreement was usually found, but some major differences were shown to occur when phase-averaged velocity and turbulent quantities were examined in aerated regions. Indeed, the resolution was mainly chosen to provide acceptable computation times with a sufficient accuracy for the free surface description considering the macro structures. This often resulted in missing the correct breaking point, which is very sensitive to any imperfections in wave generation [START_REF] Christensen | Large eddy simulation of spilling and plunging breakers[END_REF]. It was also found that the results of plunging breakers showed better agreement than spilling breakers compared to experimental measurements.

Many experimental works drew attention to the generation and the importance of air entrainment during the wave breaking process [START_REF] Miller | Role of vortices in surf zone predictions : sedimentation and wave forces[END_REF][START_REF] Nadaoka | Laboratory measurements of velocity field structure in the surf zone by ldv[END_REF][START_REF] Battjes | Surf-zone dynamics[END_REF][START_REF] Lamarre | Air entrainment and dissipation in breaking waves[END_REF][START_REF] Hwung | Energy dissipation and air bubbles mixing inside surf zone[END_REF][START_REF] Chanson | Unsteady air bubble entrainment and detrainment at a plunging breaker: dominant time scales and similarity of water level variations[END_REF]. The velocity field under broken waves is characterized by the existence of very active turbulence associated with air entrainment, which is responsible for wave energy damping in the surf zone. Lin and Hwung [START_REF] Lin | External and internal flow fields of plunging breakers[END_REF] showed that the main mechanism driving the motion in the bubble region was the vortex system generated during the jet-splash cycles. Experiments proved that the eddies contained a large quantity of air bubbles which enhanced the upwelling of sediment. Chanson and Lee [START_REF] Chanson | Plunging jet characteristics of plunging breakers[END_REF] measured the rate of energy dissipation to be increased with the bubble penetration depth. Huang et al. [START_REF] Huang | Turbulence and energy dissipations of surf-zone spilling breakers[END_REF] recently presented experimental measurements of spilling breaking waves. Significant turbulent dissipation was shown to occur initially in the roller region at the frontal wave crest and then to spread to the entire crest region after the establishment of a turbulent bore. The small-scale interactions thus have important effects on large-scale behavior [START_REF] Bradford | Progress in the observation and modeling of turbulent multi-phase flows[END_REF].

More recent numerical studies proved the importance of air entrainment for turbulence generation in numerical simulations of breaking waves. Christensen et al. [START_REF] Christensen | Vertical variation of the flow across the surf zone[END_REF] highlighted that since the mixture of air and water in the roller region had, on average, a smaller density than that of the water, the turbulence produced in the roller would have difficulties in penetrating the underlying fluid. Therefore, a large part of the production and dissipation took place in the roller before it was diffused downward, which explained the overestimation on the turbulence in the surf zone by numerical models so far. Moreover, the de-entrainment of air bubbles from the water after wave breaking may have released some wave energies into the air and may have contributed significantly to the wave energy dissipation process. Hieu and Tanimoto [START_REF] Hieu | Verification of a vof-based twophase flow model for wave breaking and wave-structure interactions[END_REF] improved their numerical model [START_REF] Hieu | Numerical simulation of breaking waves using a two-phase flow model[END_REF] showing better comparisons with experimental data concerning free surface elevations. The numerical results were also compared with those from Zhao et al. [START_REF] Zhao | Numerical simulation of breaking waves by multi-scale turbulence model[END_REF] and a higher accuracy was shown. They pointed out the major role of air entrainment and the necessity to take it into account to get accurate numerical results. Lubin et al. [START_REF] Lubin | Three-dimensional large eddy simulation of air entrainment under plunging breaking waves[END_REF] discussed the results obtained from simulating three-dimensional plunging breaking waves by solving the Navier-Stokes equations, in air and water. The splash-up mechanism was carefully detailed and vortex generation and air entrainment processes were described. The behavior of the gas pockets was analyzed and its impact on energy dissipation was shown. Recently, Wang et al. [START_REF] Wang | A coupled level set and volume-of-fluid method for sharp interface simulation of plunging breaking waves[END_REF][START_REF] Wang | Simulation of spilling breaking waves using a two phase flow cfd model[END_REF] improved the numerical methods for describing the free surface, dedicated to breaking waves. Iafrati [START_REF] Iafrati | Numerical study of the effects of the breaking intensity on wave breaking flows[END_REF] numerically analyzed the role played by the breaking intensity on free surface dynamics, air entrainment and general hydrodynamics. Clear improvements were made to the general understanding of the complicated flow. [START_REF] Kimmoun | A piv investigation on laboratory surf-zone breaking waves over a sloping beach[END_REF] recently experimented surf-zone breaking waves.

Kimmoun and Branger

Particle Image Velocimetry (PIV) experimental techniques were improved to be able to calculate velocities and void fractions in the aerated regions. Detailed pictures showed that a short spilling event occurred at the crest of the waves, before degenerating into strongly plunging breaker. Numerical works

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

usually showed better agreements when simulating plunging breaking waves than the spilling case compared with available experimental data. Fine mesh grid resolutions and appropriate numerical methods are required to accurately describe the length scales of the interface deformations experimentally measured (plunging jet, white foam, etc.).

Based on a previous numerical study, Lubin et al. [25] presented the results obtained for the LES of 2D and 3D regular waves shoaling and breaking over a sloping beach, compared with the experimental results of Kimmoun et al. [START_REF] Kimmoun | Laboratory piv measurements of wave breaking on a beach[END_REF]. A spilling/plunging breaking event was expected to occur according to the experimental measurements, but the numerical results showed discrepancies, due to the coarse mesh grid resolution. So, considering the new experimental results from Kimmoun and Branger [START_REF] Kimmoun | A piv investigation on laboratory surf-zone breaking waves over a sloping beach[END_REF] and the discrepancies highlighted by Lubin et al. [START_REF] Lubin | Large eddy simulation of regular waves breaking over a sloping beach[END_REF], 2D numerical simulation of the phenomenon has been performed to show the ability of the numerical tool to finally capture the spilling initiation process and air entrainment during the early stage of the phenomenon. The effect of air has not been studied yet in detail in most of the cited two-and three-dimensional numerical studies found in the literature. Taking the mixture of air and water into account is crucial and still remains one of the challenges of the coming years. The aim of this paper is to simulate this unsteady two-phase wave breaking motion using a LES method to gain further understanding of the complicated features of the flow, including wave overturning, occurrence of splash-ups and air entrainment. The limitations in the description of the air entrainment will be highlighted and discussed. Great care will be taken to highlight the importance of the mesh grid size. The paper is organized as follows. In section 2, we introduce the numerical model. Numerical wave generation procedure is detailed. Some preliminary numerical results and applications are provided to demonstrate the ability of the numerical tool to deal with breaking waves. We then provide, in section 3, a discussion of the numerical results obtained by simulating regular waves breaking over a sloping beach. This section is devoted to the qualitative comparison of our numerical results with the experimental observations [START_REF] Kimmoun | A piv investigation on laboratory surf-zone breaking waves over a sloping beach[END_REF]. The case study is carefully described and boundary conditions are given. We aim to accurately describe the free surface behavior, as we will focus on capturing and describing the spilling phase experimented by Kimmoun and Branger [START_REF] Kimmoun | A piv investigation on laboratory surf-zone breaking waves over a sloping beach[END_REF]. Finally, in section 4, a short critique of our numerical results is discussed. We then outline the perspectives and future works.

Description of the numerical model and validation

We solve the Navier-Stokes equations in air and water, coupled with a subgrid scale turbulence model (LES). The numerical tool is well suited to dealing with strong interface deformations occurring during wave breaking, for example, and with turbulence modeling in the presence of a free surface in a more general way. Solving the Navier-Stokes equations in an air/water configuration is still a real challenge, especially when dealing with strong interface deformations and tearing, turbulence and free surface interactions.

Governing equations

∇ • u = 0 (1)
ρ ∂u ∂t + u • ∇u = -∇p + ρg + ∇ • (µ + µ t ) ∇u + ∇ T u + F (2)
and

∂C ∂t + u • ∇C = 0 ( 3
)
where u is the velocity, C the phase function, t the time, p the pressure, g the gravity vector, ρ the density, µ the dynamic viscosity, µ t the turbulent viscosity and F the surface tension volume force.

To deal with solid obstacles within the numerical domain, it is possible to use multi-bloc domains, but it is often much simpler to consider the numerical domain as a unique porous medium [START_REF] Arquis | Sur les conditions hydrodynamiques au voisinage d'une interface milieu fluide -milieu poreux : application à la convection naturelle[END_REF][START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF][START_REF] Khadra | Fictitious domain approach for numerical modelling of navier-stokes equations[END_REF]. The permeability coefficient K defines the ability of a porous medium to let pass the fluids more or less freely through it. If this permeability coefficient is large (K → +∞), the medium is equivalent to a fluid. If it is zero, we can model an impermeable solid. A real porous medium is modeled with intermediate values of K. To take this coefficient K into account in our system of equations, we thus add an extra term in the right-hand side of the momentum equations (Eq. 2),

- µ K u.
The magnitudes of the physical characteristics of the fluids depend on the local phase. They are defined according to C in a continuous manner as:

ρ = Cρ 1 + (1 -C)ρ 0 µ = Cµ 1 + (1 -C)µ 0 (4)
where ρ 0 , ρ 1 , µ 0 and µ 1 are the densities and viscosities of fluid 0 and 1 respectively.

Based on the review of Lubin and Caltagirone [START_REF] Lubin | Large eddy simulation of the hydrodynamics generated by breaking waves[END_REF], we find that the most widely used subgrid scale model is the Smagorinsky model. However, it has been proved to be much too dissipative [START_REF] Sagaut | Large Eddy Simulation for incompressible flows -An introduction[END_REF]. In spite of its negative aspects, its simplicity is still widely appreciated. In Eq. 2), the turbulent viscosity µ t is calculated with the Mixed Scale model [START_REF] Sagaut | Large Eddy Simulation for incompressible flows -An introduction[END_REF], which has proved its accuracy for geophysical flows [START_REF] Lubin | Large eddy simulation of plunging breaking waves[END_REF][START_REF] Helluy | Numerical simulations of wave breaking[END_REF][START_REF] Lubin | Three-dimensional large eddy simulation of air entrainment under plunging breaking waves[END_REF][START_REF] Lubin | Numerical simulation of a weak breaking tidal bore[END_REF][START_REF] Lubin | Large eddy simulation of turbulence generated by a weak breaking tidal bore[END_REF]. The model exhibits a triple dependency on the large and small structures of the resolved field as a function of the cut-off length. The eddy viscosity µ t is calculated as follows (Eq. 5):

µ t (x, t) = ρC M ∆ 1+α |S| α 2 (q 2 c (x, t)) 1-α 2 ( 5
)
where S is the resolved deformation rate tensor and ∆ is the cut-off length of the filter. C M is the model constant chosen as C M = 0.06, α is a parameter with a value that varies between 0 and 1. Generally, and in the following, α is taken to be equal to 0.5. The quantity q c represents the kinetic energy of the test field extracted from the resolved velocity field through the application of a test filter associated to the cut-off lengthscale ∆ > ∆. We choose ∆ = 2∆, because it is the value which is most used and seems to give the best results.

This subgrid kinetic energy is assumed to be equal to the kinetic energy at cut-off q 2 c , evaluated in real space as (Eq. 6):

q 2 c (x, t) = 1 2 u(x, t) ′ u(x, t) ′ (6)
where the test field velocity (u ′ ) can be evaluated thanks to an explicit test filter applied to the resolved scales, noted ( .). This explicit discrete filtering operation is a linear combination of the neighboring values [START_REF] Sagaut | Large Eddy Simulation for incompressible flows -An introduction[END_REF][START_REF] Sagaut | Discrete filters for large eddy simulation[END_REF]. This test field velocity represents the high frequency part of the resolved velocity field.

The use of this model does not require a complementary wall model, as the Smagorinsky model does, because the eddy viscosity vanishes as the kinetic energy tends to zero at cut-off.

Model (Eqs. 1-4) describes all the hydrodynamic and geometrical processes involved in the motion of multiphase media.

Numerical methods

Time discretization is implicit and the equations are discretized on a staggered grid by means of the finite volume method. A dual grid, or underlying grid [START_REF] Rudman | A volume-tracking method for incompressible multifluid flows with large density variations[END_REF], is used to gain improved accuracy for the interface description, the mesh grid size being divided by two in each direction for interface tracking.

This technique also avoids the interpolations of the physical characteristics on the staggered grids, since the color function is defined on each point where viscosities and densities are needed.

The main difficulty concerning the numerical treatment of the Navier-Stokes equations (Eqs. 1-2) is the coupling between pressure and velocity and the fulfillment of the incompressibility condition. Following the work of Fortin and Glowinski [START_REF] Fortin | Méthodes de lagrangien augmenté. Application à la résolution numérique de problèmes aux limites[END_REF], the Navier-Stokes system is formulated as a velocity / pressure minimization -maximization problem requiring the computation of a saddle point (u, p) associated with the augmented Lagrangian of the problem. The pressure is considered as a Lagrange multiplier and the incompressibility constraint is introduced implicitly into the momentum equations (Eqs. 2). Then, the saddle-point (u, p) is computed using an iterative Uzawa algorithm [START_REF] Uzawa | Iterative method for concave programming[END_REF]. Parameters of the augmented Lagrangian method are calculated automatically according to the fluids and flow characteristics following [START_REF] Vincent | An adaptative augmented lagrangian method for three-dimensional multimaterial flows[END_REF]. The space derivatives of the inertial term are discretized by a hybrid Upwind-Centered scheme [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF] and the viscous term is approximated by a second order centered scheme. The MUMPS direct solver is used to solve the linear systems [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF][START_REF] Amestoy | Hybrid scheduling for the parallel solution of linear systems[END_REF].

Interface tracking is achieved by a Volume Of Fluid method (VOF), which is able to handle interface reconnections without interface reconstruction.

Lin and Liu [START_REF] Lin | Free surface tracking methods and their applications to wave hydrodynamics[END_REF] gave a complete overview and discussion of the different numerical techniques that have been used for interface tracking in numerical simulations of breaking waves. Given that Eq. 3 is hyperbolic and C is discontinuous, the explicit Total Variation Decreasing (TVD) Lax-Wendroff (LW) scheme of LeVeque [START_REF] Leveque | Numerical methods for conservation laws[END_REF] is used to directly solve the interface evolutions without the reconstruction of C. When the small scale structures of interface

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

are large compared with the grid size, the LW-TVD method is accurate and involves controlled numerical diffusion across the interface on three grid points.

The numerical code has already been extensively verified and validated through numerous test cases including mesh refinement analysis [START_REF] Lubin | Large eddy simulation of plunging breaking waves[END_REF][START_REF] Lubin | On the navier-stokes equations simulation of the head-on collision between two surface solitary waves[END_REF][START_REF] Lubin | Three-dimensional large eddy simulation of air entrainment under plunging breaking waves[END_REF][START_REF] Lubin | Large eddy simulation of turbulence generated by a weak breaking tidal bore[END_REF]. The accuracy of the numerical schemes and the conservation laws of mass and energy in the computational domain have been accurately verified.

For more details, the numerical methods have already been fully described in Lubin et al. [START_REF] Lubin | Three-dimensional large eddy simulation of air entrainment under plunging breaking waves[END_REF].

Internal wave maker

Prior to the simulation of the laboratory tests, an effort has been made

to implement and validate the procedure of regular and irregular wave generation developed by Lin and Liu [START_REF] Lin | Internal wave-maker for navier-stokes equations models[END_REF]. The method consists in introducing an internal mass source function in the continuity equation (Eq. 1) for a chosen group of cells defining the source region:

∇ • u = S(x, y, t) in Ω (7)
where S(x, y, t) is calculated thanks to any chosen analytical wave solution.

2.4. Validation: a solitary wave propagating and breaking over a submerged obstacle Lubin et al. [START_REF] Lubin | Three-dimensional large eddy simulation of air entrainment under plunging breaking waves[END_REF] validated the numerical tool considering some dam-break and rising bubbles configurations. In this section, the propagation and the overturning of a stable solitary wave over a submerged reef, consisting of a

ACCEPTED MANUSCRIPT

rectangular step-like obstacle, is investigated (Fig. 1). The reference case is taken from the work presented by Yasuda et al. [START_REF] Yasuda | Kinematics of overturning solitary waves and their relations to breaker types[END_REF]. This test case shows the ability of our numerical model to simulate overturning waves accurately, as we will compare the wave profiles at breaking. With respect to the positions of the three last gauges (referred to as the experimental values (Fig. 3). We can note that our numerical model gives very satisfactory results. The general trend is followed with a reasonable accuracy, the relative errors being presented in Table 1.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

We verify that our numerical model accurately reproduces the propagation and the interaction between the solitary wave and the obstacle.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

it is not our concern in this paper.

As recorded by the gauge standing over the reef's upstream corner (Fig. are approximately of the same length. However, we can observe that the computed jet looks "thicker" than the experimental one. This tendency is

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

generally observed when using VOF methods for simulating breaking waves [START_REF] Chen | Two-dimensional navier-stokes simulation of breaking waves[END_REF][START_REF] Lubin | Large eddy simulation of plunging breaking waves[END_REF][START_REF] Wang | A coupled level set and volume-of-fluid method for sharp interface simulation of plunging breaking waves[END_REF].

Unfortunately, Yasuda et al. [START_REF] Yasuda | Kinematics of overturning solitary waves and their relations to breaker types[END_REF] presented no experimental pictures giving information about the abscissa of the jet impact or the splash-up behavior. However, our numerical results fit very well with the published results of Yasuda et al. [START_REF] Yasuda | Kinematics of overturning solitary waves and their relations to breaker types[END_REF] considering the wave breaking initiation phenomenon.

The ability of the numerical tool to predict accurately the breaking point initiation is thus proved.

Large Eddy Simulations of breaking waves

Based on the numerical methods detailed in section 2, Lubin et al. [START_REF] Lubin | Large eddy simulation of regular waves breaking over a sloping beach[END_REF] presented the results obtained for the LES of regular 2D and 3D waves shoaling and breaking over a sloping beach, compared with the experimental results from Kimmoun et al. [START_REF] Kimmoun | Laboratory piv measurements of wave breaking on a beach[END_REF]. The main observed differences were that the first short spilling event was missed and the dislocation of the gas pockets into small bubbles could not be simulated, even though, in the numerical results, the gas pockets corresponded to some air-water mixing zones observed in the experimental pictures. In order to overcome the observed discrepancies, mainly due to the coarse mesh grid resolution, new numerical simulation was performed considering the new experimental study from Kimmoun and

Branger [START_REF] Kimmoun | A piv investigation on laboratory surf-zone breaking waves over a sloping beach[END_REF]. The mesh grid distribution was improved to be able to capture the spilling phase.

Description of the experimental configuration

The experiments were performed in the École Centrale wave tank in Mar-

Initial and boundary conditions

The computational domain is 15 m long and 1 m high (Fig. 5). The sloping beach starts at x = 3.5 m, the source function being located at

x S = 3.5 m and z S = 0.3675 m. The center of the source region is at d/2, right above the toe of the sloping beach to save computation time.

The numerical beach is considered as an impermeable solid obstacle, the permeability coefficient K being initialized at zero (Eq. 2).

522 000 mesh grid points are used to discretize the numerical domain, with nonuniform grids in both directions (∆x min ≃ 1.10 -3 m and ∆z min ≃ 2.5 × 10 -3 m). These values have to be divided by two, in both directions, for the free surface description thanks to the dual grid. Two-phase flow simulations and turbulence modeling require fine mesh cells to be very accurate, even if LES is supposed to save some mesh grid points. Moreover, when simulating two-phase flows, the interface can become smaller than the mesh grid size (droplets or bubbles). These small inclusions are thus "lost" and their contribution to the flow can lead to a wrong description of the flow. In our study, the major physical parameters leading the choice of our mesh grid distribution are the length scales of the free surface deformations, especially the overturning jet and the entrained gas pockets.

The time step is chosen to ensure a Courant-Friedrichs-Levy condition less than 1, necessary for the explicit advection of the free surface. The calculation is made with the densities and the viscosities of air and water (ρ a = 1.1768 kg.m -3 and ρ w = 1000 kg.m -3 , µ a = 1.85 × 10 -5 kg.m -1 .s -1

and µ w = 1×10 -3 kg.m -1 .s -1 ). Table 2 presents the former and new physical and numerical parameters, used in the successive studies. The mesh grid sizes have been divided by two considering the dual grid used for the interface capture (see section 2.2).

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

We adapted Fenton's [START_REF] Fenton | A fifth-order stokes theory for steady waves[END_REF] analytical developments for the fifth-order Stokes wave theory, corresponding to the experimental conditions, to the source function method (Eq. 7):

S(x, y, t) = i j 2c Ak B ij ǫ i cos i π 2 -σt -p s ( 8
)
where A is the area of the source region. c, k and σ are the celerity, wave number and angular frequency respectively. ǫ = kH 2 is the wave steepness and B ij are dimensionless coefficients. The detailed expressions for the dispersion relationship, wave celerity and B ij coefficients are given by Fenton [START_REF] Fenton | A fifth-order stokes theory for steady waves[END_REF], Sobey et al. [START_REF] Sobey | Application of stokes, cnoidal, and fourier wave theories[END_REF] or Fenton and McKee [START_REF] Fenton | On calculating the lengths of water[END_REF]. p s is the phase shift calculated to have S(x, y, 0) = 0. The source region is 0.06 m wide and 0.0735 m high. The area and the location of the source function have been designed

applying the rules described by Lin and Liu [START_REF] Lin | Internal wave-maker for navier-stokes equations models[END_REF]. The height of the source region is recommended to be in the range of 1/10 -1/2 of water depth, while the width of the source region is suggested to be less than 5% of the wavelength. The distance between the center of the source region and the still water level is advised to be in the range of 1/3 -1/2 of the water depth. The method has been extensively verified and validated compared with analytical profiles to ensure accurate wave generation.

At the beginning of the numerical simulation, the water is at rest, and, as the time increments, the source function is calculated to generate the regular waves. In Fig. 6, we show the velocity field around the source function, acting like a pump. Two trains of surface gravity waves are thus generated, as the free surface responds to a pressure increment defined within the source region cells.

ACCEPTED MANUSCRIPT

a. S(x, z, t) > 0 b. S(x, z, t) < 0

Results and discussion

Experimentally, the white cap was about 1 -mm thick as shown by Kimmoun and Branger [START_REF] Kimmoun | A piv investigation on laboratory surf-zone breaking waves over a sloping beach[END_REF] in their Figs. 4 -windows 2 and 3. So, numerically, the mesh grid refinement and the distribution were chosen to be able to capture this feature. In Figs. 7, 8, 10 and 11, we present the free surface evolution with the numerical velocity field in the water, corresponding to the 10th breaking wave. In Figs. 12, we present some snapshots of the free surface evolution with the numerical velocity fields in both media, air and water, corresponding to the same 10th breaking wave.

General flow description

Spilling breakers are observed to start as a small zone of bubbles and droplets on the forward side of the crest [START_REF] Duncan | Spilling breakers[END_REF][START_REF] Kimmoun | A piv investigation on laboratory surf-zone breaking waves over a sloping beach[END_REF]. This small region then grows by spreading downslope, most of the forward face becoming a turbulent flow region. Duncan [START_REF] Duncan | Spilling breakers[END_REF] reported that, for the long wavelengths considered here, spilling breakers can be initiated by a small jet at the crest of the wave, creating a small turbulent patch of fluid well above the mean water level.

Figs. 7 present the initiation of the spilling phase. Once the front face of the crest steepens and becomes vertical (Fig. 7 a), a thin jet of water is indeed observed to be projected (Fig. 7 b). In our numerical results, the spilling phase then starts as a very weak plunging breaking wave, with a small tongue of water thrown from the crest developing and free-falling down forward into a characteristic overturning motion. This is in accordance with reports from Duncan [START_REF] Duncan | Spilling breakers[END_REF] and Kimmoun and Branger [START_REF] Kimmoun | A piv investigation on laboratory surf-zone breaking waves over a sloping beach[END_REF]. Miller [START_REF] Miller | Role of vortices in surf zone predictions : sedimentation and wave forces[END_REF] also presented pictures showing a spilling breaking wave starting with a 5 -mm high impinging jet at the crest of the wave. Once the jet is ejected from the wave crest and plunges down, it hits the water at the plunge point, located very near the crest of the wave. The plunging jet closes over a small gas pocket (Fig. 8 c). The resulting splash is directed down the wave leading to a spilling breaker (Fig. 8 d), where white foam, consisting of a turbulent air/water mixing, should appear at the wave crest and spill down the front face of the propagating wave [START_REF] Duncan | Spilling breakers[END_REF][START_REF] Kimmoun | A piv investigation on laboratory surf-zone breaking waves over a sloping beach[END_REF].

However, compared with the 1 -mm thick layer of foam initiating the experimental spilling wave from Kimmoun and Branger [START_REF] Kimmoun | A piv investigation on laboratory surf-zone breaking waves over a sloping beach[END_REF], we find that the gas pockets entrapped are about 5 -mm thick. Fig. 9 presents the jet appearing at the wave crest and the mesh grid density used to capture the generation of the jet. The jet is clearly very small, about 1 -cm high, but it is probably thicker than experimentally. This is probably due to the same tendency detailed in section 2.4, where the plunging jet was also thicker than experimentally measured.

Then, the spilling breaker transitions into a strong plunging breaker (Fig. 10 a). A jet is ejected farther towards the lower part of the face of the steepening wave. The plunging breaking wave is then responsible for the generation of larger jet-splash cycles, this, in turn, being responsible for the generation of a sequence of large-scale coherent vortices (Fig. 10 b).

A high velocity region is located at the breaking crest of the wave. High

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

due the jet-splash cycles, during the bore propagation. The free surface is distorted and very dynamic. The bore running up the beach at some point meets the flow running down from the previous broken wave.

Nevertheless, discrepancies again appear. Experimentally, air cavities are observed to be quickly fragmented into large plumes of bubbles. The dislocation of the gas pockets into small bubbles cannot be simulated, even if, in the numerical results, the gas pockets correspond to some air-water mixing zones observed in the experimental pictures. This is due to the mesh grid resolution, which is still too coarse to be able to capture this flow feature.

Indeed, the order of magnitude of bubble radii is usually 10 -4 m [62], whereas our mesh grid resolution is ∆x min ≃ 5.10 -4 m and ∆z min ≃ 1.25 × 10 -3 m thanks to the dual grid. We are able to track the largest gas pockets and bubbles greater than 1 mm, as illustrated in the figures presented in this paper, where a large variety of inclusion length scales can be seen. Turbulence is associated with air entrainment, which is responsible for wave energy damping in the surf zone. In the experiments, it appears that the entrained air bubbles are contained mostly in the large structures and diffused towards the bottom due to the eddies. The rate of energy dissipation is increased with the bubble penetration depth and strong vertical motion is induced by the rising air bubbles. These mechanisms are mostly 3D, which cannot be taken into account in a 2D numerical simulation. Moreover, numerical diffusion due to the LW-TVD method [START_REF] Leveque | Numerical methods for conservation laws[END_REF] must also be responsible for the discrepancies. The free surface is diffused over three mesh grid cells, so any small inclusions (droplets or bubbles) are described less precisely.

We can also see strong flow dynamics in the air. We can see that the Rotating and fractioning air pockets can be clearly seen. A large range of inclusions can be described by the numerical simulation, interacting with complex free surface behavior and chaotic hydrodynamics.

ACCEPTED MANUSCRIPT

Quantitative comparisons

Table 3 compares the locations of the spilling breaking point, the plunging occurrence and the first splash-up generation, from the numerical and experimental results. The shoreline is located at x S = 14.075 m (see Fig. The general trend of the flow is found to be qualitatively well described.

Experimental [START_REF] Kimmoun | A piv investigation on laboratory surf-zone breaking waves over a sloping beach[END_REF]). Distances are given from the shore. Figure 14 compares the numerical and experimental wave profiles and velocity fields for the main events. Despite the discrepancies in the abcissae, a general agreement can be seen in the velocity magnitudes and free surface description. Figure 15 presents the numerical and experimental wave profiles when the broken wave runs-up along the beach. Air entrainment can be seen in the experimental pictures with the gray areas. As already mentioned, the gas pockets observed in the numerical pictures correspond to some air-water mixing zones presented in the experimental pictures, but the coarse mesh grid size does not allow the dislocation of the pockets into bubbles. Splash-ups cycles are well described, with horizontal eddies generation [START_REF] Lubin | Three-dimensional large eddy simulation of air entrainment under plunging breaking waves[END_REF][START_REF] Kimmoun | A piv investigation on laboratory surf-zone breaking waves over a sloping beach[END_REF].

The general trend of the flow dynamics is found to be correctly simulated. This remark raises an issue that should be addressed: how many wave cycles should be simulated to be compared with experimental results? For example, Ting and Kirby [START_REF] Ting | Observation of undertow and turbulence in a laboratory surf zone[END_REF][START_REF] Ting | Dynamics of surf-zone turbulence in a strong plunging breaker[END_REF][START_REF] Ting | Dynamics of surf-zone turbulence in a spilling breaker[END_REF] indicated that waves were generated for a minimum of 20 mn before data were taken, which gives 250 to 600 periodic waves (considering plunging and spilling breakers, respectively). Lin

and Liu [START_REF] Lin | A numerical study of breaking waves in the surf zone[END_REF] simulated about ten periodic waves considering spilling breakers [START_REF] Ting | Dynamics of surf-zone turbulence in a strong plunging breaker[END_REF]. The computed results for mean velocities and free surface profiles indicated that the waves computed in the surf zone have nearly, but not completely, reached the quasi-steady state, i.e. there was very a small difference between two successive wave cycles. Lin and Liu [START_REF] Lin | Discussion of vertical variation of the flow across the surf zone[END_REF] confirmed that, with the improvement of the wave generation mechanism and absorbing boundary conditions, the RANS models could run a long time until the waves reached the quasi-steady state, during which both the mean water level and time-averaged mean flow field could be accurately computed. However, the turbulence closure models did not give completely satisfactory results for turbulence simulation in the surf zone, the errors ranging from 25% to 100% for the turbulence calculation in the spilling breaker. With LES, Christensen [START_REF] Christensen | Large eddy simulation of spilling and plunging breakers[END_REF] computed the experimental configurations of Ting and Kirby [START_REF] Ting | Dynamics of surf-zone turbulence in a strong plunging breaker[END_REF][START_REF] Ting | Dynamics of surf-zone turbulence in a spilling breaker[END_REF], calculating average quantities over the last five wave periods considering 16 to 20

wave periods (plunging and spilling breakers, respectively). Wang et al. [START_REF] Wang | Simulation of spilling breaking waves using a two phase flow cfd model[END_REF] simulated up to ten spilling waves [START_REF] Ting | Dynamics of surf-zone turbulence in a strong plunging breaker[END_REF], calculating average quantities over the last two wave periods. Hieu and Tanimoto [START_REF] Hieu | Verification of a vof-based twophase flow model for wave breaking and wave-structure interactions[END_REF] set a computation time of 50 s, checking that time profiles of free surface elevations were almost stable after 30 s. When compared with to experimental results [START_REF] Ting | Dynamics of surf-zone turbulence in a strong plunging breaker[END_REF][START_REF] Ting | Dynamics of surf-zone turbulence in a spilling breaker[END_REF], better agreement was shown for plunging breakers. Nadoaka et al. [START_REF] Nadaoka | Structure of the turbulent flow field under breaking waves in the surf zone[END_REF] discussed the fact that the breaking point is sensitive to any imperfection in wave generation or to the effect of the previous broken wave. Therefore, the different wave generation procedures used between experiments and numerical simulations may also have contributed to the discrepancies observed in the comparisons presented in this paper. High order moments, such as skewness and flatness of the velocity signal, probably have to be estimated to check that steady state is reached in numerical simulations, to ensure accurate turbulent quantities.

It has therefore been proved that our simulation of fifteen wave periods was insufficient to get a stable and accurate breaking point. Indeed, we observed that the last three waves broke as plunging breakers. So more wave periods would be needed to clearly address this issue. Thus, we did not calculate average quantities, considering the expected discrepancies which could have resulted from the comparisons. Considering these simulations as a validation step, our numerical model gives very satisfactory and encouraging results for this 2D configuration.

Conclusions and perspectives

We focused on describing the spilling phase of the experiments detailed by Kimmoun and Branger [START_REF] Kimmoun | A piv investigation on laboratory surf-zone breaking waves over a sloping beach[END_REF]. The numerical results presented in this paper concern instantaneous quantities, simulating 2D regular waves breaking over a sloping beach. Our model was found to be reliable to describe correctly the complicated two-phase flow interactions that happen when waves break.

The breaking process, in terms of wave overturning and splash-up occurrence, is in accordance with the general observations given in the literature. Air entrainment is described, which is important as it plays a great role in the energy dissipation process. The utility of the numerical approach is to provide a complete and accurate description of free surface and velocity evolutions in both air and water media during the breaking of the waves, which must lead to the understanding of the generation processes of energy dissipation and turbulent flow structures. Nevertheless, wave breaking is a 3D two-phase turbulent problem, so the 2D numerical results presented here consisted in a first attempt.

A major limitation of our numerical methods has been illustrated in the results presented in this paper. The accuracy of the whole physical process description is closely linked to the mesh grid size. Although LES allows the use of coarser mesh grids to simulate turbulent flows, on the condition that an appropriate subgrid model is implemented, an accurate free surface description requires fine mesh grids. LES of breaking waves involve large numerical domains, to include the generation of regular waves propagating towards sloping beaches, and fine mesh grid resolutions to be able to describe a large variety of physical processes (overturning jet, air entrainment, gas pocket dislocations and coalescences). Long CPU time for the calculations is also required to simulate a sufficient number of wave periods to reach a steady state in the surf zone to perform correct averaging over time. Parallel computing will enhance access to a better level of description of the turbulent behavior of the entrained and rising air bubbles, providing refined mesh grids are used to ensure an accurate free surface description [START_REF] Lubin | Large eddy simulation of turbulence generated by a weak breaking tidal bore[END_REF].

With mesh refinment and 3D simulations, we could expect a better description of the air entrainment and energy dissipation. Considering the difficulty to capture all the small bubbles and droplets encountered in the wave breaking problem, these small interface structures can also be considered as subgrid interfaces. Mimicking the LES approach for turbulence, an appropriate model of inclusions smaller than the mesh grid size could then be proposed [START_REF] Shi | Modeling quiescent phase transport of air bubbles induced by breaking waves[END_REF].

  An incompressible multiphase phase flow between non-miscible fluids can be described by the Navier-Stokes equations in their multiphase form. In the single fluid formulation of the problem [27], a phase function C, or "color" A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT function, is used to locate the different fluids standing C = 0 in the outer medium, C = 1 in the considered medium. Intermediate values of C indicate the proportion of the medium (water) in the control volume around each node of the mesh. The interface between two media is repaired by the discontinuity of C between 0 and 1. In practice, C = 0.5 is used to characterize this surface. The governing equations for the Large Eddy Simulation (LES) of an incompressible fluid flow are classically derived by applying a convolution filter to the unsteady Navier-Stokes equations. The resulting set of equations reads (Eqs. 1-3):

  The mean water depth is d = 0.31 m and the height of the obstacle is D = 0.263 m. The initial amplitude and celerity are H = 0.131 m and c = 2.072 m.s -1 , respectively. The incident solitary wave is initialized in the numerical domain and the crest of the wave is set at x = 2 m, the internal wave maker is not used in this test case. The reef face is located at x = 4 m.The velocity field and the free surface profile are initially calculated with the theoretical third-order solitary wave solution[START_REF] Lee | Measurements of velocities in solitary waves[END_REF][START_REF] Lubin | Propagation of solitary waves in constant depths over horizontal beds[END_REF]. The numerical domain is 8 m long and 0.6 m high, discretized into 1200 × 200 nonuniform grids in the x-direction, with ∆x min ≃ 4.10 -3 m in the vicinity of the submerged obstacle. Uniform grid spacing ∆z min ≃ 3.10 -3 m is used. The solitary wave propagates towards the right side of the numerical domain. Yasuda et al. [49] used a fully nonlinear potential theory model to study the internal velocity and acceleration fields and their relationship to breaker type. They verified the accuracy of their numerical model with some experimental data. Four capacitance wave gauges recorded the free surface elevations at different locations in the canal. The first gauge was placed upstream of the reef face to check the incident waveheight. The second gauge was located above the upstream corner of the reef. The last two gauges measured the free surface at 0.515 m and 1.020 m away from the tip of the obstacle.

Figure 1 :

 1 Figure 1: Initial conditions for the propagation and the breaking of a solitary wave over a submerged rectangular reef. t = 0 s, C > 0.5.

Figure 2 :

 2 Figure 2: Solitary wave propagating over the submerged reef: initiation of the fission of the wave. The black line represents the reef profile. C > 0.5.

Figure 3 :

 3 Figure 3: Comparisons of the non-dimensional free surface elevations, η/d, plotted versus the non-dimensional time, t g/d, recorded at each location P2, P3 and P4. Black line: present numerical results; red •: experimental results [49].

3 )

 3 , the free surface elevation rises from the initial non-dimensional value η/d = 0.424 to η/d = 0.498. The wave loses symmetry aspect and starts steepening till the front face of the crest becomes vertical (Fig. 4 a). We then have a plunging breaking wave (Fig. 4 b). A jet of liquid is about to be projected from the crest of the wave. It free-falls down forward in a characteristic overturning motion. It was observed that several gas pockets were entrapped. In order to compare our results with the work of Yasuda et al. [49], we use the same non-dimensional values. Yasuda's reef face is located at x/d = 32, our reef being at x/d ≃ 12.9. So, we show in figures (4 (a) and (b)) our free surface profiles translated to the same reference, with η/d the non-dimensional free surface elevation. Figs. (4 a and b) present the non-dimensional free surface profiles of the breaking wave at two different locations. At first glance, the comparisons are not to our advantage, but we voluntary zoomed the graphs, which enlarges the discrepancies. If we consider the abscissa of the vertical front face of the wave (Fig. 4 a), the relative error between our numerical result and the experimental measurement is ≃ -1.5 %. At the instant of jet fall initiation (Fig. 4 b), if we again consider the abscissa of the vertical face of the wave, under the overhanging jet, the relative error between the numerical result and the experimental measurement is ≃ -1.2 %. The numerical and experimental jets of water

Figure 4 :

 4 Figure 4: Comparison of the free surface profiles (a) prior to the jet ejection and (b) at jet initiation. Black : present LES (C = 0.5); black •: experimental results [49]. The variables indicated on both axis are those shown in Figure 4 by Yasuda et al. [49].

Figure 5 :

 5 Figure 5: Numerical domain configuration. The small rectangular box is the wave generator, located right above the toe of the sloping beach. The dashed line shows the initial water depth. The gray box on the left side of the numerical domain is a sponge layer. The slanted line shows the sloping beach.

Figure 6 :

 6 Figure 6: Velocity fields generated by the presence of the source function.

Figure 10 :

 10 Figure 10: Wave in the process of breaking with velocity field in water. Only one vector over two is shown. Tenth breaking wave transition to strong plunging breaker with splash-ups generation and air entrainment (continuation of Figs. 7 and 8). The free surface profile corresponds to C = 0.5.

5 ).

 5 The wave height at breaking is H b = 13 cm, compared with the experimental value H b = 14 cm. The numerical breaking point is located at x b ≃ 358 cm away from the shoreline, compared with the experimental value x b ≃ 250 cm. The occurrence of breaking was numerically detected when the horizontal velocity component reached a magnitude higher than the wave celerity, matching with a vertical topography of the front face of the wave (see Fig.10 a). The sequence of the main events can be qualitatively well described by the numerical tool, as the same distance separates approximately each event.

Figure 13 compares

 13 Figure 13 compares the numerical and experimental results for the crestenvelopes (locations of the wave maxima) and the trough-envelopes (locations of the wave minima). Arrows indicate the locations of the main events (beginning of breaking, beginning of plunging and beginning of the first splash-up).

3. 3 . 3 .Figure 14 :Figure 15 :

 331415 Figure 14: Comparison for the main breaking events (spilling sb, plunging pb and first splash-up s1). Right column: numerical results; left column: experiments.

Table 1 :

 1 Maximum non-dimensional free surface elevation values, η/d, recorded at the three considered locations. Numerical values are compared with the experimental measurements [49].

	Gauges Numerical Experimental Error (%)
	P2	0.500	0.505	-1.0
	P3	0.348	0.344	1.2
	P4	0.358	0.354	1.1

P2, P3 and P4 by Yasuda et al.

[START_REF] Yasuda | Kinematics of overturning solitary waves and their relations to breaker types[END_REF]

), we plot our non-dimensional free surface elevations η/d versus non-dimensional time t g/d, compared with

Table 2 :

 2 Physical and numerical parameters used in the successive studies.

	Experimental studies	Kimmoun et al. [26] Kimmoun and Branger [1]
	d(m)	0.735	0.705
	T (s)	1.3	1.275
	a(m)	0.07	0.057
	L(m)	2.5	2.4
	H b (m)	0.137	0.14
	x b (m)	2.65	2.5
	Numerical comparisons	Lubin et al. [25]	Present study
	Sizes of domains	20m × 1.2m	15m × 1m
	Mesh grid points	520 000	522 000
	∆x min (m)	5.10 -3 m	≃ 5.10 -4 m
	∆z min (m)	≃ 1.25 × 10 -3 m	≃ 1.25 × 10 -3 m

Table 3 :

 3 Experimental and numerical values for the approximate locations of the beginning of the main events: beginning of spilling breaking sb, beginning of plunging breaking pb and beginning of first splash-up s1 (see Figs. 12 and 13 from Kimmoun and Branger

			study Numerical simulation
	H b (cm)	14	13
	sb(cm)	≃ -250	≃ -358
	pb(cm)	≃ -210	≃ -280
	s1(cm)	≃ -170	≃ -250

ACCEPTED MANUSCRIPT

Acknowledgements

The authors wish to thank the Aquitaine Regional Council for the financial support towards a 256-processor cluster investment, located in the TREFLE laboratory. This work was granted access to the HPC resources of CINES under allocation 2009-c2009026104 made by GENCI (Grand Equipement National de Calcul Intensif). The authors would like to acknowledge the financial and scientific support of the French INSU -CNRS (Institut National des Sciences de l'Univers -Centre National de la Recherche Scientifique) program IDAO ("Interactions et Dynamique de l'Atmosphère et de l'Océan").

This work was carried out under the project "Wave-induced circulation in the surf zone" coordinated by Dr. Hervé Michallet (LEGI -Grenoble).

A sketch of a wave breaking event is displayed by Kimmoun and Branger [START_REF] Kimmoun | A piv investigation on laboratory surf-zone breaking waves over a sloping beach[END_REF]. The wave starts breaking showing a brief spilling phase, the white cap was observed to be about 1 -mm high. Then a jet of liquid is rapidly ejected from the wave crest and the overturning wave front curls forward. A first splash-up is generated when the jet of liquid hits the front face of the wave.

We can then see a large amount of air entrained with foam and bubbles. Some other splash-ups are then generated. A roller propagates towards the shoreline, with a great air-water mixing area. It can be seen that the bubbles are generated in the upper part of the water column, and advected towards the bottom with a slight slanting axis. The volume of the entrained bubbles decreases gradually till the wave crosses the shoreline and runs up before coming back. This is in agreement with the general description of Peregrine [START_REF] Peregrine | Breaking waves on beaches[END_REF], for example. More details of the experiments are given by Kimmoun and Branger [START_REF] Kimmoun | A piv investigation on laboratory surf-zone breaking waves over a sloping beach[END_REF]. b andc), as observed by Bonmarin [START_REF] Bonmarin | Geometric properties of deep-water breaking waves[END_REF], Miller [START_REF] Miller | Role of vortices in surf zone predictions : sedimentation and wave forces[END_REF] or Sakaï et al. [START_REF] Sakai | Vortex formation in plunging breaker[END_REF]. This process is observed to repeat, each successive splash-up being weaker than the preceding one. Large volumes of air and water are put into rotation, the successive rebounds causing more entrainment of air and energy dissipation. Miller [START_REF] Miller | Role of vortices in surf zone predictions : sedimentation and wave forces[END_REF] indicated the formation of many air bubbles during the jet-splash cycles which illustrated that the vortices in plunging breakers significantly affect bottom flow. Large spinning gas pockets are observed in our numerical results, dissipating wave energy. Most of the rotating structures are entrained to the bottom before rising again to the surface. And once the wave moves into a bore propagating towards the shoreline (Figs. 11 d to f), we have very similar results for the free surface description detailed by Kimmoun and Branger [START_REF] Kimmoun | A piv investigation on laboratory surf-zone breaking waves over a sloping beach[END_REF]. High velocities are located near the free surface,

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT