Flow‐dependent stochastic coupling for climate models with high ocean‐to‐atmosphere resolution ratio

This study introduces a new flow‐dependent distribution sampling (FDDS) scheme for air–sea coupling. The FDDS scheme is implemented in a climate model and used to improve the simulated mean and variability of atmospheric and oceanic surface fields and thus air–sea fluxes. Most coupled circulation models use higher resolutions in the sea ice and ocean compared to the atmospheric model component, thereby explicitly simulating the atmospheric subgrid‐scale at the interface. However, the commonly applied averaging of surface fields and air–sea fluxes tends to smooth fine‐scale structures, such as oceanic fronts. The stochastic FDDS scheme samples the resolved spatial ocean (and sea ice) subgrid distribution that is usually not visible to a coarser‐resolution atmospheric model. Randomly drawn nodal ocean values are passed to the corresponding atmospheric boxes for the calculation of surface fluxes, aiming to enhance surface flux variability. The resulting surface field perturbations of the FDDS scheme are based on resolved dynamics, displaying pronounced seasonality with realistic magnitude. The AWI Climate Model is used to test the scheme on interannual time‐scales. Our set‐up features a high ocean‐to‐atmosphere resolution ratio in the Tropics, with grid‐point ratios of about 60:1. Compared to the default deterministic averaging, changes are largest in the Tropics leading to an improved spatial distribution of precipitation with bias reductions of up to 50%. Enhanced sea‐surface temperature variability in boreal winter further improves the seasonal phase locking of temperature anomalies associated with the El Niño–Southern Oscillation. Mean 2m temperature, sea ice thickness and concentration react with a contrasting dipole pattern between hemispheres but a joint increase of monthly and interannual variability. This first approach to implement a flow‐dependent stochastic coupling scheme shows considerable benefits for simulations of global climate, and various extensions and modifications of the scheme are possible.


INTRODUCTION
Current state-of-the-art climate and also weather and seasonal forecast models represent the different climate subsystems, such as atmosphere, ocean, sea ice, and land surface, through a variety of modelling approaches. Different numerical methods are used to discretize these subsystems, which often require computational grids of a certain geometry. This often produces a climate model that is constructed from independent model components that are coupled at their respective interfaces (Valcke, 2013). The typical scales of instabilities in the ocean and atmosphere differ by an order of magnitude. As a result, ocean components in current climate models participating in the fifth phase of the Coupled Model Intercomparison Project (CMIP5) apply numerical grids that are more highly resolved than their atmospheric counterparts (Taylor et al., 2012). To give an example, the (meridional) resolution in the ocean around the equatorial wave guide is often on the order of 1/3 • (≈37 km) up to 1/4 • (≈28 km) (e.g. Delworth et al., 2012), while it may be 1 • -2 • or coarser in the atmosphere. While day-to-day weather patterns are resolved at a resolution of about 1 • (e.g. Hewitt et al., 2017), the ocean resolution in the upcoming CMIP6 DECK (Diagnosis, Evaluation, and Characterization of Klima) simulations will at most be eddy-permitting at around 1/4 • resolution globally (a tentative list as of April 2019 can be found at http://rawgit.com/WCRP-CMIP/CMIP6_CVs/master/src/ CMIP6_source_id.html). This resolution allows the simulation of Tropical Instability Wave (TIW) activity in the tropical ocean (Graham, 2014;Rackow et al., 2018).
Looking more closely into the coupling procedure between atmosphere and ocean, deterministic coupling schemes in current climate models with different resolution in each component are commonly based on spatially averaged higher-resolution ocean/ice fields communicated to a coarser atmospheric grid, with ocean-to-atmosphere grid-point ratios from 4:1 up to 100:1 (see link above). This averaging procedure, however, leads to a loss in available spatial information and variability through the smoothing of surface gradients. Ultimately, spatial surface flux variability received by the atmosphere and communicated back to the ocean model is reduced. A smaller number of climate models computes the air-sea fluxes directly on the finer oceanic grid, but even then the fluxes will be averaged in the coarser atmospheric model component for further computations. In any case, potentially available surface flux variability on the oceanic subgrid-scale cannot be resolved in these coupled models. Hasselmann (1976) showed that the representation of high-frequency variability by stochastic noise in a simplified coupled system affects the model's time-scales. Frankignoul and Hasselmann (1977) further showed that long time-scale sea-surface temperature anomalies are generated via the oceanic surface layer, which integrates shorter time-scale atmospheric noise. Oceanic eddies and internal ocean dynamics are another process driving surface field and consequently surface flux variability (e.g. Frenger et al., 2013). In recent years, Williams (2012) investigated the impact of increased surface flux variability in a low-resolution coupled climate model. Williams (2012) applied a stochastic approach, where he augmented the surface fluxes communicated between the atmospheric and oceanic model components by (multiplicative) white noise. His results showed a systematic impact on the oceanic mean state and coupled El Niño-Southern Oscillation (ENSO) variability, although the applied noise was symmetric. Generally, symmetric perturbations in a nonlinear system can cause nonlinear rectifications in models (see e.g. Juricke et al., 2012;Juricke and Jung, 2014;Weisheimer et al., 2014;Strommen et al., 2019) in addition to changes in their temporal time-scales and in their amplitudes of variability (Palmer, 2012;Cooper and Zanna, 2015;Berner et al., 2017;Cooper, 2017;Davini et al., 2017;Juricke et al., 2017;Juricke et al., 2018;Leutbecher et al., 2017).
In the case of coupled climate models, subgrid-scale surface field information is in fact available, as part of it is resolved by the oceanic model component. Similar to coarse-graining methods Bessac et al., 2019;Christensen, 2019), this information can be used in a stochastic coupling scheme to better represent the statistical properties of the unresolved scales in the atmospheric model component. Using this information can improve the representation of spatial variability of surface fluxes. While previously surface fluxes and atmospheric forcing have been perturbed to analyse the impact of high-frequency noise on the coupled climate mean state (e.g. Williams, 2012;Berner et al., 2017;Christensen et al., 2017), surface field variations have only been investigated as a consequence of increased resolution, e.g. in the context of North Atlantic storm tracks (Small et al., 2013) and over the oceanic front associated with the Gulf Stream current (Minobe et al., 2008;Siqueira and Kirtman, 2016).
This study aims to incorporate some of the available subgrid surface field information into the air-sea exchange through a new coupling scheme based on "flow-dependent distribution sampling" (FDDS). By this, the new scheme potentially accounts for surface field uncertainty and also improves the representation of surface flux variability. The stochastic FDDS scheme is implemented in the Alfred Wegener Institute Climate Model (AWI-CM/ECHAM6-FESOM: Sidorenko et al., 2015;Rackow et al., 2018;2019) and compared to the model's standard deterministic coupling method. The AWI-CM in the configuration analysed by Rackow et al. (2018) is especially suited to determine effects of the FDDS scheme on tropical and global climate as it applies (isotropically) increased spatial ocean resolution of less than 25 km in a band along the Equator, while keeping mostly ≈100 km elsewhere, and a coarser resolution of ≈200 km in the atmosphere.
The article is organized as follows: in section 2, we will introduce the model, explain the old deterministic coupling scheme, motivate and describe the new stochastic FDDS scheme, and detail the experimental set-up. The effects of the FDDS coupling on the modelled atmospheric and oceanic mean fields and their variability will be investigated in section 3, with a focus on interannual time-scales. Although one would not necessarily expect that a first test case of a new coupling scheme leads to model improvements only -in particular without re-tuning the model -we can already show considerable bias reductions (e.g. for central Pacific precipitation and ENSO phase locking). Section 4 concludes with a summary of the results and discusses implications and future perspectives.

Model set-up
For the experiments of this study, AWI-CM is applied in the "REF" configuration documented by Sidorenko et al. (2015) and Rackow et al. (2018). AWI-CM is a coupled global configuration of the Finite Element Sea Ice-Ocean Model (FESOM: Timmermann et al., 2009;Sidorenko et al., 2011;Wang et al., 2014), developed at the Alfred Wegener Institute (AWI), Helmholtz Centre for Polar and Marine Research, and the state-of-the-art atmospheric model ECHAM6 developed by the Max-Planck-Institute for Meteorology in Hamburg (Stevens et al., 2013). In the configuration of this study, FESOM and ECHAM6 use a time step of 30 and 10 min, respectively, and are coupled every 6 h through exchanges of ocean surface fields and atmospheric fluxes (figure 2 in Sidorenko et al., 2015) via the OASIS3-MCT coupler (Valcke, 2013).

Ocean-to-atmosphere resolution ratio
ECHAM6 is a spectral atmospheric model and is applied in T63L47 resolution, corresponding to a horizontal resolution of 1.875 • (≈200 km at the Equator) and 47 levels in the vertical. FESOM's REF-grid has a nominal resolution of ≈100-150 km in the open ocean and is gradually refined along coastlines and toward the Equator, down to less than 25 km in the equatorial wave guide (see figure 1 in Sidorenko et al. (2015) and Rackow et al. (2018) for visualizations of the grid resolution). Thus, the number of grid points along the Equator (as well as around the coast of Greenland) is strongly increased compared to ECHAM6, with ocean-to-atmosphere grid-point ratios of 60:1 (more than 40:1) (see Figure 1 and the schematic in Figure 2).

Motivation of stochasticity in the coupling procedure
The classical approach of subgrid parametrizations is based on the assumption that a relationship between the average effect of unresolved subgrid processes on the resolved mean flow can be captured by a theoretical or heuristic approximation. For this assumption to be valid, it is necessary to have a sufficient amount of subgrid processes available in each model grid box. A common example of such an averaging procedure is the effect of single air molecules on the mean diffusion in a box. It is not necessary to know each molecule's exact movement to model the mean diffusion, since the number of molecules is large. Similarly in Atmospheric Sciences -although with less drastic magnitudes -we do not need to know each cloud within e.g. a 2 • × 2 • atmospheric grid box to approximate the average cloud effect on the mean flow. However, this averaging assumption starts to break down once the "subgrid" process actually has an extent of similar scale as the grid box. Some of the larger features can be modelled by the grid resolution now, while the smaller ones still need to be parametrized, and averaging over them might still be a valid approach (see Palmer (2012), especially Figure 3).
The ocean grid resolution applied with FESOM in the Tropics explicitly resolves TIWs. With deterministic ocean-to-atmosphere coupling, an averaging of all the ocean cells in an atmospheric grid box tends to smooth out most of these wave structures. However, the scales of these waves are often not much smaller than the actual atmospheric grid boxes. As a consequence, we often find only a few distinct subgrid surface structures within an atmospheric grid box. For example, there might be just one oceanic SST front ( Figure 2, top right), or -in the high latitudes -an ice edge location that simply separates the grid box into both a continuous open ocean and an ice-covered part. As the response of the entire atmosphere in such situations has a nonlinear F I G U R E 2 Schematic of the stochastic FDDS coupling scheme for climate models with high ocean-to-atmosphere resolution ratio, and time series of sea-surface temperature for a single atmospheric grid box in the equatorial Pacific. (Top) In our study, the ocean-to-atmosphere grid point ratio is about 60:1 in the equatorial region. The default coupling scheme computes an average over all ocean nodes inside an atmospheric box (orange box), whereas the stochastic FDDS coupling scheme is based on drawing of random ocean grid points (blue point) and their respective nodal values. (Bottom left) Six-hourly SST time series at (126.6 • W, 2.8 • N) (orange grid box in top figure) over 1 year with FDDS coupling (blue). During this simulation with active FDDS coupling, the averaged SST was written out as well (orange). Differences due to the nonlinear response of the system to the changed surface coupling are therefore not present. (Bottom right) Difference between the instantaneous time series for this specific atmospheric grid box, and the corresponding histogram of the 6 h differences. Again, these differences are taken from the same stochastic simulation where the averaged SST was written out as well. When compared to the averaged SST from the same run, the SST perturbations implied by the FDDS scheme are physics-based, they vary in magnitude over the year (year-round standard deviation of 1.0), and they have zero mean. Note that separate stochastic and deterministic runs would diverge quickly due to the chaotic nature of the climate system dependence on the ocean state, averaging is no longer a very adequate representation of the surface conditions: the surface grid-box mean is not a state that in general actually occurs. This leads to an underrepresentation of surface variability on spatial scales similar to the atmospheric grid-box size.
A possible solution is the use of stochastic parametrizations. These are designed to not just represent the mean effect of the subgrid processes on the mean flow, but also include higher-order moments such as subgrid variability. It is important to note that because a stochastic parametrization usually samples subgrid conditions from a given distribution function of possible states randomly, it is not designed to capture at each instance in time the actual (unknown) subgrid conditions. Instead, it tries to simulate in a temporally or spatially averaged sense the probability of possible subgrid states. Therefore, in one time instance it could sample the effect of just one large ocean eddy within an atmospheric box, and at a later time the effect of several such eddies, while still describing similar resolved mean flow conditions.
The advantage is two-fold. For one, the system experiences subgrid forcing with more realistic higher-order moments, which enables the atmosphere to react to the individually sampled states. These states themselves are more realistic representations of possible subgrid states than just an average. The second aspect is that stochastic parametrizations can be used as an uncertainty estimate for the flow. Since we cannot know the real state of the subgrid, we can use stochasticity (b,c) on the coarser atmospheric grid, all from the same STOCH experiment. Panel (b) shows averaged SSTs that are usually applied with deterministic coupling (but which were not applied in STOCH) while panel (c) gives the randomly drawn values. (d) Instantaneous differences between the different coupling approaches, which are strongest along the front associated with the Pacific cold tongue. The differences do not include the nonlinear response of the system since they are from the same experiment to sample possible states and their effects on the mean flow, providing us with some bounds on the estimated accuracy of our simulation. The latter generally makes use of ensemble simulations to generate probabilistic model predictions that incorporate the uncertainty of unresolved processes.
The AWI-CM configuration applied here has the advantage that the subgrid structure and therefore the respective process distribution is locally available since it is resolved by the ocean model (especially in the Tropics). The only unbiased deterministic way to map higher-resolution ocean data to a coarser atmospheric grid is spatial averaging of the surface fields (or of the locally calculated fluxes). However, from the atmospheric perspective, we can use the subgrid distribution of the unresolved surface structures to supply a stochastic parametrization with a temporally and spatially adaptable (flow-dependent) distribution function that is resolved and supplied by the ocean model. This is an advantage compared to other stochastic parametrizations that commonly assume a given, more or less static, distribution function derived from theory, previous high-resolution simulations (e.g. for the stochastically perturbed parametrization tendencies (SPPT) scheme: Leutbecher et al., 2017), or observations. We thereby sample the available, resolved subgrid ocean variability to enhance surface flux variability in the simulation without making further assumptions about temporal or spatial variations in the amplitude of the stochastic parametrization forcing. In summary, the stochastic FDDS scheme presented here is an extension to earlier approaches that used predefined distribution functions (employed, for example, by Williams (2012)).

Deterministic vs. stochastic ocean-to-atmosphere coupling
The standard deterministic coupling of AWI-CM from the ocean to the atmosphere is implemented as follows. The ocean surface fields (sea-surface temperature (SST), sea ice thickness (SIT) and concentration (SIC), snow thickness (SNT)) are averaged over the 6 h interval preceding the coupling. Note that the SST in FESOM is a "skin SST" at the actual ocean surface. Values on the atmospheric grid are determined via spatial averages over the nodal values of the ocean grid points contained in the respective atmospheric grid boxes. This is done whenever there is a reasonable ocean-to-atmosphere ratio (at least three ocean grid points in the atmospheric box: Sidorenko et al., 2015); otherwise, i.e. for similar resolution in the ocean and atmosphere, values on the atmospheric grid are determined via bilinear interpolation. This form of averaging dampens the spatial variability of ocean cells within an atmospheric grid box, especially if the ratio of ocean cells per atmospheric box is large (see Figure 1). Other models also use bilinear re-gridding for the ocean-to-atmosphere exchange, e.g. Centre National de Recherches Météorologiques CNRM-CM5, or even nearest-neighbour interpolation as done in Centro Euro-Mediterraneo sui Cambiamenti Climatici Earth System Model (CMCC-ESM: Valcke, 2013).

Construction of subgrid distributions
For the stochastic FDDS coupling, fields are still temporally averaged for each ocean node over the 6 h interval preceding the coupling. As motivated in section 2.3, the distribution functions of the oceanic surface structures within the individual atmospheric grid cells act as a temporally evolving (atmospheric) subgrid distribution. A discrete frequency distribution of the physical nodal values at the ocean grid points i can now be constructed from the ocean cell area A i , divided by the total ocean area within the atmospheric cell (ΣA i ). In other words, nodal values of smaller ocean cells have smaller probability of being drawn than nodes with larger cells. In principle, the finite-element representation of the ocean surface fields represents a continuous field that could have been sampled as well. However, we decided to use the discrete nodal representation since it provides physical values of slightly higher magnitude (both positive and negative), thereby accounting for part of the missing subgrid ocean variability that is still not resolved. To summarize the approach, for every atmospheric grid box containing at least 6 ocean nodes (which is mainly along the Equator, the coastlines and in the Arctic regions; see Figure 1), the FDDS coupling scheme randomly selects a node i based on the local probability distribution and uses its nodal value for the entire grid box, as illustrated in Figure 2 (top). The random draw is independent for the surface temperature (SST i ) and sea ice concentration (SIC j ). The sea ice thickness and snow thickness, however, are related to SIC through the ratios SIT/SIC (actual ice thickness) and SNT/SIC (actual snow thickness), whose values are requested by ECHAM6. These ratios thus need to be within a realistic range. To achieve this, the snow and ice thicknesses are consistently taken from the same node j that was drawn for the determination of SIC j , i.e. SNT j and SIT j . Enabling a separate choice for SST (e.g. SST i with the possibility of i ≠ j) in atmospheric boxes with both ice-free and ice-covered parts has only an effect on the ice-free part. ECHAM6 does not use SSTs below ice-covered parts to compute fluxes as it assumes a fixed freezing temperature under the ice. Therefore, unrealistically increased SSTs under the ice are of no consequence, while the ice-free part may potentially experience freezing temperatures which were previously assigned to an ice-covered ocean node. This may considerably change instantaneous fluxes over ice-free parts of partly ice-covered atmospheric boxes. The separate choice of ocean and ice grid points to provide surface fields for the atmospheric model presents an upper limit of potential ocean grid variations within an atmospheric box, i.e. can be viewed as a relatively strong while still realistic perturbation.

Flux conservation
By design, a global conservation of air-sea fluxes is guaranteed, since we only perturb the surface fields that enter the flux computations, not the fluxes themselves. The net fluxes of heat and freshwater are therefore still conserved between the ocean and atmosphere compartments (Sidorenko et al., 2015) within the individual simulations, so that, for example, losses of heat in one compartment are gained in the other compartment, resulting in an energetically consistent exchange. In the AWI-CM version used in this study, the fluxes are remapped from the atmospheric grid to the ocean grid with an in-house algorithm that conserves the net flux by globally distributing the residual flux on the target grid, proportionally to the pattern of the original flux (Sidorenko et al., 2015). This is similar to the OASIS3 global conservation option GLBPOS (Valcke, 2013). Thus, the FDDS scheme does not suffer from artificial mass or energy sources that could result from perturbing the fluxes directly, and the law of large numbers does not have to be invoked (as by Williams (2012)).

Characteristics of subgrid distribution sampling
From an energetics perspective, it is arguably much less critical to identically conserve fields like SST themselves between the ocean and atmosphere as long as these are within realistic ranges. As a key advantage of the FDDS scheme, the magnitude of the resulting SST (and ice/snow) perturbations is based on resolved dynamics, varies seasonally, and does not rely on fixed, predefined values (see Figure 2, bottom panels). The FDDS coupling scheme will be active in regions where the ocean is higher resolved than the atmosphere, and will increase temporal variability of the coupling fields and surface fluxes in regions of high spatial variability, e.g. along fronts, in the Tropics where TIW activity is high, or along the sea-ice edge. An example of this can be seen in Figure 3 where the effect of the two different coupling schemes -deterministic vs. stochastic -is highlighted given a certain simulated ocean SST field in the tropical Pacific. The difference between the two schemes is most dominant along the boundary of the Pacific cold tongue and along off-equatorial TIW filaments where temperature gradients are strong ( Figure 3d). The active FDDS coupling results in less smoothed surface fields not unlike the original SST on the native ocean grid (compare Figure 3c and Figure 3a). It is worth noting, however, that the resulting fields on the atmospheric grid are not vastly different to the SST fields with deterministic coupling (Figure 3b), so that we do not expect changes in the spatial gradients to negatively affect model stability.

Experimental set-up
We perform two experiments with AWI-CM. The "STOCH" experiment uses the FDDS coupling scheme and branches, after year 450, from a 1500-year present-day control simulation (Rackow et al., 2018). Since our focus is on interannual time-scales, a set of 9-year simulations is performed, initialized from 15 starting dates, resulting in a 135-year dataset.
The 9-year simulations are all started at the 1st of January and chosen 10 years apart in order to better sample decadal variability.  (Sidorenko et al., 2015;Rackow et al., 2018). Also, the FDDS scheme is not leading to a significant change in the top of the atmosphere total radiation or surface heat flux balances (see Figure S9 in Appendix S1 for global differences in these fields), or the model drift (not shown).
In section 3, model-means are simply presented as annual means over the 135-year STOCH and REF experiments. However, whenever we refer to interannual (and monthly) variability, this is computed as the square root of mean variances of all 15 linearly detrended 9-year simulations, resulting in a standard deviation (std). The same is true for the computation of histograms inside tropical boxes, which are computed separately for the fifteen 9-year simulations. Since the seasonal cycle could be affected by the FDDS scheme, we deliberately decided to present monthly variability without removing the seasonal cycle. However, while reduced in amplitude, the results are similar with the seasonal cycle removed (not shown).

(Near-)surface temperature
Figure 4a-c shows the differences in 2m air temperature mean state compared to REF, as well as interannual and monthly variability changes between STOCH and REF, averaged over the 15 9-year time slices (see also the supporting information, Figure S1 in Appendix S1, for the absolute fields of REF).
The impact of the FDDS coupling on the mean state is most pronounced in mid-to high-latitudes, with the largest signature over ocean areas. The Northern Hemisphere (NH) shows a clear increase in mean 2m temperature by more than 0.2 K in large areas over the North Pacific, the northeast of North America, the north-eastern North Atlantic, and over the Barents and Kara Seas in the Arctic. The Southern Hemisphere (SH) on the other hand exhibits temperature reductions of similar amplitude in the seas surrounding the Antarctic continent, most pronounced in the Atlantic and Indian Ocean sectors of the Southern Ocean. This remarkable dipole of NH increases and SH reductions has different origins. In the north, it is to a large degree related to the stochastic sampling of sea ice conditions near the ice edge by FDDS. The stochastic sampling increases turbulent heat fluxes from the ocean to the atmosphere, leading to an atmospheric warming (see heat flux changes in Figure S10 in Appendix S1). This will be discussed in more detail in section 3.5. For the North Pacific and North Atlantic, the enhanced sea-surface temperature variability simulated by the FDDS scheme along the coastal Kuroshio and Gulf Stream currents also plays a role in increasing heat fluxes to the atmosphere and a consequential increase in atmospheric temperatures. The SH patterns, however, are largely related to circulation changes and teleconnections associated with the impacts in the tropical Pacific and Indian Ocean. The FDDS coupling scheme does not act directly in these remote regions, since the ocean model resolution is comparable to the atmospheric resolution, as opposed to the Arctic (see section 2 and Figure 1).
Independent of the considered reference period in ERA-Interim, the SH is simulated much warmer than observed in REF, while the NH tends to be too cool in the model -in particular in the 1990-2016 and 2000-2016 periods ( Figure S2 in Appendix S1). Although it is difficult to compare a pre-industrial control run with observations from a transient climate, the dipole of warm anomalies in the NH and cooler anomalies in the SH after the introduction of FDDS in the STOCH experiment therefore tends to reduce this inter-hemispheric bias. While the mean 2m temperature change in the Tropics is small and not significant, there is a strong change in tropical variability, which is even connected to variability changes in the Southern Ocean. This is true for both the mean interannual ( Figure 4b) and monthly (Figure 4c) standard deviation. The observed widening of the eastern tropical Pacific variability tends to improve the simulation, especially for interannual variability (compare REF temperature variability biases shown in Figure S3 in Appendix S1).
The SST perturbations generated by the FDDS scheme enable a better communication of TIW variability to the atmosphere, especially with respect to anomalously high temperatures (≳27.5 • C). This is reflected by the change of the monthly 2m temperature distribution shown in Figure 4d,e. While the mean temperature in the central tropical Pacific region (grey box in Figure 4a) does not change significantly, the right tail for large temperature values is considerably extended by the stochastic perturbations, and temperatures lower than ∼27 • C are becoming less likely accordingly.

Precipitation
A direct consequence of the enhanced variability and increased probability of high tropical SSTs and 2m temperatures is the increased triggering of convective precipitation events. In order to trigger convection in the TNT (Tiedtke-Nordeng-Tiedtke) convection scheme, Stevens et al. (2013) note that ECHAM6 involves a prognostic equation for the temperature variance in the planetary boundary layer, while the predecessor ECHAM5 uses a fixed value of 0.5 • C for "the temperature excess of the triggering plume", with overall similar simulation quality between both approaches (Stevens et al., 2013). Evidently, the perturbations with the FDDS scheme are strong enough to affect the precipitation, as shown in Figure 4f,g for a tropical box (dashed blue box in Figure 4a) downstream of the temperature box.
Here, both the mean and variability are considerably changed by the stochastic SST perturbations (see also the supporting information, Figure S4 in Appendix S1, for the absolute fields of REF). Mean precipitation is enhanced from 1.84 to 2.09 mm⋅day −1 , and so is the probability of high monthly mean precipitation, e.g. of events of 10 mm⋅day −1 or higher. It should be noted again that these precipitation changes can mostly be ascribed to the changes in parametrized convective A different sign when compared to the bias pattern translates to a bias reduction. Stippling indicates significance according to a rank sum test (an F-test) at the 5% level for mean (variability). Similar to the procedure for the model results, a linear detrending has been performed on four 9-year slices of GPCP data prior to the computation of standard deviations, computed as the square root of mean variances (averaged over the set of 9-year slices) precipitation, not to large-scale precipitation, which is more dominant outside the Tropics (not shown). Figure 5 shows the global patterns of these total precipitation changes in mean and variability. The changes of increased mean and variability with STOCH in the central tropical Pacific have a very similar pattern, which resembles the precipitation bias patterns when comparing the REF simulation to available observational data (Global Precipitation Climatology Project, GPCP: Adler et al., 2003). For the mean (Figure 5b), the increase in precipitation along the Equator is accompanied by a reduction in precipitation to the north (south) of around 5-10 • N (5-10 • S), which could be related to subsidence off the Equator (Gill, 1980). Both the reduction and the increase of precipitation compensate about 10% of the double ITCZ (intertropical convergence zone) mean precipitation bias around the Equator, westward of the Date Line (Lin, 2007). In the Atlantic, however, the changes mostly lead to a slight bias increase, and it is unclear whether this is a response to the changes in the Pacific, or whether it has a local cause. It should be noted that these biases are diagnosed by comparing a long climate simulation with CO 2 forcing held at 1990 conditions to a relatively short period of observations  in a transient climate.
The increase in tropical variability for both temperature (compare Figure 4 to Figure S3 in Appendix S1) and precipitation ( Figure 5) coincides with regions where the  Figure 5c,e to the changes via the FDDS coupling in Figure 5d,f). For temperature, the increased interannual variability pattern is broader and extends further into the extratropics. Quite remarkably, for precipitation, the increase in variability on interannual and inter-monthly time-scales compensates up to 50% of the variability bias in the Pacific when compared to GPCP observations. Again, there is a slight bias increase in the eastern tropical Atlantic, although much weaker than the improvements in the Pacific and only for monthly precipitation variability (Figure 5f). Regarding area-averaged changes in precipitation biases -both for the temporal mean and variability -STOCH leads to an improvement throughout the Tropics, where changes are most dominant (see Table 1). Especially noteworthy are the improvements in the tropical central Pacific between 5 • S-5 • N and 150 • W-150 • E where STOCH reduces RMSE by 8% for the mean and 29% for both monthly and interannual standard deviation. Including the entire global tropical band between 30 • S and 30 • N reduces these improvements to around 1% for both the mean and monthly standard deviation, and 3% for interannual standard deviation. This reduction, however, is simply a result of hardly any changes -neither substantial degradations nor improvements -between REF and STOCH outside of the tropical Pacific.

El Niño-Southern Oscillation (ENSO)
Interestingly, the monthly temperature (precipitation) variability biases in the Tropics can have the opposite sign in other models (e.g. in Community Climate System Model CCSM4: Berner et al., 2018, their figure 2) since current CMIP5 models strongly differ in the strength of the simulated variability (Kim and Yu, 2012;Liu et al., 2013). Moreover, including stochasticity via the atmospheric SPPT scheme in CCSM4 leads to a reduction of the simulated ENSO amplitude, to a level that is similar to observations Berner et al., 2018). On the one hand, this is very much in line with our higher-resolution experiments using AWI-CM: the 0.25 • tropical resolution of the REF configuration applied in this study, which allows us to resolve TIWs, was shown to also reduce the ENSO amplitude when compared to a lower-resolution (1 • ) simulation (Rackow et al., 2018), settling likewise at a more realistic level. On the other hand, the FDDS coupling scheme in this study apparently slightly reinforces the simulated ENSO ( Figure S8 in Appendix S1), which suggests that once a realistic ENSO amplitude is reached, this delicate balance may again be easily perturbed and the response to added stochasticity will be model-specific. While the Niño3.4 distribution of temperature anomalies had a reasonable representation already in the REF experiment when compared to the observed distribution ( Figure S8b in Appendix S1), probabilities in the right tail of the Niño3.4 temperature distribution increase due to FDDS in STOCH compared to REF. This is more realistic according to the observed Niño3.4 index ( Figure S8b, lower panel in Appendix S1). However, while the left tail had been slightly underestimated in REF, strong cold anomalies are now overestimated in the STOCH experiment. Due to higher probabilities in STOCH than in REF for stronger (both cold and warm) anomalies, the probability of encountering small temperature anomalies close to zero decreases accordingly, which is more realistic and a better fit to the observed peak of the distribution ( Figure S8b, upper panel in Appendix S1).
Possible changes in the mean ENSO are not the dominant contributor to the observed changes between REF and STOCH: the mean response of the 2m temperature (Figure 4a, STOCH-REF) does not resemble the ENSO pattern, only the variability pattern does (Figure 4b), which implies that our results are not due to a pure change in El Niño vs. La Niña occurrences (i.e. sampling) between the STOCH and REF However, since ENSO is the main mechanism for tropical interannual variability, the increases in 2m temperature variability due to the FDDS scheme appear to project onto the ENSO pattern.
To determine whether the changed ENSO variability after introducing the seasonally varying FDDS is more realistic, we compared the phase locking of Niño3.4 temperature anomalies to the seasonal cycle in the two experiments and in observations ( Figure 6). While the observed Niño3.4 variability tends to peak in boreal winter, a minimum occurs in boreal spring. The STOCH simulation shows stronger seasonal variations of temperature anomalies, with strongest anomalies during boreal winter and weakest SST anomalies during June-August. While the magnitude of the anomalies in boreal winter is a remarkable good fit to the observed magnitude, the timing of the minimum (broad minimum in June-August) is somewhat better than in REF (clear local minimum in August), but it is still off when compared to the observed phase locking where the weakest anomalies occur earlier in boreal spring. A second smaller local maximum in spring is an issue in both STOCH and REF and the added stochasticity in STOCH slightly worsens the spring variability.
However, the general shape of the seasonal ENSO phase locking is much better represented with STOCH compared to the more uniformly distributed and generally underestimated variability with REF. This shows another seasonally varying impact of the FDDS scheme as already mentioned in the context of Figure 2 (bottom panels).

3.4
Identified mechanism involving large-scale teleconnections Figure 7 illustrates and summarizes the changes in tropical Pacific variability and mean state caused by the FDDS coupling, for a variety of variables. Two-metre temperature variability increase causes convective variability in the vertical that spreads throughout the entire tropical troposphere (Figure 7a, left). This increased temperature variability and triggering of convective events causes a change in the tropical Pacific mean precipitation (Figure 7a, right) and affects the circulation in the free troposphere. A stationary Rossby wave response is triggered (Figure 7a, right) with circulation anomalies in the extratropics. The cyclonic (anticyclonic) anomalies coincide with positive (negative) large-scale precipitation anomalies along the West Coast of North America ( Figure 5b). Therefore, while the tropical precipitation changes are mostly due to changes in convective precipitation, the changes in the midlatitudes are associated with changes in large-scale precipitation and coincide with the circulation anomalies (see circulation schematic in Figure 7a, top right).
For the ocean component, the tropical precipitation increase causes an increased freshwater flux into the Pacific Ocean, freshening of the surface layer, an increase of stratification, and a shoaling of the mixed layer (Figure 7b, centre and right). Large-scale patterns of precipitation increase (decrease), surface freshening (salinization), and mixed-layer shoaling (deepening) are strongly correlated in the tropical Pacific up to about 30 • N/S. Beyond these latitudes, SST changes are getting more dominant and the relationship between freshening and shoaling is less clear. Generally, SST changes reflect the same patterns observed for 2m temperature, with much stronger signals in the extratropics and in mid-to high-latitudes. In the Tropics, the pattern is weaker but has a similar east (lower) to west (higher) dipole pattern that is also observed for the interannual temperature variability changes. Changes in sea surface salinity (SSS) variability in the Tropics are highly correlated to changes in precipitation variability (not shown), emphasising the direct impact the convective triggering through surface temperature variability has on mean stratification and its monthly to interannual variations. F I G U R E 7 Schematic of the identified mechanism. (a) In STOCH, the increased variability of the 2m temperature in the tropical Pacific extends into the free troposphere and is associated with increased convective precipitation. This triggers a Rossby wave train, as seen in circulation anomalies at 850 hPa height. The increased cyclonic (anticyclonic) activity coincides with positive (negative) large-scale precipitation anomalies. Stippling indicates significance in variability (F-test) and in the mean (rank sum test) at the 5% level. A level of 10% is used for stream function differences. (b) The upper-ocean salinity changes in the tropical Pacific (middle) are a direct response to the precipitation changes, and the fresher surface results in shallower mixed layer depths (right). Significant changes according to a rank sum test at the 10% level are outlined with black contours

Sea ice
While most of the changes discussed so far are caused by the changes in SST variability, the Arctic region also experiences an increase in surface field variability for SIT, SIC and SNT. This variability is most dominant for SSTs and sea ice concentration along the ice edge, but is also large in regions of strong SIT gradients moving away from the multi-year ice located along the coasts of Greenland and the Arctic Archipelago. Ocean-to-atmosphere heat fluxes are strongly influenced by ice cover and respond highly nonlinearly to changes between thick and thin ice as well as open ocean to ice-covered areas. Along the ice edge or in areas of strong ice thickness gradients the FDDS scheme samples "open ocean vs. ice covered states" and "thin ice vs. thick ice states". This tends to increase variability in turbulent heat fluxes. Due to the nonlinearity in the coupling, the mean heat flux to the atmosphere is also increased: thin ice (or low concentration) has a much larger positive impact on heat fluxes than thick ice (and high concentration) has on further decoupling the ocean from the atmosphere. Consequently, the Arctic experiences a warming (see the 2m temperature change in Figure 4a) due to increased turbulent heat fluxes along the ice edge (see Figure S10c in Appendix S1) and thus a corresponding general decrease in SIT (Figure 8a). For concentration, these changes are strongly focused on the ice edge (see Figure S5a compared to Figure S6a in Appendix S1). A consequence of the reduced ice thickness and concentration along the ice edge is that thick ice can also be more easily transported out of the central Arctic, contributing to the corresponding reduction in ice thickness close to the Pole (Figure 8a). Variability of SIT on the other hand is increased especially in the central Arctic, a direct consequence of the increased variability in surface fluxes (compare sea ice variability in Figure 8b and Figure 8c to the heat flux variability in Figure S10d in Appendix S1). Only along the ice edge, i.e. in marginal ice zones where both thickness and concentration values are low, is the effect a reduction of variability, in line with the reduced thickness and concentration.
For the Antarctic, as mentioned before, the changes are not directly caused by the grid resolution ( Figure 1). Teleconnections from the Tropics are the most likely candidate for the observed changes. In general, the Antarctic experiences an increase of SIC and SIT (Figure 8d) that is in line with changes in 10m winds (not shown). Especially in the Weddell and Ross Seas, anomalous wind is advecting cold air from the continent, leading to sea ice production. In addition to that, wind changes tend to spread out the ice further into the open Southern Ocean, especially in the eastern Weddell Sea. A local reduction of thickness is observed close to the Antarctic Peninsula and in a dipole pattern in the Ross Sea, somewhat coherent with circulation changes.
Variability patterns in Figure 8e,f are mostly confined to the Weddell Sea and the ice edge, generally in line with changes in the mean. As Antarctic sea ice is much thinner than in the Arctic ( Figure S7 in Appendix S1) and the anomalies are created by teleconnections, the amplitude of the changes is comparatively small but nevertheless significant. Interestingly, while the direct changes in surface field variability lead to a reduction of sea ice in the NH, the teleconnection effects from the Tropics lead to an increase in the SH. Since AWI-CM tends to overestimate (underestimate) ice thickness in the North (South), this inter-hemispheric difference generally seems to improve the simulations, although once again one should keep in mind that the simulations are using 1990 CO 2 forcing, and the model experiences considerable drift (and commitment warming) during the first few hundred years before the initialization of these experiments. As sea ice is a prime indicator for climatic changes, improvements or model degradations should be discussed with caution.

SUMMARY AND CONCLUSIONS
We investigate the impact of a new flow-dependent stochastic coupling scheme on the simulated climate of a global climate model (AWI-CM). The new scheme is based on the fact that the ocean resolution in current CMIP5 and upcoming CMIP6 climate models is typically larger than the atmospheric resolution. Taking AWI-CM as an example, the resolution of the FESOM ocean model grid is locally refined near coastlines, in the Tropics, and in the Arctic. Consequently, the grid-point ratio between FESOM and the atmospheric model ECHAM6 can reach values of up to 60:1. In such regions, the previous deterministic coupling scheme communicated spatially (weighted-)averaged surface fields of SST, SIC, SIT and SNT for each corresponding atmospheric grid box, thereby smoothing the surface field and losing part of the atmospheric subgrid information that is resolved by the ocean model. The stochastic FDDS scheme, on the other hand, samples from the resolved subgrid distribution, with probabilities weighted by the respective ocean cell areas. First, the idea of the scheme is to enhance and improve the simulation of surface field variability. Consequently, variations in atmospheric surface fluxes based on these fields are increased. Secondly, the scheme can be interpreted to account for coupling uncertainty with respect to surface fields. An averaging procedure as applied in the deterministic case tends to spatially smooth out surface fields and does not represent the actual state of the ocean as the atmosphere would see it with a comparable resolution. Consequently, the deterministic scheme is overconfident in the averaged surface field coupling procedure and does not provide any estimates of uncertainty in those fields. The FDDS scheme presented in this study on the other hand can subsample possible ocean states. The subsampling is followed by individual dynamical nonlinear responses of the system. The result is a climate model with increased surface field variations and, in the context of ensemble forecasts, an ensemble simulation which includes a measure of coupling uncertainty associated with a simple surface field averaging.
The effect of increased flux variations on climatic scales is investigated in a set of 9-year (interannual) integrations. The results show that especially the SST sampling in the tropical Pacific has a large impact on precipitation variability and its mean, increasing the amount of large precipitation events. The double ITCZ bias is strongly reduced for mean precipitation (around 10%) in our coupled model, and an underestimation of precipitation variability in the tropical Pacific is alleviated by a remarkable 50%. These bias reductions are based on a comparison of a present-day climate simulation (CO 2 levels fixed at 1990 conditions) with present-day precipitation data (Adler et al., 2003), and therefore conclusions need to take into account that different climate regimes might be compared. In the tropical western Atlantic the FDDS scheme leads to a slight precipitation bias increase. The cause of this could be studied in a future simulation where the FDDS coupling is only active in the Pacific, in order to determine a local or remote cause. However, the FDDS coupling generally improves long-standing model biases that have also been observed in other models (Lin, 2007).
While 2m temperature mean changes due to the FDDS scheme are mostly confined to the northern (increase) and deviation, and (c) monthly standard deviation of SIT for the Arctic. For every starting year, the 9-year time series were detrended before computing the variance, and then the root of the mean variances (averaged over the set of initializations) was computed to estimate the standard deviation (std); monthly std includes the seasonal cycle. Black stippling indicates significant changes at the 5% level according to a rank sum test (F-test) for mean (variability). (d-f) Same as (a-c), but for the Antarctic southern (decrease) mid-to high-latitudes -coinciding with respective changes in sea ice -variability changes on monthly and interannual time-scales are most dominant in the tropical Pacific, in line with the precipitation changes. Especially in the tropical Pacific, an increased surface temperature variability and higher probability of SSTs ≳27.5 • C lead to enhanced triggering of convection and convective precipitation with the FDDS scheme, mostly for extreme precipitation events ≳10 mm⋅day −1 . The diabatic heating causes anomalous Rossby waves emitted from the tropical Pacific (Gill, 1980), which propagate to midlatitudes and cause changes in large-scale precipitation off the West Coast of North America. Changes in Antarctic sea ice (increase) and 2m temperature (decrease) are related to the circulation changes originating in the Tropics, as ocean-to-atmosphere resolution ratios over the Southern Ocean are close to one, and the FDDS coupling scheme is therefore not active in this area. In the Arctic, however, the stochastic sampling of SST, SIC, SIT and SNT leads to both decreased mean SIT and SIC -but increased variability.
These results show that the stochastic FDDS scheme, while acting on relatively short time-scales (6 h coupling), has large impacts not only on monthly and interannual variability in the atmosphere and ocean but also on the mean state. The prime example is the tropical Pacific, which plays a crucial role in the global climate, where the locally increased sea surface state variability triggers teleconnections that affect the large-scale circulation. The study highlights the necessity for further investigations of an accurate representation of surface field (and air-sea flux) variability and coupling. However, there is room for improvements. While the scheme is able to sample the spatial sea surface variations within one atmospheric grid box, temporal variations are still ignored during the 6 h coupling steps, especially since some of the actual temporal variation statistics are resolved more accurately by the atmospheric model with its relatively smaller time step. Changing the coupling step to a more frequent exchange of fluxes and fields (e.g. 1 h, as is now more common in climate models) could necessitate a retuning of the FDDS scheme so that the values of surface fields do not change too rapidly. This could encompass the incorporation of some form of temporal correlation between the perturbations, which are currently drawn independently in time.
In the current configuration, some variations along the sea ice edge might be overly amplified in a physical sense: for example, the SIC distribution within one atmospheric grid box near the ice edge might become bimodal and the perturbed values might therefore switch from very large to very low (or vice versa) between two coupling steps. Similarly, across an oceanic front the distribution of ocean cells within an atmospheric box might also have bimodal characteristics. By switching the randomly chosen ocean cell at each 6 h coupling time step, these spatial variations can lead to very strong temporal variability of surface fields that might be deemed overly strong. A potential way forward is to include further temporal correlation into the stochastic scheme, e.g. by only considering neighbouring ocean cells in the stochastic selection process for the next coupling time step. This random walk would prevent too strong surface field changes within 6 h in most instances. Some of the above-mentioned ideas and aspects will be investigated in future studies. All of these aspects will need to respect local conservation laws for energy and mass. This will be more difficult to achieve when adding perturbations to actual fluxes, where previous studies have shown that applying the law of large numbers generally ensures a more or less balanced zero average mean energy or mass injection by the perturbations (e.g. Williams, 2012).
In general, the stochastic FDDS coupling scheme introduced in this article is of relevance to other models with high ocean-to-atmosphere resolution ratio, including regular CMIP models. The scheme falls into the category of process-oriented schemes, which are currently experiencing a lot of interest in the community (Leutbecher et al., 2017;Zadra et al., 2017). The discrete approach presented here does not rely on a continuous (finite-element) representation of the ocean surface and is therefore generally applicable to other models. Those comprise current-generation CMIP5 models with a telescoping of grid boxes in the Tropics (e.g. Delworth et al., 2006) or next-generation CMIP6 models with ocean resolutions of typically 0.25 • and coarser atmospheric resolution. Those models generally exhibit a more or less pronounced double ITCZ bias, one of the most prominent biases of current climate models (Lin, 2007). As the improvement in the precipitation mean (variability) bias in this study is on the order of 10% (50%), the FDDS coupling introduces a new avenue as to how the important double ITCZ bias and potentially also other biases could be tackled. In particular, the seasonally varying FDDS also considerably improves the phase locking of Niño3.4 SST anomalies to the seasonal cycle. Numerical simulations on grids that apply high ocean-to-atmosphere resolution ratios in other areas of the globe where mesoscale activity is high, e.g. along western boundary currents or in the Southern Ocean (Sein et al., 2017), could also benefit from the FDDS scheme, as the local improvements in the northern high-latitudes suggest. Whether the results presented here depend on the atmospheric model's convection scheme and how these results generally translate to other climate models will be investigated in the future. Furthermore, the results of this study need to be tested in the context of other coupling strategies (e.g. with higher coupling frequencies, or accounting for surface ocean currents in flux computations). However, the implementation of the new scheme can be easily added as an option to couplers such as OASIS3 (Valcke, 2013) or to the coupling interface of any ocean model. Except for the random number generator to select the ocean cells, there are no additional computational costs involved, which makes the scheme computationally cheap and readily applicable to other model set-ups.