Seasonal variability of the . . . numerical investigation based on transfer operators

Abstract. The detection of regions in the ocean that are coherent over an extended period of time is a fundamental problem in many oceanic applications. For instance such regions are important for studying the transport of marine species and for the distribution of nutrients. In this study we demonstrate the efficacy of transfer operators in detecting and analysing such structures. We focus first on the detection of the Weddell and Ross Gyre for the four seasons spanning December 2003–November 2004 within the 3-D oceanic domain south of 30° S, and show distinct seasonal differences in both the three-dimensional structure and the persistence of the gyres. Further, we demonstrate a new technique based on the discretised transfer operators to calculate the mean residence time of water within parts of the gyres and determine pathways of water leaving and entering the gyres.


Introduction
The detection of regions in the ocean that are coherent over an extended period of time is a fundamental problem in many oceanic applications.Such coherent structures, for instance, have considerable impact on the movement of heat around the planet; they are important for the transport of marine species and for the distribution of nutrients, and they can play major roles in the mixing and redistribution of water masses.
In this study we focus on two such structures, namely the subpolar gyres of the Southern Hemisphere.These gyres are crucial for the physical and biogeochemical characteristics of the region.For instance, water whose properties are modified in the Weddell Gyre, makes up approximately 70% of the Antarctic Bottom Water that will form the abyssal wa-Correspondence to: M. Dellnitz (dellnitz@math.upb.de)ters in the Southern Ocean (Carmack, 1977).The Weddell Gyre also exports large quantities of carbon dioxide, formed through the remineralisation of falling organic material, into the deep ocean.As such, it plays a significant role in global climate changes (Hoppema, 2004).Similar processes are likely occurring in the Ross Gyre.
Coherent structures at the ocean surface can often be observed by satellite altimetry (see e.g., Fu, 2006).To assess internal flow fields, indirect estimates can be made based on the pressure fields within the ocean, derived from temperature and salinity measurements.However such observations are generally sparse in both space and time making it hard to identify in detail coherent structures in the deep ocean.In general, to resolve the full three-dimensional flow field, and subsequently deep coherent structures, we must resort to ocean general circulation models that numerically solve the equations of motion.
Given a vector field that describes the ocean flow we can use techniques developed in the dynamical systems context to detect and analyse the regions of interest.There are essentially two approaches, the geometric and the probabilistic approach, see Froyland and Padberg (2009) for a detailed comparison.The geometric approach aims at detecting barriers to transport such as finite-time invariant manifolds via finite-time Lyapunov exponents or related concepts.This technique has been very successfully employed for analysing transport in many applications (Shadden et al., 2005) and was introduced in Haller and Yuan (2000) and Haller (2001).However, in a previous study it was found to perform poorly in the context of the Southern Ocean (Froyland et al., 2007).In this contribution we will use the probabilistic approach, which attempts to detect coherent regions directly via a transfer operator.The two previous investigations of the Southern Ocean (Froyland et al., 2007(Froyland et al., , 2008) ) have demonstrated the efficacy of this method.
Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union.
M. Dellnitz et al.: Seasonal variability of the subpolar gyres in the Southern Ocean This report attempts to expand the spatial scales used in previous studies to cover the entire Southern Ocean in order to identify coherent structures in the Weddell and Ross Seas.In particular, we extend the domain of interest from our previous studies of the surface circulation (Froyland et al., 2007), and from the surface to a depth of 500 m (Froyland et al., 2008), to include the full water column and study the seasonal variance of the gyres.In addition, we use the technique to estimate residence times for water elements originating from within the gyre.We also show that our technique provides an efficient method for examining pathways for transport.Such water mass tracking is commonly done using Lagrangian particle or passive tracer experiments.

Input data and non-autonomous flow model
Our input data consists of monthly averaged output of the ORCA025 global ocean model (Barnier et al., 2006).We restrict ourselves to the period December 2003-November 2004.The ORCA025 model is integrated at 1/4 • resolution and so is classed as eddy permitting, meaning that it can resolve some of the larger mesoscale oceanic eddies that make up an important component of ocean circulation.The output consists of global 3-D discretised velocity fields and some other important ocean properties like sea surface height or mixed layer depth.In the Southern Ocean, the model grid follows a Mercator projection.
For our investigation of the Southern Ocean we work on a subset X of a solid annulus ] with S 1 parameterised in degrees from −180 • to 180 • .X is formed by the removal of continents, islands, and sea floor topography.The extension of the domain in the depth direction compared to our previous investigations (Froyland et al., 2007(Froyland et al., , 2008) ) forces this restriction in latitude to −48 • S for numerical reasons.As the focus is on coherent structures in the Ross and Weddell Seas this constraint should not effect our results.
Considered as a non-autonomous dynamical system, the ocean flow may be described by (x,t,τ ) → φ (x,t;τ ), where φ : X × R × R → X and φ (x,t;τ ) is the terminal point in X of a trajectory beginning in x ∈ X at time t and flowing for τ time units.The domain X is the oceanic domain that is covered by the output of the ORCA025 model and our domain of interest X is a subset of X.A trajectory x(t) := φ (x 0 ,t 0 ;t) is a solution to the non-autonomous ODE with initial condition x(t 0 ) = x 0 .Here x 0 ∈ X and f : X × R → R 3 is smooth.The vector field f is approximated by the output of the ORCA025 model.For our computations we concentrate on trajectories for which the initial and terminal points are in X.However, some of these may pass into X \X and then reenter X.
3 Almost-invariant sets, coherent structures, and transfer operators Let µ denote the three-dimensional volume measure, normalised so that µ(X) = 1.The measure µ is invariant under φ, i.e. for each τ ≥ 0, µ(φ (A,t;−τ )) = µ(A).We will say that a set A ⊂ X is φ-invariant over [t,t + τ ] if A = φ (A,t + s;−s) for all 0 ≤ s ≤ τ .Coherent structures obey an approximate invariance principle over short periods of time.We shall call a set The ratio in Eq. ( 2) is the proportion of the set A that remains in A under the flow from time t to time t + τ .In the sequel, we refer to ρ t,τ (A) as the coherence ratio of a set A. Clearly, the closer this ratio is to unity, the closer the set A is to being invariant.In order to discover coherent structures in the flow φ, we seek dominant almost-invariant sets, i.e. we seek sets for which the ratio in Eq. ( 2) is close to 1.The notion of almost-invariance arose as a means of discovering dominant macroscopic structures in general dynamical systems (Dellnitz and Junge, 1999) and has been refined and applied in a variety of settings, e.g.Dellnitz et al. (2005); Froyland and Dellnitz (2003); Froyland (2005).The framework introduced above can describe structures that are almost-invariant and fixed in the state space.To handle time dependent coherent regions i.e. coherent structures that move in state space, extended transfer operator frameworks are needed, such as those recently developed in Froyland et al. (2009); Santitissadeekorn et al. (2009).In the present work, we are focusing on spatially fixed gyres, and the time-independent framework is appropriate.In order to locate these almost-invariant sets we introduce a transfer operator describing flows for short periods.
We define the Perron-Frobenius operator or transfer operator P t,τ : L 1 (X,m) → L 1 (X,m) by where m denotes the normalised Lebesgue measure on X.
If there is a φ-invariant set A ⊂ X over [t,t + τ ], then P t,τ χ A = χ A , where χ A is the characteristic function and defined by χ A (x) = 1 if x ∈ A and 0 otherwise.Thus χ A is an eigenfunction of P t,τ with eigenvalue 1. Sets A that are almost-invariant correspond to eigenfunctions of P t,τ with real eigenvalues very close to 1 (Dellnitz and Junge, 1999;Froyland and Dellnitz, 2003).
To access these eigenfunctions numerically, we construct a finite-dimensional Galerkin approximation of P t,τ based on a fine box-partition {B 1 ,...,B n } of X.Following Ulam's approach (Ulam, 1960) we form the transition matrix The matrix P t,τ is stochastic, i.e. n j =1 P t,τ,i,j = 1∀i ∈ {1,...,n}.The entry P t,τ,i,j may be interpreted as the probability that a point selected uniformly at random in B i at time t will be in B j at time t + τ .In our study we fix both starting and flow time depending on the season we want to investigate.

Oceanic domain and discretisation
To calculate the entries of the transition matrix P t,τ each partition element B i ,i = 1,...,n, is filled with N uniformly distributed test points y i,l ∈ B i ,l = 1,...,N.For each i = 1,...,n we calculate φ y i,l ,t;τ ,l = 1,...,N, by numerical integration and approximate P t,τ by where # denotes the cardinality of a set.The boxdiscretisation of X and the construction of P t,τ is carried out efficiently using the software package GAIO (Dellnitz et al., 2001).

Model interpolation and trajectory integration
The ORCA025 model velocity fields for each month are given at a resolution of 1/4 degrees of longitude and latitude, and 46 non-uniform depth layers.The velocity field f (x,t) for non-grid points is produced by linear interpolation independently in each direction.The calculation of φ y i,l ,t;τ is carried out using a standard Runge-Kutta approach.The flow time τ and the box-size for each application is chosen in such a way that a vast bulk of test points leave their starting box within time τ .

Computation of ML mixing
The mixed layer (ML) refers to the near surface portion of the water column that is subject to vigorous mixing.This results from the action of the winds, waves and from surface cooling which can lead to deep convection particularly during winter.Mixing in the ML homogenises properties leading to near constant temperature and salinity, and consequently density.This layer extends from the surface down to the mixed layer depth (MLD).The monthly mean MLD is a diagnostic output of the ORCA025 model.It is defined as the depth at which the density exceeds the surface density by 0.3 kg•m −3 ; see Griffies et al. (2008).The shallowest simulated MLD in 2004 is −400 m in winter and −90 m in summer.The ML mixing is not captured by the time-averaged velocity field f and consequently must be parametrised during the tracking of Lagrangian particles.We make a standard assumption that the ML is well mixed.In order to si-mulate this with our sample trajectories, we proceed as follows.We write y ∈ X as y = (Lon(y),Lat(y),Dep(y)), and let MLD(t,Lon(y),Lat(y)) denote the mixed layer depth for the water column located at (Lon(y),Lat(y)) at time t.With this notation we can simulate the ML mixing for a test point y i,l at starting time t 0 , using the following algorithm: 1. Integrate the test point y i,l for one month (τ = 1) to produce y := φ y i,l ,t 0 ;τ .
3. Integrate y for a further month, set t 0 = t 0 + τ and go back to step 2.

Eigenvalue and eigenfunction calculation
Let A = i∈I B i , with I ⊂ {1,...,n}.Define where X = n i=1 B i and V (B i ) denotes the physical volume of a box in the oceanic domain.Then it is straightforward to show that (Froyland and Dellnitz, 2003) In fact, in the limit as n → ∞ and the diameter of the boxes {B i } n i=1 approaches zero, one obtains equality.We transform the matrix P t,τ into a "time symmetric" matrix R t,τ via R t,τ,i,j := P t,τ,i,j + Pt,τ,i,j /2, where Pt,τ,i,j := p j P t,τ,j,i p i is the transition matrix of the inverse Markov Chain.The matrix R t,τ is stochastic, has only real eigenvalues (Brémaud, 1999), and satisfies important maximisation properties (Froyland, 2005) related to almost-invariance.Following Froyland (2005) we use the right eigenvectors v (k)  of R t,τ to detect almost-invariant sets.That is, we define sets for suitable c ∈ R, such that min ρ t,τ A + c ,ρ t,τ A − c is maximised.This heuristic approach is described in detail in Froyland and Padberg (2009) (see also Froyland and Dellnitz, 2003).
The matrix R t,τ is typically very sparse and we are interested only in dominating spectral values close to 1, which may be efficiently computed by Lanczos iteration methods.
For our numerical investigations we approximate the oceanic domain by 141 426 boxes.This is done by first creating a box partitioning of the whole domain X and then deleting all boxes which are not in the oceanic domain, i.e. we removed all continents and islands.For this we choose M points z l,i l=1,...,M for every box B i and set As in our preceding investigations the individual boxes have a size of 1.4 • in both longitude and latitude.The vertical extent of boxes is non-uniform corresponding to the 46 nonuniform depth levels of the ORCA025 model grid, varying from 6.2 m near the surface to 250 m at the bottom of the ocean.
In order to obtain results for the four seasons in the Southern Hemisphere we compute transition matrices with t 0 equal to the first day of December, March, June, and September for the summer, autumn, winter and spring seasons, respectively.The flow time in each season is chosen to be τ = 90 days, compared to 60 days in our previous studies.This longer flow time is not only necessary in order to be able to study a whole season but also to guarantee that, as the box size has been increased, most of the particles will leave their initial box within time τ .For each box we use N = 512 uniformly distributed test points and a 3-day time step for the Runge-Kutta integration to compute (5).The integration step size is chosen in such a way that the vast bulk of trajectories will flow only to a neighbouring grid set in the ORCA025 velocity grid.The simulation of the ML mixing is done as described in Sect.4.3.
For each season we calculate the 40 largest eigenvalues and the corresponding eigenvectors.Unlike in our previous studies (Froyland et al., 2007(Froyland et al., , 2008) ) the Ross and Weddell Sea gyres are typically not captured by a single eigenvector.This may be due to substructures within the gyres.For the numerical investigation of the gyres we consider the union of such substructures.As almost all of the eigenvalues we use range from 0.98 to 0.99, we will have similar coherence values for the joined substructures.The eigenvalue rankings are given in Table 1.As the Weddell and Ross Sea features occupy a relatively small fraction of the Southern Ocean, the eigenvalues corresponding to these structures appear some way down the spectrum.
An important quantity for describing the coherence of the identified structures is the mass flux of water through the boundary of the structures.A box B is on the boundary of a structure when it has a neighbour box outside the structure.Let A= i∈I B i be a coherent structure and J ⊂I the index set of boxes on the boundary of the structure.The flux flux B j through a boundary box B j , with j ∈J , can be cal- Summer v (11)  0.9265 v (7)  0.9235 v (16)   Autumn v (11)  0.9112 v (5) 0.9179 v (17)   Winter v (13) 0.9106 v (4) 0.8865 v (39)  v (12)   Spring v (13) 0.9190 v (2) 0.9042 v (13)   culated as follows.We define the flux through a boundary box B j out of A by flux out B j := V B j k∈{1,...,n}\I and the flux through the box B j from outside A by This leads to flux B j = flux out B j − flux in B j .
For visualisation we define flux B j also for inner boxes B j of a coherent structure A by setting flux B j := 0. Obviously the flux of a box B j is bounded above by V B j and below by −V B j , respectively.

Coherent structures in the four seasons
In this section we present the results of our computations.Table 1 lists the eigenvectors that we have manually identified as coherent structures within the Weddell and Ross Seas, respectively.Other eigenvectors, not listed in Table 1, represent coherent features outside these regions.The values ρ t,τ (A W ) and ρ t,τ A R are the coherence ratios of the union of the identified substructures in each sea and season.The coherence ratio ρ t,τ A W = 0.9265 for the Weddell Sea in summer means for instance, that 92.65% of the water mass is retained in the identified structure after three months of flow.The transfer operator approach is able to detect coherent structures in the Weddell and Ross Seas during each season; these structures can be identified with the respective subpolar gyres.
Figure 1 shows a three-dimensional view of the gyres detected in summer.The colour scale indicates fluxes at the boundary of the gyres and gives an indication of regions where water enters or escapes from the gyre.Figure 2 shows the boxes of the detected gyres on the surface.The contouring indicates, for example, that most of the water escaping from the Weddell Gyre does so at the poleward limits of the structure.Considerable differences are evident in the shape of the gyres over the four seasons.For example, between autumn and winter, the extent of the Ross Gyre shrinks considerably and the Weddell Gyre's northern boundary contracts southwards while its eastern boundary extends considerably further east.The seasonal changes in structure of the Weddell Gyre can be seen in Fig. 3, which shows a zonal section along −64 • S.During spring and summer the Weddell Gyre extends from the surface to depths exceeding 4000 m, while in autumn and winter the main structure of the gyre is subsurface and extends to much shallower depths.These seasonal changes are also manifest in the volume of the gyres in the different seasons (Table 2), which suggests greater seasonal variability for the Weddell than the Ross Gyre.The colour scale in Fig. 3 indicates the mean residence time of water in the Weddell Gyre which is explained in detail in the next section.
Finally we remark that it would require an analysis spanning multiple years in order to test how robust these seasonal differences are to interannual variability.As such our analysis demonstrates a proof of concept and future work would need to consider longer time frames.

Transfer operator based analysis of coherent structures
To further illustrate the strength and variability of the transfer operator approach we analyse the summer Weddell Gyre in more detail.
The transfer operator approach can also be used to calculate the variability in the residence time of water within the Weddell Gyre A W = i∈I B i for each season respectively.For example for the summer season we define a transition matrix P := P Summer t,τ  sents the dynamics over one year from the start of summer.We may assume that the restriction P| A W of the transition matrix to the set A W satisfies lim k→∞ P| A W k → 0. Thus, the mean residence time t i of a particle in B i is given by the solution of the linear equation where Id denotes the identity matrix; see e.g.Froyland (2001).Figure 3a shows a zonal section along −64 • S latitude which bisects the core of the coherent structure that we identify as the summer Weddell Gyre.The colouring indicates the average time that a particle originating at different locations of the gyre will remain within the gyre.The gyre's core lies at −1250 m, with water from this region remaining within the structure for an average of about nine years.Towards the ocean surface the residence time drops off sharply as water upwells in this region, driven by surface divergence and is advected away from the gyre.The pathway of water exiting or entering the gyre can be further investigated by sequentially applying the transition matrix for each season to the Weddell Gyre or sub-regions of it.By repeating this multiple times we can simulate the spreading out of water from any given region over multiple years.Again, this repeated application means that we are neglecting year-to-year variations in the flow field.Once   W of the final state of an initial density after 50 years forward evolution by 50 iterations of the transition matrix.Boxes are coloured according to the logarithm of the normalised density.The initial region A = i∈I B i is assigned an initial density q = (q 1 ,...,q n ) with q i = 1 |I | if i ∈ I and q i = 0 otherwise.the transition matrices have been calculated, this approach has the advantage of being numerically efficient for detecting such circulation pathways.The time-stepping of large numbers of Lagrangian particles forward (or backwards) in time can be effectively achieved through a number of matrixvector multiplications.Only the initial calculation of the transition matrices needs significant computation time.Similar techniques to identify pathways of water are described in Khatiwala (2007); Khatiwala et al. (2005).They use a discretisation of the advection-diffusion operator, which leads to matrix equations describing the transport of passive tracers in the ocean.In low-diffusion settings, a drawback of using discretisations of advection-diffusion operators directly is the addition of artificial diffusion from the discretisation; this is noted in Khatiwala et al. (2005).Our approach of computing discretisations of transfer operators representing flow over short times reduces this artificial diffusion.
We examine two examples of watermass tracking using transfer matrices.
For an initial region situated within the gyre between −65 • S and −60 • S and −30 • E and −20 • E and a depth ranging from 0 m to −2712 m, the tracer primarily leaves the gyre near the surface and spreads northwards via Ekman transport, and eastwards swept along within the Antarctic Circumpolar Current (Figs.Using the transfer operator approach we can also reverse the iteration process.A similar approach is commonly used in oceanographic Lagrangian studies, where particles are backtracked from a specified region to identify the origins of water in that region.To demonstrate this we first calculate the inverse process of the Markov Chain induced by the transition matrices (4) for each season.For a given transition matrix P t,τ the inverse Markov Chain is induced by Pt,τ as defined in (9).In our example, we start with an initial region outside the gyre between −55 • S and −50 • S and −5 • E and 5 • E and a depth ranging from −2500 m to −1500 m.The water arrives, at this region, via three distinct routes.The first route is through the Drake Passage, the second originates from the eastern side of South America and the third pathway originates from the Weddell Gyre (see Fig. 6 and supplementary animation BackwardEvolution.gif, see http://www.nonlin-processes-geophys.net/16/ 655/2009/npg-16-655-2009-supplement.zip).This last route has been previously identified as an important pathway for transporting large quantities of natural CO 2 generated within the gyre, into the deep ocean (Hoppema, 2004).

Conclusions
In this investigation we extended previous work by studying the full 3-D domain of the Southern Ocean.In particular we focused on the detection of the subpolar gyres whose structure extends from the surface to abyssal depths.We applied www.nonlin-processes-geophys.net/16/655/2009/ Nonlin.Processes Geophys., 16, 655-664, 2009 well known transfer operator techniques on this large-scale domain and showed that coherent structures, identified with the Ross and Weddell Gyres, are present during each season.Significant changes in the structure of the gyres are evident during the different seasons, although the robustness of this result to interannual variability has not yet been assessed.
We have also demonstrated the ability to use the transfer operators to investigate circulation pathways.Lagrangian particle studies (see e.g., Döös, 1995) and passive tracer experiments (see e.g., Sen Gupta andEngland, 2004, 2007) are commonly used to investigate ventilation pathways and timescales.The transfer operator technique provides a new and efficient means for conducting such experiments.We demonstrate both forward tracking of tracer away from the central Weddell Gyre, and the backtracking of tracer to investigate the origins of water in a given interior region.Thus, once our transfer operators have been computed, we can use it to analyse the underlying dynamical system via a number of different approaches.

Fig. 1 .
Fig. 1.Detected gyres in the Weddell and Ross Seas in summer.Colouring shows the volume flux of water through boundary boxes over summer.

Fig. 2 .
Fig. 2. Boxes on the surface of the Southern Ocean of the detected gyres during (a) summer, (b) autumn, (c) winter, and (d) spring.Colouring shows the volume flux of water through boundary boxes over a season.

Fig. 3 .Fig. 4 .
Fig. 3. Section of Weddell Gyre during (a) summer, (b) autumn, (c) winter, and (d) spring along −64 • S. Boxes are coloured according to the mean residence time of water in the Weddell Gyre in years.
4 and 5 and supplementary animation ForwardEvolutionDepth.gif and ForwardEvolutionLon.gif, see http://www.nonlin-processes-geophys.net/16/655/2009/ npg-16-655-2009-supplement.zip).A portion of the tracer can be seen exiting across the northern edge of the domain in the eastern Pacific as part of the eastern limb of the subtropical gyre.

Fig. 5 .
Fig. 5. Horizontal surface section of the final state of an initial density after 50 years forward evolution by 50 iterations of the transition matrix.Boxes are coloured according to the logarithm of the normalised density.The initial region A and initial density are as in Fig. 4.

Fig. 6 .
Fig. 6.Horizontal section at −1561.8 m depth of the final state of an initial density after 10 years backward evolution by 10 iterations of the reverse time transition matrix.Boxes are coloured according to the logarithm of the normalised density.The initial region A and initial density are as in Fig. 4.

Table 1 .
Eigenvectors that identify the coherent structures with the corresponding coherence ratio for each season.

Table 2 .
Volume of the detected coherent structures for each season.