Forced and intrinsic variability in the response to increased wind stress of an idealized Southern Ocean

We use ensemble runs of a three-layer, quasigeostrophic idealized Southern Ocean model to explore the roles of forced and intrinsic variability in response to a linear increase of wind stress imposed over a 30-year period. We find no increase of eastward circumpolar volume transport in response to the increased wind stress. A large part of the resulting time series can be explained by a response in which the eddy kinetic energy is linearly proportional to the wind stress with a possible time lag, but no statistically significant lag is found. However, this simple relationship is not the whole story: several intrinsic timescales also influence the response. We find an e-folding timescale for growth of small perturbations of 1-2 weeks. The energy budget for intrinsic variability at periods shorter than a year is dominated by exchange between kinetic and potential energy. At longer timescales, we find an intrinsic mode with period in the region of 15 years, which is dominated by changes in potential energy and frictional dissipation in a manner consistent with that seen by Hogg and Blundell [2006]. A similar mode influences the response to changing wind stress. This influence, robust to perturbations, is different from the supposed linear relationship between wind stress and eddy kinetic energy, and persists for 5-10 years in this model, suggestive of a forced oscillatory mode with period of around 15 years. If present in the real ocean, such a mode would imply a degree of predictability of Southern Ocean dynamics on multi-year timescales.


Introduction
Observational reanalysis shows that eastward wind stress over the Southern Ocean has been increasing by approximately 0.01 Pa (or 10%) per decade over the period 1979-2010 [Swart and Fyfe, 2012; although certain caveats about sampling should be remembered: Kent et al., 2013]. Despite the strengthening winds, observations show that there is little change in baroclinic eastward volume transport through Drake Passage when sampled in the upper 3000 m or shallower [Cunningham et al., 2003;B€ oning et al., 2008;Firing et al., 2011] within this multidecadal period. There is more uncertainty about the steadiness of the deep reference-level velocity due to poorly sampled near-bottom currents [Cunningham et al., 2003;Firing et al., 2011]. Observations of the barotropic component of transport rely on bottom pressure recorders which cannot be used for more than a few years, therefore preventing detection of longer-term trends in transport. However, a review of historical measurements and modeling of the Drake Passage transport by Meredith et al. [2011] suggests that the baroclinic component accounts for about 70% of the total, and that the variability due to the barotropic component is probably only important on time scales of a few years or less. A recent model study [Hughes et al., 2014] supports the suggestion of Olbers and Lettmann [2007] that barotropic effects cease to dominate at periods longer than about 5 years.
Consideration of the change in the momentum budget due to strengthening winds has led to the promotion of the ''eddy saturation'' mechanism [Straub, 1993]. This theory explains why the mean meridional slope of isopycnals does not appear to increase even though the Ekman circulation acts to steepen them, because mesoscale eddies are responsible for maintaining a critical slope and therefore a steady baroclinic transport. There are many recent studies which have examined eddy saturation in a variety of contexts [e.g., Hallberg and Gnanadesikan, 2001;Tansley and Marshall, 2001;Hallberg and Gnanadesikan, 2006;Hogg and Blundell, 2006;Meredith and Hogg, 2006;Munday et al., 2013], including either sudden and gradual perturbations to forcing, locally or remotely to the Southern Ocean. Hogg et al. [2008] investigate the role of Ekman layer dynamics in the fast response to local, sudden changes in wind forcing, in contrast to the response of the eddy saturation mechanism. Straub [1993] and Nadeau and Straub [2012] highlight eddy saturation as one component of the transport through Drake Passage (the channel component, where Sverdrup dynamics do not hold). They also discuss the role of the subtropical gyre transport, north of the main Antarctic Circumpolar Current (ACC).
In this study, we focus only on the channel-like dynamics, as has been done previously [e.g., Straub, 1993;Meredith and Hogg, 2006;Hogg and Blundell, 2006;Hogg et al., 2008], but we also extensively consider nonlinear intrinsic variability. Since the eddy saturation mechanism assumes a central role for mesoscale eddies and since the mesoscale dynamics are significantly more nonlinear than dynamics at large scales (Rossby number 0.1 versus 0.01, respectively), and indeed are comparably nonlinear to the dynamics of atmospheric weather, ensemble forecasting is needed to distinguish between forced and intrinsic variability.
With these runs, we can test the robustness of the relationship between wind stress and eddy kinetic energy (EKE) suggested by the observational study of Meredith and Hogg [2006] and learn about the time scales and predictability of intrinsic variability. We note that the forced response of the system may exhibit a linear or nonlinear relationship to the forcing and that even a linear relationship may depend on nonlinear dynamics, for example, for any such relationship including EKE.
These experiments cover a total of 660 years at 5 km resolution using a three layer, quasi-geostrophic channel model of the Southern Ocean (Q-GCM) [Hogg et al., 2003], and include multidecadal simulations. Such experiments remain too expensive for primitive equation modeling at this spatial resolution, so quasigeostropic modeling is required to fill this niche. Many aspects of the observed circulation are reproduced (jets, eddies, topographic steering, realistic ACC volume transport) and the nonlinear dynamics give us the best estimate of the robustness of eddy saturation over interannual to decadal time scales. In contrast to previous Southern Ocean studies using Q-GCM [Hogg and Blundell, 2006;Meredith and Hogg, 2006;Hogg et al., 2008;Meredith et al., 2012] our experiments have double the lateral spatial resolution, enabling more realistic values of the baroclinic Rossby radii (22 and 12 km) in this study.
We seek to separate the forced response of the Southern Ocean to trends in wind stress from the nonlinear intrinsic variability. We will test whether the nonlinear, eddying Southern Ocean exhibits sensitive dependence on external forcing and on initial conditions, or whether it may be approximated well by a linear response without such sensitivities. For an analogous example, we do not know whether the dynamical response and adjustment time scale of 2-3 years for EKE observed by Meredith and Hogg [2006] would be similar if the predicted time series of the Southern Annular Mode (SAM) was otherwise identical, but scaled by 95% or 105% (i.e., if the wind stress forcing was marginally different in strength). Also, we do not know whether an identical SAM prediction but with slightly different ocean conditions would lead to the same response in EKE and eastward volume transport.
Our emphasis is on the response to multidecadal trends in wind stress forcing, rather than to interannual changes in the SAM, as in Meredith and Hogg [2006]. However, from a dynamical as well as a statistical perspective, it is interesting to test whether there is an obvious common time scale of adjustment, or of EKE lagging perturbations to wind forcing. Since baroclinic Rossby waves are too slow to propagate westward against the strong eastward Antarctic Circumpolar Current, there is no clear horizontal adjustment time scale other than an advective one, for processes confined to the Southern Ocean. Disentangling the relevant time scales of adjustment from the nonlinear intrinsic variability is highly relevant to improving our understanding of the Southern Ocean climate response and is the main aim of this study. For the global adjustment problem, in which the global thermocline depth north of the Southern Ocean is required to equilibrate to changes in Southern Ocean forcing, adjustment time scales in the region of 100-1000 years have been proposed [Allison et al., 2011;Jones et al., 2011;Samelson, 2011]. An important caveat on our study, using quasi-geostrophic equations in a limited latitude range domain, is that such global adjustment cannot be realistically simulated, so an equilibrium may be reached artificially rapidly. Our results can only be expected to translate to the real ocean to the extent that the processes depend on density gradients within the Southern Ocean, rather than on the absolute thermocline depth to the north or south.
Under the assumption of a perfect model, we therefore explore the sensitivity of the response of the total eastward volume transport (transport hereafter) and domain-mean, depth-integrated EKE density to: (i) changes in the strength of 30 year linear wind stress trends imposed in addition to background seasonal wind stress forcing (using a control plus 10 member boundary condition perturbation ensemble (BCPE)) and (ii) for selected members of the BCPE we examine the robustness of the response with respect to the intrinsic chaotic variability (using an initial condition perturbation ensemble (ICPE) for each selection).
In section 2, we describe the model and experiments, then we examine the results from the BCPE integrations in section 3, and the associated energy balance in section 4. In section 5, we look at the ICPE integrations to distinguish between forced intrinsic variations, with a discussion and summary in section 6.

Model and Experiments
2.1. The Quasi-Geostrophic Coupled Model (Q-GCM) Run in Ocean-Only Mode Q-GCM [Hogg et al., 2003] was developed to examine the role of the ocean mesoscale on climate processes and integrated climate metrics. After more than 10 years since its creation, primitive equation coupled ocean-atmosphere climate models that resolve ocean eddies are only now beginning to be able to simulate interannual to multidecadal periods using perturbed-parameter ensembles. As already highlighted in Hogg and Blundell [2006], there is strong evidence that the ACC cannot be represented fully by quasi-geostrophic models, but that such models do enable the simulation of nonlinear effects and eddies in an efficient manner, thereby permitting numerical experiments across a wide parameter range.
It is not our aim in this paper to simulate the exact details of particular branches of the ACC, but rather to examine the dynamics of a flow with similar characteristic spatial and temporal scales. This includes topographically steered, multiple zonal jets and a well-resolved eddy field which is able to mix tracers and momentum between large and small scales. Therefore, we use high-resolution topography but do not refer subsequently to particular geographic features. Following Hogg and Blundell [2006], the topography is obtained by averaging onto the model grid the Earth Topography 2-Minute (ETOPO2) data set [U.S. Department of Commerce, National Oceanic and Atmospheric Administration, National Geophysical Data Center, 2001], which for the latitude range we consider is essentially the data of Smith and Sandwell [1997]. At the central latitude of the domain, 55 S, the ETOPO2 grid equates to 2.1 km and contains a good representation of fracture zones, topographic slopes, and rough topography, even when averaged onto the 5 km model grid. Unlike Hogg and Blundell [2006] we use the full range of circumpolar longitudes in this study. Topographic height variation greater than ðd3lowest layer thicknessÞ 860 m, where d5 bY f0 (where b is the meridional gradient of the Coriolis parameter, Y is the meridional domain size and f 0 is the domain-mean Coriolis parameter; see Table 1), is truncated, replacing extreme values by shallow or abyssal plateaus. This avoids violation of the quasi-geostrophic assumptions. The truncated topography, which represents variations about a reference depth of 4000 m, is shown in Figure 1a.
The model governing equations are the same as those in Hogg and Blundell [2006], but certain parameters are modified to reflect the higher spatial resolution of 5 km and the more realistic stratification used to obtain more realistic baroclinic deformation radii. For a full set of the parameters used in these experiments see Table 1. For more details of the model code (version 1.4.1) and the detailed formulation and notation conventions, please refer to the appropriate user guide (version 1.4.0; see www.q-gcm.org/downloads.html; hereafter referred to as the ''Q-GCM user guide''). The wind stress parameterization depends on the 10 m wind and not on ocean surface currents. The choice of layer thicknesses and reduced gravities was not straightforward and it required careful consideration of the interaction between vertical normal modes and surface and bottom forcing, and between themselves, as outlined by Flierl [1978]. As predicted by the extension of the analysis of Flierl [1978] to three layers, the solution for the final parameter set can only satisfy a subset of all the desired constraints. Using a background profile of N 2 ðzÞ, the square of the buoyancy frequency, derived by averaging WOA05 climatology Antonov et al., 2006] over the model domain, a particular set of reduced gravity values and unperturbed layer thicknesses was found using shooting methods and the same type of vertical eigenmode constraints as in Flierl [1978]; this set was found to have good properties such as giving a sensible transport and large-scale potential vorticity distribution. Time step and dissipative parameters were then tuned, as is the convention, to get model behavior that was judged acceptable. A snapshot of potential vorticity (PV) highlights the presence of well-resolved mesoscale eddies and multiple, quasi-zonal PV fronts, as observed in the Southern Ocean and associated with jets (Figures 1b-1d). It is clear that all three layers are influenced by topographic interactions, but the topographic forcing of PV is clearest in the deepest layer ( Figure 1d). Even in the shallowest layer ( Figure  1b), the PV has regions where topographic steering causes a significantly nonzonal flow.

Spin-Up and Ensemble Experiments
The model includes a fixed depth (100 m) mixed layer with varying temperature, which requires some assumptions regarding heat fluxes to be made. In order to focus on wind stress forcing, we chose to minimize buoyancy effects by initializing the model with a mixed layer temperature profile in radiative equilibrium as described in appendices A and D of the Q-GCM user guide, which produced a reasonable N-S temperature range of 19.1 K. In the temperature calculation, the northern boundary of the mixed layer is formulated as an open boundary, with temperature beyond the boundary fixed at the initial state of radiative equilibrium. The northern boundary mixed layer temperature is the only boundary parameter which is specified and there is no relaxation of any parameters in the interior. As a result of the presence of such a mixed layer, thermodynamic balance leads to a requirement for local mass fluxes between layers 1 and 2 which are related to the mixed layer temperature and the Ekman pumping. Although these local mass fluxes integrate to zero over the layer, they lead to a thermodynamic sink of energy at the interface, given by qg where g 1 is the upward displacement of the interface between layers 1 and 2 from its resting height, and e 1 is the upwelling entrainment flux from layer 2 to layer 1 (from E.10 in the Q-GCM user guide). We include this term in our energy calculations for completeness, though it represents only a small correction to the power input from wind stress. The surface net heat fluxes are set to zero everywhere at all times.
The lateral boundary conditions on pressure, p, at the northern and southern boundaries are written as more general partial slip conditions, expressing a linear relationship between tangential stress and tangential velocity, following Haidvogel et al. [1992]. In Q-GCM notation they are @ 2 p @n 2 52 a D @p @n , and @ 4 p @n 4 52 a D @ 3 p @n 3 , where n is a coordinate normal to the boundary, and D is the grid spacing. We use a value of a 5 5, equivalent to a 5 100 in Haidvogel et al. [1992], which nondimensionalized the coefficient a in a different way. This choice means that the boundary condition is very close to no slip. See Appendix C of the Q-GCM user guide for details of the finite-difference formulation, and the different choices for nondimensionalization.
After the spin-up the zonal wind stress forcing is of the general form s x ðy; tÞ5s seas ðy; tÞ1ð11FcðtÞÞs 0 ðyÞ; (1) described below. It is constructed using the first five temporal Fourier coefficients (mean, two annual (sine and cosine), and two semiannual) provided by the scatterometer analysis of Risien and Chelton [2008] from September 1999 to August 2007, with the following modifications: the Fourier coefficients of the observed wind stress, s x ðx; y; tÞ, which are functions of (x, y), are zonally averaged; the averaged Fourier coefficients each then have their value within 275 km from the boundaries set equal to the value at 275 km from the respective boundary; finally, they are boxcar-smoothed over a meridional length scale of 150 km (30 points). The choice of zero meridional gradient in s x near the boundaries is to ensure that Sverdrup balance does not lead to a need for balancing terms at the boundary. The southern limit of valid data from Risien and Chelton [2008] is 65:625 S and the southern boundary is chosen to be 275 km south of this. [2008] contains an error that unfortunately was not spotted until after the completion of our main experiments. The semiannual harmonics that were downloaded from their website were used to reconstruct a wind stress time series as suggested by the equations in their paper. However, their error meant that the time dependence included variations at 1 and 3 cycles per year instead of 1 and 2. The annual cycle is thus correct, but the semiannual cycle has been shifted to three cycles per year. It was not possible to repeat all of the experiments with this correction, but a four member perturbed initial condition experiment comparing corrected to uncorrected integrations over years 100-150 suggests that this has no significant effect on our results.

Risien and Chelton
The time-mean of the modified wind stress for the control (F 5 0) is s 0 ðyÞ, shown in Figure 2a  To further study the robustness of the response to external forcing with respect to nonlinear intrinsic variability, three of the BCPE members, F 5 0, 0.25, and 0.5, are examined in more detail, each with an ICPE. The initial condition perturbation is generated through changes to the surface dynamics over one time step, by slightly perturbing the eastward wind stress during the reconstruction from its Fourier components. ð2pf a t n Þ, where f a is one cycle per year and t n is time in years. At the single time step of perturbation, these arguments are replaced with ð2pf a t n Þ3ð11aÞ, where a 2 f25; 24; 23; 22; 21; 1; 2; 3; 4; 5g310 26 in order to generate each of the three, 11 member ICPEs. The ICPEs are initiated at year 162, for reasons discussed in the next section, and run for 6 years.
The control run, F 5 0, is started from rest and has an initial phase over the first 10-20 years where the potential and kinetic energy grows rapidly and equilibrates. By year 30, we consider the system to be reasonably spun-up ( Figure 3a). After spin-up Figure 3 shows that the potential and kinetic energy time series are dominated by different characteristic periods. There is a change in the variability in Figure 3a after year 100, when the steady wind stress forcing changes to include seasonal variability. After this point, the potential energy time series shows clear seasonal variations, and kinetic energy increases over about 10-20 years, by about 5310 3 J m 22 .

Hogg and Blundell
[2006] find similar long-period variability, in their case of order 15-20 years. For a particular parameter regime they find an almost pure form of oscillation, which allows them to identify the associated dynamics in some detail. It consists of a cycle in which a more-zonal and shallower flow accelerates until its expression at depth leads to a topographic instability, which prompts a burst of eddy formation, which in turn enhances form stress and dissipation. Thus, the flow becomes less zonal and decays, until its deep expression becomes small enough for the cycle to start again. We will show in later sections that the dynamics of the potential energy oscillations we see in Figure 3 are similar to those seen by Hogg and Blundell [2006], and at a roughly similar period.
In primitive equation simulations, it was found [Hughes et al., 2014] that transport fluctuations were dominated by a barotropic response to winds on time scales shorter than 5-10 years, with baroclinic changes becoming dominant at longer time scales, so 15-20 years seems to be a time scale associated with baroclinic processes.

Domain-Integrated Response to Strengthening Wind Stress
Eddy saturation is usually understood to mean a situation in which a change in wind stress is balanced by a change in eddy activity, with little change in the baroclinic transport of the channel flow. As a general measure of the intensity of eddy activity, it is usual to use EKE.
Defining eddies as the difference from a temporal mean result in the artificial inclusion of long-period changes in the mean flow in the definition of eddies. An alternative is to define eddies as variations in a particular frequency band, but this involves introducing an artificial time scale into the analysis. To avoid this, we use a definition of eddies which is based on length scale. We used a 2-D, order-10 Butterworth filter to define EKE at each instant by blocking horizontal velocity length scales larger than 500 km. The velocity, u is partitioned into high-pass and low-pass components, u5u hp 1u lp . This means that there is a separation between the scale of gyres and jets, and the scale of mesoscale vortices and Rossby waves. An artificial length scale is thus introduced, but this is of less concern as our focus here is on time scales. The convolution between the filter and the velocity is periodic in x (to match the cyclic channel boundary condition) and mirrored at the meridional boundaries. A high-order filter is used to generate a sharp cutoff and the EKE is not significantly sensitive to the choice of cutoff length scale. This measure of EKE, since it does not have any temporal cutoff, also includes stationary features such as standing topographic eddies. Figure 4 shows the partition of the total KE into EKE and MKE, for a snapshot and for time series of the domainmean values. See the figure caption for definitions of EKE and MKE in terms of filtered velocity. Both the local and the domain-mean, depth-integrated EKE is very similar to KE. The MKE is generally an order of Figure 4. Illustration of the partitioning of total kinetic energy for the control, via the spatial filtering described in the text, of the velocity into ''eddy'' and ''mean'' components: u5u hp 1u lp . Snapshots at year 180 are shown for (a) the square root of depth-integrated total kinetic energy density, KE local 5 q 2 R 3 k51 ðu k Þ 2 , where k refers to the layer index; (b) the square root of depth-integrated eddy kinetic energy density, EKE local 5 q 2 R 3 k51 ðu hp k Þ 2 ; and (c) the square root of depth-integrated mean kinetic energy density, MKE local 5 q 2 R 3 k51 ðu lp k Þ 2 . Units for (a and c) are ðJm 22 Þ 0:5 . (d) The time series of the domain-mean, depth-integrated energy density terms, KE (red; as in Figure 3c), EKE (blue), and MKE (green), shows that over years 150-180, MKE is typically 14% of EKE and with 23% of its standard deviation. EKE and KE have a correlation of 0.96, with EKE having a mean which is 79% of that of KE, and a standard deviation which is 73% of that of KE. Note that KE is not simply the sum of EKE and MKE. magnitude weaker than the KE or EKE in terms of its mean and variability, although jets and topographic eddies are visible in Figure 4c.
For the BCPE, the EKE and transport for years 150-180 are shown in Figure 5. Generally, the EKE increases with eastward wind stress and transport remains fairly steady: this is the basic signature of eddy saturation. However, there are more subtle aspects contained in the response: these details are what we want to highlight. In Figure 5a, the EKE for the control run (F 5 0) contains variability on many time scales but remains 3.6-4.7 3 10 4 J m 22 . For F 5 0.25 the EKE appears to increase over the 30 years that the wind stress trend is applied, reaching up to 4.9 3 10 4 J m 22 , and also contains variability on many time scales. A similar response is present for F 5 0.5, but the EKE now reaches a larger value 5.6 3 10 4 J m 22 and has several large local maxima, the first of which occurs at around year 163.
For transport in Figure 5b, all members of the BCPE exhibit strong variability on interannual time scales and shorter, ranging from 60 to 160 Sv.
The transport response has a strong seasonal component. Figures 5c and 5d show the same fields as in Figures 5a and 5b but with a 1 year running boxcar mean applied. This serves to highlight the presence of small interannual variations in transport, with hints of an anticorrelation with EKE and of a small decrease in transport which may or may not be a response to the increase in wind stress.

The Linear, Forced Response of EKE
If we assume a linear relationship between EKE and wind stress, s, with a time lag: EKEðF; tÞ5EKE c 1Aðs lag 2s c Þ; ( 2) where, using the y-mean of the nonseasonal component of (1), s lag 5s 0 ðyÞ y Fcðt2t lag Þ effectively contains the area-mean zonal wind stress trend. EKE c , A, and s c are constants and we can apply linear regression to the results from the BCPE to find the best fit as a function of the time lag, t lag , for the full range of F.  Figure 5a, shows the pattern of EKE with which we wish to perform a lagged regression against wind stress. We find that the maximum correlation, approximately 0.72, is when t lag 50 ( Figure  6b) and that the correlation monotonically decreases with increasing t lag . It is surprising that the 2-3 year time scale of perturbation to the SAM followed by response of surface EKE found by Meredith and Hogg [2006] does not appear, but the small changes in correlation over the first few years are insufficient to distinguish between a lag of zero and of 2 or 3 years. In addition, the fit may be influenced by the rise in EKE over the first few years, in all 11 members, which may result from the starting point at year 150 being a point of relatively low KE in the control run (see Figure 3). The clear seasonal-interannual variability in EKE does not hide the basic relationship that over interannual to decadal time scales, EKE is proportional to wind stress with negligible lag. The best fit model is shown in Figure 6c  Thus, although the linear response to forcing becomes apparent when using strong trends over 30 years, other forms of variation remain as a significant confounding factor. There are hints, such as the grouping of positive anomalies around year 160 for F of 0.3 and under, and around year 163 for F of 0.4 and over, that this residual variability is not purely stochastic. We will return to this question in section 5.

Energy Balance
Using the BCPE, we can examine the basic character of the domain-mean forced-dissipative response to various wind stress trends over years 150-180 in terms of energetics. The following energy equation is similar to that in Hogg and Blundell [2006], but the Q-GCM user guide contains a full derivation: where W5q Ð Ð u 1 s x dA is the work done by the wind; B ML 52q Ð Ð g 0 1 g 1 e 1 dA is the thermodynamic source of energy at the interface (as described in section 2.2); the rates of change of potential and kinetic energy The total frictional dissipation is due to a combination of linear bottom drag, D5 qd Ek jf0j 2 Ð Ð ðu 2 3 1v 2 3 ÞdA, and the lateral viscous dissipation, V5R 3 k51 A 4 qH k Ð Ð ðv k r 4 H v k 1u k r 4 H u k ÞdA. The left-hand side of (3) is the power input, the first two terms on the right-hand side are the rate of change of total energy and the last two terms are dissipation by friction. terms have been partitioned into variability at the frequencies introduced by the harmonic wind stress forcing (Figure 7, left plots) and the residual variability, including trend (Figure 7, right plots). The temporal harmonic analysis has been performed over years 130-180, at each value of F. An earlier version of this paper calculated the harmonic analysis over years 150-180, and there is negligible difference. All of the left plots of Figure 7 are plotted with the same scale for direct intercomparison. Figure 7a shows the harmonic component of the work done by the wind, W. At these forcing periods, W balances d dt ðPEÞ (Figure 7f) remarkably well for all values of F. There is also very little variation in either term as F is varied, suggesting a separation in the response to forcing through harmonic and trend components of wind stress. The sum W1B ML 2D2V (Figure 7e) is a slightly better match to d dt ðPEÞ (Figure 7f), with the wind-driven buoyancy loss from the mixed layer to layer 1 (Figure 7b) being the next most important term. The power input by the wind acts to increase potential energy by steepening isopycnals on the seasonal forcing scales and this is not extracted significantly by conversion to kinetic energy (Figure 7g) or through friction (Figures 7c and 7d) on these particular time scales. The weak response of d dt ðKEÞ on the seasonal forcing scales might reflect the effect of baroclinic instability in the generation of EKE, which dominates KE, at time scales of this order. Baroclinic instability complicates the relationship between forcing and response over the nonlinear growth phase. However, this does not imply that there may not be a rapid response of KE to seasonal wind stress forcing, simply that this response is not phase locked to the forcing.
Over time scales different from the seasonal wind forcing (Figure 7, right plots), there is a different story. The work done by the wind (Figure 7a) increases in time and with F, in a similar pattern to the EKE in Figure  5a. An estimate of observed wind work from the zonal integral in Figure 1 of Hughes and Wilson [2008] for this latitude range is 8310 23 W m 22 , so the model is in the right regime. For the Q-GCM integrations 2. For short time scales of much less than a year, but distinct from the forcing periods, there is some variability in wind work ( Figure 7a) and lateral viscous dissipation (Figure 7d), but the strongest signal is in d dt ðPEÞ (Figure 7f) and d dt ðKEÞ (Figure 7g), a balance which is consistent with equipartition of energy expected for motions on length scales centered on the relevant Rossby radius [Gill, 1982, section 7.5]. The rapid variability in energy and its insensitivity to the strength of the wind stress trend is interesting and suggests that the slow adjustment between wind work and frictional dissipation (which does increase with wind stress) might be mediated by rapid, energetic processes which are insensitive to F. It is worth adding that, although friction may be of little importance in the momentum budget of a Southern Ocean channel with topography, as examined extensively by Wolff et al. [1991], it is important energetically.
Figure 8 (left) shows power spectra for terms in (3) averaged across all 11 members of the BCPE after detrending, so does not distinguish variation with F but rather gives an overview which relates to the intrinsic variability and any nonlinear forced response. At periods shorter than about a year, the main balance is Figure 8. (left) Power spectra of the time series of terms in the energy balance of the model run, calculated as the average of the Fourier transforms of the 11 primary model runs, after detrending. Terms shown are the rate of change of (blue) potential energy, (red) kinetic energy, (pink) total energy (kinetic plus potential), together with (black) the power input from wind (including a potential energy input to the mixed layer), frictional power loss (green), and the residual imbalance (orange). (right) Time series (after subtracting a fitted seasonal cycle) of the run with the largest trend in wind stress. Colors are as before, except that orange represents the difference between power input and frictional loss (which should be balanced by dE/dt if the residual can be neglected). The orange horizontal line represents the zero offset applied to its time series.
Journal of Geophysical Research: Oceans 10.1002/2014JC010315 an exchange between PE and KE (blue and red, with pink being the sum and much smaller). Note that the spectral slopes change character at around 20-30 cycles per year, or 12-18 days, where baroclinic instability becomes important. This is similar to the baroclinic growth period diagnosed from (temporally and spatially coarser) hydrographic climatology by Williams et al. [2007], which has a fastest time-mean exponential growth at periods between 15 and 25 days. The spectral peaks of the seasonal wind stress forcing are clearly defined in the wind work or power input (black) and d dt ðPEÞ (blue). At annual to decadal periods, there is a balance between d dt ðPEÞ (blue) and friction (D 1 V) (green). Note that d dt ðKEÞ (red) plays a reduced role at these periods, a response which would be expected for motion on scales larger than the Rossby radius.
The residual is typically 2 orders of magnitude smaller than the dominant term, except at the highest frequencies, where the fact that some terms are instantaneous and others represent averages over the sampling interval leads to an imbalance. Lines representing power laws of f 24=3 and f 23 are included for illustration only.
The right-hand plots of Figure 8 illustrate this dominant balance in terms of one particular time series (the case F 5 0.5; trend is not removed here, but seasonal harmonics are). Blue is d dt ðPEÞ, with lots of highfrequency variability, and pink is rate of change of total energy, d dt ðEÞ, where E5KE1PE, showing that the high-frequency variability is an exchange between PE and KE. The power input, W1B ML (black), does not have much time dependence apart from a trend, though it does have some. Power input (black) minus frictional loss (green) is shown in orange, which is balanced by d dt ðEÞ (pink), plus the residual (not shown). These spectral plots support statements (i) and (ii) above. The multiyear periods apparent in Figure 7e lie near the long-period limit of resolution of spectral analysis, making it difficult to draw meaningful conclusions about these time scales. The spectrum of d dt ðKEÞ is generally white at periods longer than about 10 cycles per year (cpy), but does show a noticeable drop between about 0.3 and 0.1 cpy (3 and 10 year periods). In contrast, the spectrum of d dt ðPEÞ follows a power law proportional to about f 24=3 at periods longer than about a year, which singles out no particular time scale. However, we can see that such behavior does not continue to the longest time scales by looking at the character of random time series with similar spectra (Figure 3).
Spectral analysis shows no significant peaks apart from the seasonal cycle, but the characteristic time series of the rates of change of kinetic and potential energy are evidently quite different (Figure 8, right). For the well-resolved part of the spectrum at periods longer than about 10 cycles per year (Figure 8, left), the relatively flat (white) power spectrum of d dt ðKEÞ implies a red KE spectrum proportional to f 22 where f is frequency (since time differentiation represents a multiplication by f in frequency or f 2 in power). The power spectrum of the rate of change of PE over periods longer than 1 year is steeper, approximately proportional to f 24=3 (implying an f 210=3 spectrum for the PE).
The spectra (see section 4) suggest a change of behavior in KE for periods longer than about 5 years, and we can clearly see that the red f 22 behavior does not continue to much longer periods by generating a random red noise time series with the same standard deviation as the observed KE time series (Figure 3c, green). It is only after reducing the low-frequency content by applying a running box-car smoother and subtracting the smoothed curve from the original (then renormalizing to retain the same standard deviation) that the random time series (Figure 3c, black) come to resemble the actual time series. For the kinetic energy, the appropriate smoothing is 15 years or less. In the case of the potential energy, the spectrum gives no sign of a change in behavior at long periods, but comparison with a random f 210=3 time series (Figure 3b, green) again shows that this behavior cannot continue indefinitely. In this case, subtracting a curve with smoothing of around 15 years seems to work best (Figure 3b, black). This tells us that, although we do not have sufficient information to clearly see the shape of the spectrum at the longest periods, there must be a flattening of the spectral slope at periods longer than about 15-30 years in order for the time series to be qualitatively as observed in the model.

Testing the Influence of Intrinsic Variability
Given the potential importance of nonlinear dynamics in the eddy saturation mechanism, and the amplitude of the nonlinear residual EKE response seen in Figure 6d, it is unclear how representative the responses to wind stress trends in Figure 5 may be. Small perturbations to the initial conditions may lead to a different response and Figure 5 represents just one set of realizations, each starting from the same initial Journal of Geophysical Research: Oceans 10.1002/2014JC010315 condition at year 150, but with different wind stress trends over the following 30 years. We test the robustness of Figure 5 by examining three, 11 member ICPEs for F 5 0; 0.25 and 0.5 during the time window of years 162-168 (bounded by pairs of red, blue and black circles, respectively, in Figures 5a and 5b). This period is chosen to capture the interesting local EKE maximum at year 163 and subsequent minimum in year 165.8 in the F 5 0.5 case. Figure 9 shows the EKE and transport for the individual ICPE members and control simulations and also the ensemble means (thick lines) for this 6 year window, comprising a total of 198 years of simulation. It would have been desirable to examine a period longer than 6 years, and to initialize the ICPE at several different times, but computational constraints were prohibitive. The rate of spread of the initial condition perturbations tells us how sensitive the system is at that time and for the particular background initial condition and boundary conditions. As well as examining the individual time series (Figures 9a and  9c) we also measure the ensemble spread (Figures 9b  and 9d). It takes some time for the small perturbations to produce a spread in the ensemble, with a visible spread first emerging after about 6 months and being well established by 1 year from the time of perturbation. We investigate this period in more detail in Figure  10, where the ensemble range is plotted on a logarithmic scale. For both EKE and transport, the almost linear initial growth represents an exponential growth in the effect of the perturbation, with e-folding time scales between 1 and 2 weeks, which continues for several months before starting to saturate. The time scale coincides with the region of steepest spectral slope for d dt ðKEÞ and d dt ðPEÞ, as identified in section 4. This initial growth phase can be considered complete by the start of year 163, a year after introduction of the perturbation.
Focusing first on the EKE in Figures 9a and 9b, it is immediately clear that both stochastic and deterministic aspects of the time series are significant on longer time scales. The ensemble spread is comparable to the separation of ensemble mean trajectories, demonstrating that a single model run (or a single realization, which is all we can ever have in the observations) has the potential to be quite misleading concerning underlying causes of variability. However, it seems that the particular realizations plotted in Figures 5 and 6 are not misleading in their overall behavior. The ensemble mean EKE for F 5 0.25 is consistently above that for F 5 0, and that for F 5 0.5 starts higher still, growing to a maximum in year 163, then decaying until it is indistinguishable from the F 5 0.25 case by year 165.8, just as we saw in the single realization. In fact, every single realization of the F 5 0.5 case has an annual average EKE which reduces from a higher than average value in year 163 to lower than average in year 165.8. We must therefore conclude that, despite the short e-folding time scale for perturbations to grow, there is some aspect of the system which has a multiyear memory.
Looking back to Figure 6a, the impression is given that, for larger wind stress trends, the EKE response is not simply linearly related to the forcing, but includes the excitation of an oscillation with period about 15 years. Whether this is a genuine oscillation, or simply a nonlinear response to the forcing, is unclear, but the time scale is consistent with those roughly identified from KE and PE diagnostics in section 4, and Figure 9 shows that this behavior is robust to the introduction of perturbations.
The response in transport, Figures 9c and 9d, is somewhat less clear, as the ensemble spread is much larger than the spread between ensemble means, and the annual cycle is much more important. However, the ensemble mean curves do suggest a hint of anticorrelation between EKE and transport. From year 163 onward, the F 5 0 case (red) always has the highest transport, and F 5 0.5 (black) starts with the lowest, but gradually rises until it is indistinguishable from F 5 0.25 (blue), a relationship which is the inverse of that found for the EKE ensemble means.
The ICPEs thus demonstrate that, while a single run can be quite misleading, the impression from Figure 6 of a linear, forced response to wind stress plus long time scale perturbations is accurate.
Furthermore, in this model the long time scale perturbations contain a substantial component which is deterministic; if replicated in the ocean, this component could be considered to be a potentially predictable, though nonlinear, response to the increasing wind stress.

Summary and Conclusions
We have explored the sensitivity of the response of an idealized, eddy-resolving model of the Southern Ocean to external forcing and intrinsic variability. We find that several significant forced and intrinsic time scales emerge. The eddy saturation response of EKE to increased eastward wind stress forcing is present at lowest order, with changes in EKE being linearly proportional to wind stress, but at a time lag indistinguishable from zero. This raises the question of how typical or atypical was the 2-3 year lag found by Meredith and Hogg [2006], and how one should fully account for intrinsic variability within the observed system or within numerical models, particularly where nonlinear eddy processes may play a central role.
A forced response in EKE, not linearly related to the zonally uniform wind stress, apparently oscillatory and with period of around 15 years emerges (most obviously for the larger wind stress scenarios) above the intrinsic variability (which saturates at around 6 months). The existence of this time scale gap supports potential predictability of Southern Ocean EKE over interdecadal periods. Energetics shows this mode to be associated with long-period variability in PE, which is balanced by fluctuations in the rate of viscous dissipation. KE plays a negligible direct role at these time scales, but the modulations in dissipation rate seem to be associated with modulation in EKE: there is a greater rate of dissipation when EKE is higher, as would be expected. Both lateral viscosity and bottom drag are responsible for these dissipation changes, in a ratio of about 2:1.
On the seasonal time scales of the wind forcing, there is an exchange between wind work and potential energy, where the wind simply forces the slope of isopycnals. EKE plays a negligible role in this forced Journal of Geophysical Research: Oceans 10.1002/2014JC010315 response, showing that the EKE response to wind stress is not instantaneous. For time scales shorter than a year, but distinct from those of the seasonal forcing, there is a rapid exchange between potential and kinetic energy, most likely associated with baroclinic instability and associated short-wavelength processes.
The intrinsic variability is characterized by exponential growth of perturbations with e-folding scale of 1-2 weeks, manifest through the eddy field and associated with the dominant baroclinic instability. Perturbations continue to grow until about 4-6 months, after which time the ensemble spread of EKE and transport saturate. For the three choices of external forcing that we were limited to, and for these two metrics, there does not appear to be strong variation in this behavior. Another consideration is how optimal the ICPE perturbations are. This depends on the metric in question, the time scales being studied and the size limitations of the ensemble. There is a hint that the chosen perturbations were nonoptimal for transport in the first few weeks because there was a decrease in ensemble range (Figure 10b). However, optimal perturbation strategies for ocean models are still being actively developed [Tziperman and Ioannou, 2002;Zanna and Tziperman, 2008;S evellec and Fedorov, 2010;Zanna et al., 2011], and the choice of a simple method seems sensible in the absence of a clear alternative. The MESO-CLIP project between the UK and partners in France and Australia (see http://projects.noc.ac.uk/meso-clip) is currently comparing simple and optimal perturbation schemes.
The mode of energy exchange seen in the BCPE experiments is the same as that described by Hogg and Blundell [2006], which they attributed to a cycle in which PE builds up in a more zonal ACC, until interaction with topography breaks the zonality, causing enhanced instability and a burst of eddy energy, which enhances dissipation and thus reduces the PE. It seems that this mode is also playing an important role in the channel flow response to changed wind stress, although the relationship between energy and zonality is less simple. Over years 162-168 of the ICPE, maps of depth-integrated, ensemble-and time-mean EKE and transport stream function (not shown) highlight two regions near topography, each covering about a quarter of the domain, where increased wind stress leads to either increased local EKE and spin-up of standing recirculations in the sense of causing a more nonzonal flow, or decreased local EKE and a more zonal flow. There is growing evidence for the importance of localized regions of standing waves, associated with topography, in exerting a strong influence over baroclinic instability, jet structure and transport [e.g., Thompson et al., 2010;Thompson and Naveira Garabato, 2014;Abernathey and Cessi, 2014]. The presence of blocking by such recirculations in at least one part of the domain may explain the apparent slight reduction of eastward transport. Our experiments do suggest a linear relationship between wind stress and nonzonality on the longest time scales. In addition to an eddy response that is linearly related to wind stress, as emphasized by Meredith and Hogg [2006], there is significant stochastic variability which can cloud interpretation of observations, and potentially a forced response, not linearly related to wind stress, with a natural time scale of about 15 years, with the dynamic balance first identified by Hogg and Blundell [2006].