Sylvie Van Iseghem

Emmanuelle Quillérou

Cécile Brigaudeau

Claire Macher

Olivier Guyader

Fabienne Daurès

Ensuring representative economic data: survey data-collection methods in France for implementing the Common Fisheries Policy

Keywords: data collection, economic indicators, fishery-dependent information, panel data, sampling, small vessels, statistical precision

Since 2001, Ifremer has implemented an economic data collection programme (EDCP) within the Data Collection Framework of the EU. It aims to obtain economic data from a sample of vessels representative of the entire French fishing fleet. This paper presents the strategies used for vessel sampling selection in the French EDCP and its implementation over several consecutive years. The approach is illustrated by the sampling plan for the fleet in the North Sea Channel Atlantic region. We show that the EDCP allows precise economic indicators such as gross revenue or fuel costs to be estimated for the whole fishing fleet, including small vessels (<10 m), and consequently, it facilitates sound scientific advice regarding the Common Fisheries Policy. The apparent underestimation of landings by small vessels in official statistics is highlighted.

Introduction

In 1995, the Code of Conduct for Responsible Fisheries was adopted by the UN's Food and Agriculture Organization (FAO, 1995). It emphasized that the development of fisherymanagement plans required appropriate and reliable data on all aspects of a fishery. In particular, the Code stressed that "in order to ensure the sustainable management of fisheries and to enable social and economic objectives to be achieved, sufficient knowledge of social, economic and institutional factors should be developed through data gathering, analysis and research" (FAO, 1995, p. 12). Based on these considerations, the European Council (EC) of the European Union (EU) decided in 2000 to establish a Data Collection Framework (DCF; EC, 2008) to evaluate the situation of the fishing sector under the Common Fisheries Policy (CFP). The DCF specifies the type of data to be collected and the economic indicators to be computed. It also specifies the sample representativeness to be achieved within different fishing areas and fleet segments, defined as the combination of a particular fishing method and a vessel length category.

Today, the key role of economics in the development of fishery-management plans has become well established within European and international institutions as well as in the scientific community [START_REF] Pascoe | Economics, fisheries, and the marine environment[END_REF]. Many studies have indeed shown the importance of economics in the ecosystem approach to fisheries, e.g. [START_REF] Garcia | Ecosystem approach to fisheries: a review of implementation guidelines[END_REF] and [START_REF] Daurès | Fishing fleet typology, economic dependence, and species landing profiles of the French fleets in the Bay of Biscay, 2000-2006[END_REF]. A good economic assessment of fisheries, however, depends on the quality, reliability, and completeness of the fishery-dependent data used. This consideration has been mostly neglected in the literature so far, with few studies explicitly assessing the quality of data [START_REF] Chuenpagdee | Bottom-up, global estimates of small-scale fisheries catches[END_REF]. [START_REF] Mettling | Rapport d'audit sur la situation financière des navires de pêche artisanale et des organismes d'intervention[END_REF] and Le [START_REF] Pape | The influence of vessel size and fishing strategy on the fishing effort for multispecies in north-western France[END_REF] noted the underestimation in French official statistics of the landings by small vessels (here defined as being <10 m long). As small vessels constitute 59% of the French fleet (excluding overseas fleets) in terms of number of vessels, and 25% of the total engine power, this underestimation is an important problem to be addressed. This paper presents the two-stage methodology developed by Ifremer from 2001 for its economic data collection programme (EDCP) to guarantee the representativeness of sampling data across the whole fishing fleet. The EDCP economic database has, for instance, been used to estimate the economic dependence and species-landing profiles of French fleets in the Bay of Biscay [START_REF] Daurès | Fishing fleet typology, economic dependence, and species landing profiles of the French fleets in the Bay of Biscay, 2000-2006[END_REF]. It also plays a key role in the impact assessment of management measures [START_REF] Macher | A cost-benefit analysis of improving trawl selectivity in the case of discards: the Nephrops norvegicus fishery in the Bay of Biscay[END_REF] and other ongoing work in partnership with fishers and their representatives, government administrations, biologists, sociologists, and fishery economists. The first stage optimizes the sampling plan for a given year across fleet segments and distributes samples across maritime registration districts (MRDs). This sampling design is suitable whether the economic heterogeneity within a fleet segment is higher (e.g. small vessels) or lower (e.g. longer vessels). The second stage aims to reduce non-responses to questionnaires, through differential grouping of interviewees in MRDs. A network of permanent observers who collect data for various purposes, of which economic analysis is just one, also contributes to increase response [START_REF] Leblond | The Fisheries Information System of Ifremer: a multidisciplinary monitoring network and an integrated approach for the assessment of French fisheries, including small-scale fisheries[END_REF]. Results are presented in terms of the quality of the survey design (bias and precision), following the recommendations of STECF (2009b), non-response rates, and economicindicator estimates such as gross revenue and fuel costs. We also provide evidence of the underestimation of small-vessel landings in official statistics, which has rarely been reported in the literature.

Fishing fleet and sampling strategies

Ifremer collects economic data for French vessels <40 m long operating in three supraregions defined in the DCF regulation: (i) the North Sea Channel Atlantic (NSCA), i.e. ICES Divisions IIIa, IV, VIId, and V-XIV; (ii) the Mediterranean, and (iii) other regions. The sampling strategies are identical in all three cases, so we focus on the NSCA for illustration purposes (Figure 1). In 2007, 3721 French vessels were fishing in that area, most of them small (50% of the total fleet in numbers, 25% of the total engine power).

The division of the NSCA fleet into segments follows European Council regulation (EC, 2008), and is derived from both official data supplied by the French Ministry of Fisheries and a complementary survey of fishing activities carried out by Ifremer [START_REF] Berthou | From fleet census to sampling schemes: an original collection of data on fishing activity for the assessment of the French fisheries[END_REF]. A fleet segment is defined by the vessel type and size, with small vessels (see above) being in a single size category. The heterogeneity of fishing revenues within a fleet segment is assessed by computing the coefficient of variation (CV) across all vessels of a parameter of interest, P. To demonstrate the approach, P is taken to be a proxy for fishing revenue, which is calculated by multiplying vessel capacity (length of vessel × crew size) by the intensity of fishing, i.e. the number of months in which the vessel was active for at least one day [START_REF] Daurès | Re-assessing the French small-scale coastal fisheries from fleet activity to economic performance[END_REF]. Details of the size of each NSCA fleet segment and the variation in the proxy fishing revenue within it (measured by the CVs) are presented in Table 1. The heterogeneity in proxy fishing revenue is high, because the CVs vary between 0.18 and 0.75 across fleet segments (considering those with >10 vessels). The CV is highest for small vessels using pots, and lowest for dredgers between 18 and 24 m long. More generally, the heterogeneity in proxy fishing revenue is higher for smaller vessels than for larger ones, with an overall CV of 0.7 for the small-vessel segments.

The NSCA fleet thus comprises small vessels which are very heterogeneous, and larger vessels which are more homogeneous in terms of their fishing revenues.

Two-stage sampling plan to obtain representative data

The sampling design is probability sampling (STECF, 2009b), in which data are collected from a sample of vessels selected randomly. The plan is in two stages: (i) computation of the minimum sample size to be achieved for each fleet segment to meet the requirements for precision of sample-based estimates; and (ii) the allocation of these samples across MRDs to preclude geographic bias. The MRDs correspond to the vessels' administrative attachment, so vessels from the same MRD can actually fish in different areas.

In the DCF regulation, the precision L (measured by the CV of the estimator, or CVe) to be achieved for the fleet estimate of the parameter of interest P is a CVe of 0.25 or better. Using standard sampling theory, the ideal sample-allocation plan, i.e. providing the best precision for the least cost, can be determined from Neyman allocation estimates for each fleet segment [START_REF] Ardilly | Les Techniques de Sondage[END_REF]. The minimum sample size n i for segment i is computed from the required precision (L = 0.25), the segment size N i , and the population CV of the parameter of interest (CV i ) as

2 2 4 1 1 i i i i CV L N N n   . ( 1
)
The budget allocated for the programme, however, limits the total sample size to a maximum of 600 vessels (about 17% of the NSCA fleet), which allows a precision L of the implemented sampling plan between 0.20 and 0.25 for all fleet segments to be achieved.

The n i vessels to be sampled in fleet segment i are allocated between MRDs to ensure the geographic representativeness of the collected data. For this, vessels are selected by systematic random sampling [START_REF] Tillé | Théorie des sondages[END_REF]) in five steps.

(i) All vessels of the fleet segment are ordered by maritime registration district (from north to south), and then by vessel length within each district. (ii) n i statistical units (SU) of size m i (= N i n i ) are identified from this ordered population;

each SU defines a group of vessels which is considered statistically equivalent. (iii) The first vessel of the sample is identified from a starting point chosen at random inside the first SU. (iv) Subsequent vessels are thereafter selected at regular intervals of size m i , so n i vessels are selected, one in each SU. (v) If the owner does not wish to take part in the survey, that vessel is substituted by another that belongs to the same SU. For example, if n i = 12 and N i = 100 in step (ii), m i = 8 is calculated. If now in step (iii), the 7th vessel belonging to the first SU is selected by chance among vessels 1-8, then in step (iv), the next one to be selected is the 7 + 8 = 15th vessel, then the 23rd vessel etc., moving down the list of vessels.

Organizing data collection to reduce non-responses

The EDCP is part of Ifremer's Fishery Information System (FIS; [START_REF] Leblond | The Fisheries Information System of Ifremer: a multidisciplinary monitoring network and an integrated approach for the assessment of French fisheries, including small-scale fisheries[END_REF] and is included in the French national programme implemented under the DCF (Ministry of Fisheries, 2008). The collection of economic data is made through a network of Ifremer observers covering the French MRDs. The observers are also involved in other work, such as monitoring fishing-vessel activities and biologically sampling discards. The sampling plan is organized so that each observer has a portfolio of some 30 economic surveys to carry out in a specified maritime district. The observers are responsible for both arranging the interviews and entering the information collected into the database.

Interviews take place between February and June each year to collect economic data for the previous year, i.e. the year of surveyed activity. To ensure data consistency and completeness over the surveyed year, the interviewee needs to have full knowledge over the relevant period. That is why vessel owners are interviewed personally. In order to reduce the number of non-responses, three vessel categories are distinguished. The first, referred to as the panel population, are vessels that were successfully interviewed the year before. A panel vessel cannot be interviewed for more than four consecutive years. The second category contains the non-interviewable vessels, those that cannot be interviewed either because they were sold during the previous year (so that data are incomplete over the year of activity), or they were approached the year before and refused to join the panel. The third category is the structural population, which includes all other fishing vessels. The panel and structural subpopulations are sampled independently as described above, and the global sample is the combination of the two. The number of panel vessels approached is equal to the minimum between the number of panel vessels N ip in segment i and n ipM = 0.7n i . If N ip < n ipM , then all panel vessels are approached. Subsequently, the number of structural vessels approached, n is , is n i -n ip . For example, if n i = 10, N i = 40, and N ip = 12, then n ip = min(12,10×0.7) = 7. The global sample will be a combination of seven panel vessels and three structural vessels, independently identified as described above.

The first economic questionnaire survey was carried out in 2001, in consultation and collaboration with economists from the University of Western Brittany (see [START_REF] Boncoeur | Enquête économique sur la pêche professionnelle française en Manche[END_REF], and has continued with few changes over the years. The questions are based on a monetary approach that also includes social indicators. The questionnaire consists of nine main sections in logical order; vessel owners are surveyed on their fishing activities, their related earnings and costs, and more specifically on the characteristics of their fishing vessel. The latter includes the type of hull, engine, fishing gear, electronic equipment, and the related costs. Technical documents (in French) describing the methodology, from sampling optimization to the questionnaire, and the software used to record the information collected, are available on the Ifremer website http://www.ifremer.fr/sih/.

Data collected from the economic survey are validated vessel by vessel and variable by variable in several steps, including tests of quality, consistency, and continuity between variables. All economic variables are collected in accord with the definitions in Appendix VI of EC (2008).

Results of the EDCP

As it is impractical to present all results across the 34 fleet segments, we focus here on five chosen to represent the diversity of the NSCA fleet in terms of activity and vessel length. Three are small-vessel segments, dredgers (SSDredgers), drift and/or fixed netters (SSNetters), and vessels using pots and/or traps (SSPots). The other two are demersal trawlers and/or seiners 10-12 m long (referred to as MediumDemersalTrawlers, or MDTrawlers), and 12-18 m long (referred to as LargeDemersalTrawlers, or LDTrawlers).

From 2000 to 2008, the NSCA fleet decreased by 18%, from 4200 to 3500 vessels (Figure 2a). The small-vessel fleet declined less, by some 13% (Figure 2b). The NSCA sample varied from 425 vessels in 2008, the smallest sample, to 580 vessels in 2002, the largest sample (Figure 2a). The sample rate of the NSCA fleet showed no trend during the period, however, the budget restriction (600 vessels, or 17% of the total fleet) never being realized; the sample rates achieved were between 13 and 15% (Figure 2a). Sampling rates of small vessels and the total fleet were similar (Figure 2a and 2b).

Reliability of the sample

The total sample size in 2008 was 533 vessels, i.e. 14.4% of the fleet (Table 1). The largest sample size was that for small vessels using hooks (59 vessels, sampling rate 19%). The sampling-rate range for all segments with >10 vessels was 8-30%. The range of precision for the proxy fishing revenue for small vessels was 2-12%, well within the DCF requirement. However, the precision was >30% for two small fleet segments: vessels using hooks from 18 to 24 m long, and drift and fixed netters from 24 to 40 m long. Three other fleet segments had a precision between 25 and 30%, again not in accord with the regulation.

As economic assessments are based on samples and not a census, it is important to look for any bias in the sampling design (STECF, 2009b). We examined sample representivity using Welch's t-test, an adaptation of the Student's t-test intended for use with two samples with possibly unequal variance. The Welch's t-test indicates whether the mean characteristic of the sample is the same as the population mean. For it to be applicable, the sample has to be either normally distributed or >30 vessels [START_REF] Tufféry | Data Mining et statistique décisionnelle[END_REF]. This test was performed for the five fleet segments of the NSCA fleet for the variables vessel length, vessel power, and vessel age. The null hypothesis that the means of two groups (population/sample) were equal at a risk of 5% was tested, and the results demonstrated that sample and population means were not significantly different for any of the five fleet segments selected and for all vessel characteristics tested (results not shown).

The spatial distribution across MRDs of the sample collected for the five fleet segments analysed here was close to that of the population, though with small differences (Figure 3).

Substitute vessels and the number of non-responses

During the EDCP in 2008 and based on 2007 data, 938 fishers were contacted, resulting in 533 interviews. In all, 405 fishers did not wish to participate, and another 330 were substituted by the procedure explained above. No substitution was possible for the remaining 75 cases. Ultimately, therefore, there were 75 non-responses, 13% of the successful interviews. The panel sample consisted of 251 vessels, around half the total. Nine fishers in the panel population, having been interviewed the year before and who had agreed to be on the panel then, declined to be interviewed again in 2008.

Economic performance of selected NSCA fleets

Here we show how the EDCP sample information improves estimation of some economic indicators. It is not, however, a complete economic assessment of the NSCA fleet, because that has already been reported in STECF (2009a). The comparison of gross revenue obtained from the EDCP data and with the official statistics [START_REF] Ofimer | Key figures for the fisheries and aquaculture sector in France[END_REF] for the five fleet segments selected showed that the relative percentage difference between the two was negligible for LDTrawlers, and 15% for MDTrawlers, but much larger for the segments with small vessels (30-46%; Table 2).

We evaluated the economic performance for the period 2001-2009 using the following four indicators: gross revenue, fuel cost, operational costs (i.e. the sum of all costs other than fuel and crew related to fishing effort), and employment in terms of Full Time Equivalent (FTE). The latter is calculated as the average monthly crew size in a given year. Apart from the gross revenue presented above, the other parameters are not available from official data sources. Table 3 shows estimates of these indicators in 2007 as averages per vessel for the five fleet segments. The gross revenue ranged from €65 068 for SSDredgers to €454 356 for LDTrawlers. Fuel costs spanned the range €4 270 for SSNetters to €99 535 for LDTrawlers. Crew size varied between 1.5 FTE for SSNetters to 3.6 FTE for LDTrawlers.

To examine productivity, two indicators of effort (days at sea, and hours of engine running), and one capacity indicator, the average engine power per vessel, were considered. They provide estimates of four productivity ratios (R1-R4), listed in Table 3. The segments with small vessels had lower revenue productivity than the others. Revenue productivity in terms of gross revenue per crew size times days at sea (R1) was similar between the three segments with small vessels. The SSDredgers had the lowest revenue productivity in relation to engine power (R2). Energy productivity was relatively low for the segments with small vessels (Table 3).

Discussion and conclusions

We have presented the sampling methodology designed by Ifremer for the French EDCP. It provides economic data representative of the entire French fleet in terms of gross fishing revenues, across all fleet segments (including small vessels) and MRDs. The sampling strategy used in the EDCP is probability sampling. Errors in survey design were evaluated in terms of bias and precision following the recommendations of STECF (2009b). The sampling plan is representative of the overall fleet. There was no bias in terms of vessel characteristics of the sample, but a precision of economic indicators greater than 25% for most fleet segments. These results meet the requirements of the EU Data Collection Framework regulation.

The low level of final non-responses (13%) is the result of the data collection being based on a network of Ifremer observers who collect data for various purposes, of which economic analysis is just one. It also stems from the provision for substitute vessels in the sampling plan. Nevertheless, the variability in sample sizes through years was notable, and mainly attributable to changes in the non-response rate caused by factors outside Ifremer's control, such as fishery-management policy developments, that are difficult to address.

Comparison of the EDCP data (gross revenues) with official statistics revealed differences in landing values of small vessels of up to 46%. Given the representativeness of the sampling, it can be assumed that the official statistics are underestimates for these vessels. This justifies using the EDCP to estimate the economic importance of small-vessel segments of the fleet better. This is explained by the fact that official data are derived mainly from the records of auctions, which are the only commercial outlet for large vessels, whereas small vessels have alternative outlets such as direct sale to customers.

The economic indicators revealed considerable differences among the five fleet segments selected, in particular between small and large vessels; the former have greater heterogeneity and lower productivity. This result underlines the importance of assessing the economic performances of small and large vessels separately. Sumaïla et al. (2008) reported that the energy ratio (R4) varies from one fishery to another, and can be as much as 60%.

The NSCA ratios agree well with these observations, with the maximum being 71% for the MDTrawlers. The energy ratio was less for small vessels, especially for SSPots, for which it was just 32%.

The sample size achieved corresponds to 15% of the fishing fleet, i.e. 2% less than the sampling plan. The samples collected were often not sufficient to be useful for regional assessments, but the use of additional data sources such as accounting information is one way of increasing the sample size and, consequently, improving the precision of the economic indicators derived.

There is a great diversity of sampling approaches for collecting economic data in Europe, limiting the comparability and homogenization of data between Member States (STECF, 2009b). However, the effectiveness of different sampling schemes for economic data is currently being evaluated across Europe. Three types of data-collection scheme have been identified by the Joint Working Group on Economic Affairs (SGECA): (i) Type A is a census, with data collected from the entire population, (ii) Type B a probability-sample survey, in which data are collected from a random sample of the population (examples are the French EDCP described here, and a similar scheme in Italy; Italian Ministry of Agricultural and Forestry Policy, 2009), and (iii) Type C, a non-probability sample survey, in which data are collected from a sample of population members that is not selected randomly, e.g. Netherlands Directorate of Fisheries (2008).

Even for a given type of sample collection, strategies can be based on different optimization algorithms. We used a Neyman allocation across fleet segments based on the level of heterogeneity of revenues. The Italian scheme uses the Bethel methodology, a generalization of the Neyman allocation and based on the level of heterogeneity for a bundle of variables [START_REF] Sabatella | Manual of fisheries sampling surveys: methodologies for estimation of socio-economic indicators in the Mediterranean Sea[END_REF]. The Bethel methodology allows costs to be minimized for a given level of precision for the bundle of variables selected. Future research in relation to the French EDCP needs to investigate the sample representativeness for variables other than revenue. This would allow more direct comparisons between the French and Italian data collection programmes.

We only studied landings made in NSCA MRDs and future research could include those made elsewhere. Our spatial analysis is restricted to the MRDs of the vessels, but a wider analysis has been made in the Bay of Biscay [START_REF] Daurès | Fishing fleet typology, economic dependence, and species landing profiles of the French fleets in the Bay of Biscay, 2000-2006[END_REF], and other studies by fishing area are ongoing within Ifremer.

The interest of this work lies not just in the results of the sampling strategies, but in how the methodology itself was developed; these strategies could be applied elsewhere. The index of variability is the coefficient of variation, CV (standard deviation divided by mean), of the proxy fishing revenue, estimated as the length of the vessel × crew size × fishing intensity (months) on the whole population. The sample CV of the gross revenue (CVe) is based on 2007 data. The relative percentage difference between the gross revenue derived from raised EDCP samples (ledcp) and official landings (lo) is calculated as 100(ledcp -lo)/ledcp. Small-vessel segments (<10 m): dredgers (SSDredgers), drift and/or fixed netters (SSNetters), and vessels using pots and/or traps (SSPots); medium-sized demersal trawlers and/or seiners (10-12 m; MDTrawlers); large demersal trawlers (12-18 m; LDTrawlers). Table 3. Indicators of economic performance of selected vessels of the French North Sea Channel Atlantic fleet for 2007.

Tables

Figures

Figure 1 .

 1 Figure 1. Maritime registration district (MRD) of the French North Sea Channel Atlantic fleet, and the number of fishing vessels registered by size category (<10 m and 10 m) in 2008.

Figure 2 .

 2 Figure 2. Changes in the number of vessels in the French fleet operating in the North Sea Channel Atlantic region and the number of samples by year between 2000 and 2008. (a) All vessels, (b) small vessels (<10 m).

Figure 3 .

 3 Figure 3. Spatial distribution as the percentage of vessels by maritime district for five fleet segments for the population and the sample in 2007. Small-vessel segments (<10 m): dredgers (SSDredgers), drift and/or fixed netters (SSNetters), and vessels using pots and/or traps (SSPots); medium-sized demersal trawlers and/or seiners (10-12 m; MDTrawlers); large demersal trawlers (12-18 m; LDTrawlers).

Table 1 .

 1 Characteristics of the French fleet operating in the North Sea Channel Atlantic region in 2007 by length category, and details of the sampling plan.

			Number of vessels			Index of variability CV (%)				Sample size					Sample rate (%)		Sample CV of gross revenue (Cve -%)
	Vessel
		<10
		m
	Beam trawlers			6		2		0	0	17	0	13	0			2						33				4
	Demersal trawlers and/or demersal seiners	118	193	219	233	84	13	44	35	34	30	44	30	10	32	21	25	5		8	17	10	11	6	7	1	2	1	15
	Pelagic trawlers		4	20	40	3	4	0	26	45	28	19	48		1	4	4				25	20	10			20	16	7
	Pelagic seiners		2	25	2	1	22	0	24	30	54	0	16		1	5					50	20				12	6
	Dredgers	82	77	97	11			51	42	34	18	0	0	15	13	10	2			18	17	10	18		6	5	4	5
	Vessel using other active gears	251	8					69	32	0	0	0	0	35						14					5
	Vessels using polyvalent "active" gears only	31	54	50	3			57	44	28	92	0	0	8	12	11				26	22	22			12	5	2
	Vessels using hooks	303	45	8	8	10		59	36	49	65	29	0	59	14	3	3	1		19	31	38	38	10	2	3	20	35	30
	Drift and/or fixed netters	335	196	94	40	19	1	69	36	32	44	31	0	42	19	16	6	1		13	10	17	15	5	4	2	2	11	36
	Vessels using pots and/or traps	308	65	7	11			75	34	44	20	0	0	42	17	2	7			14	26	29	64		5	2	28	1
	Vessels using other passive gears	106	6					60	31	0	0	0	0	29						27					4
	Vessels using polyvalent "passive" gears only	97	10	3				66	41	14	0	0	0	22	2					23	20				6	27
	Vessels using active and passive gears	93	92	13	1			51	35	32	0	0	0	15	14	3				16	15	23			6	3	11
	Inactive vessels	152	21	8	12	2		0	0	0	0	0	0
	Total by number	1876	773	550	361	121	40	70	41	36	36	54	66	277	125	77	47	7	0	15	16	14	39	17	0
	Percentage	50	19	15	8	3	1

Table 2 .

 2 Comparison of gross-revenue, sample-based estimates and official landings statistics[START_REF] Ofimer | Key figures for the fisheries and aquaculture sector in France[END_REF] by fleet segment for 2007.

	Parameter	SSDredgers	SSNetters	SSPots	MDTrawlers LDTrawlers
	Number of vessels in fleet	82	335	308	193	219
	Number of vessels in sample	15	42	42	31	21
	Economic indicator, whole fleet (€)
	Official landings	3 731 889	17 389 799	17 684 290	39 449 883 100 519 643
	EDCP gross revenue	5 335 563	25 510 528	33 010 802	46 354 796	99 504 055
	Relative difference	30%	32%	46%	15%	-1%

Acknowledgements

We especially thank the fishers who agreed to participate in the EDCP. From Ifremer, we thank the team of Ifremer observers involved in the programme, E. Rostiaux, P, Raulin, J. Quinquis, D. Miossec, H. Barone, P. Raguenes, and G. Morandeau, Michèle Jézéquel for IT tool design and support, Mathilde Pitel for mapping the surveyed area, and S. Demanèche, E. Leblond, and P. Berthou for their contributions to the EDCP. We are also extremely grateful to V. Trenkel, the guest editor, to D. N. MacLennan, the journal editor, and to two anonymous reviewers for their helpful suggestions on the submitted version of this paper.


C Soufflet

F Lott

F Damiens

Trapped mountain waves with critical level just below the surface

Keywords: Mountain waves, Trapped waves, Kelvin-Helmholtz, Instabilities, Critical level

Trapped mountain waves with critical level just below the surface

C Soufflet, F. Lott, F

Introduction

Trapped mountains waves generally occur when the vertical profiles of the large scale flow favor the low level confinement of stationary gravity waves. They are sometimes associated with rotors [START_REF] Doyle | The dynamics of mountain-waveinduced rotors[END_REF]] and can produce dangerous weather situations near the ground [START_REF] Keller | Lee waves associated with a commercial jetliner accident at denver international airport[END_REF]]. They also contribute to the low level mountain wave drag [Lott(1998), [START_REF] Georgelin | On the transfer of momentum by trapped lee waves. case of the iop3 of pyrex[END_REF], [START_REF] Teixeira | Drag produced by trapped lee waves and propagating mountain waves in a two-layer atmosphere[END_REF], ] and still today, it is well established that low level mountain drags improve weather forecast and climate models [START_REF] Sandu | Impacts of parameterized orographic drag on the northern hemisphere winter circulation, journal of advances in modeling earth systems[END_REF], [START_REF] Pithan | Missing orographic drag leads to climate model biases in jet streams, blocking and storm tracks[END_REF], ]. For this reason, there is some debate on which source of low-level mountain drag should be better represented, the results in [START_REF] Tsiringakis | Small-scale orographic gravity wave drag in stable boundary layers and its impact on synoptic systems and near-surface meteorology[END_REF]] suggesting that the trapped lee waves could contribute as much as the blocked flow mountain drag [START_REF] Lott | A new subgrid-scale orographic drag parametrization: Its formulation and testing[END_REF], ] or the turbulent subgrid scale orographic drag [1, ]. According to the conventional theory, the low level confinement favoring lee waves can have two dynamical mechanisms. The first attributes the trapped lee waves to free oscillations along a low level density discontinuity. This mechanism has been re-examined recently by [Vosper(2004)] and [START_REF] Sachsperger | The amplitude of lee waves on the boundary-layer inversion[END_REF]], the first extending the theory to situations where the discontinuity results from a boundary-layer inversion underneath a continuously-stratified atmosphere the latter quantifying analytically the amplitude of the waves in the nonlinear context.

The second mechanism, first identified by [Scorer(1949)], involves the trapping of vertically propagating internal gravity waves between the ground and a turning point aloft. In this theory the Scorer parameter

S = N 2 U 2 - 1 U d 2 U dz 2 , (1)
plays a central role, N, U and z being the background buoyancy frequency, the background wind and the vertical coordinate respectively. If S(z) decreases with height, due to a diminution of N with altitude for instance, some gravity waves that propagate vertically near above the surface become evanescent aloft. For a given mode with horizontal wave number k, this transition occurs at the turning level where S(z k ) = k 2 , and at this level the corresponding wave is reflected downward. The onset of trapped waves then strongly depends to what happen to the reflected waves when they return to the surface. In the initial theory of Scorer (1949), the flow is inviscid and the downward waves are entirely reflected at the surface. In this case a small and discrete number of harmonics survive the multiple interferences between the turning point and the surface, these harmonics dominate the lee waves field. Although the cases where the variations of S(z) are due variations in N(z) 2 are the most studied [START_REF]Dynamically-driven winds[END_REF]Jackson, Mayr and Vosper, [START_REF] Markowski | Mesoscale meteorology in midlatitudes[END_REF]]), it is clear from Eq. ( 1) that the variations in S(z) can also be due to the variation of the wind with altitude (examples of trapped lee waves due to increase of the jet stream with altitude in the low troposphere can be found in [START_REF] Vosper | A climatology of lee waves over the uk derived from model forecasts[END_REF], [START_REF] Teixeira | Drag produced by trapped lee waves and propagating mountain waves in a two-layer atmosphere[END_REF]]). Nevertheless and still according to Eq. ( 1), the variation of N and U should not be treated on the same footing.

To illustrate the difference we can recall that when U(z) becomes very small and even changes sign at a given altitude, S(z) becomes infinite. We are in the presence of a critical level where the waves are often absorbed in the linear case [START_REF] Booker | The critical layer for internal gravity waves in a shear flow[END_REF], ], rather than being entirely reflected as occurs at turning points. The differences become even more important in the nonlinear cases, the mountain waves systematically breaking near that level, the resulting gravity wave mean flow interactions sometimes yielding mountain waves amplification, downslope winds and foehn [START_REF] Clark | Critical level reflection and the resonant growth of nonlinear mountain waves[END_REF], [START_REF] Nappo | Wave exchange between the ground surface and a boundary-layer critical level[END_REF], ].

If we now return to the effect of the surface on the lee waves, a major limitation of the inviscid theory is that it neglects the boundary layer effect. For this reason, many numerical studies have included the effect of the boundary layer and shown that the surface friction and the treatment of the turbulent dissipation affect mountain waves ( [START_REF] Bougeault | Parameterization of orography-induced turbulence in a mesobeta-scale model[END_REF]], [START_REF] Richard | The role of surface friction in downslope windstorms[END_REF]]. More recently [START_REF] Smith | A theory of gravity wave absorption by a boundary layer[END_REF]], [START_REF] Hills | The dissipation of trapped lee waves. part ii: The relative importance of the boundary layer and the stratosphere[END_REF]], and [Teixeira (2017)] returned to the basic mechanisms and included some boundary layer absorption in the theory. More specifically, [START_REF] Smith | A theory of gravity wave absorption by a boundary layer[END_REF]] shows that flows that are more unstable near the surface absorb less gravity waves and favor the onset of lee waves. They also noticed that in the presence of a boundary layer the incident flow near the surface can become very small, causing a near critical level situation where the waves can be efficiently absorbed, if the theory of [START_REF] Booker | The critical layer for internal gravity waves in a shear flow[END_REF]] where the critical level is located within the flow applies.

Realizing that the absorption of gravity waves (GWs) by critical levels located at the surface had never been studied, [Lott(2007)] (hereinafter L07) solved this problem theoretically in the viscous case. This paper shows that the wave absorption at the "surface" critical level can be substantial and increases when the surface flow stability increases. This paper also shows that pure reflection only occurs in the inviscid limit and when the surface Richardson number, J = N 2 /( ∂u ∂z ) 2 < 0.25. To translate how this result affect trapped lee-waves, L07 applied these results to the analysis of the stationary disturbances that can exist in flows with constant stratification and wind that varies like a tanh profile above the surface. This profile was chosen because it has zero wind at the surface, near constant shear immediately aloft and a very smooth transition to constant wind in the far field: it is stationary enough in the viscous case to make the analysis of stationary solutions relevant (at least for reasonably large Reynolds numbers). In such profile, L07 confirmed that GWs absorption increases with J and also that pure trapped modes exist in the inviscid limit when J < 0.25: these neutral modes correspond to neutral modes of KH instability found by Drazin (1958) for the same profile in the unbounded case. Latter, these results were confirmed in [Lott(2016)], (hereinafter L16), where an explicit mountain forcing is introduced. The result in L16 also reveals consequences of the near surface critical level dynamics that could not be anticipated from the analysis in L07. It shows that in the presence of a near surface critical levels, the mountain forcing easily triggers strong downslope winds and foehn. Interestingly, downslope windstorms in this case do not result from upper level wave breaking, as often mentions the literature [Durran(1990), ]. In this paper we will question how these downslope winds will modify the background flow in which the trapped waves develop.

One of the more severe limitation of L16 is that it considers free-slip boundary conditions but forces the incident flow to satisfy a no-slip boundary condition far upstream: it does not use the viscous solutions presented in L07 or solutions affected by a turbulent boundary layer as in [START_REF] Belcher | Form and wave drag due to stably stratified turbulent flow over low ridges[END_REF]]. Although these simplifications are essentially made to simplify the theory they are in part justified by the fact that in the atmosphere the low level wind shears are often due to horizontal gradients in temperature, they do not solely result from boundary layer dynamics: they can extend well above the boundary layer and there, the inviscid dynamics can be applied.

The first purpose of this paper is to describe further the dynamical nature of the trapped lee waves with surface critical level found in L16 and by comparison with a more commonly accepted theory. The second is to analyze if these trapped lee waves still occur in the nonlinear context and with realistic boundary layers. For the first purpose we will compare the nature of trapped modes due to slow low level winds with and without surface critical levels using the L16 the-ory. For the second, we will follow [START_REF] Damiens | An adiabatic foehn mechanism[END_REF] Damiens, Lott, Millet and Plougonven] (hereinafter D18), who (i) adapted the theory to higher mountains than in L16, (ii) used in the theory boundary layer depths that gives result comparable to WRF, and (iii) used the WRF model in various idealized configurations and boundary layer specifications. The plan of this paper is as follows. In section 2, the formalism in L16 is extended to higher mountains, it is also applied to flows with and without surface critical levels. This will highlight how dynamically different can be situations where the changes in the background wind only affect the Scorer parameter, to the situation where they both affect the Scorer parameter and induce a near surface critical level. Section 3 then compares the theory with surface critical level and the WRF simulations for small and medium height mountains respectively.

Theory 2.1 Continuous variations in U (z)

To analyze the mountain waves produced by a stably stratified shear flow when the incident wind is null at z = 0, we follow L16 and consider the background flow profiles

U(z) = U ∞ tanh(z/d), N 2 (z) = const, (2)
incident on a 2-dimensional mountain modeled by the Witch of Agnesi profile:

h(x) = H 1 + x 2 2L 2 .
(3)

In (2), U(z) is the background horizontal wind, N(z) the Brunt Vaisala frequency, z the altitude, d is the vertical scale of the shear and U ∞ the incident wind maximum amplitude. In (3) H is the maximum mountain height, L its characteristic horizontal length, and x the horizontal coordinate. When scaling time by N -1 and distances by U ∞ /N the 2-dimensional non-rotating linear dynamics can be expressed in term of a non dimensional vertical velocity, w(x, z) of the form,

w(x, z) = +∞ -∞ f (k) ŵc (k, z)e ikx dk, (4)
where overbars denote dimensionless variables, and ŵc (k, z) is a canonical monochromatic solution of "unit" amplitude in the far field which satisfies the dissipative Taylor Goldstein equation,

d 2 ŵc dz 2 +    1 U -i z k √ J 2 - U zz U -i z k √ J -k 2    ŵc = 0. (5)
In (4) the amplitude term f (k) is obtained by numerical inversion of the non-linear free-slip boundary condition:

w(x, h(x)) = U(h) + u(x, h(x)) dh(x) dx , (6)
where u(x, z) is the horizontal wind disturbance. In ( 5)-( 6)

U(z) = tanh z √ J , and
h(x) = H N 1 + x 2 /2F 2 r (7) where J = N 2 d 2 U 2 ∞ , H N = HN U ∞ , and F r = LN U ∞ (8)
are the surface and minimum Richardson number, the non-dimensional mountain height, and the Froude number respectively. Still in (5) the dissipative vertical scale for each harmonics,

z k = z b 2 1 F r k + F r k . ( 9
)
has been expressed in terms of a global scale z b , with the first term in parenthesis resulting from Rayleigh drag and Newtonian cooling of coefficient z b /F r / √ J/2 and the second term from a viscous dissipation which acts in the x-direction only of coefficient z b F r / √ J/2 (for a more thorough discussion of these terms see discussion of Eq. 4 in L16). They are introduced to regularize the critical level dynamics for all the harmonics, i.e. the longer and the shorter ones respectively. Note that z k is written differently than in L16 to make clear that when the formula is applied to the dominant wave number F r -1 , z k in ( 9) is comparable to the dissipative scale z b 1 . To construct ŵc , L16 uses for each k an exact solution of the inviscid version of (5), ŵinv (k, z), which is based on hypergeometric functions (see [START_REF] Lott | A transition from kelvin-helmholtz instabilities to propagating wave instabilities[END_REF]]) and (34) in L16, and which asymptotic behaviors are

ŵinv (k, z >> 1) ≈ e -m z , (10) ŵinv
(k, z << 1) ≈ ŵmat (k, z) = a 1 (k)z 1/2-iµ + a 2 (k)z 1/2+iµ (11)
where,

m = |k 2 -1|, µ = |J - 1 4 |, (12)
and where a 1 (k) and a 2 (k) have analytical forms once the "unit" amplitude condition (10) is satisfied. Also, when k 2 < 1, m is changed in -isign(k)m, where the sign is to ensure upward group speed, and when J < 0.25, µ is changed in iµ. Near the surface, L16 also uses the asymptotic solution of the damped Taylor-Goldstein Equation ( 5) valid when z << 1:

ŵsrf (k, z) = a 1 (k)(z -iz k ) 1/2-iµ + a 2 (k)(z -iz k ) 1/2+iµ , (13)
where a 1 and a 2 are the same as in (11) to ensure convergence toward the "matching" function ŵmat given in ( 11) and when z → ∞. Then, following [START_REF] Bender | Advanced mathematical methods for scientists and engineers[END_REF]] we will approximate ŵc by its uniform estimate:

ŵc ≈ ŵinv + ŵsrf -ŵmat . ( 14
)
The only difference with L16, is that we use this approximation in the inversion of (6),

+∞ -∞ f (k) ŵc (k, h) -ûc (k, h) dh dx e ikx dk = U dh dx (15)
rather than using the surface solution in (13) to express ŵc and ûc in the terms between brackets in (15). This formally permits to consider more elevated ridges than in L16.

Discontinuous variations in U (z)

To contrast our results to the more common case where trapped lee waves are only due to changes in the Scorer parameter and where the ground is not a critical level we also apply our theory to a low level shear zone represented by a two layers flow,

U (z < √ J) = U 1 = 0.5 and U(z > √ J) = U 2 = 1. ( 16
)
This profile is given in the non-dimensional formalism of subsection (1.1), see also Fig. 1a. For this profile the dissipative Taylor Goldstein equation takes the form,

d 2 ŵc dz 2 + 1 (U j -iz k ) 2 -k 2 ŵc = 0 ( 17
)
where the j = 1, 2 index indicates the layer of interest. Also, the continuity of pressure and the uniqueness of the material displacement at the altitude z = √ J require that

(U j -iz k ) d ŵc dz j=2 j=1
= 0, and ŵc

U j j=2 j=1 = 0 (18)
respectively. If we then define vertical "wavenumbers" as

m j = ǫ k 2 - 1 (U j -iz k ) 2 , ( 19
)
where the sign ǫ = ±1 is chosen so that ℜ {m j } > 0, the canonical solution writes, ŵc = e -m 2 (z-

√ J) and ŵc = 1 4 + m 2 m 1 e -m 1 (z- √ J) + 1 4 - m 2 m 1 e +m 1 (z- √ J) , (20)
for z > √ J and z < √ J, respectively (see also [START_REF] Teixeira | Drag produced by trapped lee waves and propagating mountain waves in a two-layer atmosphere[END_REF]]).

Inviscid result

To show the fundamental dynamical differences between the two profiles we next search the neutral solutions that exist in both the inviscid case and when the surface is flat. When the variations in U(z) are continuous, L07 shows that the dynamics is intimately related to the absorptive properties of the near-surface critical level. To expose them in a concise way we can use the near surface inviscid solution in (11) and express the wave momentum flux

F z = ûc ŵ * c + û * c ŵc 4 = µ 2k (|a 1 | 2 -|a 2 | 2 ) when J > 0.25, µ 2k (a 1 a * 2 + a * 1 a 2 ) when J < 0.25. ( 21
)
L07 showed that for fixed a 1 , the inviscid solution (11) can be matched to dissipative solutions that satisfy the surface boundary condition, a matching that determines a 2 . It also shows that |a 2 | < |a 1 | when J > 0.25, which implies F z > 0 according to (21): the upward waves dominate the downward propagating ones. The pure surface reflection which is at the basis of the classical trapped lee wave theory [Scorer(1949), ] cannot be satisfied. L07 also shows that the ratio |a 2 /a 1 | tend to decrease when J increases, which in lee waves theory translate into the fact that the wave field decay more rapidly with downstream distance when J increases [START_REF] Smith | A theory of gravity wave absorption by a boundary layer[END_REF], ]. L07 also shows that these results still hold in the inviscid limit. When J < 0.25 L07 also shows that for fixed |a 1 | one can have |a 2 | = 0 in the inviscid limit. In this case F z = 0 according to (21) which means that the reflection is total and trapped lee waves can exist. To find them, the near surface analysis in (11) is nevertheless not sufficient. We have to find if there are indeed solutions of the global inviscid problem that satisfies a 2 = 0 near z = 0. L07 also shows that such solution do exist, they correspond to half of the branch of the neutral modes of KH instabilities found by [Drazin(1958)] (see his Eq. 25). In our definition of the non dimensional parameters it corresponds to the modes for which

k 2 = 1 - √ 1 -4J 2J , (22)
This dispersion curve is shown in Fig. 1b. When the variations in U(z) are discontinuous, the dynamics is almost as in [START_REF] Teixeira | Drag produced by trapped lee waves and propagating mountain waves in a two-layer atmosphere[END_REF]], and the onset of trapped modes correspond to wave numbers for which ŵc (z = 0) = 0, a condition that imposes,

4 tanh(m 1 √ J) = - m 1 m 2 . ( 23
)
In the inviscid case, the only wave numbers that satisfy this condition are those for which 1 < k < 2 and when:

4 tan √ J 4 -k 2 = - 4 -k 2 k -1 . (24)
The corresponding dispersion curves are shown in Fig. 1b. It shows that pure trapped waves occur for much "deeper" low level shear layer, the onset of trapped waves being now conditional to J > π 2 /12 ≈ 0.822 rather than J < 0.25 in the continuous case. We also see that when J increases the number of trapped modes increases, this number being given by the resolution of (24) when k = 1, e.g. by the number n satisfying

π 2 + nπ < √ 3J < 3π 2 + nπ. (25)

Forced dissipative results

To include an explicit forcing, we next proceed numerically as in L16, discretize in the spectral space the Fourier integral in (4), taking for ŵc (14) in the continuous case or (20) in the stairway case. In both, cases, f (k) is determined through the numerical inversion of the nonlinear boundary condition (15). In the discrete space this inversion is a matrix inversion that only converges with increasing resolution when there is dissipation and because the surface boundary condition becomes singular in the inviscid case when there are resonant modes. In this section, we take a domain of length 500, spanned by 2048 equally spaced points, the boundary layer parameter z b = 0.01, and the non-dimensional mountain height H N = 0.1. We also take F r = 2, which means that waves with wavenumbers around F r -1 = 0.5 will be forced. This value is sufficiently near the lower bound of the dispersion curves in Fig. 1b to ensure that some trapped waves will be well excited. The two panels in Figs. 1c) and 1d) show the vertical velocities around the center of the shear layer in z = √ J/2, for the continuous and discontinuous profiles respectively. As expected from L16, one sees that in the presence of a surface critical level, trapped lee waves develop more easily when the flow is unstable, i.e. for small J's, whereas it is the other way round when the background flow is discontinuous. In the later case, trapped waves develop more easily when the parameter J becomes larger than one, consistent with the inviscid results for the stairway case in section 2.3). These results corroborate that reasoning in terms of trapping only can be misleading. If the surface critical level dynamics is neglected, increasing the depth of the shear increases the depth of the low level wave duct: more and more modes can be trapped within it. When including reasoning in terms of shear flow stability, the surface critical level dynamics becomes central, and, when the depth of the shear zone decreases, the surface critical level absorption decreases, which favors trapped lee waves.

3 Fully nonlinear simulations

Experimental setup

To test if the relation between flow stability and trapped mountain waves still hold in the fully non linear context, we next follow D18 and use the WRF model in the 2D mountain flow configuration available on line [START_REF] Skamarock | A description of the advanced research wrf version 2[END_REF], ]. In all the simulations presented we consider a mountain of length L = 1.1 km which height typically varies between 50 m and 500 m. To capture well the gravity waves dynamics the model domain is 80 km long and 9 km high, and both directions are discretized by 800 x 450 equally spaced points respectively. These yield an horizontal resolution of ∆x = 100 m and a vertical resolution of ∆z = 20 m, both are near an order of magnitude smaller than the mountain dimension. These rather high resolutions force to use a time step of ∆t = 0.5 s, they are necessary to guarantee that the near surface critical level dynamics is well captured, and we verified that our results are not much changed when these resolutions and domain sizes are changed by a factor 2. To allow gravity waves to propagate away from the domain without lateral and upper boundary reflections we use open lateral boundary conditions and introduce a 6km high damping layer below the domain lid.

To impose a constant buoyancy frequency, and because the WRF model is fully compressible we take an isothermal atmosphere with T 0 = 288 K. This results in a constant Brunt-Vaisala frequency N, with N 2 = g θ 0 dθ 0 dz = κg 2 RT 0 = 3.32 10 -4 s -2 , and where κ, R and g have their usual earth values, and θ 0 is the background potential temperature. The background wind is defined by ( 7) with U ∞ = 10 m s -1 . In this setup the froude number value F r = 2, as in section 2 and will not vary, the sensitivity to the surface Richardson number J will be analyzed by changing the value of the shear layer depth d, and that to the mountain height by varying the mountain height H N . In most case, and except when specified the results are shown after 4 h of integration, that is much longer than the advective time-scale L U∞ ≈ 100 s, and sufficiently long for the trapped lee waves to develop downstream.

As the theory in Section 2 is based on an free-slip treatment of the lower boundary, we will first present WRF simulations where the constant diffusion coefficients are put to zero and where the surface boundary condition is free-slip. Nevertheless, this does not necessarily make the simulations inviscid because irreversible or diabatic processes certainly occur near the surface : it happens that the WRF dynamical core is stable and dissipative enough to handle these processes (see discussion in D18). When we return to the theory, it is a priory difficult to say how these irreversible processes can affect the wave dynamics. Nevertheless, as the theory has a free parameter to represent dissipation, the boundary layer scale z b , we can vary it and try to identify a value for which there is a match between the nonlinear simulations and the theory. In D18, it was found that a dependance in √ J helped the comparison with WRF and we take here

z b / √ J = 0.005. ( 26
)
We well keep this value in all the following experiments with the theoretical model.

To extent even further the significance of our results, we have also used WRF with two boundary layer schemes that have been extremely validated by the community, the Yonsei University Planetary Boundary Layer scheme (hereinafter YSU BLYR, [START_REF] Hong | A new vertical diffusion package with an explicit treatment of entrainment processes[END_REF]]) and the Mellor Yamada scheme (hereinafter MY BLYR, [START_REF] Mellor | Development of a turbulence closure model for geophysical fluid problems[END_REF]]). Although the two schemes have fundamental differences, the YSU BLYR is based on the nonlocal diffusion concept of [START_REF] Troen | A simple model of the atmospheric boundary layer; sensitivity to surface evaporation[END_REF]] whereas the MY calculate local diffusion coefficients with amplitude controlled by the amplitude of a turbulent kinetic energy evaluated at each time step, we do not expect them to behave very differently in all the dry and initially convectively stable cases we consider. The intention here is more to test robustness using well known schemes, and although the two schemes may well regularize the convectively unstable regions produced by the waves quite differently. In these two configurations where the interaction with the surface is more thoroughly taken into account we force the potential temperature of the surface to equals the background potential temperature, θ s (x) = T 0 e κh(x)/H 0 , where the characteristic height H 0 = RT 0 /g, whereas the roughness length has its default value z 0 = 0.1 m.

Low mountain

The left panels in Fig. 2 show the vertical velocity fields produced by the theoretical model for three different values of J, and the right panels show the same fields evaluated with the free-slip WRF simulations. In all panels, the results are shown using non-dimensional variables to ease comparison with theory and the mountain maximum height H N = 0.1. The first striking result is that the theory and WRF compare very well qualitatively and quantitatively. This is maybe the most important result of this work: the theory is based on a linear approximation that is often questionned in the presence of critical levels [START_REF] Lott | A transition from kelvin-helmholtz instabilities to propagating wave instabilities[END_REF], [START_REF] Dörnbrack | A note on the application of linear wave theory at a critical level[END_REF]].

If we return to the fields themselves, we see that for J = 3 in Figs. 2a and2b, the vertical velocity field is largely dominated by vertically propagating mountain waves, the sign of the constant phase lines slope being opposite to that of the incident wind. In both models, the wave field substantially extents downstream of the ridge, which is characteristic of non-hydrostatic effect when the background flow is uniform [Queney(1947), ] and the Froude number near F r = 1. Here the interpretation holds above the shear layer, i.e. when z > √ J. Still for J = 3, note that in WRF, the wave amplitude is a little smaller than in the theory, but this difference could be sorted out by moderately increasing the boundary layer depth z b in the latter (not shown).

If we now look at smaller values of J, we see that for J = 0.75 the theory and WRF in Figs. 2c and 2d still predict comparable fields, trapped lee waves start to appear, although they are still quite small in both models. Note that now, it is the theory that predicts smaller amplitude wave field than WRF, which could suggest that we should now decrease the boundary layer depth z b in the theory. As this amplitude factor changes from one experiment to another this indicates that the relation between z b and J we have taken here and in D18, should eventually be improved or made dependent on the Froude number or on H N . Nevertheless as this does not affect our interpretation of the results we prefer to keep the relation (26).

In both models, and for even smaller values of J, (i.e, J = 0.16) we see again that WRF and theory give comparable results, but now both present a well defined train of lee waves, extending downstream down to more than 50 km. These trains of lee waves are characterized by almost vertical phase lines, the vertical velocity fields being essentially confined within the shear layer (see U(z) in the left part of each panel). A difference is that in the theory, the train of lee waves seems to dissipate faster than in WRF, again this can be easily sorted out by decreasing further z b in the theory, an issue that is not critical for our interpretation of the results. A more interesting aspect concern the horizontal wavelength. We see that in both models it is near equals λ ≈ 5, which corresponds to k ≈ 1.05, whereas the corresponding KH mode in ( 22) has k ≈ 1.10 for J = 0.16. If we remember that the Drazin solutions are confined to the shear layer and have only one node for the vertical velocity, the inflection point of the background wind (here z = 0), it is clear that the wave field downstream in Figs. 2e and2f has structure of a neutral mode of KH instability.

As the central objective of this paper is focused on the trapped lee-waves development, we next choose a representation that emphasize them. For this purpose, the Fig. 3 shows the horizontal variations of the vertical velocity at the altitude z = 1. The panels in the left column are for the theory (Fig. 3a), and are consistent with the results in L16, which establishes that more unstable flows favor trapped lee waves. In L16 nevertheless dissipations are almost absent, which makes that the lee waves for J < 0.25 are almost pure trapped lee waves, their amplitude do not decay downstream, the mountain exciting a neutral mode of KH instability. Also, in L16, the boundary layer depth does not increase with J, which here absorbs even further the waves when J increase: this makes that in Fig. 3a) almost no waves are present when J > 1.

A striking result is that WRF reproduces well this dependence with J, in the free-slip case which we already discussed but confirm here more systematically. What is more interesting is that with more sophisticated boundary layer parameterization schemes in Figs. 3c and3d the transition from more unstable flows with more trapped waves to more stable flows with less trapped waves still hold. As said in the model description, our choice was to use well tuned and stable schemes to confirm our findings beyond the free slip assumption, rather than to discuss with details the differences. What can be said here is that imposing a surface drag via similarity theory, what fundamentally do both the MY-BLYR and the YSU-BLYR schemes, does not fundamentally change the absorptive properties of the waves near the surface compared to the inviscid case, and hence the downstream development of the trapped lee-waves. A reason is probably that our background winds are small near the surface, and the roughness length is also quite small, two reasons that make the frictional drag small and the results quite consistent with the free-slip cases.

Medium height mountain

To test if some of our results hold for higher mountains, we next present experiments with H N = 0.7, that is H = 404 m in dimensional units. The results for the horizontal profiles of the vertical velocity in Fig 4 confirm that some aspects of the theory are still relevant. The most remarkable one is that trapped waves do occur in most experiments, they all have horizontal wavenumber around k = 1, their spatial structure also resemble to KH neutral modes (not shown). Although there is an overall tendency for the simulations with small J to produce trapped waves propagating over longer distances downstream, some remarkable differences start to occur.

The first is that for small J and except when the MY BLYR scheme is used, the first minimum in vertical velocity located immediately downstream of the ridge becomes substantially larger in amplitude compared to the oscillations associated with the lee wave. This follows that enhanced non-linearities increase downslope winds, a behavior already analyzed in D18. We will not discuss it further here, essentially because this effect is present in the theory in Fig. 4 for J = 0.16, J = 0.23, and J = 0.3, as well as in the free slip case, and to a lesser extent when using the YSU BLYR. Note that this enhancement of downslope winds is almost absent with MY-BLYR, indicating that turbulent parameterizations affect seriously the downslope winds. As we shall see with more details below this impacts in return the development of the trapped lee waves: at small J and with MY-BLYR in Fig. 4d), the downslope winds are not as intense as with the other BLYR, which seems to favor the downstream development of the lee waves.

The second difference is more directly related to the lee wave downstream extension. In some cases it can almost disappear locally at small J (for instance in the free-slip case in Fig. 4b when J = 0.16 between 1 < x < 30). Conversely lee waves can develop downstream more substantially than in the theory, as occur quite systematically in WRF when J > 1. To test if these changes challenge our interpretation of trapped lee waves in terms of stability, we next analyze with more details the cases with the two extreme values J = 0.16 and J = 5.

The results above suggest that when H N approaches 1 the interpretation in term of upstream flow surface stability becomes incomplete. Interestingly for such values of H N , we also know that strong downslope winds and low level wave breaking can occur. These certainly change the large-scale flow downstream and hence the local properties of the flow in which the trapped waves develop. To substantiate this relation between local flow stability and trapped waves, we next use the fact that the lee wave fields shown in Fig. 4 have rather well defined horizontal wavelength: we can make a separation between the "large-scale flow" and the waves by averaging fields over horizontal distances that compare with the horizontal wavelength. In the following, as the trapped waves have wavelength that are always near and below 3 km, we will apply a top hat horizontal filter with a fixed length of 3 km to extract the large scale, but we found little sensitivity when moderately increasing the horizontal length of the filter.

For J = 0.16 the top panels in Fig. 5 show the temporal evolution of the altitude of the "unfiltered" isentrope θ = 291 K, and which upstream altitude is almost equal to H. The time evolution in the free slip case in Fig. 5a shows the development of a moderate Foehn downslope, the region where the isentrope altitude is below the mountain top extents with time to reach near x = 20 km after t = 240 min. This region of moderate Foehn does not end by an abrupt hydraulic jump, consistent with the fact that strong Foehn are favored when J is large (see L16, D18). In term of trapped waves, what is interesting is that the oscillations with wavelength below 3 km are stationary, as expected for mountain waves, but begins downstream the Foehn region. The panel below in Fig. 5d shows the raw fields of potential temperature altitude at t = 240 min, superimposed onto the "large scale" Richardson number,

R if = g θ f ∂θ f ∂z / ∂u f ∂z 2 (27)
where θ f and u f are the filtered field of potential temperature and horizontal wind respectively.

Here we see that the region where the trapped waves are present are also characterized by regions where R if is small near the surface, and which are quite far downstream of the ridge.

In Figs. 5b one sees that an effect of the YSU-BLYR layer parameterization is to reduce the Foehn intensity, the descent of the θ = 291 K surface downstream is less pronounced than in the inviscid case, and the flow near the surface still has R if < 0.25, as indicates the green curves and at least until the x = 20 km. Above this region substantial small scale oscillations do occur, whereas beyond x = 20 km, the surface flow is stabilized and the isentropes no longer present small scale oscillations. This correspondence between surface flow stability and trapped waves is also found when the MY BLYR scheme is used (see , in this case also the scheme reduces the foehn amplitude, the flow stays unstable on the lee-side of the ridge, it even stays so down to at least x = 40 km and trapped waves are present almost everywhere in this zone.

In all cases, we see that the correspondence between trapped waves and the near surface flow stability continues to hold.

The differences between the theory and the WRF are even more pronounced when the initial background flow is stable as show the simulations with J = 5 in Fig. 6. First, in the free slip case, the θ = 291 K surface in Fig 5a presents a pronounced foehn immediately downstream of the ridge, followed by an abrupt transition beyond which trapped lee-waves develop. Again the trapped lee waves field develops in time, but is stationary: the first pronounced crest at around x = 10 km for instance stays almost at the same place, and this is true for the other crests. The fields of large scale flow stability and isentropic surface in Fig 5d show that on the lee side, the dynamics has not much modified the surface flow stability (the Richardson number near z = 0 is still much larger than 1), but produced in the mountain wake a large unstable zone that extent with time. We see in Fig. 5d that this unstable zone is also where the trapped wave signal is the most pronounced. According to this result, it seems reasonable to suggest that the trapped waves are again related to an instability triggered by the mountain, as suggests the stationarity of the lee wave.

To support this, one needs to verify that where the flow is unstable the large-scale wind is also small, and following that (i) KH instabilities have a phase speed near the large-scale wind value in the region where they are produced, and (ii) mountains force stationnary waves. This is verified in Fig. 7a, which shows profiles from filtered filelds at x = 13 km. In it, we see that in the region of instability, the large scale wind is quite small, the near unstable modes, if they exist, are likely to have small phase speed and to be excited by an orographic forcing. The results from other boundary layer parameterizations in Figs. 6b-6e and Figs. 6c-6f somehow corroborate these results. In all of them, trapped waves develop whereas the surface stability stays strong. They seem associated to the unstable regions located well above the surface in the wake of the mountain, and in all cases the large scale winds in the region where the Richardson number is small are quite small (Figs. 7b and7c). A marked difference is that between the free slip case and the cases with BLYR parameterization the horizontal wavelength of the waves is quite different, an other is that the intensity of the foehn is much less pronounced with the BLYR parameterizations than in the free slip case.

Conclusion

Using theory and WRF simulations with various boundary layer parameterization, this paper has shown that the onset of trapped lee waves is strongly linked to the stability property of the background flow. This complement our conventional view of trapped lee waves where the dynamics is mainly explained in term of variations of the Scorer parameter. We believe that this view becomes mandatory when the variations in the Scorer parameter are due to variations of the incident wind rather than of stratification, and when the incident wind is quite small near the surface. For small mountains, it happens that the theory in L16 explains well this dependence. Small surface Richardson numbers result in large gravity wave surface reflections and favor the trapped wave development. For medium height mountain, the nonlinear dynamics in the wake of the mountain modify the large scale flow, and therefore the background properties in which the trapped waves develop. For small J, non-linearities yield in an increase of the surface Richardson number near the surface and the trapped waves are attenuated. For large J, strong downslope windstorms happen and yield mixing that destabilizes the flow in the wake of the mountain. In this case, we attribute the onset of trapped waves to the triggering of small phase speed unstable modes.

Despite the fact that there is some coherency in our results concerning the links between trapped waves and instabilities, it is important to note that the boundary layer parameterizations have strong quantitative impacts. These impacts are extremely evident in the medium height mountain simulation (section 3.3). In this case and when the background flow is unstable near the surface (J = 0.16), the dynamics stabilizes the surface flow in the free slip case, whereas it does not affect it much when the MY-BLYR scheme is used. This impacts the location of the trapped waves. For very stable cases J = 5, the dynamics produce strong foehn in the free-slip case, whereas the Foehn is less intense with the YSU and MY boundary layer parameterizations. Again the location of the trapped waves is affected in return. These results clearly illustrate that the boundary layer parameterizations have strong impacts on mountain wave dynamics, at least when the incident wind present strong shears in the low troposphere. In essence, the results obtained illustrate again the significance of low level wind shears and stability on mountain flow dynamics, something often noticed in the context of trapped lee waves [START_REF] Reinecke | Initial-condition sensitivities and the predictability of downslope winds[END_REF], [START_REF] Georgelin | On the transfer of momentum by trapped lee waves. case of the iop3 of pyrex[END_REF]]. In this context it is interesting to note that there are many field experiments where upstream soundings show low level shears [START_REF] Doyle | An intercomparison of t-rex mountain-wave simulations and implications for mesoscale predictability[END_REF], [START_REF] Sheridan | Influence of lee waves on the near surface flow downwind of the pennines[END_REF]] generally related to advancing fronts [START_REF] Lothon | Experimental study of five foehn events during the mesoscale alpine programme: From synoptic scale to turbulence[END_REF]], when the low level shear is associated with thermal wind balance.

Figure 1 :Figure 2 :Figure 3 :Figure 4 :Figure 5 :Figure 7 :

 123457 Figure 1: (a) Mountain profile and background wind configurations used in the theory: U in tanh (black solid), U stairway (grey dashed); (b) Dispersion curves (22) and (24) (same color convention); c) vertical velocity at z = √ J/2 from the dissipative model when U is in tanh: d) Same as c) but when U is stairway.

For consistency, note also that the boundary layer depth z B used in L16 is related to the one here by z B = 5z b .

Acknowledgements:

This work was supported by the ANR/JPI-Climate/Belmont Forum project GOTHAM (ANR-15-JCLI-0004-01).

acknowledgements

This publication was partly funded by the ANR/JPI-Climate/Belmont Forum project GOTHAM (ANR-15-JCLI-0004-01).