Estimating sea-surface temperature transport fields from stochastically-forced fluctuations

In this article we introduce a method to quantify the transport of sea surface temperature (SST) from SST fluctuations. Previous studies have estimated the advective transport of SST from time-lag correlation of SST anomalies. However this approach does not consider diffusive SST transport or relaxation to atmospheric temperatures. To quantify the transport more completely we use a response function (Greenʼs function) which solves the SST continuity equation for an impulsive forcing. The response function is estimated from SST anomalies using a fluctuation-dissipation approach. An assumption of spatial locality in the linear operator used to produce the response significantly improves the accuracy of the method. Using 100 years of data from a stochastically-forced prototypical SST model, the method estimates the SST transport response function to within 10% error. Decomposing the linear operator into symmetric, anti-symmetric, and divergent operators enables estimates of the modelʼs spatially dependent velocity vector, diffusivity tensor, and relaxation rate which converge at the same rate as the response function.


Introduction
Accurate estimates of sea-surface temperature (SST) transport are fundamental to predicting climate change. One method to obtain such estimates is to analyze the propagation of anomalously warm and cold patches of water (SST anomalies) in surface currents. Advective transport times are often identified by the time-lag of greatest correlation between time series of SST anomalies at remote locations [1][2][3]. The disadvantage of this approach is that the modal time of the two-point correlation function does not necessarily occur at the advective time in flows where SST anomalies are advected by currents, diffused by eddies, and relaxed to atmospheric temperatures.
A more complete quantification of the transport is given by the response function (Greenʼs function) which solves the SST continuity equation for an impulsive forcing. Instead of a single advective time, the Greenʼs function gives the distribution of travel times that arise from the complex interactions of advection, diffusion and relaxation. This approach has provided a foundation for quantifying the transport of a variety of tracers in geophysical reservoirs [4][5][6]. If SST arises from a stochastic atmospheric forcing [7] and is advected and diffused over certain temporal and spatial scales [8], then the transport response function is related to correlation functions by the covariance structure of the atmospheric forcing [9].
An alternative method of estimating transport from SST anomalies may be possible using the fluctuation-dissipation theorem. Under certain conditions, natural fluctuations of a system can be used to construct a linear operator that estimates the systemʼs response to an external forcing [10,11]. It has been suggested that the climate has properties which approximately satisfy the required conditions [12,13] and the theorem has recently been used to estimate response functions in global circulation models [14,15]. If the system can be modeled by a set of stochastically-forced linear differential equations, then estimating the linear operator amounts to an inverse modeling problem [16,17]. The linear stochastic model has shown to be a good approximation of SST anomaly dynamics in the tropical Pacific [18], North Atlantic [19], and recently on a global scale [20].
This article combines and builds upon the body of work outlined above. We use a fluctuation-dissipation approach including the inverse linear method to estimate the timeaverage SST transport response from SST anomalies. We make a spatial locality assumption in the calculation of the linear operator that significantly improves the accuracy of the response function and enables estimates of individual field parameters in the transport equation. In section 2 the method is described for a generic physical system. In section 3 we apply the method to a prototypical numerical model of SST transport. A summary of the results and future work is provided in section 4.

Time average of a continuous model
The method we present applies to physical systems that can be usefully modeled by stochastically-forced linear partial differential equations, Here the field variable c and the stochastic source term q vary over space x and time t. The linear operator  contains spatial derivatives and coefficients that also depend on space and time.
Let each term in (1) be the sum of a time-mean (denoted by¯) and fluctuation (denoted by ′ ) component, We are interested in statistically-steady solutions to (1) so we assume thatc,  andq are well defined and disregard transient behavior associated with the initial condition. Inserting (2) into (1) and subtracting the time-mean leaves an equation for fluctuations in the field variable, where Here we have created a time-varying forcing term ′ f which sums fluctuations in the source term with fluctuations in the operator acting on the field. Note that this forcing has zero time mean because  ′ = ′¯= q c 0 and   ′ ′ − ′ ′ = c c 0.

Spatial discretization
A spatial discretization of (3) is Here ′ t c ( ) and ′ t f ( ) are × M 1 vectors of field and forcing fluctuations, respectively, at M spatial locations. The × M M matrix B is a discretization of the time-average continuous operator  . The system is stable if all eigenvalues of B have negative real parts. It has been shown that (5) provides an accurate model of SST anomaly dynamics for characterization and predictive purposes [18,19].
The × M M response (Greenʼs) function t G( ) is given by (5) with an impulsive forcing at each location, where δ t I ( ) is an × M M identity matrix multiplying a Dirac delta function. The solution of (6) is = t G( ) e t B , and the solution of (5) is t where we have assumed causality Note that the frequency content of ′ c in (7) results from the smoothing of ′ f by G. If the forcingʼs temporal autocorrelation is negligible (white in time), which has been found to be a good approximation for the forcing of linearized SST dynamics at monthly timescales [18,19], then the timescales of variability in ′ t c ( ) results from the distribution of timescales in G. In this case the timescales in G determine the length of a sample of ′ c that is required to accurately estimate the systemʼs equilibrium statistics.

Response estimation-global method
Now consider that we possess a sufficiently long timeseries of ′ t c ( ) and wish to produce estimates (denoted byˆ) of the discrete operator B and response function G(t). Using the fluctuation-dissipation approach by way of the inverse modeling method [16], we multiply (5) by ′ c T and then average over time, As discussed above we assume the temporal autocorrelation of the forcing is negligible. This implies that the covariance between ′ c and ′ f at zero time lag is also negligible, and a least squares minimization of ′ ′ f c T provides an estimate for the discrete operator, The corresponding estimated response function iŝ =t G( ) e . (10) t B We refer to (9) as the global method because ′ ′ c c T contains all × M M covariances in the system. In practice the covariance matrices of high dimensional systems often have high condition numbers and cannot be inverted accurately, so measures are taken to reduce the dimensionality of the system. Empirical orthogonal function (EOF) decomposition is a common strategy for dimensionality reduction (e.g. [19]), but our interest is in quantifying the transport in physical space instead of EOF space, so we take a different approach.

Response estimation-local method
The local method assumes spatial locality in the linear operator in order to reduce the dimension of the inverted covariance matrix in (9). This approach assumes that the dominant physical processes in  can be expressed as local spatial derivatives, such as advection and diffusion. In addition we assume that the spatial scale of these process acting over the timestep of interest is smaller than the size of the stencil used in the localization. We are interested in monthly SST propagation, so the size of the stencil must be larger than the characteristic advective-diffusive propagation distance of an SST anomaly in one month.
If the discretization of  ′ c at each grid point m is confined to a local spatial stencil of N neighbors surrounding m, then B is a sparse matrix with sparse rows b m for = m M 1, 2, 3 ,..., , .
Let the × N 1 row vectorb m be the stencil-only subset of non-zero elements in b m . Following steps (5)-(9), the estimate of B is (12). This essentially removes the noisy, non-neighboring covariances from the calculation. Assuming the conditions of locality mentioned above are appropriate, this will improve the estimates of B and t G( ) for a sample of ′ t c ( ) of a given length, compared to the global method.

Decomposing a centered stencil
For a centered stencil such as the one to the right, we can separateB into components that act symmetrically, anti-symmetrically, and divergently about the central point. In contrast to the common symmetric-antisymmetric matrix decomposition, this method pairs matrix elements on the same rowB mn andˆ′ B mn , where for each stencil neighbor n of center cell m, the opposing neighbor is ′ n .
Lettingˆ=ˆ+B S A, the operator that acts symmetricallyŜ iŝ =ˆ+ˆ′ ( ) mn mn mn and the operator that acts anti-symmetricallyÂ iŝ mn mn mn for all neighbors n of all grid cells m. If sums along the rows ofB are not zero, thenB contains a divergent partĴ which is diagonal, mm n mn This decomposition is useful whenŜ,Â, andĴ, can be used to distinguish local physical processes in  . In the next section we useŜ to estimate the SST diffusion,Â to estimate SST advection, andĴ to estimate the rate at which SST relaxes to atmospheric conditions.

Model description
We now apply the method described above to a model of SST transport. The model is a stochastically-forced two-dimensional heat transport equation in a Cartesian plane.
which is the time-average of a velocity vector = u , and relaxation rate r t x ( , ). The diffusivity is a continuous parameterization of all molecular and turbulent mixing processes occurring below a finite spatial scale Δ, the velocity represents the average advective current over Δ, and the relaxation rate is the timescale for an SST anomaly of size Δ to relax to the atmospheric temperature. The structure of the flow is illustrated in figure 1. The values of the flow parameters are chosen to represent observed SST transport [1].
The continuous equations are spatially discretized to satisfy (5). The rectangular grid   The modelʼs forcing is generated randomly with a defined covariance structure and SST fluctuations are produced by integrating (5) numerically with a time step Δt = 1 month. The forcing has zero temporal decorrelation time (white in time) and its spatial decorrelation distance is tuned to produce SST anomalies characteristic of those seen in the Gulf Stream region of the Hadley Centre Sea Ice and Sea Surface Temperature dataset [21]. The time series below the schematic in figure 1 shows 10 years of SST fluctuations for model locations a b , and c. The anomaly signal is obtained by removing the average 12-month seasonal cycle from the SST fluctuation time series.
The transport between locations a b , and c is quantified by elements G t ( ) ab , G t ( ) bc , and G t ( ) ac of the true response function t G( ) as shown (by the solid lines) in figure 2. The → a b response arrives sooner and has a larger magnitude than the → a c response because a is nearer to b than c. However it is not the case (as might naively be assumed) that the time-lag of peak response results only from the advection, the width results from the diffusion, and the magnitude results from the relaxation. Instead these properties of the response function arise from a complex interaction of the individual mechanisms.

Estimating the transport response and flow parameters
We now estimate the modelʼs SST transport response function and flow parameters using the SST anomaly time series. We estimate the discretized operatorB using the local method (12) where ′ t c ( ) is a finite length time series of monthly SST anomalies. The estimated responsê t G( ) is calculated by insertingB into (10). Figure 2 shows elements of the estimated responsê   Figure 3 shows estimates of the flow parameters using 100 and 1,000 years of SST model output. The true value of each flow parameter (right column) has a unique spatial pattern which is visible in the estimates made with 100 years of data (left column). With 1000 years (center column), the errors in the estimates is 3-5%. Figure 4 shows the convergence of κ¯ū G, , xx andr as the length ℓ of the SST time series increases from 10 to 10,000 years. The error decreases as ℓ 1/ for each quantity. For a given ℓ, estimates of the velocity and relaxation rate are more accurate than the diffusivity and less accurate than the response. The difference in accuracy of the flow parameters appears to be related to their relative magnitudes in the transport operator (the magnitude of the average velocity in the antisymmetric operator is 4 times greater than the average diffusivity in the symmetric operator). We expect the response to be more accurate than the flow parameters because the calculation of the response (10) averages random errors in the parameter fields. A more thorough analysis of these differences is an area for future investigation. Figure 4 also shows the error in the response function estimated by the global method (9) which does not assume locality in the linear operator. This method cannot reliably produce a response estimate with less than 100 years because the condition number of ′ ′ c c T is too large. Between 100 and 10,000 years the error in the global methodʼs response decreases only marginally. The reason is that even when ′ ′ c c T can be inverted, it is dominated by the noisy covariances of non-neighboring locations. Inverting the covariance matrices in the local method (12) is more accurate because non-neighboring covariances have been removed. A local assumption is also required to estimate the velocity, diffusivity, and relaxation rates, so they do not appear for the global method in figure 4.

Summary and future work
The goal of this study is to demonstrate a method for estimating SST transport information from a field of SST anomaly time series. Using 100 years of output from a stochastically-forced SST model ( figure 1) the method estimates the SST transport response function to within 10% error (figure 2). A key aspect of our approach is assuming spatial locality in the linear operator (12) used to construct the response. This significantly reduces noise that otherwise prohibits the required inversion of the anomaly covariance matrix. The locality assumption also enables estimates of the modelʼs spatially dependent velocity vector, diffusivity tensor, and relaxation rates (figure 3) which converge at the same rate as the response function ( figure 4).
We intend to use this method to estimate response functions in global circulation models and climatological datasets. We hope the method will provide independent estimates of SST transport in the North Atlantic to compare with previous studies. The transport of other oceanographic properties such as salinity or ocean color can be estimated using the same approach. More generally, the method presented here may be used to estimate atmospheric or other climatic response functions when the spatial locality assumption applies. , and relaxation rate (r). The true value of each flow parameter (right column) has a unique spatial pattern. This pattern is visible in the estimates made with 100 years of SST data (left columnn). After 1000 years (center column) the average error in the flow parameter estimates is 3-5%. xx and the base 10 logarithm ofĜ. Using the local method (12) the error decreases as ℓ 1/ . The locality assumption also makes it possible to estimate the velocity vector, diffusivity tensor, and relaxation rates which converge at the same rate. The global method (9) cannot reliably estimate a response with < ℓ 100 years because the condition number of the covariance matrix is too large. For < < ℓ 100 10, 000 years, the error in the global methodʼs response function decreases only marginally due to noisy covariances between non-neighboring grid locations.