S Close

T Penduff

S Speich

J.-M Molines

A means of estimating the intrinsic and atmospherically-forced contributions to sea surface height variability applied to altimetric observations

Drawing on a 50-member ocean ensemble hindcast, the magnitude and characteristic temporal and spatial scales of intrinsic and forced sea surface height (SSH) variability are evaluated over a 37-year period. The intrinsic and forced contributions derived from the ensemble are found to have similar temporal spectra, but different characteristic spatial scales. These results suggest that, with an appropriate choice of cutoff scales, simple spatial filtering can be used to estimate the forced and intrinsic contributions given either a single model run, or an observational data set. The method is tested using a single member drawn from the ensemble, before being applied to the observed altimetric record. Two sample applications with relevance to large-scale climate are used to illustrate the method's potential utility. Firstly, the long-term trends calculated from the total and recreated forced components using the altimetric record are compared and local differences highlighted. Second, the recreated forced SSH is shown to covary with the North Atlantic Oscillation at seasonal time scales in regions where no such influence can be found using the original SSH signal. Some limitations and uses for which the method may prove unsuitable are also briefly considered.

Introduction

The variability and long-term tendency of sea surface height (SSH) is of wide interest as a climatic indicator, often being interpreted as a metric of the integrated heat content of the water column [START_REF] Church | Climate Change 2013: The Physical Science Basis, Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, chap. Sea Level Change[END_REF], and the advent of satellite altimetry has ensured that SSH is one of the properties for which the longest, continuous, global observational records are available. Its variability reflects the redistribution of oceanic mass, heat and salt, and SSH thus further provides insight into ocean dynamics across all resolved scales.

In the presence of mesoscale processes, 1 the variability of SSH, and indeed that of other dynamically-influenced quantities, at any given location can be considered to comprise two components: a deterministic contribution, occurring in response to an applied forcing, and a non-deterministic part, arising predominantly from small-scale (i.e. mesoscale and smaller), turbulent processes. In contrast to the hypothetically predictable, deterministic, atmospherically-forced response, this non-deterministic contribution is inherently random.

Although intrinsic variability is often associated with the mesoscale, small-scale processes are not its sole means of generation; indeed, numerous studies have demonstrated the existence of decadal-scale oscillatory behaviour associated with large-scale baroclinic instability in numerical models (e.g. Colin de Verdière and [START_REF] De Verdière | Baroclinic instability: an oceanic wavemaker for interdecadal variability[END_REF][START_REF] Pierini | A Kuroshio extension system model study: decadal chaotic self-sustained oscillations[END_REF][START_REF] Huck | Multidecadal variability of the overturning circulation in presence of eddy turbulence[END_REF]. Where intrinsic variability has its genesis in mesoscale processes, the inchoate signal may evolve to yield coherent behaviour over much larger spatial and temporal scales. In the particular case of SSH, [START_REF] Sérazin | Intrinsic Variability of Sea Level from Global Ocean Simulations: Spatiotemporal Scales[END_REF][START_REF] Sérazin | Quantifying uncertainties on regional sea level change induced by multidecadal intrinsic oceanic variability[END_REF] suggest that intrinsic variability can contribute non-negligibly to the total variance at horizontal spatial scales > 12°and at up to multi-decadal time scales; in a similar vein, [START_REF] Ocaña | Stochastic secular trends in sea level rise[END_REF] analyse long-term tendencies and suggest that the hypothesis that recent trends in SSH may be due to stochastic processes cannot, at present, be rejected. Using a model ensemble, [START_REF] Llovel | Contributions of Atmospheric Forcing and Chaotic Ocean Variability to Regional Sea Level Trends Over 1993-2015[END_REF] suggest further support for these results with respect to regional SSH trends over the period of the satellite record. In a state-estimate-based analysis, [START_REF] Forget | The partition of regional sea level variability[END_REF] also find that mesoscale (< 3°spatial scale, classified therein as intrinsic) variability dominates over large-scale variability at time scales up to 3 years, and contributes significantly to the total variability at time scales longer than this.

Whilst the link between SSH and the large-scale circulation is known to be complex [START_REF] Kanzow | Basinwide integrated volume transports in an eddy-filled ocean[END_REF][START_REF] Hirschi | Sea surface height signals as indicators for oceanic meridional mass transports[END_REF][START_REF] Clément | Vertical structure of eddies and Rossby waves, and their effect on the Atlantic meridional overturning circulation at 26.5[END_REF], considering integrated metrics, a number of studies have also shown that a substantial fraction of the variability of the Atlantic meridional overturning circulation is intrinsic at multiple locations and time scales [START_REF] Hirschi | Chaotic variability of the meridional overturning circulation on subannual to interannual timescales[END_REF][START_REF] Thomas | Eddy-induced variability of the meridional overturning circulation in a model of the North Atlantic[END_REF][START_REF] Grégorio | Intrinsic Variability of the Atlantic Meridional Overturning Circulation at Interannual-to-Multidecadal Time Scales[END_REF][START_REF] Leroux | Intrinsic and Atmospherically Forced Variability of the AMOC: Insights from a Large-Ensemble Ocean Hindcast[END_REF]. Further studies have similarly inferred an important role of intrinsic variability in determining multiannual-to decadal-scale changes in the Southern Ocean circulation [START_REF] O'kane | Decadal variability in an OGCM Southern Ocean: Intrinsic modes, forced modes and metastable states[END_REF][START_REF] Wilson | Forced and intrinsic variability in the response to increased wind stress of an idealized Southern Ocean[END_REF]. These studies thus suggest that intrinsic variability can have an impact at climatic temporal and spatial scales, and suggest further motivation to quantify its effects.

Although intrinsic variability is presumed to exist in observational data as well as in eddy-permitting numerical model output, the estimation of its contribution to the total variability of a given quantity necessitates the use of a(n eddy-permitting) numerical model, since some means must be provided by which either the forced or intrinsic component can be isolated. Previous studies have generally addressed this problem by using one of two methodologies: either a realistic atmospherically-forced simulation may be compared with one subject to only climatological forcing (e.g. [START_REF] Penduff | Sea level expression of intrinsic and forced ocean variabilities at interannual time scales[END_REF], or, alternatively, and similarly to approaches long used in numerical weather prediction, an ensemble of simulations with perturbed initial conditions may be employed (e.g. [START_REF] Hirschi | Chaotic variability of the meridional overturning circulation on subannual to interannual timescales[END_REF][START_REF] Roberts | Atmosphere drives recent interannual variability of the Atlantic meridional overturning circulation at 26.5[END_REF][START_REF] Bessières | Development of a probabilistic ocean modelling system based on NEMO 3.5: application at eddying resolution[END_REF]. The former approach permits the comparison of statistical metrics in the temporal dimension (e.g. the temporal standard deviations of the climatological and time-varying forcing runs may be compared at a given spatial point). However, the forced and intrinsic signals cannot themselves be derived from the model outputs using this particular approach, since the intrinsic variability, being inherently random, differs from simulation to simulation. In contrast, the latter, ensemble, approach introduces an additional dimension (i.e. that of the multiple simulations) and permits both an estimation of the forced (ensemble mean) and intrinsic (residual after subtraction of the ensemble mean from each member) signals, as well as an estimation of the statistical properties in a temporal, as well as in a spatial sense. For example, where the former approach permits comparison of the temporal-dimension standard deviation at a given spatial point, thus providing a metric of the overall difference in variance at the specified location, using the ensemble approach, the standard deviation of the component ensemble members may be calculated for the same point at each instant in time, yielding a time-varying estimate of the intrinsic contribution. The estimate of the atmospherically-forced contribution to the total variability that emerges from this methodology represents the response of the ocean to both anthropogenic variability and also the internal variability of the atmosphere at all resolvable time scales.

In this work, we examine the statistical properties of the forced and intrinsic contributions to total SSH variability, drawing on an ensemble of ocean hindcasts, and, based on these results, present a method for estimating the contributions of the two components given a single timevarying experiment (either model, or observation-based). In Section 2, the ensemble experiment and methodology used in its analysis are presented. Some statistical properties are presented in Section 3: directly exploiting the ensemble dimension, we calculate the spatial, temporal and spatiotemporal spectral coherence of the forced and intrinsic components with respect to the total signal. Motivated by these results, in Section 4, a simple method of estimating the two contributions is proposed, verified, and some applications discussed. We conclude with a discussion and summary in Sections 5 and 6.

Data and methods

Model data

The data used in the analysis are outputs from the OCCIPUT (OceaniC Chaos -ImPacts strUcture predicTability project) experiment, the motivation and methodology behind which are described more fully in [START_REF] Penduff | Ensembles of eddying ocean simulations for climate, CLIVAR Exchanges[END_REF] and [START_REF] Bessières | Development of a probabilistic ocean modelling system based on NEMO 3.5: application at eddying resolution[END_REF] respectively. Briefly, the data set comprises a 50-member ensemble of global ocean-sea ice hindcasts, covering the period 1960-2015, and driven with historical atmospheric forcing using the Drakkar Forcing Set DFS5.2 [START_REF] De Verdière | Baroclinic instability: an oceanic wavemaker for interdecadal variability[END_REF]. The model configuration is based on NEMO (Nucleus for European Modelling of the Ocean) v.3.5 [START_REF] Madec | NEMO ocean engine, Note du pôle de modélisation 27[END_REF], with 0.25°h orizontal resolution and 75 vertical levels. The outputs used in this work are successive 5-day averages. The ensemble is generated by first producing a single simulation, which is spun-up over a 21-year period; at the end of this period, the 50 ensemble members are started from this simulation, and a weak, time-varying stochastic perturbation with locally-dependent amplitude [START_REF] Brankart | Impact of uncertainties in the horizontal density gradient upon low resolution global ocean modelling[END_REF][START_REF] Brankart | A generic approach to explicit simulation of uncertainty in the NEMO ocean model[END_REF] applied in the density equation in each member over a 1-year period to generate the ensemble dispersion. The perturbation is switched off after this time, and the 50 members then experience identical forcing over the 56 year integration, but depart from slightly differing initial states. A complete description of the methodology and associated technical aspects can be found in [START_REF] Bessières | Development of a probabilistic ocean modelling system based on NEMO 3.5: application at eddying resolution[END_REF]. Although the model experiment begins in 1960, in this study only the period 1979-2015 will be analysed. This decision is made in order to obtain a period over which the atmospheric forcing is consistent, since prior to the satellite era a number of quantities (notably radiative and freshwater fluxes, see Table 1 in [START_REF] De Verdière | Baroclinic instability: an oceanic wavemaker for interdecadal variability[END_REF] are represented by a climatological cycle; in preliminary analyses, it was found that this may have an effect on the ensemble spread, and this early period is thus excluded from this study, yielding an effective spin up time for each member of 41 years. Prior work [START_REF] Bessières | Development of a probabilistic ocean modelling system based on NEMO 3.5: application at eddying resolution[END_REF] has found that the ensemble spread saturates rapidly, within a few months or years of the stochastic parameterisation being switched on, and so the spread is approximately stable (but subject to fluctuations on sub-up to multiannual timescales) over the period considered here.

Observational data

In addition to the model data described above, SSH observations derived from altimetry are also employed. The data set comprises the 0.25°Merged-Gridded Sea Level Anomalies and Absolute Dynamic Topography heights and currents in Delayed-time (DT-MSLA-H Global) product distributed by the Copernicus Marine Environment Monitoring Service. The product combines data from all existing altimeter missions to yield a consistent and homogeneous record, and the data set is available from 1993 onwards. In this work, we employ all years for which data were available in all months at the time of analysis, i.e. 1993-2017.

Preprocessing of model data

For models that conserve volume, rather than mass, the long-term signal in the simulated SSH is not physically meaningful [START_REF] Greatbatch | A note on the representation of steric sea level in models that conserve volume rather than mass[END_REF]. In a first step, the global mean SSH signal is thus subtracted from each member at each time step; as detailed in Greatbatch, 1994, this corrects for the spurious extra terms introduced by the use of the Boussinesq approximation in volume-conserving models, and the solution obtained according to volume-conserving kinematics is then exactly equivalent to that obtained from mass-conserving kinematics. In a second step, to remove long-term behaviour with periods longer than that resolved by the 56-year run (including any effect of model drift), the data are detrended prior to analysis. This is performed here by first applying a LOWESS (LOcally WEighted Scatterplot Smoothing) filter [START_REF] Cleveland | Robust locally weighted regression and smoothing scatterplots[END_REF] with 45% span to the yearly mean data, where the span is determined by optimisation of the mean average error on a sample set of 200 points, selected with even spacing across the model grid. A 5th order spline is then fitted to this smoothed form; the spline is constrained to pass through all points, and this function used to estimate and remove the non-linear trend at each time step for the 5-daily SSH data. The residual mean state is approximately zero after application of this procedure. Finally, after detrending, the mean seasonal cycle, calculated over 1979-2015, is removed from each member.

Definition of the forced and intrinsic contributions

The atmospherically-forced signal is estimated as the mean of all ensemble members at every time step (yielding a 37-year time series at each spatial point). The intrinsic signal is then defined as the difference between the original signal of each member and this ensemble mean (yielding a 37-year time series at each spatial point for each of the 50 members of the ensemble). By construction, the average of the intrinsic signal is thus zero both at each time step and also over the total population. Although this implies that the intrinsic signal and forced signal are uncorrelated when considering the full ensemble, note that any given member of the ensemble considered in isolation is not bound by this constraint, and it is thus possible for non-negligible correlation to exist between the forced and intrinsic variability of any given member.

Notation

As noted above, the use of an ensemble introduces an additional dimension that will be exploited in this study. Means and variances may thus be calculated in the ensemble dimension, as well as in the temporal dimension. In the following, we will broadly follow the notation used in [START_REF] Leroux | Intrinsic and Atmospherically Forced Variability of the AMOC: Insights from a Large-Ensemble Ocean Hindcast[END_REF], and the temporal mean will thus be denoted using the overbar symbol (x ¯), and the ensemble mean using brackets (〈 〉 x ). The temporal variance will be denoted σ x ( )

2

, and the ensemble variance ∊ x ( )

2

. Additionally to the terminology introduced in Leroux et al. ( 2018), we will also refer to the population variance, s x ( )

2

, which is calculated over both the temporal and ensemble dimensions, and where the ensemble is taken to represent the population of all possible states (although we acknowledge that, being a finite-sized sample, this is not strictly the case).

Where a quantity relates to an individual member, the subscript m will be used. The SSH for a given member will thus be written η t ( )

m
, for example, where t ( ) indicates that η m is a function of time. The ensemble population (i.e. the set of all members of the ensemble) will be indicated by the use of bold, thus the ensemble population of SSH will be written η m t ( , ), where m indicates that η is a function of ensemble member (in addition to time, t, as for each member).

In discussing the total, forced and intrinsic contributions to the variability, the subscripts T, F and I will be used respectively. ∊ η m t ( ( , ))

T 2 thus denotes the ensemble variance of the total (forced + intrinsic: = + η η η T F I ) signal for the ensemble population, for example, and σ η t ( ( ))

I 2 m denotes the temporal variance of the intrinsic contribution for a given member m.

When considering the dispersion of the ensemble (i.e. the spread) in isolation, we will generally use the standard deviation, to obtain a quantity with the same units as the original signal. In contrast, when comparing variability of the forced and intrinsic contributions either against each other or against the variability of the total signal, where the comparison involves the calculation of the variability of one quantity as a fraction of the other, we will employ rather the variance, since the use of the standard deviation in this case would imply that the resultant fraction would not vary linearly.

The ensemble variance of a time-varying quantity remains a function of time. The ensemble variances of either the total SSH or of the intrinsic contribution, ∊ η m t ( ( , ))

T 2 or ∊ η m t ( ( , )) I 2
, are thus timevarying (and are equivalent). In contrast, the forced contribution is, by definition, the same for all members, and thus it is only possible to consider the temporal variance for this quantity:

η σ t ( ( )) F 2
. To compare the intrinsic and forced variances, the temporal variance of the forced contribution will generally be compared with the population variance of the intrinsic contribution, yielding a single value at each gridpoint, permitting the calculation of a simple fraction describing the relative magnitudes of the two.

Reliability of the forced estimate

In cases where, at a given instant, the ensemble spread is very large, and the forced signal relatively weak, the forced signal might be considered to be ill-defined where the ensemble mean has an amplitude that is smaller than the ensemble standard error at a given time and location. To evaluate the reliability of the forced estimate, we thus calculate

⎜ ⎟ ⎛ ⎝ ∊ ⎞ ⎠ - η η t m t n | ( )| ( ( , )) F I 1
for each time, t, where = n 50 is the number of members comprising the ensemble. The percentage of time that this quantity is greater than 1 is then evaluated, to yield a metric of the accuracy of the estimate of the forced contribution (shown in the supplementary material, Fig. SI1). With the exception of very few regions, which correspond with those associated with a strong intrinsic variability (discussed below), the forced contribution is clearly distinguishable from the intrinsic spread over 75% of the time in all other areas, and thus considered to be meaningfully defined outside of these areas.

Magnitudes and characteristic scales of forced and intrinsic SSH variability

Standard deviations of the forced and intrinsic variability

In a first step, the standard deviations of the forced and intrinsic variability are calculated. Both quantities are computed as the square root of the variances of their respective data sets; for the forced contribution this is equivalent to simply calculating the temporal standard deviation, η σ t ( ( )) F , whilst for the intrinsic contribution the ensemble dimension is also included in the calculation: η s m t ( ( , ))

I

. The total variance is then simply the sum of the variances of these two contributions, since the total population forced-intrinsic covariance is zero by construction, and the total standard deviation, η s m t ( ( , )) T , the square root of this quantity. Maps of the total, forced and intrinsic signals' standard deviations are shown in Fig. 1.

The patterns shown for the total signal and intrinsic component (Fig. 1a/c) are very similar to those presented by the previous studies of [START_REF] Penduff | Sea level expression of intrinsic and forced ocean variabilities at interannual time scales[END_REF]Penduff et al. ( , 2015)), in which pairs of atmospherically-forced and climatological model runs were compared: in both cases, the variance is elevated in the major current systems, with local maxima in the Agulhas region, the Brazil-Falklands Confluence region (hereafter BFC), the Gulf of Mexico, the Gulf Stream, the Kuroshio and at localised points along the path of the Antarctic Circumpolar Current (hereafter ACC) in the Pacific sector. The spatial patterns of the total and intrinsic signals show a very strong resemblance, being characterised primarily by the elevated variances found in these regions.

The variability of the forced contribution (Fig. 1b) has a greatly reduced range relative to that of the total and intrinsic contributions2 . The patterns of enhanced variability emerging from this pointwise analysis are generally associated with much greater spatial scales than in the intrinsic case (e.g. spanning the tropical Pacific, the Indian Ocean or the Bellingshausen Basin), consistent with the large-scale nature of atmospherically-forced modes of variability. Although a number of smaller scale patterns can also be discerned (e.g. in the Agulhas region, the Gulf of Mexico and along the length of the Gulf Stream), these are generally co-located with regions of intense intrinsic variability, and coincide with regions where the estimate of forced variability obtained from the model is not well-defined (shown by white contours in Fig. 1b).

Relative magnitudes of forced and intrinsic variability

Having evaluated the standard deviations of the forced and intrinsic contributions, the relative magnitudes of the forced and intrinsic SSH variability are now compared. The ratio of the forced to the total variance is defined as

= η η σ t s m t RF ( ( )) ( ( , )) F T 2 2
and is shown in Fig. 2 for all time scales (Fig. 2a), subannual time scales (Fig. 2b) and interannual and greater time scales (Fig. 2c). The temporal filtering is performed using a Kaiser-Bessel filter with a cutoff period of 1 year.

Although the forced contribution dominates over the intrinsic over large regions, particularly at interannual and greater time scales (Fig. 2c), there nevertheless exist local regions where the SSH variability is almost entirely intrinsic. This is particularly the case for strongly eddying regions, such as the Agulhas Return Current region, the BFC, the Gulf of Mexico Loop Current, the core of the Gulf Stream and local regions along the path of the ACC, where the variability remains almost entirely intrinsic even at interannual time scales.

Strong similarity can be noted between RF calculated over all time scales (Fig. 2a) and that calculated over subannual time scales (Fig. 2b), highlighting the strong role of high frequency variability in setting the total variance. Notably, at subannual time scales, large regions of the Atlantic sector of the Southern Ocean are dominated by the intrinsic contribution. Basin scale means for the major basins (the Arctic excluded) are detailed in Table 1. There is a strong contrast between the high latitude Southern Ocean, which has the lowest fractional contribution of intrinsic variability, and the northern sector, containing the ACC, which has the highest. Considering the Southern Ocean, the sea ice variability in the model is reasonable compared to observations; however, whilst the ensemble spread of certain other surface quantities (notably sea surface temperature) exhibits a notable seasonal cycle, no such distinct seasonal cycle is found in the SSH over most of this region. This may potentially indicate a non-negligible role of the subsurface water column in determining the SSH evolution in the Southern Ocean; however, an alternative hypothesis, given the greater relative effect of meltwater injection on salinity than on temperature (e.g. [START_REF] Close | Control of Mode and Intermediate Water Mass Properties in Drake Passage by the Amundsen Sea Low[END_REF] and the dominance of salinity in the equation of state in the high latitude Southern Ocean, could be that the time scale necessary to mix the meltwater input sufficiently for the salinity to exhibit notable divergence amongst the members (and thus for any such divergence to emerge as a signal in the SSH) is longer than the time scale associated with the seasonal cycle of the sea ice. Of the remaining regions, the South Atlantic has a slightly higher contribution of intrinsic variability on average compared to the rest; this originates from the strongly intrinsic zone lying along the southern boundary of the domain, suggesting the possible incursion of this signal into the South Atlantic.

At interannual and greater time scales (Fig. 2c), RF is, globally, much weaker, and the areas of the regions identified above as being strongly dominated by intrinsic variability are notably reduced. Certain regions nevertheless retain their strongly intrinsic character; referring to Fig. 2c, the core of the ACC and, in particular, the Agulhas Return Current and BFC stand out in this regard. These results thus hint that certain localised regions may have a highly intrinsic character across a wide range of time scales; this dependence of RF on frequency will be explored further below.

Temporal and spatial scales of the forced and intrinsic variability

Temporal coherence

The analyses of the absolute and relative values of the variance presented above describe the total SSH variance at any given point (evaluated in a temporal, and temporal/ensemble sense). The patterns revealed by the analyses suggest both a certain spatial coherence (through the large-scale patterns evident in the point-wise analyses), and also a temporal dependency of the magnitudes of the forced and intrinsic contributions. To understand the relative importance of the temporal and spatial scales at play, the spectral coherence is calculated. For two signals, x and y, the coherence is defined as:

× CSD PSD PSD | | x y x y , 2
where CSD indicates the cross-spectral density between the two signals, and PSD the power spectral density of each given signal. In all of the analyses discussed below, the associated phase:

I R ⎜ ⎟ = ⎛ ⎝ - ⎞ ⎠ ϕ m CSD e CSD arctan ( ) ( ) x y x y , ,
is found to be close to zero, and is hence neither shown nor discussed further.

We begin by considering the temporal coherence between the forced and total SSH signals. Maps of the coherence for a selection of frequency ranges are shown in Fig. 3. The coherence between the intrinsic and total signals is essentially simply 1 minus these values, since the forced-intrinsic coherence is weak and has a diffuse spatial distribution; in the interests of brevity, the intrinsic-total and forced-intrinsic coherences are thus not shown. In all cases shown below, the power and cross spectral densities are estimated using Welch's method, with a 50% overlap between segments and employing a Hann window as the windowing function. Spectral coherence is then calculated using the averaged periodograms. Full details of the number of segments used and significance levels, calculated using Monte-Carlo simulation, are given in the supplementary information (Section SI2.2). Whilst spectral coherence is used throughout in the following analyses, RF has also been calculated for all cases, and found to show consistent results. Where the forced signal is strongly coherent with the total signal, it thus also has a stronger magnitude than the intrinsic signal in the spatiotemporal range under consideration (and vice versa). Hence, in the interests of brevity, only the coherence will be shown here.

At very high frequencies, Fig. 3 (top row) suggests that large regions of all basins are predominantly atmospherically-forced, with little intrinsic contribution. Nevertheless, the rapidity with which the areas of the intrinsic regions spread is notable even between frequencies with 10-(Fig. 3a) and 20-day periods (Fig. 3b), where the coherence between the total and forced contributions diminishes markedly in the equatorial tropics in all three basins, and the localised regions that show very low coherence at 10-day periods increase rapidly in area over these time scales. The short time scales associated with the zonally large-scale changes in the equatorial tropics suggest tropical instability waves to be a likely cause (e.g. [START_REF] Willett | Eddies and Tropical Instability Waves in the eastern tropical Pacific: a review[END_REF][START_REF] Vialard | Sensitivity of Pacific Ocean Tropical Instability Waves to Initial Conditions[END_REF].

At frequencies with periods greater than one month and up to one year (Fig. 3, middle row), the patterns of coherence between the forced and total signals are largely established and show comparatively little evolution relative to the sub-monthly time scales. The primary changes during these periods occur in the subtropics, particularly in the Pacific and Indian Oceans, where the 30-60°N/S ranges become more intrinsically-dominated at lower frequencies. In contrast, the equatorial region begins to transition back towards atmospherically-forced behaviour over these time scales, with the predominantly forced equatorial band widening between the 1-3 and 6-12 month periods (Fig. 3d/f). Globally, the coherence between the forced and total signals is at its weakest in the 6-12 month period range.

At interannual and greater time scales (Fig. 3, bottom row), the coherence between the forced and total signals begins to increase. Although the intrinsic contribution is still important at 1-2 year time scales (Fig. 3g), the total signal begins to show greater coherence with the forced contribution in the 2-4 year range (Fig. 3h), and the intrinsic contribution is greatly restricted at time scales greater than 4 years (Fig. 3i), consistent with the results shown previously in Fig. 2. Finally, considering the results obtained over all frequencies, over certain regions, essentially corresponding to those that show weak coherence (< ca. 0.3) at 10 day periods (Fig. 3a), the coherence between the forced and total signals remains low at all frequencies; this highlights the dominance of the intrinsic contribution at all time scales in these localised areas.

Pre-empting the spatiotemporal analysis presented below, the global minimum in temporal coherence between the forced and total signals noted in the 6-12 month range in Fig. 3 is investigated further by briefly analysing the evolution of the spatial distribution of the intrinsic power, averaging the spatiotemporal intrinsic PSD, and the forced-intrinsic CSD over the frequency ranges shown in Fig. 3 (not shown). In all basins, the emergence of the peak in power in the mesoscale wavelength range begins in the 1-3 month time period range, becoming well-established in the 3-6 month range. The maximum power is then attained in the 6-12 month range, consistent with the global synopsis of Fig. 3f, before dying down over longer time periods. There is some evidence that the spatial scales associated with this peak in energy also evolve over time, where the peak in mesoscale energy shifts from ~300 km to ~450 km, stabilising at this higher value, between the 1-3 and 6-12 month analyses, which could represent the (spatial) growth of the intrinsic perturbations over time. Meanwhile, the CSD show peaks in energy at wavelengths ca. 250-300 km, with a secondary maximum around the 875 km wavelength. The amplitude of the CSD is large relative to the amplitude of the forced PSD within the wavelength range < ~700 km, and large relative to the amplitude of the intrinsic PSD within the wavelength range > ~1000 km. The forms of the CSD are compatible with the notion that some energy transfer from the intrinsic to the forced component may take place at length scales ca. 700-1000 km, and thus that interactions leading to the synchronisation of the intrinsic component may occur beyond certain characteristic length scales (consistent, for example, with the hypothesis of barotropisation outlined in [START_REF] Scott | Direct Evidence of an Oceanic Inverse Kinetic Energy Cascade from Satellite Altimetry[END_REF]). However, dedicated analyses of spectral flux would be necessary to verify these ideas robustly, and have not been performed here.

The wide-reaching nature of the intrinsic influence on time scales up to and including the multiannual found here suggests that it is likely to be difficult to separate the forced and intrinsic contributions based on a purely temporal criterion. Indeed, it is apparently only at very high (and, in certain regions, low) frequencies that, at a global scale, the total signal could be reasonably approximated to be representative of the forced contribution. In a second step, the spectral coherences are thus evaluated in a spatial sense, with the aims both of better understanding their spatial imprints and also of investigating the potential separability of the forced and intrinsic contributions.

Spatial coherence

The spatial spectral coherences are calculated over the regions defined in Table 1, and shown in Fig. 4. The spectra are estimated using a Hann-windowed 2D Fourier transform (as detailed above), and thus exhibit central symmetry (the number of segments and 95% significance levels are detailed in the supplementary information, Section SI2.2). A general pattern, in which the forced variability shows relatively strong coherence with the total variability at both small and large scales, but weak coherence over the mesoscale range, is common to all the analyses. Nevertheless, there are some clear differences in the form of the spectra amongst the various basins: notably, the North and South Atlantic, North Pacific and Indian Oceans show stronger anisotropy in the distribution than the South Pacific and, particularly, Southern Ocean, where the coherence distribution is almost perfectly symmetrical. The more asymmetric distributions show weak forced-total spectral coherence (and, logically, strong intrinsic-total spectral coherence; not shown) along the positive-positive/negative-negative (i.e. northwestward/southeastward) zonal-meridional wavenumber diagonal. This behaviour is consistent with a strong impact of the westward propagation of eddies on the variance in the mesoscale wavenumber range (where the direction of propagation has been diagnosed using the 3D PSD of the intrinsic component, used in the spatiotemporal analyses presented in the following section). It is hypothesized that the relative impact of this mechanism is reduced in the presence of the eastwardflowing ACC, leading to the more symmetric distributions evident in the Southern Ocean and South Pacific; inversely, in the South Atlantic, the strong asymmetry of Fig. 4b highlights the dominant role of Agulhas Rings in setting the variance in the mesoscale wavelength range in the South Atlantic.

Spatiotemporal coherence

Finally, we undertake a full 3D spatiotemporal analysis of the coherence, shown in Fig. 5. Although, as previously shown in Fig. 4, the zonal-/meridional-wavenumber spectra are not perfectly isotropic, the radially-averaged wavenumber-frequency coherence is shown here to simplify the presentation of the results.

The broad form of the spatiotemporal coherence between the forced and total signals is rather consistent across all regions examined here. The forced signal is coherent with the total signal at both very short and also long wavelengths, and, in most regions, also at very high frequencies (< ca. 1 month) over a wide range of wavelengths (Fig. 5). The region of weakest coherence (and also the strongest coherence between the intrinsic and total signals, not shown) is centred approximately between wavelengths of 128-512 km, and 1 month to 4 year periods. The Indian Ocean (Fig. 5e) appears to show weaker forced/ total coherence at high frequencies than the other regions, with the characteristic intrinsic range being shifted towards higher frequencies relative to the other basins.

The spatiotemporal coherences can be interpreted in the context of the separate spatial and temporal analyses, where many of the most salient features of the spatiotemporal representation were previously noted. The strong influence of the forced contribution on the total signal at sub-monthly frequencies is again found here, where the localised strongly intrinsic regions identified in Fig. 3a-b nevertheless yield a significant coherence between the intrinsic and total signals over their characteristic length scales within this frequency range.

Consistent characteristic length scales are found in both the 2D and 3D analyses, where the 0.5 coherence contours found at wavelengths ca. 64-128 km in Fig. 5 correspond to wavenumbers of ~0.016-0.008 km -1 in Fig. 4, and the contours associated with longer wavelengths, ca. 1024 km, in Fig. 5 correspond to wavenumbers of ~0.001 km -1 in Fig. 5. Although the increase in coherence between the forced and total signals at long time scales noted previously in Fig. 3i can also be identified in the spatiotemporal analysis, it is notable that, in all cases, following the establishment of the intrinsic signal at frequencies with periods ca. 1 month, the mesoscale wavelength range remains dominated by the intrinsic contribution even up to the longest frequencies resolved in the analysis (here, periods of ~14 years).

Table 1

Basin mean values of RF. Areas used to calculate the means are as follows: North Atlantic (8-77°W,5-65°N), South Atlantic (70°W-24°E,35-5°S), North Pacific (104°E-107°W,5-60°N), South Pacific (132°E-75°W,57-5°S), Indian Ocean (24-103°E,35°S-28°N), Southern Ocean (circumpolar within the latitude ranges specified in the table). 1. The thick white lines show the coherence = 0.5 contours. For ease of identification of (an) isotropy, note that each axis is drawn with equally-scaled spacing in the x-and y directions; elongated shapes thus indicate anisotropy.

Representativity of a given model run

The spatiotemporal properties of the forced and intrinsic variability indicate that the two signals have largely distinct signatures in wavenumber space, hinting that it may be possible to separate the two components, and thus potentially to estimate the two signals in the absence of an ensemble. To understand whether such an estimate is likely to yield useful information, it is necessary also to understand to what extent any given simulation (here, member of the ensemble) can be considered representative of the full population. This question can be considered to comprise two parts: firstly, whether or not the total temporal variance of any given ensemble member is comparable to the mean of the total temporal variances of all ensemble members and, second, whether the partitioning of the forced and intrinsic variability of any given ensemble member is likely to be typical of the ensemble as a whole.

To evaluate these two subquestions, an analysis of the ensemble spread of the total temporal variance of any given member, and of the contributions made to the temporal variability of any individual member by the forced-intrinsic covariance and the intrinsic variance (detailed fully in the supplementary information, Section SI2.1) is performed. These analyses indicate that the combined effects of 1. The thick white lines show the coherence = 0.5 contours.

differences in the intrinsic variance and/or in the forced-intrinsic covariability for any given individual member could lead, in certain localised regions (shown in Fig. SI2), to the obtention of results that differ from the mean ensemble state due to the choice of member; however, globally the differences in variance due to the choice of member appear to be small. Based on this analysis, we proceed with the assumption that any given member (or the observational record) is likely, globally, to provide a representative sample of the ensemble. The above analyses of the temporal and spatial partitioning and coherence of the variance suggest that it may be possible to make a meaningful estimation of the atmospherically-forced and intrinsic contributions to the total variability, given only a single data set, by a simple filtering operation, based on the characteristic length scales. The motivation for this is clear in the context of the observational record, which, to create an analogy with the model-based experiment context, can be considered as the equivalent of having only a single model run; no ensemble or twin experiment is available here, yet the obtention of an estimate of the atmospherically-forced oceanic response would clearly be of relevance to, for example, our understanding of ongoing climate change. Pragmatically, given the very high computational expense of performing an ensemble experiment, such a method may also be of use in understanding the results obtained from eddy-permitting/resolving model simulations, not only in the scientific analysis of a given model run, but also for model development purposes, where a quantification of the intrinsic variability may be useful in the assessment of whether a given introduced change has made a significant difference to the forced contribution of the simulated variability of a given quantity, for example.

Methodology

The scales to be used in the separation of the intrinsic and forced components are chosen by optimisation of the power spectral densities of the forced and intrinsic components, weighted by their respective coherences with the total signal so that wavenumber/frequency ranges over which the total signal contains useful information about each component are accorded higher priority. Since the numerical model is not a perfect representation of reality, we expect that, for example, certain dynamical features may be displaced with respect to observations; because of this, we aim to choose scales that work acceptably well over a global range, rather than attempting to compile a more nuanced, targeted set of regional scales that may give better results here, but be maladapted to other data sets.

To determine the scales to be used, we seek a cutoff with the constraint that the variance retained in the weighted forced component should be maximised and the variance lost in the weighted intrinsic component minimised at length scales greater than the cut off, using the area-weighted combined PSD (and corresponding re-calculated coherence) of all regions evaluated in Fig. 5. This yields an optimal limit of 855 km, conserving 89% of all intrinsic variance, and 70% of the forced variance. Considering very short length-scales, since Figs. 4 and 5 also suggest the forced-total coherence to be larger than the intrinsictotal coherence for most regions and timescales, we also search for a lower cutoff limit for the intrinsic range. By minimising the loss of weighted intrinsic variance whilst maximising the gain of weighted forced variance, a cutoff is determined at ca. 115 km. Since in both the forced and intrinsic cases the power is concentrated at longer lengthscales, this has relatively little effect on the overall partitioning of variance achieved using only the previous cutoff limit, increasing the forced variance preserved by 2%, and decreasing the intrinsic variance preserved by 1%. However, despite the relatively slight changes in total variance that this second cutoff affords, we nevertheless choose to implement the separation as a bandpass, rather than a simple high-/ lowpass filter, with the aim of achieving a better representation across a range of scales.

The filter is implemented as a bandpass Kaiser-Bessel filter (with parameter = α 4.5). Two configurations are tested for ten members (members #1-10) over a ten-year sample period: in one, the filtering is performed at fixed length scales (cutoffs of 115 and 855 km, as above), and in the second the filtering is performed with the cutoff expressed in degrees, based on the conversion of the appropriate lengths from km to degrees in the midlatitudes, at 45°N/S (yielding equivalent cutoffs of 1.5°and ~10.5°). Evaluation of the mean absolute error (MAE) of the residual (relative to the true ensemble mean) reveals that the difference between the estimates based on fixed length scale or fixed spacing in degrees is negligible, with a very slightly smaller error for the spacing expressed in degrees (global MAE of 0.014 m cf. 0.015 m over the tenyear period; the distribution in both cases is similar to that over the full run, shown and discussed below). Since this configuration is also significantly less computationally expensive, we thus choose to implement the filter using the cutoffs expressed in degrees. The estimate of the forced contribution is thus created by the application of this bandpass filter, and the intrinsic contribution estimated as the difference between the total and the estimated forced contribution for a randomly chosen member of the ensemble (member #28). Finally, the same filtering procedure is applied to the altimetric global sea level anomaly product [START_REF] Close | Estimate of the atmosphericallyforced contribution to sea surface height variability based on altimetric observations[END_REF]. The estimates of the forced/intrinsic contributions thus created will be referred to as the recreated fields hereinafter.

Consistency of model outputs and observational data

The model configuration has a spatial resolution of 0.25°, consistent with the grid spacing used in the production of the available gridded altimetric data. However, a number of small-scale processes are parameterised in the model, and it is thus possible that these parameterised representations (and, indirectly, the model resolution) could influence the distribution of variance across the various spatial scales, with potential implications for the results obtained here. The influence might be expected to be greatest at small scales (within the lower range of the resolved mesoscale), since larger spatial scales should be fully resolved by the model. Nevertheless, this does not exclude the possibility of differences at larger scales, either due to an imperfect representation of the underlying physics in the model, or due to interactions between the unresolved smaller-, and larger-scale variability.

The possibility of a disparity in variance amongst the scales between the model and observations has been evaluated by comparing the radially-integrated wavelength-frequency PSD for the six regions considered in Fig. 5 for both the total model SSH signal and altimetric data. Whilst there is no systematic bias at the smallest scales (below ca. 130 km), the observations have notably more energy than the model over the wavelength range ~130-500 km in all regions except the Pacific (not shown). Since this range falls within that classified as intrinsic in the method proposed here, it is hypothesized that the lower limit for the intrinsic classification, set here at 115 km, is likely to be valid. The possibility that higher resolution models/ observations may potentially yield a different distribution of variance is thus acknowledged, but, after comparison of the variance distributions of the two data sets, it is believed that the methodology adopted here based on the 0.25°model setup is suitable for use on the altimetric product.

Reliability of the forced estimate

As noted in Section 3.2 and shown in Fig. SI1, the forced signal derived from the ensemble, defined here as the ensemble mean, is generally well-defined, with the exception of certain localised regions, corresponding to areas where the SSH variability is almost entirely intrinsic. Both RF and RF m (evaluated for any given member) are very strongly correlated with the percentage of time that η F is well-defined (Spearman's rank correlation, = r 0.98). This indicates that in regions where the intrinsic variability dominates strongly over the forced contribution, the forced signal is generally not well-defined; the strength of the correlation relationship further suggests that a cutoff threshold for one quantity can be skillfully estimated from the other. Our motivation for determining this equivalent cutoff is that the calculation of the percentage of time that η F is well-defined requires knowledge of the ensemble spread; this quantity thus cannot be estimated in the absence of an ensemble, and we wish to define an equivalent value for application to a quantity that is equally representative, and can be directly calculated from the recreated fields.

Motivated by the above, the relationship between the percentage of time that η F is well-defined and RF m is fitted using a cubic spline. This function is then used to evaluate a cutoff threshold in RF m (which can be calculated directly for our recreated fields) equivalent to a value of 75% for the percentage of time that η F is well-defined. An equivalent value is derived for the observational data based on the relationship between the model and observed values of RF. This cutoff (RF = 0.18) will be used to identify regions in which we suggest that the recreated forced contribution should not be considered to be accurate, and the cutoff for the recreated intrinsic contribution to be considered as inaccurate is then taken to be 0.82 (assuming a symmetric response). Further evaluation of the validity of these values may be found in the supplementary material (Section SI2.3), where the forced and intrinsic spectra of regions where η F is ill-defined (based on the above criterion that the signal is well-defined less than 75% of the time) are examined.

Assessment of the performance of the method

Comparison of original and recreated standard deviations

As a first overview of the representativity of the recreated fields, the temporal standard deviation is calculated for both the fields recreated from the model and also the altimetric record. The results are shown in Fig. 6. For both the forced and intrinsic fields, there is strong similarity between the population standard deviations (Fig. 1b/c) and the corresponding fields recreated using the single ensemble member (Fig. 6a/ b), with both the magnitudes of the two estimates and the overall spatial patterns being similar. The estimates of the forced and intrinsic components retrieved from the application of the filtering method to the altimetric data are also similar to the model-based estimates, with the spatial patterns and magnitudes again resembling the "true" ensemble values. The recreated forced estimate for the raw altimetric data has slightly higher variability globally than in the model-based estimate (Fig. 6a andc); this difference is simply due to the fact that the altimetric data contain a long-term trend, whereas, as detailed in Section 2.3, all such signals have been removed from the model data. When the long-term signal is removed from the raw observational data (using the same method as for the model) prior to the application of the filter, the resulting patterns and magnitudes of variability correspond much better with the model-based recreation (cf. Fig. 6a ande). Despite the good global agreement, a number of hydrographic features can be seen to be displaced in the model with respect to altimetry: the Kuroshio appears to have a shorter zonal extension, and the core of the Gulf Stream/North Atlantic Current lies further from the Newfoundland coast, for example. Nevertheless, overall, the variability reproduced appears to be consistent between the true model fields and those recreated by application of the filter, and both the magnitude and patterns of variability of the fields derived from the altimetric data also strongly resemble the model-based outputs.

Assessment of the average magnitude of the errors

The performance of the recreated forced and intrinsic estimators are evaluated by calculation of the MAE and the correlation between each estimator and its respective "true" value, obtained from the ensemble. The 95% significance level for the correlation is evaluated at each spatial point by calculation of the critical r value based on the effective number of degrees of freedom, which is estimated following [START_REF] Chelton | Effects of sampling errors in statistical estimation[END_REF] for each pair of time series. The results are shown in Fig. 7.

The MAE of the forced and intrinsic estimates are identical by definition 3 , and shown in Fig. 7a. With the exception of a few localised regions of elevated MAE, the error relative to the true value is generally small, with a global mean value of 0.014 m. The distribution of the error is such that the largest residuals occur in regions that have a strong intrinsic contribution across all time scales (cf. Fig. 3a), and also relatively high variance (cf. Fig. 1).

A second estimate of the error on the recreation of the forced contribution is provided in the form of the population spread on this estimator, where this is calculated for all 50 members of the ensemble. Due to the computational expense of performing this calculation for all 50 members, this quantity is estimated over a limited 10-year sample period, rather than over the full 37-year range. This spread should theoretically be zero for a perfect estimator, since the true forced contribution is defined as the ensemble mean. The population spread is shown in Fig. 7b, and shows strong resemblance to Fig. 7a, suggesting that the regions of higher MAE represent a difficulty in obtaining an accurate partitioning of the signal in a few, localised, highly intrinsic regions, rather than from any widespread, fundamental misrepresentation of the signal arising from the choice of partitioning. Furthermore, the largest errors are confined to regions in which the forced signal may be ill-defined (in the model and/or our estimate), based on our RF threshold criteria (shown by dashed white contours in Fig. 7).

Correlation between the original and recreated fields

The correlation between the two recreated model fields and their corresponding true ensemble outputs indicate generally good agreement between the two pairs (Fig. 7c/d). The correlation between the forced signal obtained from the model ensemble and that of the recreated field exceeds 0.9 over 50% of the global ocean, and over 74% of the global ocean when regions where the forced signal may be ill-defined are excluded. In the seasonally ice-covered regions and some coastal areas, the intrinsic contribution is not significantly correlated with the estimate at the 95% level, indicating that the recreated intrinsic field should be interpreted cautiously here. However, as for the MAE, these areas are identified effectively by our RF threshold criterion. The recreated forced variability shows generally very good agreement with the ensemble mean, with only very small regions of the Kuroshio and the Agulhas Return Current, and points in the Gulf Stream and BFC that are not significantly correlated. Again, these points lie within regions where the RF threshold criterion indicates that the forced signal may not be well-defined in the ensemble and/or the recreated field. As noted previously, these correspond to regions in which the variability is almost purely intrinsic at all time scales. 4.3. Two applications using the observational record: long-term trends in SSH and detection of the influence of the North Atlantic Oscillation 4.3.1. Estimation of the forced and intrinsic contributions to long-term trends in SSH Being a metric of the evolution of the depth-integrated ocean heat content, the rate of change of sea level is an important and widely used climatic indicator (e.g. Church et al., 2013 and references therein). Motivated by this, as a first practical application of the separation method we consider the trend in observed SSH from altimetry over the 1993-2017 period, as derived from both the raw (total) signal, and the recreated field of forced variability obtained using the method described above. The results, and the differences between the two 3 Let ηF and ηI represent the recreated forced and intrinsic contributions.

Since = + η η η T F I
and also

̂= + η η η T F I , it follows trivially that ̂- = -- η η η η ( ) F F I
I , and so the MAE is the same.

estimates of the trend, are shown in Fig. 8. Considering the large-scale spatial patterns, the two trend estimates (Fig. 8a/b) are, broadly, extremely similar. This is to some extent unsurprising, since the very definition of our recreated forced estimate is that it comprises large spatial scales, and is consistent with the modelbased results of [START_REF] Llovel | Contributions of Atmospheric Forcing and Chaotic Ocean Variability to Regional Sea Level Trends Over 1993-2015[END_REF]. Nevertheless, although the result might appear predictable, it is not inevitable, since the trend is here calculated on a pointwise basis. It might thus be predicted that the two estimates should resemble one another in regions where the forced variability dominates over the intrinsic; however, in regions where the intrinsic variability is strong, it is not given that there should be a similarity between the two signals. Indeed, it can be anticipated that ultimately this will depend on the time scales associated with the intrinsic signal at any such point, and on the magnitude of the low-frequency intrinsic variability relative to the forced trend.

Consistently with these assumptions, the regions in which the two estimates of the trend differ the most are essentially those in which the intrinsic variability is strongest. This is exemplified by the Southern Ocean, where the trends associated with the ACC comprise a series of small regions with random phase (Fig. 8a), associated with an intrinsic response (Fig. 8b), and that leave very little residual signature in the recreated forced SSH (Fig. 8c). The Gulf Stream/ North Atlantic Current and Kuroshio regions, also associated with strong eddy activity, are similarly modified. Notable differences can also be seen in the subtropics, away from these major current systems. Along the 30°S parallel, for example, the trend based on the recreated forced SSH appears weaker both in the west Pacific and Indian Oceans.

Calculation of the difference between the trends estimated using the total signal and the recreated forced SSH (Fig. 8c; equivalent to the trend of the recreated intrinsic field) reveals a striking coherence in the resulting spatial patterns, evoking the widescale prevalence of striations previously described by [START_REF] Maximenko | Observational evidence of alternating zonal jets in the world ocean[END_REF], and suggesting a weak, but quasi-global influence of low-frequency intrinsic variability. These alternating bands are strongest in the Kuroshio region (consistent with previous descriptions of decadal-scale mesoscale variability here, e.g. [START_REF] Qiu | Eddy-mean flow interaction in the decadally modulating Kuroshio Extension system[END_REF]), but are nevertheless evident across large regions of the Pacific and Atlantic Oceans, albeit with much weaker magnitude. The Gulf Stream shows similar coherent organisation, with anomaly magnitudes being comparable to the Kuroshio region. In contrast, the differences between the two trends in the ACC region are much less structured, yielding a jumble of patches of positive and negative trend with comparatively small spatial scale, associated with the random phase of the low-frequency intrinsic signal. Investigation of a selection of these patches suggests that many show moderate anticorrelation (of order -0.5) with neighbouring patches of opposing sign; coherence analyses of these examples indicate that the maximum covariability tends to occur in the inter-to multiannual time scale range. We thus hypothesize that, in the Southern Ocean, these patches may represent low-frequency meandering of the ACC.

The ratio of the difference between the total and recreated forced trends to the recreated forced trend is shown in Fig. 8d. The estimate confirms the general schema noted above, that the large-scale forced and total signal trends generally show good resemblance, but also highlights the strong effect that the low-frequency intrinsic variability can have on local estimates. This intrinsic influence leads not only to strong relative uncertainty on the size of the long-term trend in regions where the recreated forced signal indicates a weak trend (e.g. the Gulf Stream/ North Atlantic Current, or the Ross sector of the Pacific Southern Ocean), but also a much more moderate, but still discernible influence in regions where the forced trend is strong (e.g. mean values of approximately 14% and 11% in the subtropical Pacific and Indian Oceans).

Overall, these results demonstrate the highly coherent structure of the influence of intrinsic variability at long time scales, and suggest that, although the intrinsic contribution does not appear to significantly alter the large-scale patterns of variability, it can nevertheless modify them. Consistently with the model-based results of [START_REF] Llovel | Contributions of Atmospheric Forcing and Chaotic Ocean Variability to Regional Sea Level Trends Over 1993-2015[END_REF], however, it is found that, although the effect of the intrinsic variability on long-term trends can be large, it is most likely to be influential at a regional, rather than global scale.

Detection of the oceanic response to the North Atlantic Oscillation

As a second application of the method, the utility of the recreated forced field as a method of clarifying the impact of large-scale forcing in a region of strong intrinsic variability is considered; this will be elucidated here by examining the influence of the North Atlantic Oscillation (NAO) on the recreated forced SSH field obtained from the altimetric record. We consider the winter (JFM) mean of all signals, where the seasonal mean NAO index (DJF) is created from the monthly index, and both the total and recreated forced seasonal mean SSH fields are created as means of their respective 5-daily fields. This decision is made firstly because the influence of the NAO is strongest in the winter season, and, second, to highlight the relative importance of the intrinsic contribution at seasonal time scales.

The loading coefficients of the regression of the total SSH signal and the recreated forced field on to the NAO index are shown in Fig. 9a/c. The loading coefficients for the total signal are characterised by a patchwork of small regions with high values spanning the latitude band 30-50°N. In contrast, in the recreated forced estimate, these small regions are largely absent, with a weaker, but much more uniform background tripole pattern emerging from the coefficients. The extreme but spatially localised values noted in Fig. 9a suggest a strong role of mesoscale features in determining the local variability. This suggestion is confirmed in the evaluation of the local percentage of variance explained (Fig. 9b/d). Here, the region defined by the swath of positive loading coefficients in Fig. 9a/c is characterised by very weak correlation with the NAO when considering the total signal. However, when considering the recreated forced SSH, a band of increased correlation emerges across the southern flank of this central pole; in the total signal, this large-scale influence of the NAO (hypothesized to represent the intergyre gyre of [START_REF] Marshall | A Study of the Interaction of the North Atlantic Oscillation with Ocean Circulation[END_REF] is almost totally obscured by intrinsic mesoscale variability, even when considering seasonal average values. This is illustrated further in Fig. 9e/f by two time series drawn from this central pole, demonstrating both the much stronger variance of the total signal (blue lines), and that the recreated forced SSH (red lines) has a form that strongly resembles that of the NAO (r = 0.27 and 0.22 for the total signal/NAO correlations for panels e/f respectively, cf. r = 0.78 and 0.75 for the recreated forced/NAO correlations).

These two applications highlight two means by which intrinsic variability can disrupt efforts to interpret SSH data at climatic time scales. In the case of the long-term SSH trends, low-frequency intrinsic variability can strongly modify the total estimated trend in certain regions, leading to the obtention of a result that is not representative of the atmospherically-forced response. In contrast, in the case of the detection of the influence of the NAO, the intrinsic variability was found to dominate over the forced response at seasonal to interannual time scales, thus effectively obscuring the large-scale SSH response via a patchwork of strong, uncorrelated, mesoscale signals.

Discussion

The summary metrics and good performance of the separation method presented above suggest that it is possible to generalise the spatiotemporal distribution of the forced/intrinsic contributions to the total SSH variability to obtain useful information about the two contributions over large regions of the global ocean. Consideration of the temporal spectral coherence between the total and the forced and intrinsic contributions suggests that, although the forced contribution predominates at short (up to approximately 2 weeks) and, in certain regions, long (>4 yr) time scales, over the intervening time scales, the signatures of the intrinsic and forced contributions are not distinct. There is an apparent contrast between this result and the findings of Forget and Ponte (2015), who found, based on analyses of a state estimate with 1°resolution and altimetric data, that the variance associated with spatial scales of less than 3°(termed therein as noise, or intrinsic variability) exceeded the variance associated with larger scales up to three years. This discrepancy may potentially be attributable to a number of factors relating to the partitioning of variance amongst the scales; notably, our results indicate that the intrinsic contribution dominates the total signal on spatial scales much larger than the 3°c utoff used in [START_REF] Forget | The partition of regional sea level variability[END_REF] (up to ~855 km), suggesting, then, that the form of the temporal power spectrum may depend more strongly on the dominant spatial scale range than the intrinsic or forced The results presented here may have implications both for modeland observation-based studies. From a modelling perspective, analysis of the ensemble suggests that inter-member differences in both the intrinsic variance and, particularly, the forced-intrinsic covariance can lead to differences in the total variance simulated in a given run. With the exception of a few localised regions, over the 37-year period considered here, the differences between the total temporal variance among the ensemble members appears to be generally moderate. Nevertheless, locally (e.g. in the North Pacific), the differences can be substantial (discussed further in the supplementary materials, Section SI2.2). Further, considering the analyses of temporal coherence presented here (Fig. 3), it is possible that greater disparity would be observed amongst short simulations (< approx. 5 years). Application of the method outlined above may potentially prove useful to gain insight into the extent of this problem in the absence of an ensemble: by estimating the forced and intrinsic contributions and evaluating the relative magnitudes and associated time scales of the two, it should be possible to obtain a better understanding of the relative importance of the non-deterministic contribution in the total simulated variability, and thus the representativity of a given run.

In this work, the method developed for the estimation of the intrinsic and atmospherically-forced contributions to SSH variability has been applied to the altimetric record, a situation in which it is clearly not possible to obtain multiple versions of the oceanic response to a prescribed forcing. However, the statistical properties of the forced and intrinsic contributions suggest further implications for the interpretation of (particularly in situ) observations and, potentially, for observing system design. Both the pointwise analysis of RF (Fig. 2) and the spectral coherence analysis of the forced/total temporal variability on a pointwise basis (Fig. 3) highlight the important role of intrinsic variability over relatively large regions when all spatial scales are integrated. This should be borne in mind in the interpretation and attribution of the variability of time series originating from a single location (e.g. moorings, tide gauges) in locations where the intrinsic variability has been identified to make a significant contribution, even at time scales ranging up to the multiannual. The analyses presented in this work indicate that, for the case of SSH, it is not correct to assume that the impact of intrinsic variability can be removed by simple temporal filtering of a given time series in all cases.

Despite the potential utility of the method demonstrated by these applications, certain limitations are nevertheless associated with this approach. The model resolution may potentially impact the form of the spectra obtained (discussed previously in Section 4.1.3), and further differences in the ocean dynamics may potentially be found were the atmosphere to be coupled to the ocean. We thus acknowledge that the model outputs are perhaps best considered as a tool by which to learn about the processes at play, rather than a true representation of reality. This view has guided the approach employed here, and, despite the potential sources of error associated with the use of a numerical model, the good agreement between the recreated observed fields and the model output is encouraging, and suggests that the approach does have merit.

As a further limitation, we emphasize that, since the method is ultimately based on a simple separation of scales, it should be borne in mind that it will not be well suited for all purposes. In particular, for applications where the SSH is used to make inferences about ocean dynamics, it is anticipated that the recreated fields are unlikely to yield accurate results since the spatial scales of many of the features of interest (e.g. western boundary currents) are comparable to or smaller than the 1.5°lower cutoff bound, and thus simply approximated as forced by our method. We thus caution against the use of the fields produced by this method for such purposes.

Considering particularly long-term behaviour (e.g. estimation of decadal-scale variability, or the linear trend at a given point), the presence of a non-negligible intrinsic contribution to the total variability may affect the interpretation of the signal under consideration in two ways. Firstly, where the intrinsic signal is of comparable size to the forced contribution over the temporal and spatial scales considered, the magnitude of the total variability and/or estimated linear trend derived from the total signal will contain a non-negligible, random contribution from the intrinsic component. Examples of this effect are evident in the analysis of the long-term trend presented in Fig. 8: in regions that are strongly characterised by mesoscale variability (notably the major current systems, but also, for example, the southwest Pacific), although the associated power is strongest at subannual-up to 1-2 year time scales, the variance associated with these characteristic intrinsic wavelengths can nevertheless remain stronger than that of longer wavelengths at multiannual time scales (i.e. than the forced trend). These results are consistent with those obtained by [START_REF] Llovel | Contributions of Atmospheric Forcing and Chaotic Ocean Variability to Regional Sea Level Trends Over 1993-2015[END_REF] in an analysis of long-term trends based on the ensemble used here in the inference of the forced and intrinsic characteristic wavelengths. As shown in Fig. 5 (and exploited in the forced/intrinsic estimation method used here) over wavelengths of ~115-855 km, the total signal experiences a dominant influence from the intrinsic contribution at all resolvable time scales. This strong, long-term intrinsic influence over these wavelengths may present a problem in the context of attribution studies, and advocates for a careful evaluation of not the temporal, but also the dominant spatial scales that characterise the signal under consideration before interpretation of such results.

The second means by which the intrinsic signal may influence the interpretation of the variability of a given signal occurs when the intrinsic contribution has an extremely large variance. Amongst the regions with the highest intrinsic variance are those associated with bimodal states; as noted by [START_REF] Hughes | Identification of jets and mixing barriers from sea level and vorticity measurements using simple statistics[END_REF], such regions occur at the interface between zones of positive and negative skew, and may thus occur when orthogonally traversing a front or jet, being associated with the approximately equal sampling of two different hydrographic regimes at a fixed point as the jet meanders. Especially large variances due to this behaviour at a fixed point then occur when there is a strong contrast between the two hydrographic regimes that are sampled; an excellent example of this is the Agulhas region. Although the variance of the forced contribution here is, in absolute terms, not negligibly small (being equivalent to a standard deviation of approx 8 cm), the intrinsic contribution to the total variance is much larger ( ~35×) than that of the forced component (equivalent to a standard deviation of 45 cm), thus dwarfing the forced contribution. In such scenarios, detection of any change in the forced signal is likely to be, in practical terms, impossible; for example, considering the variances of the two signals above over the period of the model data availability (i.e. 37 years, 1979-2015), it would not be possible to distinguish a trend in the forced signal from noise for values less than ~2 mm/yr. However, over the shorter length of the current altimetric record (i.e. 26 years, 1993-2018), the signal would not be detectable for values less than ~5 mm/yr. For the model-length period, this equates to a total change of at least 8 cm over the 37-year period, equivalent to the standard deviation, whilst the shorter, 25-year period would implicate a total change of at least 12.5 cm. In addition to its potential to contribute significantly to the variance at any given frequency, the intrinsic signal may thus also affect the interpretation of the total variability by a second mechanism, in which the stationary, but highly variable intrinsic signal is able to render changes in the forced signal undetectable.

Summary and conclusions

The main findings of this work are that the intrinsic and atmospherically-forced contributions to the total SSH variance have, globally, distinct signatures in their respective spatial power spectra, permitting effective estimation of the two quantities by means of a simple filtering operation employing a well-adapted choice of cutoff length scales (approximately 1.5°and 10.5°). Application of the method to both a single model simulation and to the observed altimetric record yields estimates of the forced and intrinsic variability that agree well with the corresponding fields calculated using the full model ensemble. We thus suggest that spatial filtering can be used to meaningfully approximate the forced and intrinsic contributions to SSH variability given either a single model experiment, or observational data.

In contrast with the distinct spatial spectra, our analyses suggest that the temporal spectra of the two components generally overlap, and that attempts to separate the two quantities based solely on frequency would be unlikely to be effective on a global scale. Given two isolated time series, we thus find that it is not possible to separate the forced and intrinsic contributions using only a temporal criterion in all regions, even on multiannual time scales: additional spatial context may be needed. This is confirmed in our spatiotemporal analysis, which highlights that the wavelength range classified here as intrinsic remains intrinsic even at the longest time scales resolved (here, up to ~14 years). Further, this result is also manifest in the estimation of the forced variability and long-term trends from altimetry, where we note widespread structure in the difference field, consistent with the existence of previously-reported intrinsic striations and suggestive of a low-frequency intrinsic contribution to the long-term signal.

The results presented here have implications both for the interpretation of observations and for numerical modelling experiments. Whilst spatial filtering may provide an effective way to estimate the atmospherically-forced part of the signal, we note that this may present a challenge for any localised in situ observations taken in highly intrinsic regions in the absence of a wider spatial context, and suggest that the variability described by such data should be interpreted cautiously at all time scales. In a numerical modelling context, our results largely confirm those of previous studies based on twin experiments, and both provide a quantification of the magnitude likely differences in variance among runs driven by the same forcing, and suggest a potential way to estimate the atmospherically-forced contribution in the absence of an ensemble.

The findings of this work confirm previous reports that the intrinsic component of SSH variability can have an impact over a wide-range of time scales, but show that a meaningful estimate of the magnitude of this contribution can nevertheless be made. Given the computational expense of ensemble simulation and its inapplicability to the observational record, this suggests cause for optimism in ongoing efforts to disentangle the deterministic and non-deterministic parts of the SSH signal.

Fig. 1 .

 1 Fig. 1. Standard deviation of (a) total signal, (b) atmospherically-forced contribution and (c) intrinsic contribution. White contours in (b) show the regions in which the forced signal is well-defined less than 75% of the time.

Fig. 2 .

 2 Fig. 2. Ratio of forced to total variance for (a) all time scales, (b) subannual time scales and (c) interannual and greater time scales.

Fig. 3 .

 3 Fig. 3. Spectral coherence in frequency space between the forced and total SSH on a pointwise basis, for periods corresponding to: top row: (a) 10 days, (b) 20 days, (c) 30 days; middle row: (d) 1-3 months, (e) 3-6 months, (f) 6-12 months; bottom row: (g) 1-2 years, (h) 2-4 years, (i) 4 years and longer. Thin black contours show the 95% significance level (coh = 0.64, see supplementary information). The values shown are band averages, calculated using the coherence.

Fig. 4 .

 4 Fig. 4. Time-average spatial coherence between the forced contribution and total signal for: (a) North Atlantic, (b) South Atlantic, (c) North Pacific, (d) South Pacific, (e) Indian Ocean and (f) Southern Ocean, where the regions are defined as in Table1. The thick white lines show the coherence = 0.5 contours. For ease of identification of (an) isotropy, note that each axis is drawn with equally-scaled spacing in the x-and y directions; elongated shapes thus indicate anisotropy.

Fig. 5 .

 5 Fig. 5. Spectral coherence between the forced contribution and total signal in radially-averaged wavenumber-frequency space for: (a) North Atlantic, (b) South Atlantic, (c) North Pacific, (d) South Pacific, (e) Indian Ocean and (f) Southern Ocean, where the regions are defined as in Table1. The thick white lines show the coherence = 0.5 contours.

Fig. 6 .

 6 Fig. 6. Standard deviation of the recreated SSH fields for: (a) forced and (b) intrinsic components using the total SSH signal of a single member of the ensemble; (c) forced and (d) intrinsic components from altimetric data without removal of trend; (e/f) as (c/d) but with trend removed from the total observed signal using the same method as for the model data. Hatched lines indicate regions not covered by the observational data. White dashed contours indicate regions where our criteria based on RF threshold values indicates that the forced (a,c,e) and intrinsic (b,d,f) signals may be unreliable.

Fig. 7 .

 7 Fig. 7. (a) MAE of recreated field relative to true ensemble values (same for intrinsic and forced components). (b) Population standard deviation of the recreated forced contributions calculated for all 50 members over a 10-year period. Correlation between recreated and true time series for (c) forced and (d) intrinsic components. Hatched black contours indicate areas where the correlation is not significant at the 95% level after taking into account the effective number of degrees of freedom. White dashed contours indicate regions where our criteria based on RF threshold values indicates that the forced (a,b,c) and intrinsic (d) signals may be unreliable.

Fig. 8 .

 8 Fig. 8. (a) Trend in SSH over 1993-2017 from the total observed signal. (b) Trend in SSH from the recreated forced observed signal. (c) Difference between the full and recreated forced observed trends. (d) Difference between the full and recreated forced observed trends expressed as a fraction of the recreated forced trend.

Fig. 9 .

 9 Fig. 9. Loading coefficients (a) and % of variance explained (b) by regression of total signal JFM mean SSH from altimetry on to the DJF NAO index. (c/d) as (a/b) but for the recreated forced field. Black stars in (c) show the location of the two sets of time series shown in (e) and (f), where the blue lines show the total SSH at the given point, the red lines show the recreated forced SSH and the black lines show the NAO index. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

this result may be anticipated from the definition of the forced contribution as the ensemble mean

Acknowledgements

The Ssalto/Duacs altimeter products were produced and distributed by the Copernicus Marine and Environment Monitoring Service (CMEMS) (http://www.marine.copernicus.eu). Funding: This work is a contribution to the AtlantOS project, and has received funding from the European Union Horizon 2020 research and innovation program under grant agreement No 633211. This is also a contribution to the PIRATE project funded by CNES through the Ocean Surface Topography Science Team (OST-ST), and to the GLO-HR project funded by the Copernicus Marine Environment Monitoring Service (CMEMS); CMEMS is implemented by Mercator Ocean International in the framework of a delegation agreement with the European Union. This work was also supported by the French national programme LEFE/INSU. The ensemble simulation used in the study was performed as part of the OCCIPUT project, funded by the ANR through contract ANR-13-BS06-0007-01. We acknowledge that the results of this research have been achieved using the PRACE Research Infrastructure resource CURIE based in France at TGCC; some of the computations presented in this study were performed at TGCC under allocations granted by GENCI. The model data set used in this study is available upon request (contact: thierry.penduff@cnrs.fr). The estimates of the atmospherically-forced SSH variability based on altimetry have been archived in the Zenodo data repository and are freely available (Close et al., 2020). We thank the three reviewers for their constructive comments, which improved the quality of the manuscript.

Declaration of Competing Interest

The authors declare no competing interests.

Appendix A. Supplementary material

Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.pocean.2020.102314.