Influence of cropping system on Fusarium head blight and mycotoxin levels in winter wheat

Keywords: Cropping systems, mycotoxins, deoxynivalenol (DON), Fusarium, Fusarium head blight, wheat

In this study, we investigated the effects of various wheat cropping systems on Fusarium head blight severity and mycotoxin (deoxynivalenol, nivalenol, zearalenone) levels. We investigated four wheat-cropping systems: a conventional system, an integrated system, a direct drilling system and an organic system, in fields with natural contamination. The study was carried out over three one-year periods: 1999-2000, 2000-2001 and 2001-2002. The severity of Fusarium head blight and toxin levels differed considerably between years, reflecting climatic effects. In a given year, the severity of Fusarium head blight also differed between cropping systems, with the most severe disease observed in the direct drilling and conventional system plots in 2000. The level of grain contamination by mycotoxins therefore depends on both climate and cropping system. Mycotoxin levels were highest in the year with the highest disease severity: 2000. Contamination levels were highest with the direct drilling system. We were unable to rank the organic and conventional systems, because neither was consistently more contaminated than the other. Moreover, no clear relationship was found between disease severity and levels of contamination with deoxynivalenol, nivalenol or zearalenone under conditions of natural contamination.

Introduction

In recent years, the food industry has been faced with several major health crises. This has led to an increase in public concern about food safety. Several mycotoxins cause major safety problems in cereal production [START_REF] Boulc'h | Un test innovant pour contôler la qualité sanitaire des céréales : détection des champignons potentiellement producteurs d'une famille de mycotoxines, les trichothécènes[END_REF]. The ingestion of mycotoxins may result in vomiting, reproductive disturbances, leukoencephalomalacia, pulmonary oedema, impairment of the humoral and cellular immune responses, nervous disorders, myocardial hypertrophy and several cancers [START_REF] Placinta | A review of worldwide contamination of cereal grains and animal feed with Fusarium mycotoxins[END_REF]. Some of these toxins are produced before harvest (trichothecenes, fumonisins and zearalenone), whereas others are produced after harvest (ochratoxin A and aflatoxins). Several species of Fusarium are responsible for toxin production before harvest [START_REF] Maurin | Espèces de Fusarium présentes dans les cultures de céréales en France[END_REF]. These fungi, together with Microdochium, are responsible for Fusarium head blight, or scab [START_REF] Cassini | Facteurs favorables ou défavorables au développement des fusarioses et septorioses du blé[END_REF][START_REF] Ban | International symposium on wheat ingouvernement for scab resistance. Section II. Genetic analysis and gene mapping for scab resistance[END_REF].

As the deleterious effects of this disease include yield losses of up to 50% [START_REF] Atanasoff | Fusarium blight (scab) of wheat and other cereals[END_REF][START_REF] Parry | Fusarium ear blight (scab) in small grain cereals -a review[END_REF] and the production of various mycotoxins: zearalenone, deoxynivalenol and derivatives, nivalenol, T2-toxin, HT2-toxin and other toxins [START_REF] Placinta | A review of worldwide contamination of cereal grains and animal feed with Fusarium mycotoxins[END_REF].

Climate is thought to be a key determinant of Fusarium head blight [START_REF] Langseth | The influence of cultivars, location and climate on deoxynivalenol contamination of Norwegian oats 1985-1990[END_REF][START_REF] Lacey | Effects of infection time and moisture on development of ear blight and deoxynivalenol production by Fusarium spp. in wheat[END_REF]. However, in a given year, neighbouring fields may display different levels of contamination: the cropping system may therefore affect the level of mycotoxin production by Fusarium. The use of fungicides during growth of the wheat crop does not guarantee the absence of toxins in the harvested grain. The entire cropping system must be considered. Several studies have demonstrated effects of soil tillage [START_REF] Teich | Effect of cultural practices, soil phosphorus, potassium, and pH on the incidence of Fusarium head blight and deoxynivalenol levels in wheat[END_REF][START_REF] Dill-Macky | The effect of previous crop residues and tillage on Fusarium head blight of wheat[END_REF], preceding crop [START_REF] Teich | Survey of Fusarium head blight and possible effects of cultural practices in wheat fields in Lambton County in 1983[END_REF][START_REF] Dill-Macky | The effect of previous crop residues and tillage on Fusarium head blight of wheat[END_REF] and wheat genotype [START_REF] Mesterhazy | Breeding wheat for resistance to Fusarium graminearum and Fusarium culmorum[END_REF][START_REF] Teich | Effect of cultural practices, soil phosphorus, potassium, and pH on the incidence of Fusarium head blight and deoxynivalenol levels in wheat[END_REF][START_REF] Snijders | Fusarium head blight and mycotoxin contamination of wheat, a review[END_REF][START_REF] Dill-Macky | The effect of previous crop residues and tillage on Fusarium head blight of wheat[END_REF] on scab development and/or mycotoxin contamination.

However, few studies have analysed both disease severity and mycotoxin contamination.

There is also a lack of studies comparing the effects of entire cropping systems rather than aspects of crop management. The aim of this study was to investigate the effects of four cropping systems (conventional, integrated, direct drilling and organic systems) on Fusarium head blight severity and on levels of several mycotoxins in grains. We also tried to correlate both disease severity and levels of contamination with deoxynivalenol, nivalenol and zearalenone levels under conditions of natural contamination.

Materials and methods

Experimental design

were compared using "crop X -wheat -crop Y -wheat" crop rotations. We assigned two large plots (0.5 ha each) to each cropping system. Each plot was divided in half. On each half (subplot), we organised the rotation such that there was one wheat crop every year on each plot. The conventional system corresponded to intensive agriculture, in which the aim is to achieve the highest net income by maximising yield. The integrated system, in which input levels (and yields) were lower, aimed to maximise net income whilst preserving the environment. In the direct drilling system, crops had been planted without soil tillage and permanent vegetation cover (fescue) had been maintained since 1998-1999. We evaluated this system as an alternative to the integrated system, aiming to achieve the same objectives. No chemical inputs were used in the organic system. Each cropping system was defined in terms of a specific combination of technical choices, each closely related to the others. For example, the choice of cultivar depended on whether a fungicide was to be used. Thus, the cropping systems differed in several aspects: cultivar, tillage, chemical treatments, nitrogen fertilisation, time and density of sowing, and crop sequence. Furthermore, the techniques used within a given cropping system were prone to change from year to year, and depended on climate in particular. As part of a larger study, we compared the agronomic, environmental and economic aspects of the various systems (results not shown). We began to study Fusarium head blight in these experimental plots in 1999/2000. We set up additional treatments in each system to enable us to increase our understanding of the determinants of pathogen attack: the Charger cultivar (the cultivar used in the conventional system) was cropped without fungicide in a subplot of each system.

Over the three-year study period, the preceding crop and fungicide treatment changed each year for each system. In 1999-2000, the fungicide treatments used contained strobularins, whereas in 2000-2001 and 2001-2002, a mixture of strobularins and triazoles was used. Each system remained internally consistent from year to year and the technical details of each system are presented in Table 1. In 2000-2001, sowing conditions were extremely difficult for the direct drilling system due to large amounts of crop residue on the soil surface and high moisture levels in the seedbed (high rainfall during autumn). This led to poor emergence of the wheat crop in this system. We therefore took no samples from the direct drilling system in 2000-2001.

Evaluation of kernel infection by head blight

For each cropping system (and each subplot), we used four randomly selected samples of 25 heads to determine disease severity, every week between flowering and physiological maturity of the grain. Disease severity was estimated by determining the percentage of spikelets per head for which at least one third of the external surface of glumes presented symptoms of head blight.

Mycotoxin extraction and analysis

We determined the mycotoxin content of grains from plants sampled at harvest. Plant sampling was based on the 98/53/CE directive, which lays down the sampling procedure for official controls of aflatoxin level. We used this directive as there is no similar text dealing with Fusarium toxins. According to this directive, for plots with yields lower than 1 t, 10 samples of 100 g must be harvested and pooled to give a total sample of 1 kg. We collected 1 kg grain samples in this manner from each plot. Samples were dried (48h -80°C) and sent to Qualtech (Vandoeuvre-les-Nancy, France) for testing. In the first two years, we determined the levels of the following mycotoxins: zearalenone, group A and B trichothecenes

Statistical methods

We used the classic ANOVA (one-way analysis of variance) program of the SAS software suite to compare the systems. If a significant F ratio was obtained (α<0.05), we analysed differences between treatments, using Bonferroni's multiple range test. In this study, there were interactions between year and climate and between year and changes in the cropping system. We therefore chose not to analyse the year effect directly. We focused instead on cropping system effects. ANOVA was carried out with all the individual measurements from each replicate.

Results

Disease severity

The severity of Fusarium head blight differed considerably from year to year (Table 2). In 2000, severe Fusarium head blight was observed, whereas no symptoms were detected in 2001. Disease severity differed significantly between the four cropping systems (Table 2-a) in 2000. No difference between the cropping systems was observed in 2001, and in 2002 slight differences were observed between cropping systems but the range of severity observed was limited. In 2000, the direct drilling system displayed the highest level of contamination, followed by the conventional and organic systems. In 2002, the conventional system was the most heavily contaminated. These rankings, in particular that for 2002, which was less clearcut, were confirmed by analysis of variance: cropping system was found to have a significant effect on disease severity (α < 0.001).

The pattern of disease severity observed in the subsystems was similar to that for the systems as a whole (Table 2-b): contamination was greatest in 2000, no symptoms were observed in 2001, and 2002 was intermediate. Analysis of the results for 2000 showed that in two systems (conventional and integrated), the subsystem (conventional cultivar used without fungicide) crop was more heavily infected that the adapted system (appropriate cultivar for the system) crop. In contrast, the subsystem and adapted plots were similar for the direct drilling and organic systems. The direct drilling subsystem displayed the highest level of disease severity.

Thus, the factors responsible for this high level of disease were presumably not cultivar or fungicide strategy. In contrast, in the other three systems, disease severity seemed to depend on the cultivar and/or fungicide treatment used. In 2001 and 2002, the adapted cropping systems and the subsystems presented similar levels of disease. A plot effect was detected for the subsystem treatment, with plot I more heavily contaminated than plot II. This effect was significant (α ≤ 0.001) in 2000.

The accumulation of mycotoxins in grain

Very little deoxynivalenol (DON) was accumulated in the harvested grain (Table 3-a). The DON content was below the maximum levels (750 µg/kg) cited in provisional European recommendations, for most treatments. The highest levels of accumulation of this toxin were recorded in 2000. In addition to this year effect, a cropping system effect was observed: the direct drilling system (in 2000 and 2002) was systematically and significantly the most heavily contaminated with this toxin. In 2001, when no samples were available for the direct drilling system, the organic system showed significantly higher levels of DON contamination than the other systems.

The results obtained for the subsystems were very similar (Table 3-b), with the same year and cropping system effects identified. However, the effects of not applying fungicide differed according to the year. DON contamination of the conventional system plots was lower in the absence than in the presence of fungicide in 2000 and 2001, and higher in 2002. Conversely, the Charger cultivar systematically displayed higher levels of contamination than the cultivar used in the organic system. These factors had different effects on mycotoxin contamination.

The integrated and direct drilling systems (for which fungicide protection and cultivar varied simultaneously) displayed higher levels of mycotoxin contamination for the subsystems than for the adapted system, with the exception of the direct drilling system in 2000. However, the choice of cultivar and the decision not to treat had a relatively small effect in this system in this year because DON levels consistently exceeded 4000 µg/kg (which is much higher than the maximum level recommended). No significant plot effect was identified.

The fields were also contaminated with nivalenol and zearalenone. Nivalenol contamination differed significantly (α≤0.05) between cropping systems in 2001. In all the years studied, nivalenol was detected only in the organic and direct drilling systems. Zearalenone contamination differed considerably from year to year. Zearalenone levels were highest in 2000, with the direct drilling system displaying significantly higher levels of contamination than the other systems (α≤0.05) in that year. In 2001, only one field (the organic subsystem) had zearalenone contamination. Zearalenone was not detected in 2002. We detected no other mycotoxins.

The relationship between disease severity and deoxynivalenol accumulation

No relationship between disease severity and mycotoxin levels was observed, between or within years (Figure 1). For example, disease severity was scored at 0% despite DON levels of 0-600 µg/kg (years 2001 and 2002). Furthermore, DON levels were low (350 µg/kg) on plots in which disease severity score was higher (22%).

The lack of a relationship between Fusarium head blight severity and mycotoxin content was even more evident for the subsystems. No disease was observed despite DON levels of 0-recorded for plots with high levels of disease severity (63%). We also found that there was no relationship between disease severity and mycotoxin levels for nivalenol and zearalenone.

Discussion

Our results demonstrate that Fusarium head blight severity depends mainly on climatic effects. Similar results have been obtained in previous studies [START_REF] Andersen | The development of Gibberella zeae head blight of wheat[END_REF][START_REF] Langseth | The influence of cultivars, location and climate on deoxynivalenol contamination of Norwegian oats 1985-1990[END_REF]. However, in years with average or high levels of disease, cropping system affects disease severity. The direct drilling system resulted in the highest disease severity in both years for which results were available. In 2000, the direct drilling system wheat crop followed a maize crop whereas the wheat crop of the integrated system followed oilseed rape. In 2002, the preceding crop was pea in both cases. However, the major difference between the direct drilling and integrated systems was the absence of tillage in the direct drilling system: comparison of these two systems, in 2002, can therefore be used to determine the effect of ploughing in residues. Our results seem to support previous conclusions concerning the effects of land preparation [START_REF] Teich | Survey of Fusarium head blight and possible effects of cultural practices in wheat fields in Lambton County in 1983[END_REF][START_REF] Teich | Effect of cultural practices, soil phosphorus, potassium, and pH on the incidence of Fusarium head blight and deoxynivalenol levels in wheat[END_REF][START_REF] Dill-Macky | The effect of previous crop residues and tillage on Fusarium head blight of wheat[END_REF].

An analysis of agronomic traits (head density, yield, grain protein content, Septoria leaf spot severity, eyespot severity) revealed no clear reason for the plot effect observed for disease severity in 2000. Follow-up of scab severity between flowering and harvest provided a potential clue: disease severity rapidly increased between June 27 and July 3 in plot I, but not in plot II. This increase in disease severity was probably due to heavy rainfall (21 mm) on July 2nd, after 20 days without precipitation. The experimental field was surrounded by trees and plot II was more protected from rain than plot I, which may have resulted in an edge effect.

This work shows that the Fusarium mycotoxin (deoxynivalenol, nivalenol and zearalenone) content of wheat crops depends on both year and cropping system. These results are consistent with those of previous studies [START_REF] Teich | Effect of cultural practices, soil phosphorus, potassium, and pH on the incidence of Fusarium head blight and deoxynivalenol levels in wheat[END_REF][START_REF] Dill-Macky | The effect of previous crop residues and tillage on Fusarium head blight of wheat[END_REF]: the direct drilling system generally results in the highest levels of contamination, particularly in years with high levels of disease, such as 2000. Moreover, in 2000, the presence of maize immediately preceding the wheat crop in the rotation of the direct drilling system seemed to increase the production of mycotoxins [START_REF] Teich | Effect of cultural practices, soil phosphorus, potassium, and pH on the incidence of Fusarium head blight and deoxynivalenol levels in wheat[END_REF]; Dill-pattern emerged. In 2000 and 2002, the conventional system was more contaminated than the organic system, whereas the reverse was true for 2001. Higher levels of contamination in the conventional system may result from disease severity being high in that year (which was the case in 2000) or may be due to the fungicide used. Indeed in 2000, we used a mixture of strobularins whereas in 2001 and 2002, we used a mixture of strobularins and triazoles.

Strobularins belong to a family of fungicide molecules that are particularly effective against Microdochium strains but some Fusarium strains are limited effectively only by triazole molecules [START_REF] Maumené | Incidence des fongicides et de la variété sur le risque de contamination du blé tendre par le déoxynivalénol (DON) en france : premiers résultats[END_REF]. As Microdochium produces no toxins, the application of strobularins does not limit mycotoxin levels. Conversely, it might favour the development of pathogenic strains that produce toxins. Thus, the fungicide treatment used in 2000 might have favoured an increase in mycotoxin levels in the adapted conventional system whereas competition between Fusarium and Microdochium strains may have limited the production of mycotoxins in the organic systems and in the conventional subsystem. In 2000, fungicide application might have resulted in an increase in disease severity in the adapted conventional system, whereas saprophytes might have limited disease development in the adapted organic system. In contrast, in 2002, the absence of fungicide treatment resulted in higher levels of disease severity. In this case, the molecule used was probably more suitable for the infecting pathogenic strains.

The "contaminating complex" effect, operating under natural contamination conditions is probably responsible of the absence of a relationship between disease severity and DON levels. Our results, which are consistent with those of other authors [START_REF] Liu | Comparison of visual head blight ratings, seed infection levels, and deoxynivalenol production for assessment of resistance in cereals inoculated with Fusarium culmorum[END_REF], are especially important as they relate to moderate or low disease levels, which may conceal significant contamination by mycotoxins. Indeed, in contrast to studies in which crops were artificially contaminated with a single pathogenic strain of Fusarium culmorum, in which such a relationship was found [START_REF] Teich | Effect of cultural practices, soil phosphorus, potassium, and pH on the incidence of Fusarium head blight and deoxynivalenol levels in wheat[END_REF][START_REF] Snijders | Effects of head blight caused by Fusarium culmorum on toxin content and weight of wheat kernels[END_REF], several species of Fusarium and Microdochium are present in the field in natural contamination conditions. Hence, DON levels, disease severity and the relationship between these variables probably depend on the traits of the species present -capacity to produce disease, capacity to produce mycotoxins (deoxynivalenol, nivalenol, zearalenone), capacity to produce symptoms and mycotoxins-together with the time at which contamination occurs, the climatic conditions between the time of contamination and harvest, and the nutritive resources that the plant can provide to the pathogenic strains: in other words, the whole cropping system. Table 1a: Principal traits of the four adapted cropping systems over the three years. Mixture 1: Somme, Malacca, Renan, Camp Remy; Mixture 2: Somme, Malacca, Virtuose, Apache; Mixture 3: Charger, Virtuose, Renan, Paindor.

The organic system received no nitrogen fertilisation in 1999-2000 and 2000-2001. In the 2000-2001 study period, the wheat crop in the organic system was sown in January 2001. We planned to sow this crop in November, but incessant rain between October 2000 and January 2001 rendered this impossible.

Table 1b: Principal traits of the four cropping subsystems over the three years. No fungicide treatment, cultivar Charger.

The organic system received no nitrogen fertilisation in 1999-2000 and 2000-2001. In the 2000-2001 study period, the wheat crop in the organic system was sown in January 2001. We planned to sow this crop in November, but incessant rain between October 2000 and January 2001 rendered this impossible.

Table 2: Comparison of disease severity in the four cropping systems during three years.

Charger is the cultivar used in the adapted conventional system. Table 3: Comparison of deoxynivalenol accumulation (µg/kg) for four cropping systems over three years. nd: not detected (30 µg/kg is the detection limit).

(

  diacetoxyscirpenol (DAS), T-2, HT-2, nivalenol (NIV), deoxynivalenol (DON) and derivatives (3-aDON and 15-aDON), fusarenone (FUX), and fumonisins B1 and B2. In the final year, we determined levels of zearalenone, group A and B trichothecenes (diacetoxyscirpenol (DAS), monoacetoxyscirpenol (MAS), T-2, HT-2, neosolaniol, T2-triol, nivalenol (NIV), deoxynivalenol (DON) and derivatives (3-aDON and 15-aDON), and fusarenone (FUX). Zearalenone levels were determined with an immunoaffinity highperformance liquid chromatography (HPLC) method, validated by the French norm NF V03-110. Fumonisines were determined by HPLC, and trichothecenes by gas chromatographymass spectroscopy (GC-MS), validated by the French norm NF EN ISO/CEI 17025.

Figure 1 :

 1 Figure 1: Relationship between disease severity and deoxynivalenol accumulation in four cropping systems over three years. C conventional, I integrated, DD direct drilling, O organic. 00 2000, 01 2001, 02 2002 (♦) Adapted cropping system, (□) Subsystem: cultivar Charger without fungicide

Table 1

 1

	368	(b)
	Subsystem	Sowing time	Flowering time	Harvest time	Nitrogen (kg/ha)	Seed rate (seeds/m²)	Fungicide treatment	Cultivar	Previous crop
	Conventional 11/10/99 23/05/00		247	260			Oilseed
	Integrated				197	165			Rape
	Direct drilling		27/10/99 25/05/00	03/08/00	197	220	None	Charger	Maize with fescue
	Organic	10/11/99 30/05/00		0	450			Oilseed Rape
	Conventional 04/10/00 28/05/01		249	260			Oilseed
	Integrated	24/10/00 01/06/01	30/07/01	201	165	None	Charger	Rape
	Organic	17/01/01 27/06/01		0	450			Spring lupin
	Conventional 05/10/01 06/06/02 29/07/02	208	260
	Integrated	17/10/01 06/06/02		131	165
	Direct drilling		26/10/01 06/06/02	30/07/02	141	220	None	Charger	Pea
	Organic	31/10/01 06/06/02		54	450
	369

Acknowledgements

This research was funded by the DGAL and by INRA. We would like to thank P. Verger for his suggestion at the start of this study. We would also like to thank Mrs Tanis-Plant and Mrs Jerome for editorial advice and Mr Saulas, Mr Grandeau, Mrs Tanneau and Mr Mille for their contributions to the experimental work.


Dimitrios Mitsotakis
email: dimitrios.mitsotakis@vuw.ac.nz

Costas Synolakis
email: costas@usc.edu

Mark Mcguinness
email: mark.mcguinness@vuw.ac.nz

A MODIFIED GALERKIN / FINITE ELEMENT METHOD FOR THE NUMERICAL SOLUTION OF THE SERRE-GREEN-NAGHDI SYSTEM

Keywords: 2010 Mathematics Subject Classification. 76B15, 76B25, 65M08 Finite element methods, Solitary waves, Green-Naghdi system, Serre equations

A new modified Galerkin / Finite Element Method is proposed for the numerical solution of the fully nonlinear shallow water wave equations. The new numerical method allows the use of low-order Lagrange finite element spaces, despite the fact that the system contains third order spatial partial derivatives for the depth averaged velocity of the fluid. After studying the efficacy and the conservation properties of the new numerical method, we proceed with the validation of the new numerical model and boundary conditions by comparing the numerical solutions with laboratory experiments and with available theoretical asymptotic results.

Introduction

The motion of an ideal (inviscid, irrotational) fluid bounded above by a free surface and below by an impermeable bottom is governed by the full Euler equations of water wave theory, [START_REF] Whitham | Linear and nonlinear waves[END_REF]. Because of the complexity of the Euler equations, a number of simplified models describing inviscid fluid flow have been derived such as various Boussinesq type (BT) models. The Serre-Green-Naghdi (SGN) system can be considered to be a BT model that approximates the Euler equations, and models one-dimensional, two-way propagation of long waves, without any restrictive conditions on the wave height. The SGN system is a fully non-linear system of the form,

h t + (hu) x = 0 , (1a)
h + T h b u t + gh(h + b) x + huu x + Q h u + Q h b u = 0 , (1b)
where

h(x, t) . = η(x, t) -b(x) , (1c)
is the total depth of the water between the bottom b(x) and the free surface elevation η(x, t), u(x, t) is the depth averaged horizontal velocity of the fluid, and g the acceleration due to gravity. The operators T h b , Q h and Q h b depend on h, and are defined as follows:

T h b w = h h x b x + 1 2 hb xx + b 2 x w - 1 3 h 3 w x x , (1d)
Q h w = - 1 3 h 3 (ww xx -w 2 x ) x , (1e)
Q h b w = 1 2 h 2 (w 2 b xx + ww x b x ) x - 1 2 h 2 ww xx -w 2 x b x +hw 2 b x b xx +hww x b 2 x . (1f)
In dimensional and unscaled form, the independent variable, x ∈ R is a spatial variable and t ≥ 0 represents the time.

The SGN equations as derived by Seabra-Santos et.al. in [START_REF] Seabra-Santos | Numerical and experimental study of the transformation of a solitary wave over a shelf or isolated obstacle[END_REF] have also been derived in a three-dimensional form in [START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF] and in a different formulation by Green and Naghdi [START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF]. In the case of a flat bottom (i.e. b x = 0) (1) is simplified to the so-called Serre system of equations derived first by Serre [START_REF] Serre | Contribution á l' étude des écoulements permanents et variables dans les canaux[END_REF][START_REF] Serre | Contribution á l' étude des écoulements permanents et variables dans les canaux[END_REF] and re-derived later by Su and Gardner, [START_REF] Su | Korteweg-de Vries equation and generalizations. III. Derivation of the Korteweg-de Vries equation and Burgers equation[END_REF]. For these reasons the equations (1) are also known as the Serre, or Green-Nagdhi, or Su and Gardner equations. We will henceforth refer to them here as the Serre-Green-Naghdi (SGN) equations.

Under the additional assumption of small amplitude waves (i.e. the solutions are of small amplitude), the SGN system reduces to the Peregrine system, [START_REF] Peregrine | Long waves on a beach[END_REF]:

h t + (hu) x = 0 , (2a)
u t + g(h + b) x + uu x - b 2 [bu] xxt - b 2 6 u xxt = 0 . (2b)
Peregrine's system belongs to the weakly dispersive and weakly nonlinear BT systems. There are also other BT systems that are asymptotically equivalent to Peregrine's system, cf. [START_REF] Mitsotakis | Boussinesq systems in two space dimensions over a variable bottom for the generation and propagation of tsunami waves[END_REF][START_REF] Nwogu | Alternative form of Boussinesq equations for nearshore wave propagation[END_REF]. Differences between the SGN equations and BT models are explained in [START_REF] Dias | On the fully-nonlinear shallow-water generalized Serre equations[END_REF]. In the same work the inclusion of surface tension effects have been included and explained in detail.

Although equations ( 2) can be derived from the SGN system, their solutions have different properties. For example the solutions of the SGN equations are invariant under the Galilean boost, while the respective solutions of (2) are not. It is noted that the SGN equations have a Hamiltonian formulation, [START_REF] Li | Hamiltonian structure and linear stability of solitary waves of the Green-Naghdi equations[END_REF][START_REF] Israwi | Large time existence for 1D Green-Naghdi equations[END_REF]. Specifically, for a stationary bathymetry the SGN system conserves the total energy functional [START_REF] Israwi | Large time existence for 1D Green-Naghdi equations[END_REF]:

I(t) = R gη 2 + hu 2 + T h b u • u dx , (3)
in the sense that I(t) = I(0), for all t > 0. The conservation of this Hamiltonian will be used to measure accuracy and conservation properties of the proposed numerical methods.

Although both systems are known to admit solitary wave solutions propagating without change in their shape over a horizontal bottom y = -b 0 , only the solitary waves of the SGN system have known formulas in a closed form. Specifically, a solitary wave of the SGN system with amplitude A can be written in the form:

h s (x, t) = b 0 + Asech 2 [λ(x -c s t)], u s (x, t) = c 0 1 - b 0 h(x, t) , (4a)
λ = 3A 4b 2 0 (b 0 + A) and c s = c 0 1 + A b 0 , (4b)
where c s is the phase speed of the solitary wave and c 0 = √ gb 0 is the linear wave speed.

In addition to the above properties the system (1) is known to have several favourable well-posedness properties. For example, it is locally well-posed in time for smooth bottom functions in R. Specifically, if C ∞ b (I) is the space of bounded and continuously differentiable functions on the interval I and H s (I) denotes the usual Sobolev space of s-order weakly differentiable functions on I, and if b ∈ C ∞ b (R), κ > 1/2, s ≥ κ + 1 with the initial condition is (η 0 , u 0 ) ∈ H s (R) × H s+1 (R), then there is at least a maximal time T max > 0, such that the SGN equations admit a unique solution (η, u) ∈ H s (R) × H s+1 (R), cf. [START_REF] Israwi | Large time existence for 1D Green-Naghdi equations[END_REF]. We also mention that the solution will satisfy a non-vanishing depth condition when the initial data satisfy the same non-vanishing condition, i.e. there is an α > 0 such that for all t ≥ 0 η(x, t) -b(x) ≥ α > 0 .

(

) 5
For more information on the model equations, including the derivation, theory and justification, refer to [START_REF] Lannes | The water wave problem: Mathematical analysis and asymptotics[END_REF][START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF][START_REF] Barthelemy | Nonlinear shallow water theories for coastal waves[END_REF]. We consider here the SGN equations on a finite interval I = [a, b]. In addition to periodic boundary conditions, we will study physically relevant reflective boundary conditions. The linearized equations (1) about the trivial solution coincide with the linearized Peregrine system. Initial-boundary value problems for Peregrine's system have been studied in [START_REF] Fokas | Boundary value problems for Boussinesq type systems[END_REF], where the existence of solutions was proven for boundary data for u(a, t) and u(b, t). In a similar manner, the boundary conditions u(a, t) = u(b, t) = 0 for t ≥ 0 are sufficient for the system (1) to model the reflection of the waves on solid wall boundaries. These specific wall boundary conditions have been used widely to describe reflection of waves on the boundaries for various numerical models including BT systems such as the Nwogu system, [START_REF] Kazolea | A well-balanced shock-capturing hybrid finite volume-finite difference numerical scheme for extended 1D Boussinesq models[END_REF].

The numerical discretization of the SGN equations is a challenging problem due to their complicated form. Not only do the operators in front of the temporal derivatives depend on the unknown function h(x, t), but they also contain high order derivatives in nonlinear terms. Because of the presence of those high order spatial partial derivatives the solution to the system should be smooth. Recently several schemes have been proposed such as Finite Difference / Finite Volume Schemes [START_REF] Mirie | Collision between two solitary waves. Part 2. A numerical study[END_REF][START_REF] Seabra-Santos | Numerical and experimental study of the transformation of a solitary wave over a shelf or isolated obstacle[END_REF][START_REF] Antunes Do | Numerical solution of the generalized Serre equations with the MacCormack finite-difference scheme[END_REF][START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part I: Model development and analysis[END_REF][START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part II: Boundary conditions and model validation[END_REF][START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF] and Discontinuous Galerkin methods [START_REF] Li | High order well-balanced CDG-FE methods for shallow water waves by a Green-Naghdi model[END_REF][START_REF] Panda | Discontinuous Galerkin methods for solving Boussinesq-Green-Naghdi equations in resolving non-linear and dispersive surface water waves[END_REF][START_REF] Panda | Discontinuous Galerkin methods for solving Boussinesq-Green-Naghdi equations in resolving non-linear and dispersive surface water waves[END_REF]. Although these methods are very useful to study practical problems, such as the runup of nonlinear waves on slopes, they can be highly dissipative methods in the sense that they introduce numerical dissipation or dispersion, due to the approximation of the nonlinear terms by dissipative flux functions and the use of low order slope limiters. These methods appear to have good conservation properties though when high-order WENO methods are used or other non-classical numerical fluxes are used, [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF]. Some other highly accurate numerical methods have been developed for the Serre equations in the case of a horizontal bottom, such as spectral methods [START_REF] Dutykh | Finite volume and pseudo-spectral schemes for the fully nonlinear 1D Serre equations[END_REF] and Galerkin / Finite Element Method (FEM) [START_REF] Mitsotakis | On the Galerkin / finite-element method for the Serre equations[END_REF]. Although these methods appear to have satisfactory conservation properties, it is very difficult to extend them to the SGN equations especially in the case of two horizontal dimensions. For example the standard Galerkin method of [START_REF] Mitsotakis | On the Galerkin / finite-element method for the Serre equations[END_REF] requires tensor products of cubic splines in order to be consistent in two spatial dimensions in a similar manner to [START_REF] Dougalis | On some Boussinesq systems in two space dimensions: Theory and numerical analysis[END_REF].

In this study, a new modified Galerkin / Finite Element Method (FEM) is proposed for the numerical solution of the SGN equations. This method allows the use of low-order finite elements such as piecewise quadratic (P 2 ) or even piecewise linear (P 1 ) Lagrange finite elements. Two of the main advantages of this method is that it is highly accurate, and it has very good conservation properties. Other advantages are the sparsity of the resulting linear systems, the low complexity of the algorithm due to the use of low-order finite element spaces and finally its potential to be extended to the two-dimensional model equations of [START_REF] Lannes | The water wave problem: Mathematical analysis and asymptotics[END_REF].

Similar techniques that reduce the requirement of high-order finite elements (for example the use of cubic splines) have been used previously for weakly-nonlinear Boussinesq systems in [START_REF] Walkley | A finite element method for the one-dimensional extended Boussinesq equations[END_REF][START_REF] Walkley | A finite element method for the two-dimensional extended Boussinesq equations[END_REF][START_REF] Dougalis | Initial-boundary-value problems for Boussinesq systems of Bona-Smith type on a plain domain: theory and numerical analysis[END_REF] where the second derivative in the linear dispersive terms has been replaced by either the discrete Laplacian operator or the solution of an intermediate problem. In the case of the Bona-Smith type of Boussinesq systems with wall boundary conditions, the modified Galerkin method converges at an optimal rate showing great performance contrary to the suboptimal convergence rates achieved with the standard Galerkin / FEM method, [START_REF] Dougalis | Initial-boundary-value problems for Boussinesq systems of Bona-Smith type on a plain domain: theory and numerical analysis[END_REF]. Similarly to the behavior of the modified Galerkin method for the weakly-nonlinear Boussinesq equations, the proposed FEM scheme for the SGN equations can achieve optimal convergence rates depending on the choice of the trial function spaces, contrary to the suboptimal rates obtained for the standard Galerkin method for the SGN system and also for Peregrine's system as shown in [START_REF] Antonopoulos | Numerical solution of the 'classical' Boussinesq system[END_REF]. Specifically, in order to achieve optimal convergence properties one may use spaces of piecewise linear elements for the free surface elevation and piecewise quadratic elements for the horizontal velocity.

The validated numerical method is applied to study the SGN equations with reflective boundary conditions in a systematic way through a series of numerical experiments. In particular, we focus on the following issues:

• accuracy of the modified Galerkin method and invariant conservation; • reflection of solitary waves at a vertical wall; and • shoaling of solitary waves on plain or composite beaches.

The convergence properties of the new numerical method are also tested in the case of periodic boundary conditions. For more information about the behavior and the properties of the Galerkin / FEM method with periodic boundary conditions we refer to [START_REF] Mitsotakis | On the Galerkin / finite-element method for the Serre equations[END_REF].

This paper is organized as follows. Section 2 presents the fully discrete schemes for the SGN equations. In Section 3 we study the convergence, accuracy and numerical stability of the modified Galerkin method. Finally, Section 4 presents computational studies of shoaling and reflected waves validating both the choice of boundary conditions and the numerical scheme. We close the paper with conclusions in Section 5.

The numerical methods

We consider the initial-boundary value problem (IBVP) comprising System (1), subject to reflective boundary conditions:

h t + (hu) x = 0 , (6a)
h + T h b u t + gh(h + b) x + huu x + Q h u + Q h b u = 0 , (6b) u(a, t) = u(b, t) = 0, h(x, 0) = h 0 (x), u(x, 0) = u 0 (x) , (6c)
where again

T h b w = h h x b x + 1 2 hb xx + b 2 x w - 1 3 h 3 w x x , (6d)
Q h w = - 1 3 h 3 (ww xx -w 2 x ) x , (6e)
Q h b w = 1 2 h 2 (w 2 b xx + ww x b x ) x - 1 2 h 2 ww xx -w 2 x b x +hw 2 b x b xx +hww x b 2 x , (6f)
x ∈ (a, b) ⊂ R and t ∈ [0, T ]. We assume that (6) possesses a unique solution, such that h and u are sufficiently smooth and, for any t ∈ [0, T ], in suitable Sobolev spaces h(x, •) ∈ H s , u(x, •) ∈ H s+1 0 , where s ≥ 1. Here and below, • s denotes the standard norm in H s while H s 0 will denote the subspace of H s whose elements vanish at x = a and x = b. We also use the inner product in L 2 ≡ H 0 , denoted by (•, •), which is

(u, v) = b a uv dx , for all u, v ∈ L 2 .
A spatial grid of the interval [a, b] is a collection of points x i = a + i ∆x, for i = 0, 1, • • • , N , where ∆x is the grid size, and N ∈ N, such that ∆x = (b -a)/N . Let ( h, ũ) ∈ S h × S u be the corresponding spatially discretized solutions of the Galerkin / Finite Element Method (FEM) for suitable finite-dimensional spaces S h and S u . First we present the standard Galerkin / FEM semidiscretization.

2.1. The standard Galerkin method. For the standard Galerkin method, we consider the space of smooth splines

S r = ϕ ∈ C r-1 [a, b] ϕ| [xi,xi+1] ∈ P r , 0 ≤ i ≤ N -1 ,
where P r is the space of polynomials of degree r. The standard Galerkin method requires r ≥ 3. Here, we take r = 3. The trial function space for the first equation and for the solution h is chosen as S h = S r , while for the second equation and for the approximation of the depth averaged velocity of the fluid ũ is

S u = S r ∩ { ϕ ∈ C[a, b]| ϕ(a) = ϕ(b) = 0}.
To state the associated semi-discrete problem, let φ ∈ S h and ψ ∈ S u be arbitrary test functions. After taking inner products, and using integration by parts, the semi-discrete problem takes the form:

( ht , φ) + ( hũ) x , φ = 0 , (7a)
B(ũ t , ψ; h) + h g( h + b) x + ũũ x , ψ + Q(ũ, ψ; h) + Q b (ũ, ψ; h) = 0 , (7b)
where B, Q and Q b are defined for any ω, ψ ∈ S u as

B(ω, ψ; h) . = h 1 + hx bx + 1 2 hb xx + b2 x ω, ψ + 1 3 h3 w x , ψ x , (7c)
Q(ω, ψ; h) . = 1 3 h3 ωω xx -ω 2 x , ψ x , (7d)
Q b (ω, ψ; h) . = - 1 2 h2 (ω 2 bxx + ωω x bx ), ψ x - 1 2 hb x h ωω xx -ω 2 x -ω 2 bxx -ωω x bx , ψ . (7e)
This is a system of ordinary differential equations (ODEs). Given the initial conditions h(x, 0) . = h0 = P h (h 0 ) and ũ(x, 0) . = ũ0 = P u (u 0 ) where P h and P u are appropriate projections on S h and S u respectively, we assume that the system [START_REF] Antunes | On breaking waves and wave-current interaction in shallow water: A 2DH finite element model[END_REF] has a unique solution. Appropriate projections of the initial conditions could be the standard L 2 -projections on S h and S u , defined as H 0 ∈ S h and U 0 ∈ S u such that The presence of the term ω xx in Q and Q b requires the use of at least C 2 smooth splines. In particular, one may use cubic splines, which correspond to S h ≡ S 3 , i.e., r = 3. A basis for the space S h for a uniform grid of mesh-length ∆x can be formed by the functions

φ j (x) = B(x -x j /∆x) [a,b] , j = -1, 0, • • • , N + 1 where x -1 = a -∆x, x N +1 = b + ∆x and B(x) =            1 4 (x + 2) 3 , -2 ≤ x ≤ -1 , 1 4 [1 + 3(x + 1) + 3(x + 1) 2 -3(x + 1) 3 ], -1 ≤ x ≤ 0 , 1 4 [1 + 3(1 -x) + 3(1 -x) 2 -3(1 -x) 3 ], 0 ≤ x ≤ 1 , 1 4 (2 -x) 3 , 1 ≤ x ≤ 2 , 0, x ∈ R -[-2, 2] .
The basis of S u = S 3 0 can be described by the functions ψ j (x) = φ j (x) for 2 ≤ j ≤ N -2, plus four functions ψ 0 , ψ 1 , ψ N -1 , ψ N , taken as linear combinations of the φ -1 , φ 0 , φ 1 and φ [START_REF] Schultz | Spline Analysis[END_REF]. The convergence properties of the standard Galerkin method with cubic splines are very similar to those of the modified Galerkin method with S 3 elements and thus are not presented here. Some details and numerical experiments with the standard Galerkin method can be found in [START_REF] Kalisch | Mechanical balance laws for fully nonlinear and weakly dispersive water waves[END_REF].

N -1 , φ N , φ N +1 , which are such that ψ 0 (a) = ψ 1 (a) = ψ N -1 (b) = ψ N (b) = 0. For example, we take ψ 0 = φ 0 -4φ -1 and ψ 1 = φ 1 -φ -1 ,
2.2. The modified Galerkin method. In real-world applications, the use of loworder finite element methods can lead to faster computations. We derive a numerical method that does not require high-order finite element spaces. This can be done by using similar techniques to those proposed in [START_REF] Walkley | A finite element method for the one-dimensional extended Boussinesq equations[END_REF][START_REF] Walkley | A finite element method for the two-dimensional extended Boussinesq equations[END_REF][START_REF] Dougalis | Initial-boundary-value problems for Boussinesq systems of Bona-Smith type on a plain domain: theory and numerical analysis[END_REF], but applied to the nonlinear term uu xx . The resulting method is a modified Galerkin method that allows the use of Lagrange elements as low as P 1 in order, consisting of piecewise linear functions. We proceed with the derivation of the modified Galerkin method but first we introduce the notation for the Lagrange finite element spaces that are subspaces of H 1 . It is noted that we cannot use Lagrange finite elements with the standard Galerkin method due to the second order derivative.

The space P r of Lagrange finite elements is defined on the grid

x i = x 0 + i∆x, i = 0, 1, • • • , N as: P r = χ ∈ C 0 [a, b] | χ| [xi,xi+1] ∈ P r , 0 ≤ i ≤ N -1 , (8)
We will restrict the analysis below to P 1 , P 2 and P 3 Lagrange finite element spaces.

The P 1 Lagrange elements can be defined using the basis functions

χ i (x) =    x-xi-1 ∆x , if x ∈ [x i-1 , x i ] , xi+1-x ∆x , if x ∈ [x i , x i+1 ] , 0 , otherwise .
The P 2 Lagrange finite element space can be defined on the same grid x i = x 0 +i∆x for i = 0, 1, • • • , N and at the midpoints x i+1/2 = x i + ∆x/2 for i = 0, 1, • • • , N -1 using the basis functions

χ i (x) = ϕ x -x i ∆x , i = 0, 1, • • • , N , and
χ i+1/2 (x) = ψ x -x i+1/2 ∆x , i = 0, 1, • • • , N -1 , with ϕ(x) =    (1 + x)(1 + 2x), -1 ≤ x ≤ 0 , (1 -x)(1 -2x), 0 ≤ x ≤ 1 , 0 otherwise , and
ψ(x) = 1 -4x 2 , |x| ≤ 1/2 , 0, otherwise .
Then a function w ∈ P 2 is written as

w(x) = N i=0 w(x i )χ i (x) + N -1 i=0 w(x i+1/2 )χ i+1/2 (x) .
For the construction of the Lagrange basis function of the general space P r we refer to [START_REF] Ern | Theory and Practice of Finite Elements[END_REF]. We also consider the Lagrange finite element spaces P r 0 = {χ ∈ P r |χ(a) = χ(b) = 0}. These spaces are subspaces of H 1 0 and will be used to approximate the depth averaged horizontal velocity of the water. For example, we will use the spaces S h = P r and S u = P q 0 , for some integers r, q or we will use the spaces of smooth splines described in the previous section. For more information related to Lagrange finite element spaces and its approximation properties, we refer to [START_REF] Ern | Theory and Practice of Finite Elements[END_REF]. It is worth mentioning that, for the discretization of the bottom boundary, we use the L 2 projection of the bathymetry and its derivatives.

Given v ∈ H 1 0 , we define the non-linear discrete Laplacian operator ∂2 :

H 1 0 → S u such that ∂2 v, ψ = -v 2 x , ψ -(vv x , ψ x ) , for all ψ ∈ S u . (9)
In fact, the function ∂2 v ∈ S u approximates the function vv xx as if v was a smooth C 2 function and substitution in (7d) and (7e) leads to the modified Galerkin semidiscretization:

( ht , φ) + ( hũ) x , φ = 0, φ ∈ S h , (10a)
B(ũ t , ψ; h) + h g( h + b) x + ũũ x , ψ + Q(ũ, ψ; h) + Qb (ũ, ψ; h) = 0, ψ ∈ S u , (10b)
where B, Q and Qb are defined for any ω, ψ ∈ S u as

B(ω, ψ; h) . = h 1 + hx bx + 1 2 hb xx + b2 x ω, ψ + 1 3 h3 w x , ψ x , (10c)
Q(ω, ψ; h) . = 1 3 h3 ∂2 ω -ω 2 x , ψ x , (10d)
Qb (ω, ψ; h) . = - 1 2 h2 (ω 2 bxx + ωω x bx ), ψ x - 1 2 hb x h ∂2 ω -ω 2 x -ω 2 bxx -ωω x bx , ψ . (10e)
Although we assume that the bottom function has the appropriate smoothness, the absence of second order spatial derivatives of the depth integrated horizontal velocity in the semidiscrete scheme allows the use of Lagrange elements, as well as high-order elements such as cubic or quintic splines. In what follows, we test the efficiency of the modified Galerkin method using the spaces P 1 , P 2 and P 3 of Lagrange elements and the space S 3 of cubic splines with periodic and reflective boundary conditions.

Remark 2.1 (Mass lumping). In order to compute the non-linear discrete Laplacian ∂w one needs to solve the linear system obtained by the discretization of ( 9).

In the case of wall boundary conditions if, for example, ψ i denotes basis functions of S u the system can be written as Mw = f where the mass matrix is a banded matrix with entries M ij = (ψ i , ψ j ) and

f i = -(v 2 x , ψ j ) -(vv x , ψ ′ j ).
To improve the speed of the numerical method, one may apply the method of mass lumping in the formulation of the matrix M. This can be done, for example in the case of quadratic Lagrange elements, by approximating the integrals with Simpson's rule. This leads to a diagonal matrix that can be inverted trivially. All the numerical experiments with P 1 and P 2 elements have been performed with mass lumping, in addition to the standard matrix formulation with comparable results.

Remark 2.2. The choice of the discrete Laplacian is not unique. For example, when periodic boundary conditions are used, then one may consider the linear discrete Laplacian, which replaces the term u xx as proposed in [START_REF] Walkley | A finite element method for the one-dimensional extended Boussinesq equations[END_REF][START_REF] Walkley | A finite element method for the two-dimensional extended Boussinesq equations[END_REF][START_REF] Dougalis | Initial-boundary-value problems for Boussinesq systems of Bona-Smith type on a plain domain: theory and numerical analysis[END_REF].

Remark 2.3. The presence of the term b xx in the semidiscrete scheme implies the typical requirement of a smooth bottom. When a piecewise linear bottom topography is given, then the use of an appropriate projection onto the finite element space (such as the elliptic projection) or local smoothing of the bottom is required. In this paper the bathymetry is usually a piecewise linear function and thus we use the smoothing method described in [START_REF] Alenitsyn | On smoothing of non-smooth functions[END_REF].

Remark 2.4. The modified Galerkin method can be expressed in an equivalent way by introducing an additional independent variable w satisfying the equation w = uu xx , [START_REF] Walkley | A finite element method for the one-dimensional extended Boussinesq equations[END_REF]. This increase of the degrees of freedom in the model equations by one allows the approximation of the unknown function u by low-order finite element spaces such as P 1 or P 2 finite element spaces.

Temporal discretization.

In [START_REF] Mitsotakis | On the Galerkin / finite-element method for the Serre equations[END_REF], it was shown by numerical means that the standard Galerkin method for the discretization of the SGN equations with flat bottom leads to a system of ODEs that is not stiff. Also, the classical, explicit, fourstage, fourth order Runge-Kutta (RK) method, described by the following Butcher tableau:

A b τ = 0 0 0 0 1/6 1/2 0 0 0 1/3 0 1/2 0 0 1/3 0 0 1/2 0 1/6 0 1/2 1/2 1 , (11)
is able to integrate the respective semi-discrete system numerically in time without imposing restrictive stability conditions on the ratio ∆t/∆x, but only mild restrictions on the mesh length such as ∆t/∆x ≤ 2 for smooth solutions.

Concerning the time integration, the modified Galerkin method has very similar behavior to the standard Galerkin method. Upon choosing appropriate basis functions for the spaces S h and S u , the semidiscrete system (10) represents a system of ODEs. We use a uniform time-step ∆t such that ∆t = T /K for K ∈ N. The temporal grid is then t n = n∆t, where n = 0, 1, • • • , K. Given the ODE y ′ = Φ(t, y), one step of this four-stage RK scheme (with y n approximating y(t n )) is:

for i = 1 → 4 do ỹi = y n + ∆t i-1 j=1 a ij y n,j y n,i = Φ(t n,i , ỹi ), evaluated at t n,i ≡ t n + τ i ∆t end for y n+1 = y n + ∆t 4 j=1 b j y n,j ,
where a ij , τ i , b i are given in table [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF]. Applying this scheme to [START_REF] Beji | A formal derivation and numerical modelling of the improved Boussinesq equations for varying depth[END_REF] and denoting by H n and U n the fully discrete approximations in S h and S u of h(•, t n ), u(•, t n ), respectively, leads to Algorithm 1. Given bases {φ i } of S h and {ψ i } of S u , the Algorithm 1 Time-marching FEM scheme for the IBVP of the system (10)

H 0 = P{h 0 } U 0 = P{u 0 } for n = 0 → N -1 do for i = 1 → 4 do Hi = H n + ∆t i-1 j=1 a ij H n,j Ũ i = U n + ∆t i-1 j=1 a ij U n,j ∂2 U i , φ = -((U i x ) 2 , φ) -(U i U i x , φ x ) (H n,i , ψ) = -(( Hi Ũ i ) x , ψ), evaluated at t n,i ≡ t n + τ i ∆t B(U n,i , φ; Hi ) = -Hi g( Hi -b) x + Ũ i Ũ i x , φ -P( Ũ i , φ; Hi , b) = 0 where P( Ũ i , φ; Hi , b) = Q( Ũ i , φ; Hi ) + Qb ( Ũ i , φ; Hi )
end for (c) Four linear systems with the time-independent matrix (ψ i , ψ j ) . The four linear systems in (a) arise from the discretization of the equation (10a), while the four linear systems in (c) arise from the computation of the discrete Laplacian [START_REF] Barthelemy | Nonlinear shallow water theories for coastal waves[END_REF]. All of these matrices are banded and symmetric, and one can use either direct methods for banded systems or classical iterative methods for sparse systems. To approximate the inner products, we use Gauss-Legendre quadrature with 3 nodes per ∆x for P 1 elements, 5 nodes for P 2 and P 3 , and 8 nodes for S 3 elements. It is noted that most of the experiments with P 1 , P 2 and mixed elements have also been tested using adaptive time-stepping methods such as the Runge-Kutta-Fehlberg method, [START_REF] Hairer | Solving ordinary differential equations: Nonstiff problems[END_REF], to ensure that the errors introduced by the temporal integration are negligible.

H n+1 = H n + ∆t 4 j=1 b j H n,j U n+1 = U n + ∆t

Accuracy and convergence

As it was pointed out in the introduction, the 'classical' Boussinesq (cB) system is similar in structure to the SGN system, given that both systems admit the same number of boundary conditions, and the equation for the free surface is the same. We thus expect similar behavior for the convergence of the numerical method. In [START_REF] Antonopoulos | Galerkin approximations of periodic solutions of Boussinesq systems[END_REF], it was shown that the convergence of the Galerkin / FEM method for the cB system with periodic boundary conditions is optimal, while in [START_REF] Antonopoulos | Error estimates for Galerkin approximations of the "classical" Boussinesq system[END_REF][START_REF] Antonopoulos | Numerical solution of the 'classical' Boussinesq system[END_REF] it was shown that the convergence of the Galerkin / FEM method with cubic splines for the cB system is suboptimal, and that in order to achieve an optimal rate of convergence, a non-standard method should be used. Specifically, it was shown that for wall boundary conditions in the case of P 1 elements, the following error estimates hold:

max 0≤t≤T h -h ≤ C ∆x 3/2 , max 0≤t≤T u -ũ ≤ C ∆x 2 , (12)
which are suboptimal for h and optimal for u. In the case of cubic splines, the following error estimates hold:

max 0≤≤T h -h ≤ C ∆x 3.5 ln 1 ∆x , max 0≤t≤T u -ũ ≤ C ∆x 4 ln 1 ∆x , (13)
which is suboptimal in both h and u but the factor ln ∆x is not dominant and generally it has an effect on the accuracy of the method that is too small to observe. In [START_REF] Kalisch | Mechanical balance laws for fully nonlinear and weakly dispersive water waves[END_REF], it is shown, by numerical means that the standard Galerkin method with cubic splines for the SGN system with wall boundary conditions appears to have similar convergence properties. The standard Galerkin method for the initial-periodic boundary value problem for the SGN system has been studied in [START_REF] Mitsotakis | On the Galerkin / finite-element method for the Serre equations[END_REF], where it was shown that the convergence is optimal. It is worth mentioning that the approximation ( 9) can be used also with periodic boundary conditions, and the resulting modified Galerkin method has optimal convergence rates including all the advantages of the modification. Specifically, the errors in the L 2 norm for the periodic case are of O(∆x 2 ) for the P 1 elements for both h and u, O(∆x 3 ) for the P 2 elements and for the S 3 elements is of O(∆x 4 ). Because the results for the periodic problem are very similar to those of [START_REF] Mitsotakis | On the Galerkin / finite-element method for the Serre equations[END_REF] and to the modified Galerkin method for the SGN system with periodic boundary conditions, we don't present them here.

We continue with the wall boundary conditions. In order to study the convergence of the numerical method, we use the SGN equations in nondimensional but unscaled form. We start with the IBVP [START_REF] Antunes | Extended Serre equation for applications in intermediate water depths[END_REF]. Because we don't know any analytical solutions for this problem, we consider the non-homogenous problem, which has the exact solution h(x, t) = 1 + e 2t (cos(πx) + x + 2) and u(x, t) = e -tx x sin(πx) satisfying the equations (6a) and (6b) with the appropriate right-hand sides. The solution is computed using the modified Galerkin method in the interval [0, 1] and for t ∈ (0, T ] with T = 1. The finite element spaces of our choice are P 1 and P 2 Lagrange finite element spaces and the space S 3 of cubic splines. It is noted that we also tried other "exact" solutions with different right-hand sides in (6a) and (6b) to verify the computed rates of convergence. The results have been all very similar, and thus we only present the results of one case.

In order to study the accuracy and the convergence of the modified Galerkin method, several error indicators have been computed. The computed errors are normalized and defined as

E s [F ] . = F (x, T ; ∆x) -F exact (x, T ) s F exact (x, T ) s , (14)
where F = F (•; ∆x) is the computed solution, i.e., either H ≈ h(x, T ) or U ≈ u(x, T ), F exact is the corresponding exact solution and s = 0, 1, 2, ∞ correspond to the L 2 , H 1 , H 2 and L ∞ norms, respectively. The analogous rates of convergence are defined as

rate for E s [F ] . = ln(E s [F (•; ∆x k-1 )]/E s [F (•; ∆x k )]) ln(∆x k-1 /∆x k ) , (15)
where ∆x k is the grid size listed in row k in each table. To ensure that the errors incurred by the temporal integration do not affect the rates of convergence we use ∆t ≪ ∆x while we take ∆x = 1/N . Tables 1, 2 and 3 present the errors and the corresponding rates of convergence of the modified Galerkin method with P 1 finite elements. It is shown that the rate of convergence is suboptimal for the total depth h and optimal for velocity u. Specifically, Table 1 suggests that h -h = O(∆x 3/2 ) and u -ũ = O(∆x 2 ). Table 2 suggests that h-h 1 = O(∆x 1/2 ) and u -ũ 1 = O(∆x). Finally, the L ∞ estimates shown in 3 are worse than the L 2 estimates, suggesting h-h ∞ = O(∆x) and u -ũ ∞ = O(∆x 2 ). All these results are very similar to those obtained theoretically and numerically in the case of Peregrine's system [START_REF] Antonopoulos | Error estimates for Galerkin approximations of the "classical" Boussinesq system[END_REF]. 1. Spatial errors and rates of convergence for the exact solution with P 1 finite elements using the L 2 norm.

N E 0 [H] rate for E 0 [H] E 0 [U ] rate for E 0 [U ] 10 1.4661 × 10 -2 - 2.
In the case of quadratic P 2 elements, the results with the same initial conditions are similar, but the convergence rates for h are different to the expected suboptimal rates, analogous to the P 1 case. The respective convergence rates for u are optimal as expected. Tables 4, 5 and 6 present the errors and the convergence rates for the 1.9769 Table 3. Spatial errors and rates of convergence for the exact solution with P 1 finite elements using the L ∞ norm.

N E 1 [H] rate for E 1 [H] E 1 [U ] rate for E 1 [U ] 10 2.8836 × 10 -1 - 1.6945 × 10 -1 - 20
modified Galerkin method with P 2 elements. Table 4 suggests that the h -h = O(∆x 2 ) and u -ũ = O(∆x 3 ). Table 5 suggests that h -h 1 = O(∆x) and u -ũ = O(∆x 2 ). Finally, Table 6 suggests the same estimates in the L ∞ norm with L 2 norm, i.e. h -h ∞ = O(∆x 2 ) and u -ũ ∞ = O(∆x 3 ). The same behaviour was observed when we used P 3 elements, i.e. the convergence rate for h was suboptimal h -h s = O(∆x 3-s ) while the convergence for u was optimal u -ũ s = O(∆x 4-s ), for s = 0, 1. 2.0021 1.9773 × 10 -9 3.0051 Table 4. Spatial errors and rates of convergence for the exact solution with P 2 finite elements using the L 2 norm.

N E 0 [H] rate for E 0 [H] E 0 [U ] rate for E 0 [U ] 10
We also studied the errors and convergence rates of a mixed modified Galerkin method, where we considered S h = P 1 and S u = P 2 i.e. we used linear elements for the approximation of h and quadratic elements for the approximation of u. This was the only case where optimal rates of convergence for both h and u were obtained. Tables 7-9 suggest that h -h s = O(∆x 2-s ) and u -ũ s = O(∆x 2.0000 2.5810 × 10 -9 2.9994 Table 6. Spatial errors and rates of convergence for the exact solution with P 2 finite elements using the L ∞ norm.

for s = 0, 1 while in the maximum norm we found that h -h ∞ = O(∆x 2 ) and u -ũ ∞ = O(∆x 3 ). Using mixed elements of higher order does not give analogously optimal results. For example, when we took S h = P 2 and S u = P 3 the convergence was suboptimal for h and optimal for u. More specifically, the errors were h -h s = O(∆x 2-s ) and u -ũ s = O(∆x 4-s ) for s = 0, 1 while in the maximum norm we found that h -h ∞ = O(∆x 2 ) and u -ũ ∞ = O(∆x 4 ). Observe that the convergence for h achieved with the mixed P 1 -P 2 elements was the same with the mixed P 2 -P 3 elements. 2.0014 1.9751 × 10 -9 3.0223 Table 7. Spatial errors and rates of convergence for the exact solution with S h = P 1 and S u = P 2 using the L 2 norm. Finally, Tables 10, 11 and 12 present the errors and convergence rates for the modified Galerkin method with cubic splines, i.e. with S 3 elements. In this case, the results are similar to those of the standard Galerkin method for the SGN and Peregrine's system. Specifically, Table 10 8. Spatial errors and rates of convergence for the exact solution with S h = P 1 and S u = P 2 using the H 1 norm. 2.0014 1.9751 × 10 -9 3.0223 Table 9. Spatial errors and rates of convergence for the exact solution with S h = P 1 and S u = P 2 using the L ∞ norm.

N E 0 [H] rate for E 0 [H] E 0 [U ] rate for E 0 [U ] 10
N E ∞ [H] rate for E ∞ [H] E ∞ [U ] rate for E ∞ [U ] 10
3.5 which is in a good agreement with the estimate h -h = O(∆x 3.5 ln(1/∆x)) suggested for the standard Galerkin method for Peregrine's system [START_REF] Antonopoulos | Error estimates for Galerkin approximations of the "classical" Boussinesq system[END_REF]. Also Table 10 suggests a rate of convergence for h close to 4, which agrees with the estimate u -ũ = O(∆x 4 ln(1/∆x)) of [START_REF] Antonopoulos | Error estimates for Galerkin approximations of the "classical" Boussinesq system[END_REF]. 3.4679 0.7968 × 10 -12 3.9708 Table 10. Spatial errors and rates of convergence for the exact solution with S 3 finite elements using the L 2 norm.

N E 0 [H] rate for E 0 [H] E 0 [U ]
To study further the accuracy of the method we compute the conservation of the energy functional I(t). As initial condition we use a solitary wave of amplitude a = 0.2 with a smooth bottom given by the function b(x) = 1 + 0.1 sin(πx/2). Because the SGN system together with the wall boundary conditions does not conserve the energy functional, we consider the interval [-100, 100] and maximum time T = 50 where the solution remained practically zero at the endpoints of the domain. We also used a fairly coarse grid with ∆x = 0.1. Due to the dissipative properties of the explicit Runge-Kutta method, conservation of the invariant 12. Spatial errors and rates of convergence for the exact solution with S 3 finite elements using the H 2 norm. 13. Spatial errors and rates of convergence for the exact solution with S 3 finite elements using the L ∞ norm. sensitive to the size of the timestep ∆t. For example, when we used ∆x = ∆t = 0.1 with elements in S 3 , the invariant remained at the value I(t) = 0.314542, when we used ∆t = 0.01, then the invariant remained at the value I(t) = 0.31454249795 conserving the digits shown.

I(t) is N E 1 [H] rate for E 1 [H] E 1 [U ] rate
N E ∞ [H] rate for E ∞ [H] E ∞ [U ] rate for E ∞ [U ] 200 0.
Although the energy functional I(t) depends on the choice of ∆t, the value of the Hamiltonian remains practically constant for any small value of ∆t < 0.01 and the error |I(t) -I(0)| is of O(10 -12 ). Figure 1 shows the sensitivity of the invariant I(t) to ∆t for fixed ∆x = 0.1 for S 3 finite elements. Figure 2 shows the error of the energy functional as a function of ∆x for different finite element spaces. We observe that the convergence of the computed energy functional is connected to the convergence of the numerical method. Although the explicit Runge-Kutta methods we used are not conservative (in the sense that they introduce a small amount

I (t ) -I (0) ∞ ) P 1 P 1 -P 2 P 2 P 2 -P 3 P 3 Figure 2.
Sensitivity of energy conservation to the value of ∆x for small values of ∆t of numerical dissipation), and their contribution to the error in the conservation of energy is larger than the error coming from the spatial discretization of the same order, it seems that the actual error is not important and upon choosing appropriate small values of ∆t, it can be considered negligible. When we used P 1 and P 2 elements, then the error from the spatial discretisation is larger than the errors embedded by the temporal discretisation that are considered unimportant.

Numerical experiments

In this section we present a series of numerical experiments that serve as benchmarks to verify the accuracy of the modified Galerkin method, in simulations with variable bathymetry and wall boundary conditions. We tested all the numerical methods for all the numerical experiments and we report all the important differences in the numerical results. All the graphs contain the results obtained with P 1 elements.

4.1. Shoaling of solitary waves. Shoaling of solitary waves and the nonlinear mechanisms behind the shoaling of a solitary wave have been studied theoretically and experimentally by various authors, [START_REF] Grilli | Shoaling of solitary waves on plane beaches[END_REF][START_REF] Grilli | Breaking criterion and characteristics for solitary waves on slopes[END_REF][START_REF] Synolakis | Green's law and the evolution of solitary waves[END_REF][START_REF] Synolakis | Evolution of maximum amplitude of solitary waves on plane beaches[END_REF], and shoaling is closely related to the runup of solitary wave on a plain beach, [START_REF] Synolakis | The runup of solitary waves[END_REF]. In [START_REF] Synolakis | Evolution of maximum amplitude of solitary waves on plane beaches[END_REF] after reviewing the derivation of Green's law for the amplitude evolution of shoaling waves for BT systems, that is,

η max /A ∼ b(x) -1/4 , (16)
several experiments with shoaling of solitary waves are presented. In [START_REF] Synolakis | Green's law and the evolution of solitary waves[END_REF] Green's law η max /A ∼ b(x) -1 for the Shallow Water Wave equations is derived. In [START_REF] Synolakis | Evolution of maximum amplitude of solitary waves on plane beaches[END_REF][START_REF] Skjelbreia | Observation of breaking waves on beaches by use of an LDV[END_REF] four regions of shoaling were determined, where the general Green's law η max ∼ b(x) -α applies for different values of the parameter α: The zone of gradual shoaling (α < 1), the zone of rapid shoaling α ≥ 1, the zone of rapid decay α < -1 and the zone of gradual decay α ≥ -1. Green's law for Boussinesq systems is valid in the zone of gradual and rapid shoaling, while Green's law for the shallow water system derived in [START_REF] Synolakis | Green's law and the evolution of solitary waves[END_REF] applies in the zone of rapid shoaling. In this work, we consider first the shoaling of solitary waves on a plane beach of (mild) slope 1 : 35. This experiment has been proposed by Grilli et. al. [START_REF] Grilli | Shoaling of solitary waves on plane beaches[END_REF][START_REF] Grilli | Breaking criterion and characteristics for solitary waves on slopes[END_REF].

Here we use experimental data taken from [START_REF] Grilli | Shoaling of solitary waves on plane beaches[END_REF] to compare with our numerical solution. For the numerical simulation we considered the domain shown in Figure 3. For this experiment we take a uniform mesh on [-100, with meshlength ∆x = 0.1. We also translate the solitary waves so that the crest amplitude is achieved at x = -20.1171. Because of the low regularity of the bottom at x = 0, artifacts of the numerical method might appear especially in the case of S 3 elements. For this reason, we ensure the required regularity of the seafloor (required by the model equations), by approximating the bottom topography function with a quadratic polynomial near x = 0. Although the SGN system cannot model accurately the breaking of solitary waves, it appears that it models shoaling with higher accuracy than other Boussinesq models [START_REF] Filippini | On nonlinear properties of enhanced Boussinesq models[END_REF], even very close to the breaking point. For example when A = 0.2 the steep wave observed in the laboratory at the location called gauge 9 of [START_REF] Grilli | Shoaling of solitary waves on plane beaches[END_REF] is approximated by a smooth solution in the SGN system.

Next, we study the shoaling of solitary waves with normalized amplitude A = 0.1, 0.15, 0.2 and 0.25. We monitor the numerical solution on the gauges (enumerated as in [START_REF] Grilli | Shoaling of solitary waves on plane beaches[END_REF]) with number 0, 1, 3, 5, 7 and 9 located at x = -5.0, 20.96, 22.55, 23.68, 24.68, 25.91 respectively. The results with P 1 elements in the case of the shoaling of the solitary wave with A = 0.2 are shown in Figure 4. The results with the other methods are almost identical and are not shown here. In these experiments, we also monitored the invariant I(t) for values of t small enough that the waves do not interact with the boundaries. In the case of S 3 elements on a uniform grid with ∆x = 0.1, the values of the invariant up to time t = 30 are given in the Table 14. For the same experiments, the respective invariants were conserved to 4 decimal digits in the case of P 1 elements with ∆x = 0.1. Finally, we computed the relative wave height defined as H(x * )/b(x * ) A I(t) 0.10 0.104058609813 0.15 0.197139475070 0.20 0.312548348249 0.25 0.449208354485 Table 14. The conserved invariant I(t) for up to t = 30.

where H(x * ) = max x {|ζ(x, t)|} is the crest amplitude of the wave normalized by the with experimental data of [START_REF] Grilli | Shoaling of solitary waves on plane beaches[END_REF] for the shoaling of solitary waves on a plane beach of slope 1 : 35, and with Green's law. '-': Numerical solutions, '•': Experimental data, '-•': Green's law for Boussinesq systems local depth b(x * ) evaluated at the same point x * . Figure 5 presents the experimental and numerical data. We observe a very good agreement between the numerical and the experimental results. We note that the numerical results obtained by the numerical solution of the full water wave problem in [START_REF] Grilli | Shoaling of solitary waves on plane beaches[END_REF] fit the experimental data (see e.g. Fig. 4 of [START_REF] Grilli | Shoaling of solitary waves on plane beaches[END_REF]), as well as the SGN model.

It is noted that the evolution of the maximum of the solution according to Green's law for Boussinesq systems, [START_REF] Synolakis | Evolution of maximum amplitude of solitary waves on plane beaches[END_REF], represented in Figure 5 by a broken line is initially close to the numerical and experimental data, then underestimates the evolution of the maximum especially during the shoaling of large amplitude solitary waves. This indicates that the solitary waves are in between the zone of rapid and gradual shoaling since they are nearly breaking waves. On the other hand Green's law for the nonlinear shallow water system, [START_REF] Synolakis | Green's law and the evolution of solitary waves[END_REF], which is omitted here, overestimates the amplitude of the shoaling solitary wave, likely because of the simplifications inherent in the derivation of Green's law. These results verify the conclusions made in [START_REF] Synolakis | Evolution of maximum amplitude of solitary waves on plane beaches[END_REF] about the zone of gradual shoaling where a law that is similar to Green's law η max ∼ b -1/4 is valid. In our case we computed by experimentation that η max ∼ b -2/7 describes quite well the evolution of the amplitude of a shoaling solitary wave on a plane beach of slope : 35.

4.2.

Reflection of solitary waves on vertical wall. We study also the reflection of the solitary waves on a vertical wall. In general, a vertical wall can be modeled by assuming that there is no flux through the vertical wall, i.e. the horizontal velocity of the fluid on the wall is u = 0. The reflection of a solitary wave at a vertical wall is equivalent to the head-on collision of two counter-propagating solitary waves of the same shape. Because the head-on collision implies that at the center of the symmetric interaction η x = 0 then one might argue that the symmetric head-on collision is governed by different mathematical properties compared to the interactions of a solitary wave at a wall.

In practice, during the interaction of a solitary wave with a wall the derivative η x on the boundary is negligible and so the additional condition seems to be satisfied and the reflection at the wall is apparently equivalent to the head-on collision of two counter propagating waves, [START_REF] Antonopoulos | Numerical solution of the 'classical' Boussinesq system[END_REF]. (Although the boundary condition η x = 0 is not necessary, its imposition does not require any modification to the numerical method or to the finite element spaces). The interaction of solitary waves with a vertical wall and the head-on collision of the solitary waves has been studied by theoretical and numerical means in [START_REF] Mirie | Collision between two solitary waves. Part 2. A numerical study[END_REF]. In [START_REF] Mirie | Collision between two solitary waves. Part 2. A numerical study[END_REF], an asymptotic solution for the maximum runup has been derived, namely, if α is the normalized amplitude A/b of the impinging wave, then the normalized maximum runup Rmax

= R max /b is approximately Rmax ∼ 2α + 1 2 α 2 + 1 2 α 3 . (17)
In [START_REF] Mirie | Collision between two solitary waves. Part 2. A numerical study[END_REF], the SGN system with horizontal bottom has been solved numerically using a finite difference scheme. Finite difference schemes for the SGN equations have certain disadvantages, including the introduction of numerical dispersion or dissipation. In addition it is necessary to use less physical boundary conditions on the wall, such as h x = u xx = u -1 3h (h 3 u x ) x = 0, [START_REF] Mirie | Collision between two solitary waves. Part 2. A numerical study[END_REF]. A fourth-order compact finite volume scheme was used to solve the SGN with general bathymetry in [START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part I: Model development and analysis[END_REF][START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part II: Boundary conditions and model validation[END_REF]. The derivation of the boundary conditions proposed in [START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part II: Boundary conditions and model validation[END_REF] follow analogous techniques with [START_REF] Peregrine | Long waves on beaches[END_REF] imposing also more boundary conditions for both h and u that approximate the wall boundary conditions. Other useful boundary conditions such as absorbing boundary conditions were constructed in [START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part II: Boundary conditions and model validation[END_REF] and can be implemented equally well with the present numerical model.

The head-on collision of solitary waves in the full water wave problem has also been studied asymptotically, in [START_REF] Su | On head-on collisions between two solitary waves[END_REF]. Specifically, a similar formula to the one for the SGN equations has been derived, and is given by the formula 2α + 1 2 α 2 + 3 4 α 3 . Other somewhat less accurate asymptotic approximations have been derived in [START_REF] Pelinovsky | Runup of tsunami waves on a vertical wall in a basin of complex topography[END_REF]. Recent numerical experiments showed that BT models do not describe accurately the reflection of large amplitude solitary waves at a vertical wall. Specifically, in [START_REF] Cooker | Reflection of a high-amplitude solitary wave at a vertical wall[END_REF] solving the full water wave equations using boundary element methods, it was shown that large amplitude waves can achieve a higher maximum runup during a head-on collision with a vertical wall, than the predicted values of the respective asymptotic solutions of [START_REF] Su | On head-on collisions between two solitary waves[END_REF][START_REF] Mirie | Collision between two solitary waves. Part 2. A numerical study[END_REF] for solitary wave reflection at a vertical wall. These results have been verified in [START_REF] Chambarel | Head-on collision of two solitary waves and residual falling jet formation[END_REF][START_REF] Chen | An experintal study of steep solitary wave reflaction at a vertical wall[END_REF]. Additionally, the formation of a residual jet was observed during the head-on collision of two large amplitude solitary waves, which was responsible for the large values of the maximum runup height.

To verify the ability of the numerical method to model accurately the reflection of solitary waves, we considered the case of a horizontal bottom b(x) = -1 for x ∈ [-100, 0] and solitary waves with amplitudes A = 0.075, 0.1, 0.15, • • • , 0.7. This case and these waves match the experiments which are presented in [START_REF] Cooker | Reflection of a high-amplitude solitary wave at a vertical wall[END_REF] and which serve as benchmarks for the models in [START_REF] Chambarel | Head-on collision of two solitary waves and residual falling jet formation[END_REF][START_REF] Chen | An experintal study of steep solitary wave reflaction at a vertical wall[END_REF]. We also consider a uniform grid on the interval [-100, 0] with ∆x = 0.1 while we translate the solitary waves so that their crest amplitude is at x = -50. Figure 6 shows the reflection of solitary waves with A = 0.075 and 0.65.

The reflection at a vertical wall of the extreme case of solitary waves of amplitude A = 0.7 is presented in Figure 7, together with the solution of the head-on collision of two symmetric solitary waves. In this figure, we cannot observe any differences between the solutions of the reflected wave and of the colliding solitary waves (within graphical accuracy). During this interaction a jet-like structure is visible, similar to those reported in [START_REF] Chambarel | Head-on collision of two solitary waves and residual falling jet formation[END_REF][START_REF] Chen | An experintal study of steep solitary wave reflaction at a vertical wall[END_REF][START_REF] Maxworthy | Experiments on collisions between solitary waves[END_REF]. The jet formed in the case of the SGN system is a solution of the mathematical model while the respective jet computed in [START_REF] Chambarel | Head-on collision of two solitary waves and residual falling jet formation[END_REF] is not. It is noted that no wave breaking is observed during these head-on collisions. Due to the inelastic interaction, the reflected solitary wave is followed by a dispersive tail. The magnitude of the dispersive tail depends on the amplitude of the colliding solitary wave. As it can be observed in Figure 6 the reflection of large solitary waves results in the generation of large dispersive tails.

In [START_REF] Chambarel | Head-on collision of two solitary waves and residual falling jet formation[END_REF], it was reported that the formation of a residual jet for the full water wave equations starts when the normalized amplitude of two incident solitary waves is larger than 0.60. In the case of the SGN equations, this jet formation was observed in reflections of solitary waves with amplitudes larger than 0.65. A comparison of the numerical maximum runup values and the asymptotic formula [START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part II: Boundary conditions and model validation[END_REF] in Figure 8 shows a very close match between the numerical solutions and the asymptotic results. In this figure, the numerical results of both the piecewise linear and the cubic spline finite elements are presented but no differences can be observed within graphical accuracy. Moreover, the numerical results of [START_REF] Cooker | Reflection of a high-amplitude solitary wave at a vertical wall[END_REF] are compared with the numerical results obtained with the modified Galerkin method and we verify the difference in the head-on collision processes between the SGN and the Euler equations, [START_REF] Cooker | Reflection of a high-amplitude solitary wave at a vertical wall[END_REF][START_REF] Chambarel | Head-on collision of two solitary waves and residual falling jet formation[END_REF]. This can be explained by noting that the jet cannot be modeled by the SGN equations in the present form, since this artifact cannot be described by smooth function but only by a parametric curve, [START_REF] Chambarel | Head-on collision of two solitary waves and residual falling jet formation[END_REF].

Remark 4.1. In some cases it is considered that the reflection of a wave-train of more than one pulse can result in extreme runup values. Extreme wave runup on a vertical wall has been studied using periodic-boundary conditions and the head-on collision of wave-trains propagating in different directions in [START_REF] Viotti | Conditions for extreme wave runup on a vertical barrier by nonlinear dispersion[END_REF][START_REF] Carbone | Extreme wave run-up on a vertical cliff[END_REF]. We repeated several of these experiments using half of the domain (required by the periodic code and with the wall boundary conditions described in this paper). The results were identical (except for the accuracy) to those obtained in [START_REF] Viotti | Conditions for extreme wave runup on a vertical barrier by nonlinear dispersion[END_REF][START_REF] Carbone | Extreme wave run-up on a vertical cliff[END_REF] while our code remained stable during the strong interactions of the waves with the wall. This verifies that the wall boundary conditions and the new numerical method can describe accurately strong interactions and extreme wave runup on a wall. As an indication of the results obtained we report here only the case of the reflection of a wave-train with three pulses (N w = 3, λ 0 = 125d in the notation of [START_REF] Viotti | Conditions for extreme wave runup on a vertical barrier by nonlinear dispersion[END_REF]). In this case the normalized maximum runup observed was R max /a 0 ≈ 5.545. 4.3. Reflection of shoaling waves. In our final set of numerical experiments we consider the benchmark described in [START_REF] Walkley | A finite element method for the one-dimensional extended Boussinesq equations[END_REF][START_REF] Dodd | A numerical model of wave run-up, overtopping and regeneration[END_REF]. Specifically, we consider the computational domain [-100, 20] (with ∆x = 0.1) and a bottom topography consisting of a horizontal seafloor in [-100, 0] and a plane beach of slope 1 : 50 for x ∈ [0, 20]. A sketch of the computational domain is presented in Figure 9. The vertical wall is located at x = 20.

In this set of experiments we study the reflection of solitary waves at a vertical wall after the waves have climbed up the sloping beach. We consider two cases, one with a solitary wave of amplitude A = 0.07 and another with A = 0.12. The initial conditions have been translated so that the maximum of the crest is at x = -30. Figures 10 and11 show experimental data recorded at three locations, by wave gauges g 1 , g 2 , and g 3 , placed at x = 0, 16.25 and 17.75. The comparison of numerical solutions with experimental data shows a better match from the present SGN model than is obtained using any weakly nonlinear and weakly dispersive Boussinesq system, [START_REF] Filippini | On nonlinear properties of enhanced Boussinesq models[END_REF][START_REF] Walkley | A finite element method for the one-dimensional extended Boussinesq equations[END_REF]. These numerical experiments verify the ability of the present numerical scheme to approximate with high accuracy and model fully nonlinear and weakly dispersive waves with wall boundary conditions at the endpoints of the computational domain. In order to study the stability of the modified Galerkin method in more demanding situations, we consider the propagation of a solitary wave over a composite beach simulating the geometrical dimensions of Revere Beach, and its reflection by vertical wall. These experiments were conducted at the Coastal Engineering Laboratory of the U.S. Army Corps of Engineers, Vicksburg, Mississippi facility, [START_REF] Kanoglu | Long wave runup on piecewise linear topographies[END_REF] and serve as benchmarks for the reflection of nonbreaking, nearly breaking and breaking solitary waves by vertical wall. The composite beach consists of three piecewise linear segments while the bathymetry is constant away from the beach and equal to b 0 = 0.218 m. The bathymetry can be realized by the function: The sketch of the domain is presented in Figure 12.

Here, we considered three solitary waves with normalized amplitudes A/b 0 = 0.05, 0.3 and 0.7. We monitored the water depth at x locations that correspond to gauges 5, 7 and 9 in [START_REF] Kanoglu | Long wave runup on piecewise linear topographies[END_REF]. In this experiment, we used ∆x = 0.1 in [-11.77, 23.23]. The normalized maximum runup computed at the location of the vertical wall in the first case was R/b 0 = 0.122 which is very close to the experimental value measured in the laboratory R/b 0 = 0.13. Moreover, the computed solution is very close to the experimental data at the gauges and we don't present these results here. It is noted that because a very small relative amplitude wave is involved in this case, its reflection can be modelled quite accurately, even by nondispersive models. The other two cases involve a nearly breaking and a breaking wave. Due to the steepness of the wave in the third case a wave breaking mechanism should be considered in order to approximate the solution in a stable manner. In the second experiment, the solitary wave is a nearly breaking wave. Although in this case the wave becomes very steep during shoaling, the numerical maximum runup computed was 0.46 m, which is again very close to the experimentally recorded runup value 0.45 m. The solution at the wave gauges is presented in Figure 13. The results can be improved by considering wave breaking mechanisms [START_REF] Sørensen | Boussinesq-type modelling using an unstructured finite element technique[END_REF][START_REF] Antunes | On breaking waves and wave-current interaction in shallow water: A 2DH finite element model[END_REF][START_REF] Zelt | The run-up of nonbreaking and breaking solitary waves[END_REF] and Green-Naghdi models with improved dispersion characteristics such as those proposed in [START_REF] Beji | A formal derivation and numerical modelling of the improved Boussinesq equations for varying depth[END_REF][START_REF] Dias | On the fully-nonlinear shallow-water generalized Serre equations[END_REF][START_REF] Antunes | Extended Serre equation for applications in intermediate water depths[END_REF][START_REF] Antunes | Boussinesq and Serre type models with improved linear dispersion characteristics: Applications[END_REF][START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF]. It is noted that in both cases only linear terms with second order derivatives are included locally around the regions where the waves become steep. These terms can be incorporated easily by any numerical method while the possible stiffness induced by the new terms to the problem can be handled by taking smaller time-steps or by using a time integration method appropriate for stiff problems, [START_REF] Antunes | On breaking waves and wave-current interaction in shallow water: A 2DH finite element model[END_REF].

Summary and Conclusions

We present a fully discrete numerical scheme for the Serre-Green-Naghdi (SGN) system with wall boundary conditions. Semidiscretization of the model equations is based on a modification of the standard Galerkin / finite-element method that allows solutions in function spaces of low regularity. The time discretization is based on a fourth-order, four-stage, explicit Runge-Kutta method. A detailed computational study of the convergence properties of the numerical scheme shows that the method converges with similar convergence rates to those of Peregrine's system. Some of the advantages of the new numerical method are:

• the high accuracy and the very good conservation properties;

• the sparsity of the resulting linear systems;

• the low complexity of the algorithm due to the use of low-order finite element spaces; and • its potential to be extended to the two-dimensional model equations. In addition, we perform a series of very accurate numerical experiments to verify the efficacy of the numerical scheme in studies of shoaling and reflecting solitary waves. The numerical solutions are compared with available experimental data and theoretical approximations, whenever possible. The numerical model appears to be efficient and the match between numerical results, experimental data, and theoretical approximations is very satisfactory and shows better performance than other shallow water wave systems when there is no wave breaking. Wave breaking can be treated by following heuristic methodologies, such as discarding dispersive terms or adding new dissipative terms.

  dx, for all φ ∈ S h , dx, for all ψ ∈ S u .

4

  j=1 b j U n,j end for implementation of Algorithm 1 requires solving at each time step the following linear systems:(a) Four linear systems with the time-independent matrix (φ i , φ j ); (b) Four linear systems with the time-dependent matrix B(ψ i , ψ j ; h);

Figure 1 .

 1 Figure 1. Sensitivity of energy conservation to the value of ∆t for fixed ∆x = 0.1

Figure 3 .

 3 Figure 3. Sketch of the domain for the shoaling of solitary waves on a plain beach of slope 1 : 35.

Figure 4 .

 4 Figure 4. Solution on various wave gauges for the shoaling on a plane beach of slope 1 : 35 of a solitary wave with A = 0.2. Circles show experimental data, and lines show numerical solutions.

Figure 5 .

 5 Figure 5. Comparison of the computed relative wave height H/bwith experimental data of[START_REF] Grilli | Shoaling of solitary waves on plane beaches[END_REF] for the shoaling of solitary waves on a plane beach of slope 1 : 35, and with Green's law. '-': Numerical solutions, '•': Experimental data, '-•': Green's law for Boussinesq systems

Figure 6 .

 6 Figure 6. Reflection of solitary waves at a vertical wall located at x = 0.

41 Figure 7 .

 417 Figure 7. Reflection of solitary waves of normalized amplitude A = 0.7 at a vertical wall located at x = 0, plotted together with the head-on collision of two solitary waves of the same amplitude to show their equivalence.

P 1 Figure 8 .

 18 Figure 8. Maximum runup values during the reflection of solitary waves at a vertical wall located at x = 0.

Figure 9 .

 9 Figure 9. Sketch of the numerical experiment for the reflection at a vertical wall located at x = 20 of a shoaling wave over a plain beach of slope 1 : 50.

Figure 10 .

 10 Figure 10. Reflection at a vertical wall located at x = 20 of a shoaling wave over a plain beach of slope 1 : 50. Initial solitary wave amplitude A = 0.07.

Figure 11 .Figure 12 .

 1112 Figure 11. Reflection at a vertical wall located at x = 20 of a shoaling wave over a plain beach of slope 1 : 50. Initial solitary wave amplitude A = 0.12.

Figure 13 .

 13 Figure 13. Reflection of a solitary wave over a composite beach. Vertical wall located at 23.23. Initial solitary wave amplitude A/b 0 = 0.3.

Table 2 .

 2 Spatial errors and rates of convergence for the exact solution with P 1 finite elements using the H 1 norm.

	1.8855 × 10 -1 40 1.3684 × 10 -1 80 1.0135 × 10 -1 160 7.3388 × 10 -2 320 5.2503 × 10 -2 640 3.7340 × 10 -2	0.6129 0.4625 0.4331 0.4658 0.4831 0.4917	6.6788 × 10 -2 3.2009 × 10 -2 1.5859 × 10 -2 7.9074 × 10 -3 3.9492 × 10 -3 1.9736 × 10 -3	1.3433 1.0611 1.0131 1.0041 1.0016 1.0007
	N	E ∞ [H]	rate for E ∞ [H] E ∞ [U ]	rate for E ∞ [U ]
	10 3.4775 × 10 -2 20 1.1039 × 10 -2 40 5.8666 × 10 -3 80 3.4929 × 10 -3 160 1.8861 × 10 -3 320 9.7698 × 10 -4 640 4.9684 × 10 -4	1.4587 1.6554 0.9121 0.7481 0.8890 0.9490 0.9755	3.9193 × 10 -2 5.7364 × 10 -3 1.0352 × 10 -3 2.1191 × 10 -4 5.2158 × 10 -5 1.3480 × 10 -5 3.4245 × 10 -6	1.4068 2.7724 2.4701 2.2884 2.0225 1.9520

Table 5 .

 5 Spatial errors and rates of convergence for the exact solution with P 2 finite elements using the H 1 norm.

	3-s )

  suggests a rate of convergence for h close to

	N	E 1 [H]	rate for E 1 [H] E 1 [U ]	rate for E 1 [U ]
	10 6.9247 × 10 -2 20 3.3950 × 10 -2 40 1.6779 × 10 -2 80 8.3767 × 10 -3 160 4.1860 × 10 -3 320 2.0924 × 10 -3 640 1.0460 × 10 -3 Table	-1.0283 1.0168 1.0022 1.0008 1.0004 1.0002	1.5110 × 10 -2 2.6966 × 10 -3 5.7095 × 10 -4 1.3474 × 10 -4 3.3114 × 10 -5 8.3275 × 10 -6 2.0570 × 10 -6	-2.4863 2.2397 2.0831 2.0247 1.9915 2.0173

Table 11 .

 11 Spatial errors and rates of convergence for the exact solution with S 3 finite elements using the H 1 norm.

	for E 1 [U ]

Acknowledgment

The authors were supported by the Marsden Fund administered by the Royal Society of New Zealand.