The sensitivity of probabilistic convective‐scale forecasts of an extratropical cyclone to atmosphere–ocean–wave coupling

The benefits of dynamical atmosphere–ocean–wave coupling in probabilistic weather forecasts generated using convective‐scale ensemble prediction systems are to date unknown. We investigate the respective impacts of atmosphere–ocean–wave coupling, and initial condition (IC), lateral boundary condition (LBC), and stochastic physics perturbations within a convective‐scale ensemble coupled system for an extratropical cyclone case study. Towards this aim, we developed the first 18‐member, 2.2 km grid spacing ensemble regional coupled system (Ensemble‐RCS) with domain covering the British Isles and surrounding seas. Ensemble‐RCS coupled and uncoupled simulations of cyclone Ciara (February 2020) were performed. Adding stochastic perturbations to the model physics parametrizations enhances the ensemble spread of the uncoupled atmosphere‐only ensemble driven by IC and LBC perturbations, while slightly reducing (by up to 0.5 m · s −1 ) the median of the ensemble 95th percentile 10‐m wind speeds from its value of about 24 m · s −1 at peak time. A substantial proportion of this impact is attributable to Charnock parameter perturbations alone. By coupling the atmosphere‐only ensemble, with stochastic physics, to the ocean, the ensemble median and spread is mainly unaffected. However, additional coupling to waves reduces the median wind speed by 1 m · s −1 , which leads to reductions of up to 70% in strong wind strike probability, and halving of the spatial coverage of high values (>50%) of this probability. Finally, we demonstrate the usefulness of two metrics originally developed for precipitation verification – the neighbourhood‐based Fraction Skill Score (FSS) and the object‐based Structure, Amplitude, Location (SAL) – for examining the spread in convective‐scale ensemble forecasts. It is concluded that coupling has a consistent impact across the ensemble members. Remarkably, the impact of coupling to waves is found to be comparable in size to that of adding IC, LBC and stochastic physics perturbations to the uncoupled atmosphere‐only ensemble simulation, implying that the dynamical coupling to ocean and sea‐state are important aspects of model uncertainty.

stochastic physics perturbations to the uncoupled atmosphere-only ensemble simulation, implying that the dynamical coupling to ocean and sea-state are important aspects of model uncertainty.

K E Y W O R D S
atmosphere-ocean-wave coupling, convective-scale ensemble forecasts, extratropical cyclones, FSS, initial condition error, SAL, stochastic physics perturbations, surface wind speeds INTRODUCTION Extratropical cyclones can generate intense surface wind speeds and gusts as they pass over the UK (e.g., Martínez-Alvarado et al., 2014), leading to large socio-economic impacts (e.g., Hewston and Dorling, 2011). Numerical weather prediction systems can predict the synoptic-scale evolution of cyclones with reasonable skill (Frame et al., 2015) but, because the atmosphere is a chaotic system (Lorenz, 1963), small errors in the initial conditions (ICs) of forecasts and in model physics approximations grow exponentially over time, limiting the predictability of the location and strength of mesoscale cyclone features such as intense winds (Buizza et al., 2005;Bowler et al., 2008). The forecast uncertainty yielded by ICs and model errors can be estimated by generating probabilistic realizations of events using ensemble prediction systems (EPSs; Buizza et al., 1999). Compared to IC errors, the correct estimation of model errors represents a much greater challenge (Buizza et al., 2005). Candidate sources of model errors include the sea-state-independent model parametrizations of air-sea momentum, heat and moisture fluxes that control, at convective scale, the location and magnitude of strong winds near the ocean surface (Janssen, 2004;Lewis et al., 2018;Lewis et al., 2019). Because of the steep vertical gradients of the air-sea fluxes at the air-sea interface, small model errors in the parametrization of air-sea fluxes can result in large forecast errors and model biases. A number of studies (Wahle et al., 2017;Lewis et al., 2019) have shown that the direct simulation of the effect of the dynamical ocean and wave state on air-sea surface exchange coefficients in a deterministic coupled system has a great potential to improve the prediction of cyclone wind speed forecasts at or close to convective scale. However, to date, a study based on the development of a dynamically coupled convective-scale EPS has not yet been published. Here we develop such a system and, by applying it to a cyclone case-study, determine the sensitivity of convective-scale ensemble forecasts of wind speeds to ocean and wave coupling and the relative sensitivity to coupling compared to sensitivity to ICs, lateral boundary conditions (LBCs) and stochastic model physics perturbations.
Medium-range EPSs have been used operationally for 30 years and are now capable of generating accurate and reliable forecasts of the synoptic-scale development of cyclones with several days lead time (Buizza and Palmer, 1995;Houtekamer et al., 1996;Molteni et al., 1996;Toth and Kalnay, 1997;Buizza et al., 2005;Bowler et al., 2008). Medium-range EPSs focus mainly on the uncertainties associated with synoptic-scale baroclinic instability (Bowler et al., 2008), but also attempt to represent the uncertainty at the convective scales of motion (Tennant et al., 2011) either by using multiple physics parametrization schemes, by using multi-parameter schemes within a single-model approach (Bowler et al., 2009;Berner et al., 2011), or by widening the sample range of possible forecast scenarios by using a multi-model approach (Hagedron et al., 2005). To represent the interaction between the uncertainty in the atmospheric and oceanic boundary layers, the European Centre for Medium-Range Weather Forecasts (ECMWF;Bauer et al., 2020) and the Naval Research Laboratory (NRL; Holt et al., 2011) have dynamically coupled the atmospheric dynamical core of their operational medium-range EPSs with their operational ocean and wave models. Also the Met Office plans to run an operational medium-range atmosphere-ocean coupled ensemble due to the promising results obtained by the corresponding global deterministic forecast system (Vellinga et al., 2020). Despite the increase in physics complexity and in the degree of coupling between model components, medium-range EPSs cannot exploit the full benefits of dynamical integration of the atmosphere, ocean, and wave feedbacks (Holt et al., 2011), as the relatively coarse resolutions cannot provide detailed km-scale forecasts over regions of particular interest. For example, in the case of a rapidly developing cyclone, Lean and Clark (2003) established that a grid spacing of (10) km is sufficient to represent the overall frontal structure of an extratropical cyclone and the associated airflow; however, a convective-scale grid spacing of (1) km is required to resolve smaller-scale multiple slantwise circulations, the 3D structure of convection lines, and the peak cyclone surface wind speed.
At convective scale, coupling atmosphere, ocean, and wave model components in deterministic prediction systems has shown potential to reduce the biases of cyclone wind speed deterministic simulations. For example, Wahle et al. (2017), simulating a historical series of cyclones that hit the southern North Sea, demonstrated that coupling the deterministic WAM wave model (Komen et al., 1994) with the deterministic COSMO regional atmosphere model (Rockel et al., 2008) (run at 5 and 10 km grid spacing, respectively) reduced wave heights and wind speeds by up to 8% and 3%, respectively, with associated equivalent improvements in the forecast skill of wind and waves. Comparable reductions in wind speeds and wave heights, associated with skill improvements of 20% and 5%, respectively, were achieved by Lewis et al. (2019) in coupled simulations of part of the 2014 UK cyclone season. These simulations were carried out using the deterministic Met Office Regional Coupled System (RCS), termed UKC3 in Lewis et al. (2019), which couples models of the atmosphere, land-surface, shelf-sea ocean, and ocean surface waves. A further study (Gentile et al., 2021) demonstrated that the improvements in cyclone wind speed forecasts generated by the RCS are mainly due to the impact of coupling to the wave model component. Bias reductions in wind speeds associated with severe weather conditions on coupling to a wave model have also been demonstrated for cyclonic events over the Mediterranean Sea (Katsafados et al., 2016;Varlas et al., 2017;Ricchi et al., 2017Ricchi et al., ,2019. Although experimental convective-scale coupled deterministic systems have undoubtedly shown benefits for the prediction of severe weather events, it has not yet been established whether the benefits from coupling are still retained when running as an ensemble, or what the size is of the impact of coupling relative to those of the perturbations applied to generate the ensemble. When convective-scale EPSs were designed about 15 years ago, researchers mainly focused on exploring the benefits of increasing the ensemble grid spacing for forecast skill improvements and of increasing the ensemble size for better sampling the uncertainty; at the time, there were not enough computational resources to additionally implement coupling of the atmospheric dynamical model core with the ocean and the wave models (Holt et al., 2011). For example, Schellander-Gorgas et al. (2017) compared the performance of the 16-member convective-scale 2.5 km horizontal grid spacing AROME-EPS to that of the 16-member mesoscale 11 km grid spacing ALADIN-LAEF EPS. For a summer case-study over the Alpine region, they found that the AROME-EPS significantly improved the precipitation forecast skill due to better resolution of the local effects of mountain topography. A further study by Raynaud and Bouttier (2017) varying resolution and ensemble size of convective-scale AROME-EPS found that increasing the resolution only led to short-range forecast improvements, whereas increasing the ensemble size exhibited a larger impact at longer forecast ranges. In this context, the Met Office developed the short-range high-resolution Met Office Global and Regional Ensemble Prediction System (MOGREPS; Bowler et al., 2008), consisting of two related EPSs: the medium-range global 12-member MOGREPS-G with 33 km grid-spacing and the short-range, limited-area (British Isles) 12-member MOGREPS-UK with 2.2 km grid-spacing (Hagelin et al., 2017). Since 2016 the MOGREPS-UK has been inititalised by perturbing about an analysis of the Met Office's UK variable resolution (UKV) deterministic model using the MOGREPS-G perturbed ensemble members (Hagelin et al., 2017). Investigating a three-month-long summer case-study, these authors found that the convective-scale ensemble MOGREPS-UK not only improved the precipitation forecast skill compared to the deterministic UKV but also the skill of other meteorological variables of interest such as 1.5 m air temperature, wind speeds, visibility, cloud-base height, and cloud cover. Further investigation by Porson et al. (2020) using a MOGREPS-UK configuration with 18 time-lagged members revealed substantial improvements in skill for the same meteorological variables considered in Hagelin et al. (2017). However, convective-scale EPS forecasts are more challenging to verify than coarser-scale ones because of the well known "double penalty" problem: that is, that a forecast with a slightly misplaced weather feature is penalised both for not having forecast the feature in the correct place and for forecasting the feature in a place where it did not occur. To avoid this problem, new metrics have been developed such as the neighbourhood-based fraction skill score (FSS; Roberts, 2008) and the object-based structure, amplitude, location (SAL; Wernli et al., 2008) score. Dey et al. (2014) and Zschenderlein et al. (2019) introduced and successfully used a dispersion form of FSS and SAL metrics to verify convective-scale EPS forecasts of precipitation events against radar observations.
Despite the many benefits of convective-scale EPSs compared to their medium-range counterparts described above, representation of model error at convective scale remains challenging because the smaller scales of motion simulated by convective-scale EPSs lead to a faster error growth (Clark et al., 2010;Baker et al., 2014;McCabe et al., 2016;Hagelin et al., 2017). Clark et al. (2021), investigating a series of precipitation events simulated by a convective-scale EPS, demonstrated that model error impact is greatest when it arises from parametrization of boundary-layer structure. Moreover, Flack et al. (2021) further demonstrated that uncertainty from stochastic perturbations applied to the model boundary-layer parametrization is statistically comparable to IC and LBC uncertainty. These findings clearly indicate that convective-scale EPSs need a more accurate representation of boundary-layer processes to reduce forecast error and better estimate model uncertainty of the meteorological variables that are sensitive to boundary-layer parametrizations. Although dynamical coupling of atmosphere, ocean, and wave models has been shown to improve the accuracy of boundary-layer parametrizations in deterministic simulations, its potential benefits for convective-scale EPSs have not yet been explored, with the exception of a single study from Bousquet et al. (2020), which investigated the performance of the ensemble configuration of the convective-scale atmosphere model AROME-IO (Indian Ocean) coupled with a 1D ocean mixed-layer model, for eleven major storms. As shown by Bousquet et al. (2020), coupling AROME-EPS to a 1D ocean mixed-layer model rather than to a dynamical 3D ocean model led to important limitations in capturing the simulated physics, such as not reproducing the outward horizontal transport of warm waters located near the tropical cyclone core as well as underestimating the storm-core sea-surface cooling. A further limitation of the ensemble 1D ocean-coupled AROME-IO of Bousquet et al. (2020), was that the model did not incorporate coupling to a wave model, despite atmosphere-wave feedbacks having proven equivalently important to atmosphere-ocean feedbacks for simulations of idealised cyclones (Doyle, 1995;2002), tropical cyclones (Holt et al., 2011), Mediterranean cyclones (Varlas et al., 2017;Ricchi et al., 2019), and extratropical cyclones affecting the UK and the North Sea (Wahle et al., 2017;Lewis et al.,, 2018;Gentile et al., 2021).
To address this question, we present the first study, to date, on the sensitivity of convective-scale EPS simulations to dynamical atmosphere-ocean-wave coupling relative to IC, LBC and stochastic model physics perturbations. Towards this aim, we have developed the first coupled convective-scale ensemble, termed Ensemble-RCS, focused on the British Isles domain and surrounding seas. This EPS is an ensemble implementation of the RCS for an atmosphere model domain with 2.2 km horizontal grid-spacing, matching that used in the MOGREPS-UK ensemble, with the atmosphere component uncertainty driven by MOGREPS-UK ICs, LBCs, and stochastic perturbations. Coupled (ocean and wave model), partially coupled (only ocean model) and uncoupled Ensemble-RCS simulations were run from 7 to 10 February 2020 during which time cyclone Ciara crossed the UK, the most intense cyclone since storm Tini (12 February 2014; Kendon, 2020). The mean, median, standard deviation and strike probability of the different coupled and uncoupled Ensemble-RCS forecasts are analyzed and compared. The ensemble spread characteristics are also evaluated using the dispersion FSS (dFSS) and dispersion SAL (dSAL) methods. Moreover, we evaluate the robustness of the Ensemble-RCS forecasts of Ciara's intense wind speeds on atmosphere-ocean-wave coupling and on perturbations of ICs, LBCs, and model parametrizations. The atmosphere, ocean, and wave feedbacks on which we focus are the air-sea turbulent momentum and heat exchange, mediated by the way the Charnock parameter and the sea-surface temperature (SST), respectively, are simulated and communicated to the atmosphere model.
The remainder of this article is organised as follows. Section 2 describes the weather conditions associated with the case-study cyclone Ciara. The Ensemble-RCS coupled and uncoupled configurations are described in Section 3, where the neighbourhood and object-based metrics used to characterize the ensemble spread are also explained in detail. Results on the relative sensitivity of stochastic perturbations to ICs and LBCs perturbations are presented and discussed in Section 4, followed by results of the sensitivity to coupling to ocean and coupling to waves and comparison of coupled and uncoupled Ensemble-RCS configuration spread characteristics. Conclusions are given in Section 5.

CASE-STUDY: STORM CIARA
Cyclone Ciara, officially named by the Met Office on 5 February 2020, was selected as a case-study because it was the most intense and damaging cyclone to have hit the British Isles since 2014 (Kendon, 2020). After originating from a weak area of low pressure in the southeastern US on 5 February 2020, cyclone Ciara entered the Met Office's high-resolution analysis (UKV) domain at 0000 UTC on 9 February 2020. It then separated from the left exit region of the jet stream (not shown), continuing to intensify until 0900 UTC on that day, when the central mean sea-level pressure (MSLP) reached a minimum of 950 hPa. At this time, as represented by the UKV analysis field in Figure 1b, the vertical pressure structure exhibited only a slight westward tilt from comparison of the MSLP field with that of 500 hPa geopotential height (Figure 1b), as expected given that the cyclone had already undergone most of its rapid intensification stage by this time. As can be seen from the Met Office analysis at 0600 UTC on 9 February 2020 in Figure 1a, Ciara has a complex frontal structure (with a possible frontal wave) and we therefore use information from the boundary-layer depth rather than a temperature field to characterise the warm and cold sectors of the cyclone. The map of boundary-layer depth, F I G U R E 1 (a) The Met Office surface analysis at 0600 UTC on 9 February 2020. UKV analysis maps at 0900 UTC on 9 February 2020 for (b) mean sea level pressure (white solid contours, interval 4 hPa), and 500 hPa geopotential height (colour shading), (c) 10-m wind direction (white arrows) and speed (colour shading), with 95th percentile of 10-m wind speed field plotted as a contour (dashed blue lines), (d) boundary-layer depth (greyscale shading). In (d) the mean sea level pressure minimum is marked by a purple cross [Colour figure can be viewed at wileyonlinelibrary.com] Figure 1d, shows cold sector air, characterised by localised regions of greater depth (≈ 1,400 m) interspersed with shallower regions, primarily to the south and west of the pressure low centre, while the warm sector, characterised by a broad region of consistent shallower depth (≈ 600 m), lay east of the centre as well as wrapping round the centre to the north. Ciara was associated with gale-force (>17 m⋅s −1 ) near-surface wind speeds over large regions of the North Atlantic (off the coast of Ireland), the Irish Sea, English Channel, and parts of the southern North Sea (Figure 1c). Consequently, waves as high as 10 m battered the exposed coastlines of Ireland, Wales, and Cornwall, over-topping sea defences in many places (Kendon, 2020). Although over the land surface winds were reduced due to increased surface friction, the vast majority of the UK's land surface stations reported widespread gusts exceeding 30 m⋅s −1 , most likely due to the unstable conditions of the daytime boundary layer, which favoured turbulent transport of fast-flowing air from Ciara's low-level jets to the surface. These extreme weather conditions, together with the exceptionally intense precipitation brought by Ciara (totalling 50-75% of the 1981-2020 February average of rain for some UK regions), resulted in widespread flooding across the British Isles, 675,000 homes without power for several hours, £1.6 billion of damage, and tragically three fatalities (Kendon, 2020).

The ensemble regional coupled system
We developed the new 18-member Ensemble-RCS, a convective-scale regional ensemble coupled modelling framework, to integrate the features of the operational regional atmosphere ensemble forecast system MOGREPS-UK (Hagelin et al., 2017) with the atmosphere-ocean-wave coupled modelling framework RCS (Lewis et al., 2018). RCS incorporates convective-scale atmosphere, land, ocean and wave model components which can be run in coupled and uncoupled modes. The atmosphere and land component feedbacks are integrated in each coupled and uncoupled RCS mode, consistently with the latest MOGREPS-UK version (Porson et al., 2020). In this implementation, the RCS atmosphere (Met Office Unified Model, MetUM; Walters et al., 2017) and land surface (the Joint UK Land Environment Simulator, JULES; Best et al., 2011) components are run on a domain with 2.2 km horizontal grid spacing, matching the grid of the operational MOGREPS_UK system, and use the RAL1-M convective-scale science configuration (Bush et al., 2020). As described in Lewis et al. (2019), the shelf-sea ocean component uses the Atlantic Margin Model (AMM15) regional configuration of the Nucleus for European Modelling of the Ocean (NEMO; Madec 2016) coded on a grid with 1.5 km grid spacing and ocean surface waves are simulated using the corresponding regional configuration of WAVEWATCH III ® (WW3DG, 2016), defined on the AMM15 domain using spherical multi-cell grid (Li, 2011;Valiente et al., 2021) with 1.5 km grid spacing (termed AMM15-SMC). The LBCs for the regional wave model configuration are taken from a global wave model simulation, as described in Lewis et al. (2019), which accounts for sizeable waves in the domain which are generated by storms located outside the domain. The surface flux parametrizations and ensemble perturbations used by the Ensemble-RCS are described in the following two subsections, and more details can be found in Lewis et al. (2018) and Gentile et al. (2021).

Surface flux parametrizations
The turbulent air-sea momentum flux, 0 , sensible heat flux, H 0 , and moisture flux, E 0 , are parametrized assuming that the Monin-Obukhov similarity theory holds in the surface-layer scheme, yielding and where 0 is the air surface density; C D the momentum transfer coefficient; c p the specific heat capacity; C h the scalar transfer coefficient; U the log-law wind speed; z 0m and z 0h the momentum and scalar roughness lengths; z 1 the height of the bottom model layer above the surface; and Δv, ΔT, and Δq the gradients between surface and bottom model level of the velocity vector v, temperature T, and specific humidity q fields, respectively. In Equations (1)-(3), the transfer coefficients for momentum, C D , and scalar flux, C h , depend on atmospheric stability, Ψ, and momentum, z 0m , and scalar roughness length, z 0h (Bush et al., 2020) as where k = 0.4 is the von Kármán constant, Ψ m and Ψ h are the Monin-Obukhov stability function for momentum and scalar fluxes, and z is the vertical height coordinate. In neutral conditions, parametrizations described in Equations (4) and (5) model an approximately linear growth of transfer coefficients C D and C h with 10-m wind speeds. Over the sea, the momentum roughness length, z 0m , is described by the generalized Charnock formula where = 14 × 10 −6 m 2 ⋅s −1 is the kinematic viscosity of air, u * is the friction velocity, and is the Charnock parameter, which accounts for increased roughness as wave heights grow due to increasing surface stress (Charnock, 1955;Smith, 1988). The RCS computes roughness length for scalar fluxes over the sea, z 0h (sea), based on z 0m according to the surface divergence theory from Csanady (2001); full details are in Edwards (2007) and Gentile et al. (2021). The friction velocity, computed as u * = √ C D |Δv|, together with the value of the momentum roughness length z 0m , determines the wind speed, U, according to the log-law The Ensemble-RCS can be run in the uncoupled atmosphere-only mode, the partially coupled atmosphere-ocean (AO) mode and the fully coupled atmosphere-ocean-wave (AOW) mode. When the RCS runs in uncoupled atmosphere-only mode (Atm-only), the momentum, heat and moisture fluxes are parametrized by the MetUM by making the following assumptions: • The lower boundary of the atmosphere is at rest; • The SST lower boundary condition provided by the Operational Sea Surface Temperature and Sea Ice Analysis (OSTIA; Donlon et al., 2012), persists throughout the simulation; • The Charnock parameter, , in Equation (6) is set to the empirically based value, = 1.1 × 10 −2 , and persists throughout the simulation.
In the partially coupled AO mode, the MetUM atmosphere component receives hourly SST and ocean currents fields from the AMM15 ocean model to update the lower temperature and velocity boundary conditions, respectively, and in turn the scalar and momentum surface fluxes, surface pressure and velocity fields are used to drive the ocean component. As in the Atm-only mode, the Charnock parameter is a constant set to the empirically based value = 0.011. Consequently, both the Atm-only and AO modes compute the momentum roughness length, z 0m , and the momentum and scalar transfer coefficients are a function of only the atmospheric surface layer.
In the fully coupled AOW mode, the MetUM atmosphere component sends the velocity field to both the AMM15 ocean model and the AMM15-SMC wave model. In turn, AMM15-SMC sends to the MetUM a spatially and temporally varying Charnock parameter which the MetUM uses to update the aerodynamic surface-roughness estimate through Equation (6). As a result, the momentum roughness length z 0m and, therefore, the scalar and momentum transfer coefficients, C D and C h , depend on both the atmospheric surface layer flow and the underlying dynamic wave state via Equations (4) and (5)

Ensemble perturbations
The 18-member Ensemble-RCS comprises a reference simulation (without perturbations) and 17 perturbed simulations that implement perturbations of the atmospheric state, but not of the ocean or wave state. The atmosphere components of the 17 perturbed members are initialised by perturbing about the upscaled UKV analysis with the downscaled perturbations from MOGREPS-G, while the reference member is initialised from the unperturbed upscaled UKV analysis. The design and workflow of the perturbations of the ICs and LBCs of the atmospheric component of the Ensemble-RCS forecast is represented in Figure 2, along with the scheme of coupling between the ensemble atmosphere members and the deterministic ocean and wave models. A detailed explanation of the workflow is now given. The Ensemble-RCS builds the perturbed ICs for the perturbed ensemble members by adding to the UKV analysis fields (upscaled from 1.5 to 2.2 km grid spacing) the 17 difference fields between the MOGREPS-G ensemble members and the ensemble mean at analysis time (downscaled from 20 to 2.2 km grid spacing) through the Incremental Analysis Update scheme (Bloom et al., 1991). The reference member uses unperturbed ICs from the UKV analysis upscaled to 2.2 km grid spacing and uses LBCs from the unperturbed MOGREPS-G forecast downscaled to 2.2 km grid spacing. The perturbed ensemble members use the corresponding downscaled LBCs from the perturbed MOGREPS-G forecasts.
MOGREPS-G uses an Ensemble Transform Kalman Filter (ETKF; Bishop et al., 2001) to generate the ensemble perturbations of all the IC and LBC fields except the LBC SST field. The perturbations of the SST field are generated using the following perturbation approach of Tennant and Beare (2014): 1. Compute the power spectrum of the climatological average of day-to-day differences in SST; 2. Set the amplitude of each mode of the power spectrum to that of each mode of a triangular spherical harmonic expansion in the horizontal; 3. Map back to the physical space; and finally 4. Multiply the mapped spectrum by a monthly-mean day-to-day change, calibrating perturbations in the range ≈ ±2 K.
This perturbation methodology generates greater SST uncertainty in regions where day-to-day fluctuations are greater, for example, the Gulf Stream, Kuroshio Current and the Polar Front, which are also the regions above which cyclones form. Optionally, stochastic uncertainty is added to the physics parametrizations via the Random Parameter 2 scheme (RP2; McCabe et al., 2016). A total of ten parameters from the boundary-layer and microphysics parametrization schemes are perturbed in the RP2 scheme affecting boundary-layer mixing, cloud formation, precipitation and droplet settling near the surface. Unlike the operational MOGREPS-UK system, in which the parameters are updated every five minutes via a first-order autoregression model, the RCS-UK stochastic perturbations can only be applied once at the start. More precisely, each parameter P is computed as F I G U R E 2 Ensemble-RCS design and workflow: integration of the atmosphere-only regional ensemble MOGREPS-UK capabilities with the regional coupled system RCS. where is the arithmetic average (P min + P max )∕2, P min and P max are the minimum and maximum value bounds, and is the stochastic shock term sampled from a uniform distribution in the range ±(P max − P min )∕3. The P min and P max bounds are estimated by experts of the field as the values that guarantee the perturbations are physically meaningful and realistic. It is important to note that the arithmetic average is different from the model default value of P used in the unperturbed forecasts. When the Ensemble-RCS is run in the partially coupled AO mode, the atmosphere component of each of the 18 ensemble members is coupled to the AMM15 ocean model component, which runs on the 1.5 km grid and uses the OASIS3-MCT libraries (Valcke et al., 2015) to map the ocean fields on the atmosphere component 2.2 km model grid. When the Ensemble-RCS is run in the fully coupled AOW mode, the atmosphere component of each of the 18 ensemble members is coupled to both the AMM15 ocean and the AMM15-SMC wave components. Identical ocean ICs and LBCs are applied to the ocean component of all ensemble members, and thus, in the absence of additional perturbations, differences between ocean members are due solely to the differences in atmospheric forcing. Hence, to maintain the spread in SST for the atmosphere when run in coupled mode, the IC perturbations derived from MOGREPS-G (Tennant and Beare, 2014) are added back to the coupled SST after the SST field is passed from AMM15 to MetUM at each coupling time step (hourly frequency).

Ensemble regional coupled system experiments
Three Ensemble-RCS experiments (defined below) were carried out to investigate the sensitivity of Ensemble-RCS forecasts to AOW coupling, and to three types of perturbations (to IC, LBC, and model physics parametrizations) for intense wind speeds during cyclone Ciara. Table 1 summarizes the five Ensemble-RCS configurations used in the three experiments, indicating for each ensemble configuration the perturbation strategy and degree of coupling.
The aim of the first numerical experiment was twofold: first, to isolate the impact of applying the RP2 stochastic physics scheme to the 18-member ensemble created by using IC and LBC perturbations and second, to determine what fraction of the sensitivity to the stochastic perturbations is attributable to perturbation of the Charnock parameter. We do not investigate the sensitivity of Ciara's intense surface wind speeds to other parameters perturbed in the full RP2 scheme, such as the stability function and mixing-length parameters, as we expect the sensitivity to Charnock parameter perturbations to dominate over them. Extreme midlatitude cyclones, such as Ciara, are associated with a strong wind forcing at the top of the boundary layer, which leads to a strong wind shear in the boundary layer, dominating in turn the production of turbulent kinetic energy. Because of the mechanical nature of turbulence production, momentum transfer from air to sea is more sensitive to changes in the sea-state roughness z 0m (sea), via the Charnock parameter which in turn affects the sea-surface roughness and frictional dissipation at the surface (Equation (6) and more details in Lewis et al., 2019 andGentile et al., 2021), than to changes in the stability of the boundary layer.
Thus, in this experiment, three uncoupled Ensemble-RCS configurations were set up as follows.
1. The Control_ENS configuration, with IC, LBC, and RP2 stochastic perturbations. The latter is applied to the model parametrizations using P min and P max values listed in McCabe et al. (2016). 2. The Atm-crk_ENS configuration with IC and LBC perturbations and with RP2 stochastic physics perturbations applied only to the Charnock parameter. In particular, the RP2 scheme sets P min = 1 × 10 −2 and P max = 2.6 × 10 −2 for the Charnock parameter. 3. The Atm-only_ENS configuration using only IC and LBC perturbations.
The aim of the second numerical experiment was to isolate and quantify the impact of AO coupling on the Control_ENS configuration of Ensemble-RCS, and in particular the impact of dynamical SST simulated by the ocean model. Therefore, the two-way AO coupled configuration AO_ENS of Ensemble-RCS was set up by coupling each of the 18 Control_ENS members with the AMM15 deterministic ocean model.
Lastly, the aim of the third experiment was to isolate and quantify the impact of additionally coupling to the wave model on the Control_ENS configuration. The two-way AOW coupled configuration AOW_ENS of Ensemble-RCS was set up by coupling each of the 18 Con-trol_ENS members with the AMM15 deterministic ocean model and the AMM15-SMC deterministic wave model.
All the simulations were initialised at 0000 UTC on 7 February 2020, after cyclone Ciara had formed in the MOGREPS-G and well before it came into the MOGREPS-UK domain. The simulations were run for five days until 0000 UTC on 11 February 2020; however, the analysis of meteorological fields of interest does not extend beyond 0000 UTC on 10 February 2020 as after this time cyclone Ciara exited the domain.

Characterization of ensemble dispersion
The spread between the ensemble members of a given Ensemble-RCS configuration provides an estimate of the uncertainty of that Ensemble-RCS forecast. To quantify the ensemble spread between the ensemble member forecasts of cyclone Ciara's 10-m wind speed of each of the five Ensemble-RCS configurations, two dispersion metrics were employed besides the commonly used box and whisker plots: the dispersion form of the Fractions Skill Score (dFSS), which provides a neighbourhood-based perspective on the ensemble spread, and the dispersion form of the Structure, Amplitude, Location score (dSAL), which provides an object-based perspective.
FSS (Roberts, 2008) and SAL (Wernli et al., 2008) were initially developed, and have then been applied, to verify deterministic or ensemble convective-scale precipitation forecasts against radar observations (mainly for convective events). It is well known that, for such events, the precipitation field tends to be noisy, and so using well-established grid-point forecast verification metrics likely leads to the double penalty problem where a small spatial displacement of localised precipitation in a forecast is penalised twice (for being absent where the precipitation was observed and present where no precipitation was observed). The FSS and SAL metrics can be used to overcome this problem. The FSS overcomes the problem because it considers the grid point fraction exceeding the chosen percentile threshold in a neighbourhood of grid cells, whose centre corresponds to the observing site. The SAL overcomes the problem by considering separately the structure (S), amplitude (A), and location (L) components of the ensemble forecast objects, identified by the contiguous sets of grid cells exceeding the chosen percentile threshold. Thus, a displaced precipitation object with the correct amplitude and structure would only have a poor location score.
Similarly to the precipitation field, the 10-m wind speed field in convective-scale simulations also tends to be noisy which implies that the FSS and SAL metrics could be more suitable than grid-point forecast verification metrics for the assessment of the Ensemble-RCS configurations forecasts of cyclone Ciara's wind speeds. To date, this is the first study in which the FSS metric is applied to 10-m wind speeds while the SAL metric has previously been applied to 10-m wind gust speeds associated with wind storms by Zschenderlein et al. (2019). Since we want to compare ensemble member simulations against one another to determine the ensemble spread, rather than verifying a forecast against observations, the dispersion forms of FSS and SAL, termed dFSS and dSAL, are used here.
The dFSS, introduced by Dey et al. (2014) (termed dFSSmean there) and applied to ensemble precipitation forecasts, is used here to compare the similarity in wind fields between the 18 ensemble members of each Ensemble-RCS configuration. First, for each grid cell in each ensemble member forecast, the fraction of surrounding grid cells within a given size square neighbourhood whose field values exceed a specified wind speed percentile threshold is determined. Then, for each of the pairs of ensemble member forecasts, the FSS is computed as where FBS is the Fractions Brier Score, a variation of the Brier score (Brier, 1950) in which both forecast probabilities (fractions) can take any value between 0 and 1, according to where N is the number of grid cells in the domain, and M l and M n are the ensemble member forecast fractions of the (l, n) pair at grid cell j. FBS worst is the worst possible FBS (largest) that is obtained when there is no collocation of non-zero fractions. It is calculated as Second, the dFSS is computed for neighbourhoods of increasing size, by evaluating the mean of the 153 FSS scores in Equation (9). Larger values of dFSS indicate greater spread between ensemble members with possible values between zero and one.
The dSAL is the dispersion form of the object-based verification method SAL, computed here by evaluating the spread of the SAL scores for each of the 153 pairs of ensemble members. The dSAL presented here compares ensemble member pairs to evaluate ensemble dispersion. It differs from that discussed in Zschenderlein et al. (2019) which is instead based on comparing each ensemble member and the ensemble mean with truth. Following the description of SAL in Wernli et al. (2008), a percentile wind speed threshold value is chosen for each ensemble pair (l, n) to identify contiguous wind speed objects. The amplitude component for each pair (l, n) is computed as where the D operator performs the average of a forecast field over the model domain. The location component L is the sum of two parts L = L 1 + L 2 . The first term, L 1 , measures the normalized distance between the two centres of mass x(M l ) and x(M n ) of the (l, n) ensemble member pair: where d is the largest distance between two boundary points of the domain. The second term, L 2 , measures the averaged distance between the centre of mass of the system of the forecast objects and the individual forecast objects. The structure component S describes the shape and size of the forecast objects. It is computed as the difference between the weighted means of the scaled forecast field volume, V, of all objects in the ensemble member pair (l, n), normalized as in component A: The amplitude, A, and structure, S, components are scaled so that their values extend from -2 to +2, and the possible values of location component, L, range from 0 to 1.

RESULTS
In this section we present the results of the relative sensitivities of the wind speed forecast by the ensembles to stochastic perturbations, compared to ICs and LBCs perturbations (the first experiment described in Section 3.2) followed by results of the sensitivity to coupling to ocean and coupling additionally to waves (the other two described experiments). We conclude by comparing the spread characteristics of the coupled and uncoupled Ensemble-RCS forecasts.

Sensitivity to stochastic physics perturbations
This sensitivity experiment aims at determining the relative impacts of applying the RP2 stochastic physics scheme and the IC and LBC perturbations on the ensemble forecasts of cyclone Ciara's low-level wind jets. Figure 3a presents the 12-hourly time series of the distributions of the 95th percentile 10-m wind speed for the model domain obtained from each of the three AO ensembles: the control ensemble with the full set of stochastic perturbations (Control_ENS), the ensemble with stochastic perturbations to the Charnock parameter only (Atm-crk_ENS), and the ensemble with no stochastic perturbations (Atm-only_ENS). For each simulation, the 95th percentile was computed over sea-only grid points to isolate the impact of stochastic perturbations of the Charnock parameter. The median (from the 18 ensemble members) of the Control_ENS 95th percentile 10-m wind speed increases from 15.7 m⋅s −1 at 0000 UTC on 7 February 2020 (hour 0) to the peak value of 23.8 m⋅s −1 at 0900 UTC on 9 February (hour 57). This rapid increase of wind speed highlights the rapid intensification of cyclone Ciara while crossing the North Atlantic, as discussed in Section 2. At peak wind time, the maps of the spread of the Control_ENS 10-m wind speed and MSLP, measured as the standard deviation of the ensemble members with respect to the mean, are depicted in Figure 4a, b. An overlay of the spread of 10-m wind speed with the MSLP spread shows a close spatial agreement. The largest 10-m wind speed spread is localised around the area with sharp wind speed gradient, between the near-surface winds exceeding 18 m⋅s −1 associated with Ciara's densely packed isobars and the winds slower than 10 m⋅s −1 to the north . In this region of sharp gradient, the spreads in 10-m wind speed and MSLP arise from the slightly different tracks and frontal positions that each ensemble member simulates for Ciara. These differences between ensemble members in forecasts of Ciara's large-scale evolution are likely linked to the Control_ENS perturbations in ICs and LBCs. Shortly after cyclone Ciara made landfall in the UK, the Control_ENS shows a rapid decrease of the median 95th percentile 10-m wind speed (between 1200 and 1800 UTC, corresponding to hours 60 and 66 respectively, in Figure 3a) associated with Ciara's decay. The Atm-only_ENS and Atm-crk_ENS of 10-m wind speed follow the same trend as the Control_ENS F I G U R E 3 Time evolution of the spread of the 95th percentile 10-m wind speed forecast (computed at sea-only points) for the Control_ENS (black), Atm-only_ENS (red), and Atm-crk_ENS (green) ensembles. (a) 12-hourly box and whiskers plot for the three ensembles over the period 0000 UTC on 7 February (hour 0) to 0000 UTC on 11 February 2020 (hour 96). The three plotted boxes of different colours relative to the same simulation hour are represented side by side, instead of overlapping, for greater clarity. (b) Hourly time series for the three ensembles over the period 0000 UTC 9 February to 0000 UTC 10 February of the upper (solid) and lower (dashed) whiskers.
In ( in Figure 3a over the entire simulation time. However, the median speed of the Atm-only_ENS ensemble forecast members is consistently slightly higher, by up to 0.5 m⋅s −1 , than that of the Control_ENS and Atm-crk_ENS. To examine in more detail the impact on the ensemble spread of applying the stochastic perturbations to the model physics parametrizations, the hourly time series of the upper and lower whiskers of the distribution of the 95th percentile 10-m wind speeds (indicating values 3 from the mean) are shown in Figure 3b, focusing on the time period during which Ciara crossed the UK. The values of the upper whiskers for the Atm-only_ENS and Atm-crk_ENS follow those for Control_ENS values closely with only slight differences overall (values up to 0.3 m⋅s −1 smaller). In contrast, the lower-whisker values of the Atm-only_ENS and Atm-crk_ENS are more different from those of the Control_ENS. The lower-whisker values of the Atm-only_ENS are consistently above those of the other two ensembles (by up to roughly 1.5 m⋅s −1 ).
The lower-whiskers values of the Atm-crk_ENS generally lie between those of the Control_ENS and Atm-only_ENS until 1200 UTC on 9 February 2020 (hour 60); after this time, the Atm-crk_ENS values collapse to the Con-trol_ENS values, and closely follow them for the rest of the period shown. Note that the dip at 59 hr in the upper whisker of the Atm-crk_ENS is an artificial feature that arises because, unlike at the other times, the second highest value of the 95th percentile 10-m wind speed distribution is plotted rather than the first.
The results of this first sensitivity experiment demonstrate that the larger ensemble spread between the upper and lower whisker values physics parametrizations. In addition, the reduced median of the Atm-crk_ENS compared to the Atm-only_ENS can be explained by the increased surface stress simulated by the stochastic perturbations applied to the Charnock parameters. The perturbed Charnock parameters vary for the different ensemble members, but for a given ensemble member the Charnock parameter is spatially uniform over the grid cells. As a consequence of the P min and P max values used by the RP2 stochastic physics scheme for the Charnock parameter (given in Section 3.2), 16 of the 17 perturbed Charnock parameters created via Equation (8) are larger than the constant value of = 1.1 × 10 −2 used by the Atm-only_ENS (without stochastic perturbations) and Atm-crk_ENS reference (unperturbed) member, with a maximum of = 2.6 × 10 −2 . According to Equation (6), a larger Charnock parameter leads to a larger aerodynamic sea-surface momentum roughness length, a larger air-sea momentum flux, and therefore reduced 10-m wind speeds in Ciara's low-level jets. It is also important to note the close similarity between the ensemble member spread of the 10-m wind fields for the Control_ENS and the Atm-crk_ENS between 60 and 72 hr (Figure 3b), which suggests that the stochastic perturbations in the Charnock parameter are the main cause of the additional ensemble spread due to the stochastic perturbations during this time. However, the relative spread of the Control_ENS and Atm-crk_ENS may also be affected by the gradual divergence of ensemble member forecasts with time.

Sensitivity to coupling to the ocean and waves
An indication of the sensitivity of cyclone Ciara's 95th percentile 10-m wind speed to ocean and wave coupling is given in Figure 5 which shows time series of the absolute difference from the Control_ENS reference (unperturbed) member for the same members from the AO_ENS and AOW_ENS. To place the differences between the AO_ENS and AOW_ENS reference members and the Control_ENS reference member in context, the ensemble spread of the Control_ENS (spread of absolute differences of the members from the ensemble mean), arising from the perturbations to the ICs, LBCs and stochastic physics, is shown by the grey shading.
The absolute difference between the reference members of the AOW_ENS and Control_ENS is much bigger than that between the AO_ENS and Control_ENS (reaching maxima of 1.2 and 0.2 m⋅s −1 , respectively). This implies that the impact of coupling to both ocean and wave models is much greater than that of coupling to the ocean model alone. Moreover, the difference between the AO_ENS and Control_ENS reference member is found to be comparable to about the lower quartile value of the control ensemble, while the difference between the AOW_ENS and Con-trol_ENS reference member is found to be comparable to about the upper quartile value of the control ensemble. This result suggests that the impact of coupling the atmospheric model to both ocean and wave models is at least comparable in size to that of adding IC, LBC and stochastic physics perturbations to the ensemble system when the near-surface wind speeds are intense (above 20 m⋅s −1 , exceeding gale force).
The changes in behaviour of the simulated ensemble wind speeds on coupling to the ocean model, and to the F I G U R E 5 Time series of the quartile spread (light and dark grey bands) of the absolute difference in 95th percentile 10-m wind speed between the 18 Control_ENS members and the Control_ENS mean. The light grey indicates the control ensemble (Control_ENS) spread values below the lower quartile and above the upper quartile, while the dark grey represents the lower-upper quartile range. Overlaid is the absolute difference of the same quantity between the AO_ENS and Control_ENS reference members (purple) and between the AOW_ENS and Control_ENS reference members (red). Vertical green dashed lines indicate the range of hours in 9 February 2020 [Colour figure can be viewed at wileyonlinelibrary.com] ocean and wave models, are now described and further quantified in the following two subsections.

4.2.1
Sensitivity to coupling to the ocean The impact of coupling was further quantified by comparing the time series of the median, quartiles and extremes of the AO_ENS and the AOW_ENS against the Control_ENS ( Figure 6). The small sensitivity of cyclone Ciara's wind speeds to coupling to the ocean, for the reference simulations, relative to the Control_ENS spread arising from stochastic physics, IC and LBC perturbations indicated in Figure 5 is confirmed by examination of Figure 6a, b: the box plots of AO_ENS (red) and Control_ENS (black) wind speeds present negligible differences and the corresponding upper and lower whisker line plots nearly overlap. As a consequence, the spread of Ciara's wind speed forecasts in the AO_ENS can be attributed to the IC, LBC and stochastic physics perturbations with no marked effect from the coupling to ocean. A closer inspection of instantaneous differences between the AO_ENS and the Control_ENS means of 10-m wind speed field at the peak wind speed time 0900 UTC on 9 February 2020 reveals that the sensitivity of Ciara's winds to coupling to ocean is localised to three small regions in the domain (Figure 7a). The AO_ENS 10-m wind speeds exhibit reductions by up to 0.7 m⋅s −1 compared to the Control_ENS off the west coast of Ireland, off the northern coast of France, and between the English Channel and the North Sea. In contrast, in the English Channel and in the Irish Sea increases of up to 1.5 m⋅s −1 are found. Increases in the AO_ENS mean 10-m wind speed occur in those regions where this ensemble predicts higher SSTs than the Control_ENS (Figure 7b), and decreases occur where the SSTs are lower; recall that the SST values for the Control_ENS are persisted from the ICs. More specifically, SST reductions (and associated wind speed reductions) occur off the coast of Brittany and in the southern North Sea (up to 0.4 K in both regions), and off the west coast of Ireland (up to 0.7 K ). Ciara's strong low-level jets, identified by the 95th percentile 10-m wind speed contour, occur in the latter two regions at this time. A plausible physical argument for the response of the AO_ENS 10-m wind speed to SST follows. Lower SSTs likely increase the near-surface potential temperature gradient, enhancing the stability of the surface layer in the AO_ENS forecasts. Since, typically, in the intense wind speed regions of a cyclone the boundary layer is quasi-neutral, a small increase in stability does not consistently alter the overall turbulent kinetic energy budget dominated by wind shear. This lack of consistency leads to very localised and patchy reductions in 10-m wind speeds over the domain. In the regions characterised by weaker wind speeds, such as off the coast of Brittany, the contribution of wind shear to the turbulent kinetic energy budget is likely to be less dominant. Consequently, near-surface wind speeds can be more susceptible to stability increases, as described by Monin-Obukhov similarity theory, leading to a uniform reduction in 10-m wind speeds across the region. In contrast, in regions such as the Irish Sea, where coupling to the ocean increases the SSTs, the 10-m winds also increase because higher SSTs reduce boundary-layer static stability, increasing the entrainment of fast-flowing air from aloft. The impact of this mechanism is most evident far from the most intense Ciara winds.
The spread of the AO_ENS forecasts of SST at peak wind time increases poleward, as the distance from the greater reductions in SSTs increases, reaching a maximum of 1.2 K (roughly six times larger than that in the regions where the impact of coupling to ocean is more important; Figure 7c). Consequently, although the ensemble mean wind speed reduction on coupling to the ocean occurs in regions where there is a reduction in SST, and the wind speed increase occurs where there is increased SST, these regions are characterised by relatively small ensemble spread in SST. Moreover, the spread of SSTs in the AO_ENS at the initialisation and peak wind speed times (0000 UTC on 7 February and 0900 UTC on 9 February, respectively) is very similar (compare Figures 7d, c). Hence, the structure and size of the SST ensemble spread is predominantly attributable to the size and distribution of IC perturbations of SSTs, with these dominating over the dynamically evolved changes in the SST field. This F I G U R E 7 Maps of differences at 0900 UTC on 9 February 2020 between the ensemble means of AO_ENS and Control_ENS for (a) 10-m wind speeds (colour shading) and 95th percentile contour (dashed black line) of the Control_ENS 10-m wind field, and (b) SSTs (colour shading). Maps of AO_ENS SST ensemble spread (grey shading) and SST AO_ENS ensemble mean contours (solid red) at (c) 0900 UTC on 9 February 2020 (d) 0000 UTC on 7 February 2020 [Colour figure can be viewed at wileyonlinelibrary.com] result is expected because the MOGREPS-G SST perturbations (Tennant and Beare, 2014) downscaled to 2.2 km are added back to the AO_ENS SST ensemble forecasts at each coupling time step (as described in Section 3.1.2).

4.2.2
Sensitivity to coupling to waves The time series illustrated in Figure 6a, b shows that the ensemble simulations of cyclone Ciara's wind speeds are characterised only by a small sensitivity to coupling to the ocean, but by a strong sensitivity to additional coupling to waves. The AOW_ENS median of the 95th percentile 10-m wind speeds is consistently reduced by ≈1 m⋅s −1 compared to both the Control_ENS and AO_ENS ensemble median.
Also the values of the upper and lower whiskers are reduced on additional coupling to waves, particularly for the upper whisker, with a maximum reduction of ≈1.5 m⋅s −1 around 0900-1200 UTC on 9 September (especially evident after the peak wind speed time). As the Control_ENS and AO_ENS distributions are nearly identical, the downward shift of the AOW_ENS distribution of Ciara's 10-m wind speeds can be attributed to coupling to the wave model rather than to coupling to the ocean model. The spatial structure of the sensitivity of Ciara's peak wind speeds to coupling to waves is illustrated by the instantaneous differences between the means of the AOW_ENS and the Control_ENS 10-m wind speed forecasts at the peak wind time (Figure 8a). The mean values of the 10-m wind speeds are decreased over most of the domain in the AOW_ENS relative to the Control_ENS. More specifically, the wind speed is reduced by up to 1.5 m⋅s −1 in the strong-wind-speed regions extending from the central to the southern North Sea and off the west coast of Ireland (where the strong-wind region is indicated by a dashed contour outlining the 95th percentile wind speed) and by up to ≈ 1 m⋅s −1 off the north coast of Brittany. From comparison of Figures 7a and 8a, it is evident that the reductions off Brittany and Ireland on coupling to ocean and wave are about 50% larger in size and geographical extent than those on coupling to ocean only. These results indicate that, for this cyclone case, the impact on the ensemble forecasts of additionally coupling to wave is stronger and more coherent than coupling to ocean only, as might be anticipated from the deterministic coupled simulations of cyclone case-studies in Lewis et al. (2018), Wahle et al. (2017) , Lewis et al. (2019), and Gentile et al. (2021).
Regarding the dynamical evolution of the wave state at peak wind time, the snapshot of the AOW_ENS mean of the Charnock parameters (Figure 8b) shows that the momentum transfer into the sea (as inferred from the Charnock parameter) is more pronounced where the 10-m wind speed reductions on additionally coupling to waves are larger (Figure 8a). Comparison of Figure 8b, c reveals that the mean of the AOW_ENS Charnock parameter field increases by up to 1 × 10 −2 from the Control_ENS mean value of 1.4 × 10 −2 , reaching local maxima of 2.4 × 10 −2 where 10-m wind speed reductions are largest (nearing −1.5 m⋅s −1 ). Recall that, unlike in the Control_ENS, the AOW_ENS Charnock parameter field varies temporally as well as spatially. According to the momentum roughness-length parametrization (Equation (6)), the increases in the Charnock parameters on coupling to ocean and waves lead to an enhanced simulated sea-surface aerodynamic roughness. Therefore, according to Equation (4), the sea-surface stress is increased and consequently also the air-sea momentum transfer, explaining the 10-m wind speed reductions in the AOW_ENS compared to the Control_ENS.
Time series of the spread (over the ensemble members) of the Charnock parameter values are shown in Figure 8d by the lower and upper quartile values; note that the range in these values shown for the AOW_ENS arises from the range across the domain (values are constant across the domain for the Control_ENS). The Charnock parameter values in the AOW_ENS are already markedly above those in the Control_ENS even at the start of the period during which Ciara crossed the model domain. As Ciara intensified, the upper quartile of the AOW_ENS distribution of the Charnock parameter field diverged from that of the Control_ENS, growing in size and spread, with the maximum value of ≈2.1-2.2 × 10 −2 reached at peak wind time. The AOW_ENS upper quartile of Charnock parameter increased with time to a peak value at about 50 hr, similar to that of the 95th percentile distribution of wind speed shown in Figure 6. In contrast, the lower AOW_ENS quartile only increased consistently after the peak wind speed time to lie about 0.7 × 10 −2 above the corresponding Control_ENS lower quartile. Hence, the increases in the Charnock parameter on coupling to waves seen in the instantaneous data shown in Figure 8c are consistent across the simulation hours. Moreover, the upper quartile of the AOW_ENS Charnock parameters increased until a peak at 0900 UTC on 9 February 2020 when the 95th percentile 10-m wind speeds are a maximum, which suggests that the most intense near-surface wind speeds are associated with regions of strongest wave growth.
The physical relationship between wind speeds and Charnock parameter explains the similarity between the temporal trends of the upper quartile of the AOW_ENS Charnock parameters and the 95th percentile 10-m wind speeds as follows. As the value of the upper quartile of the AOW_ENS distribution of Charnock parameters increases, the ocean becomes rougher and young growing waves extract an increasingly larger amount of momentum from the overlying airflow, driving the reductions in the AOW_ENS 10-m wind speeds observed in Figure 6a, b. Space-and time-varying Charnock parameters simulated by the AOW_ENS moderate the dynamic response of Ciara's near-surface wind speeds to the evolving wave state. This dynamic response cannot be captured by the Control_ENS and AO_ENS because they use stochastically perturbed Charnock parameter values that are constant at each ocean grid cell throughout the simulation time.
Finally, the impact of coupling to waves was also assessed in terms of wind strike probability for the control and the two coupled ensemble simulations during 9 February 2020, when cyclone Ciara crossed the model domain ( Figure 9). For each member, each grid point in the domain was labelled with a binary value to identify whether the 10-m wind speed exceeded, or not, a storm wind threshold chosen as U 10 = 24 m⋅s −1 , in a 300 km radius circle having that grid point as its centre. The wind storm threshold corresponds to a Beaufort scale of 9 (strong gale), which has associated effects over the sea of high waves with their crests beginning to topple, tumble and roll over. At the end of the labelling process, the probability of the ensemble members exceeding the storm wind threshold (strike probability) was computed for each grid point. As expected, the largest Control_ENS wind strike probability features (>50%) are located in regions corresponding to the areas swathed by the cold and warm conveyor belts of cyclone Ciara (Figure 9a). The strike probabilities for the AO_ENS are generally slightly higher (by between about 10-30%) over a wide region. Consistent with the peak wind speed instantaneous difference maps and the time series of the distributions of ensemble winds, the strike probability shows a stronger sensitivity to coupling to waves and ocean than coupling to ocean alone. Although the strike probabilities reach up to 100% for all three ensembles, the AOW_ENS wind strike probabilities are consistently reduced compared to the Control_ENS, by at least 10% across most of the domain, resulting in an approximate halving of the size of regions where the strike probabilities exceed 50%. In the central North Sea, the Control_ENS wind strike probability is reduced by as much as 70%, with reductions between 50% and 70% in the southern North Sea and the English Channel. The marked reductions in strike probabilities on coupling to waves occur for all the times that Ciara is in the domain, consistent with the dynamical interpretation of the wind-wave feedback provided earlier.

Ensemble dispersion
Two convective-scale metrics were used to characterize the spread of the five ensemble forecasts of Ciara's 10-m wind speed at peak wind speed time (0900 UTC on 9 February 2020): the dFSS and dSAL metrics described in Section 3.3 applied using the 95th percentile wind speed threshold. The dFSS as a function of neighbourhood size for all the Ensemble-RCS configurations is shown in Figure 10. Larger values of dFSS imply that the ensemble members are more similar to one another and the neighbourhood size is the width of the domain over which the similarity is calculated at each grid point. The Atm-only_ENS has the highest dFSS (i.e., most similar members) starting at ∼0.55 at 2.2 km horizontal scale and reaching ∼0.75 at 68.2 km horizontal scale. The Atm-crk_ENS dFSS is less than that of the Atm-only_ENS only by ∼0.05 at the 2.2 km scale but, as the neighbourhood size increases, diverges slightly from that of Atm-only_ENS reaching ∼0.65 at 68.2 km horizontal scale (0.1 lower than the Atm-only_ENS value). The dFSS values for the remaining three ensembles are similar to one another and lower than for the Atm-only_ENS by up to 0.2. The difference between the dFSS for the Atm-only_ENS and Atm-crk_ENS indicates that Charnock perturbations increase the ensemble spread by ≈20% relative to IC and LBC perturbations alone. Adding the complete set of stochastic perturbations (yielding the Control_ENS) enhances the ensemble spread by an additional ≈20%. In contrast, coupling to ocean and wave models does not change the dFSS of the ensembles relative to the Control_ENS. Overall, these results imply that coupling preserves, at all spatial scales considered, the dFSS in 95th percentile 10-m wind speeds, unlike the addition of stochastic physics perturbations which instead increases the dFSS values.
The dispersion of the 95th percentile 10-m wind speed objects was investigated with the dSAL metric. Scatterplots showing the three dSAL components for each pair of ensemble members are shown in Figure 11, with each panel showing results for a different ensemble. For each ensemble, larger amplitude component values are generally associated with larger structure component values, so if the wind speed objects for an ensemble member have larger wind speed values (compared to the other member in the pair being considered), then they are usually also wider or flatter. Note that, because we are comparing pairs of equally likely ensemble members, the two members are arbitrarily assigned as member l or m in Equations (12)-(14) and hence also the sign of the structure and amplitude component is arbitrary as they could both have equivalently the opposite signs, but the same magnitudes, to those plotted. Higher location component values are generally associated with larger magnitudes of the amplitude and structure components, implying that larger differences in amplitude or structure between the ensemble pairs are typically also associated with larger relative shifts in the 95th percentile wind speed area in the ensemble forecast pairs, as most of intense wind speed features occur in single objects.
The spread of the values of each component across the ensemble member pairs can be considered as a further metric for the ensemble spread. Comparison of the dSAL Control_ENS results in Figure 11a against those for the Atm-only_ENS and Atm-crk_ENS in Figure 11b-d, respectively, indicate that applying stochastic physics perturbations in addition to the IC and LBC perturbations enhances the ensemble spread, particularly in the amplitude component. The ensemble spread can also be inferred from the median values of the absolute dSAL components, as given in Table 2. The median values are largest for the structure component, followed by the location and amplitude components. However, comparison of the medians for the different ensembles (Table 2) shows, consistently with the visual impression from Figure 11, that the proportional differences in the spread between the different ensembles is largest for the amplitude component. The spread in the amplitude component is smallest for the Atm-only_ENS and increases on adding the stochastic perturbations to the Charnock parameter (Atm-crk_ENS) and then further on adding the full set of stochatic perturbations (Control_ENS). As can be seen from Figure 11a, c, e, the distribution of the Control_ENS dSAL scatter values is nearly identical to that of the AO_ENS, but exhibits a nearly 20% larger dispersion in the amplitude component, and slightly less than 10% larger dispersion in the structure and location components compared to the AOW_ENS. These results imply that coupling to ocean does not substantially affect the spread of the dSAL components. However, a small reduction in the medians of the absolute structure, location, and, particularly, amplitude component can be observed when coupling also to waves, likely attributable to the reduced wind speeds in the AOW_ENS simulations.
Overall, the neighbourhood-based dFSS and object-based dSAL results agree with the grid-point-based box and whisker plots in indicating an enhancement of atmospheric ensemble spread on applying the stochastic perturbations to the Charnock parameter, and a further enhancement of spread on applying the full set of stochastic perturbations. Grid-point-, neighbourhood-, and object-based metrics also agree on the roughly equal size of the ensemble spread of the Control_ENS and AO_ENS. In contrast, when considering the AOW_ENS ensemble spread relative to those of the AO_ENS and Control_ENS, the dFSS disagrees with grid-point-based and object-based dSAL metrics. In particular, the box and whisker plots and the amplitude component of the dSAL (but less so for the structure and location components) show the AOW_ENS is slightly underspread compared to the AO_ENS and Con-trol_ENS, while the dFSS results indicate that AOW_ENS is slightly overspread compared to AO_ENS and Con-trol_ENS. This difference is likely due to the box and whisker plots and the amplitude component of the dSAL

CONCLUSIONS
Integrating regional atmosphere, ocean, and wave model components into a coupled system is being increasingly trialled by research groups and operational centres (Ricchi et al., 2017;Varlas et al., 2017;Wahle et al., 2017;Lewis et al., 2018). Direct simulation of the effect of the dynamical ocean and wave state on air-sea surface exchange coefficients is expected to better represent surface fluxes and so improve forecasts of weather systems with strong near-surface wind speeds, such as extratropical cyclones. Undoubted benefits have been shown by coupling in convection-permitting, (1) km grid-spacing deterministic models and coarser medium-range, (10) km grid-spacing EPSs (Holt et al., 2011;Ricchi et al., 2017;Varlas et al., 2017;Lewis et al., 2019). However, with the exception of Bousquet et al. (2020) who investigated the performance of a 1D ocean-mixed-layer model ensemble, there are, to date, no published studies considering whether dynamical coupling could also bring benefits to convective-scale short-range EPS. The present study investigated, for an intense cyclone case-study, the respective impacts of atmosphere-ocean-wave coupling, IC and LBC perturbations and stochastic physics perturbations within a convective-scale ensemble coupled system. For that purpose, we developed the first (to our knowledge) ensemble regional coupled system, the 18-member 2.2 km grid-spacing Ensemble-RCS focused on the British Isles and surrounding seas, by implementing the ensemble capability of the Met Office's operational ensemble system (MOGREPS) in their regional coupled system (RCS). The resulting Ensemble-RCS was run in coupled and uncoupled modes for the severe weather period of 7-10 February 2020, which saw the most intense cyclone, named Ciara, crossing the UK since storm Tini (12 February 2014). The impacts of coupling to ocean only and to both ocean and waves were assessed. Three different versions of the uncoupled system were also run to determine the impact of the stochastic physics perturbations, with these perturbations applied at initialisation time and persisting throughout the forecast. The spread between the ensemble members of a given configuration was quantified using, besides the commonly used box and whisker plots, the dispersion forms of the neighbourhood-based metric Fractions Skill Score (dFSS) and the object-based metric Structure, Amplitude, Location score (dSAL).
Applying the stochastic physics scheme to the ensemble created by using IC and LBC perturbations (Con-trol_ENS and Atm-only_ENS, respectively) enhanced the ensemble spread of Ciara's 95th percentile wind speeds, as assessed by both box and whisker plots of wind speed values and the dFSS and dSAL metrics. Moreover, the reduction by up to 0.5 m⋅s −1 in the median of the ensemble wind speed indicated that the ensemble distribution of Ciara's intense wind speeds was shifted slightly downward by the stochastic physics perturbations. Comparison of Con-trol_ENS and Atm-crk_ENS configurations showed that a substantial proportion of the impact obtained from applying the full set of stochastic physics perturbations was attributable to stochastic perturbations of the Charnock parameter alone, revealing the relative importance of this specific perturbation.
The ensemble spread was nearly unchanged on coupling to ocean (AO_ENS), indicating the localised nature of the ocean's impact on the ensemble. Small increases in the 10-m wind speeds occurred, likely due to decreased static stability in the surface layer by SST warming; similarly, small decreases in winds occurred, likely due to increased static stability by SST cooling. However, only the increases in wind speeds affected the AO_ENS strike probabilities, corresponding to enhancements of roughly 10-30%.
Additionally coupling to waves (AOW_ENS) led to a marked downward shift by 1 m⋅s −1 in the median of the ensemble distribution of Ciara's intense wind speeds, along with a consistent reduction of the ensemble mean by up to 1.5 m⋅s −1 at peak wind time. Although a comparable impact on coupling to waves has been observed by a number of studies running deterministic coupled forecasts of cyclones striking both the Mediterranean basin and the North Atlantic (Ricchi et al.,, 2017;Varlas et al., 2017;Wahle et al., 2017;Gentile et al., 2021) here we report the consistency of this impact across all the perturbed (and the unperturbed) members of the ensemble forecast of Ciara for intense 10-m wind speeds. Moreover, this impact proved to be consistent for all Ciara's simulation hours, as indicated by the observed reduction by up to 70% in value and 50% in the extent of high (exceeding 50%) wind strike probability relative to the Control_ENS. The relationship between the ensemble distribution of the AOW_ENS Charnock parameters and the 95th percentile ensemble Ciara wind speeds confirmed that the coupling to wave impact could be attributed to the response of wind speeds to the young ocean waves dynamically simulated by the wave model. In contrast, the other ensemble configurations were unable to capture this feedback because, for a given simulation member, the Charnock parameter was constant in time and for all grid cells, though it is worth noting that a higher minimum bound for the stochastic Charnock parameters in Control_ENS would have probably lessened the impact of dynamical coupling to waves. Remarkably, the impact of coupling to both ocean and wave models on cyclone Ciara's intense wind speeds is found at least comparable in size to that of adding IC, LBC and stochastic physics perturbations to the ensemble system. Although a strong sensitivity to coupling to waves has been already observed in deterministic simulations of cyclone wind speeds (Varlas et al., 2017;Wahle et al., 2017;Lewis et al.,, 2018;Gentile et al., 2021) the Ensemble-RCS results establish, for the first time, the size of this sensitivity relative to the ensemble spread from IC, LBC and stochastic physics perturbations.
The spread characteristics of the Ensemble-RCS on coupling were further examined with dFSS and dSAL metrics. The reduced ensemble wind speeds on coupling to waves (in the AOW_ENS) led only to a slight impact on the dFSS metric and the structure and location component of the dSAL metric, compared to the Control_ENS, but to a 20% reduction in the median of the absolute value of the amplitude component of the dSAL metric. However, this reduction was only ≈ 1∕2 of that from removing the full stochastic perturbations applied to the model physics parametrizations. Together with the fact that AO_ENS showed similar dSAL and dFSS values to Control_ENS, these findings corroborate that coupling preserves the ensemble spread driven by IC, LBC and stochastic physics perturbations.
This study demonstrates that coupling to waves, thereby addressing model physics errors in convective-scale EPS arising from the failure in atmosphere-only model configurations to account for a dynamic sea state in the parametrization of air-sea momentum flux, can have a consistent impact across ensemble members while preserving the ensemble spread driven by IC, LBC and stochastic physics perturbations.
Moreover, we demonstrate that object-based dSAL and neighbourhood-based dFSS dispersion metrics are useful for assessing the spread of convective-scale ensemble forecasts of cyclone wind speeds.
For future work, a broader range of weather conditions should be tested, such as a convective summer case-study, assessing a wider range of meteorological variables. Moreover, the ensemble simulations should be compared against in situ and satellite observations over a season-long case-study to determine the impact on forecast skill. This aspect is outside the scope of this initial study, focused on characterising the sensitivity of cyclone near-surface wind speeds to ocean and wave feedbacks. Further investigation should also focus on the link between the uncertainty in the atmosphere and ocean boundary layers, implementing IC perturbations in the ocean and wave models and then integrating with an atmosphere-ocean-wave coupled data assimilation system (Lea et al., 2021).