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

The method we present applies to physical systems that can be usefully modeled by stochastically-forced linear partial differential equations,

$$\begin{eqnarray}&&\frac{{\rm d}}{{\rm d}t}c({\bf x},t)-\mathcal{L}({\bf x},t)c=q({\bf x},t).\end{eqnarray} \tag{ 1 }$$

Here the field variable c and the stochastic source term q vary over space ${\bf x}$ and time t. The linear operator $\mathcal{L}$ 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$\overline{{}}$) and fluctuation (denoted by $^{\prime}$) component,

$$\begin{eqnarray}\begin{array}{rcl} c({\bf x},t) & = & \bar{c}({\bf x})+c^{\prime} ({\bf x},t) \\ \mathcal{L}({\bf x},t) & = & \bar{\mathcal{L}}({\bf x})+\mathcal{L}^{\prime} ({\bf x},t) \\ q({\bf x},t) & = & \bar{q}({\bf x})+q^{\prime} ({\bf x},t). \\ \end{array}\end{eqnarray} \tag{ 2 }$$

We are interested in statistically-steady solutions to (1) so we assume that $\bar{c}$, $\bar{\mathcal{L}}$ and $\bar{q}$ 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,

$$\begin{eqnarray}&&\frac{{\rm d}}{{\rm d}t}c^{\prime} -\bar{\mathcal{L}}c^{\prime} =f^{\prime} ({\bf x},t),\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&f^{\prime} ({\bf x},t)=\mathcal{L}^{\prime} \bar{c}+\mathcal{L}^{\prime} c^{\prime} -\overline{\mathcal{L}^{\prime} c^{\prime} }+q^{\prime} .\end{eqnarray} \tag{ 4 }$$

Here we have created a time-varying forcing term $f^{\prime}$ 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 $\overline{q^{\prime} }=\overline{\mathcal{L}^{\prime} \bar{c}}=0$ and $\overline{\mathcal{L}^{\prime} c^{\prime} }-\overline{\overline{\mathcal{L}^{\prime} c^{\prime} }}=0$.

A spatial discretization of (3) is

$$\begin{eqnarray}&&\frac{{\rm d}}{{\rm d}t}{\bf c}^{\prime} (t)-{\bf Bc}^{\prime} ={\bf f}^{\prime} (t).\end{eqnarray} \tag{ 5 }$$

Here ${\bf c}^{\prime} (t)$ and ${\bf f}^{\prime} (t)$ are $M\;\times \;1$ vectors of field and forcing fluctuations, respectively, at M spatial locations. The $M\;\times \;M$ matrix ${\bf B}$ is a discretization of the time-average continuous operator $\bar{\mathcal{L}}$. The system is stable if all eigenvalues of ${\bf 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\;\times \;M$ response (Greenʼs) function ${\bf G}(t)$ is given by (5) with an impulsive forcing at each location,

$$\begin{eqnarray}&&\frac{{\rm d}}{{\rm d}t}{\bf G}(t)-{\bf BG}={\bf I}\delta (t),\end{eqnarray} \tag{ 6 }$$

where ${\bf I}\delta (t)$ is an $M\times M$ identity matrix multiplying a Dirac delta function. The solution of (6) is ${\bf G}(t)={{{\rm e}}^{{\bf B}t}}$, and the solution of (5) is

$$\begin{eqnarray}&&{\bf c}^{\prime} (t)=\int _{-\infty }^{t}{\bf G}(t-s){\bf f}^{\prime} (s){\rm d}s,\end{eqnarray} \tag{ 7 }$$

where we have assumed causality ${\bf G}(t\lt 0)=0$.

Note that the frequency content of ${\bf c}^{\prime}$ in (7) results from the smoothing of ${\bf f}^{\prime}$ by ${\bf 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 ${\bf c}^{\prime} (t)$ results from the distribution of timescales in ${\bf G}$. In this case the timescales in ${\bf G}$ determine the length of a sample of ${\bf c}^{\prime}$ that is required to accurately estimate the systemʼs equilibrium statistics.

Now consider that we possess a sufficiently long timeseries of ${\bf c}^{\prime} (t)$ and wish to produce estimates (denoted by $\widehat{{}}$) of the discrete operator ${\bf B}$ and response function ${\bf G}$(t). Using the fluctuation-dissipation approach by way of the inverse modeling method [16], we multiply (5) by ${\bf c}{{^{\prime} }^{{\rm T}}}$ and then average over time,

$$\begin{eqnarray}&&\overline{\frac{{\rm d}}{{\rm d}t}{\bf c}^{\prime} {\bf c}{{^{\prime} }^{{\rm T}}}}-{\bf B}\overline{{\bf c}^{\prime} {\bf c}{{^{\prime} }^{{\rm T}}}}=\overline{{\bf f}^{\prime} {\bf c}{{^{\prime} }^{{\rm T}}}}.\end{eqnarray} \tag{ 8 }$$

As discussed above we assume the temporal autocorrelation of the forcing is negligible. This implies that the covariance between ${\bf c}^{\prime}$ and ${\bf f}^{\prime}$ at zero time lag is also negligible, and a least squares minimization of $\overline{{\bf f}^{\prime} {\bf c}{{^{\prime} }^{{\rm T}}}}$ provides an estimate for the discrete operator,

$$\begin{eqnarray}&&{\bf \hat{B}}={\mkern 1mu} \overline{\frac{{\rm d}}{{\rm d}t}{\bf c}^{\prime} {\bf c}{{^{\prime} }^{{\rm T}}}}{{\left( \overline{{\bf c}^{\prime} {\bf c}{{^{\prime} }^{{\rm T}}}} \right)}^{-1}}.\end{eqnarray} \tag{ 9 }$$

The corresponding estimated response function is

$$\begin{eqnarray}&&{\bf \hat{G}}(t)={{{\rm e}}^{{\bf \hat{B}}t}}.\end{eqnarray} \tag{ 10 }$$

We refer to (9) as the global method because $\overline{{\bf c}^{\prime} {\bf c}{{^{\prime} }^{{\rm T}}}}$ contains all $M\;\times \;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.

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 $\bar{\mathcal{L}}$ 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 $\bar{\mathcal{L}}c^{\prime}$ at each grid point m is confined to a local spatial stencil of N neighbors surrounding m, then ${\bf B}$ is a sparse matrix with sparse rows ${{{\bf b}}_{m}}$ for $m=1,2,3,...,M$,

$$\begin{eqnarray}&&{\bf B}=\left[ {{{\bf b}}_{1}};{{{\bf b}}_{2}};{{{\bf b}}_{3}};...;{{{\bf b}}_{M}} \right].\end{eqnarray} \tag{ 11 }$$

Let the $1\times N$ row vector ${{\tilde{{\bf b}}}_{m}}$ be the stencil-only subset of non-zero elements in ${{{\bf b}}_{m}}$. Following steps (5)–(9), the estimate of ${\bf B}$ is

$$\begin{eqnarray}&&{{\hat{\tilde{\bf b}}}_{m}}= \overline{\frac{{\rm d}}{{\rm d}t}c_{m}^{\prime }{\tilde{{\bf c}}^{\prime}}_{m}^{{\rm T}}}\left( \overline{{\tilde{{\bf c}}^{\prime }}_{m}{\tilde{{\bf c}}^{\prime}}_{m}^{{\rm T}}} \right)^{-1},\end{eqnarray} \tag{ 12 }$$

where $\tilde{{\bf c}}_{m}^{\prime }$ is the corresponding $N\times 1$ stencil-only subset of ${\bf c}^{\prime}$, and $c_{m}^{\prime }$ is element m of ${\bf c}^{\prime}$.

The local method changes the single $M\;\times \;M$ matrix inversion in (9) to several $N\;\times \;N$ matrix inversions in (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 ${\bf B}$ and ${\bf G}(t)$ for a sample of ${\bf c}^{\prime} (t)$ of a given length, compared to the global method.

For a centered stencil such as the one to the right, we can separate ${\bf \hat{B}}$ 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 row ${{\hat{B}}_{mn}}$ and ${{\hat{B}}_{mn^{\prime} }}$, where for each stencil neighbor n of center cell m, the opposing neighbor is $n^{\prime}$.

Letting ${\bf \hat{B}}={\bf \hat{S}}+{\bf \hat{A}}$, the operator that acts symmetrically ${\bf \hat{S}}$ is

$$\begin{eqnarray}&&{{\hat{S}}_{mn}}=\frac{1}{2}\left( {{{\hat{B}}}_{mn}}+{{{\hat{B}}}_{mn^{\prime} }} \right),\end{eqnarray} \tag{ 13 }$$

and the operator that acts anti-symmetrically ${\bf \hat{A}}$ is

$$\begin{eqnarray}&&{{\hat{A}}_{mn}}=\frac{1}{2}\left( {{{\hat{B}}}_{mn}}-{{{\hat{B}}}_{mn^{\prime} }} \right),\end{eqnarray} \tag{ 14 }$$

for all neighbors n of all grid cells m. If sums along the rows of ${\bf \hat{B}}$ are not zero, then ${\bf \hat{B}}$ contains a divergent part ${\bf \hat{J}}$ which is diagonal,

$$\begin{eqnarray}&&{{\hat{J}}_{mm}}=\mathop{\sum }\limits_{n}{{\hat{B}}_{mn}}.\end{eqnarray} \tag{ 15 }$$

This decomposition is useful when ${\bf \hat{S}}$, ${\bf \hat{A}}$, and ${\bf \hat{J}}$, can be used to distinguish local physical processes in $\bar{\mathcal{L}}$. In the next section we use ${\bf \hat{S}}$ to estimate the SST diffusion, ${\bf \hat{A}}$ to estimate SST advection, and ${\bf \hat{J}}$ to estimate the rate at which SST relaxes to atmospheric conditions.

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. The continuous equation satisfies (3) where $c^{\prime} ({\bf x},t)$ is the SST fluctuation in space ${\bf x}=\left[ x,y \right]$ and time t. The modelʼs SST transport operator is

$$\begin{eqnarray}&&\bar{\mathcal{L}}=\nabla \cdot {\bf \bar{K}}({\bf x})\nabla -{\bf \bar{u}}({\bf x})\cdot \nabla -\bar{r}({\bf x})\end{eqnarray} \tag{ 16 }$$

which is the time-average of a velocity vector ${\bf u}=[u({\bf x},t);v({\bf x},t)]$, symmetric diffusivity tensor ${\bf K}=\left[ {{\kappa }^{xx}}({\bf x},t),{{\kappa }^{xy}}({\bf x},t);{{\kappa }^{xy}}({\bf x},t),{{\kappa }^{yy}}({\bf x},t) \right]$, and relaxation rate $r({\bf x},t)$. 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].

Zoom In Zoom Out Reset image size

The continuous equations are spatially discretized to satisfy (5). The rectangular grid (figure 1) has a total of $M=W\;\times \;H$ cells with W horizontal cells of size $\Delta x$ and H vertical cells of size $\Delta y$. The discretization of $\bar{\mathcal{L}}c^{\prime}$ uses a central differencing stencil,

$$\begin{eqnarray}\begin{array}{rcl} \bar{\kappa }_{m}^{xx}\frac{{{\partial }^{2}}}{\partial {{x}^{2}}}c_{m}^{\prime } & \approx & \frac{\bar{\kappa }_{m}^{xx}}{\Delta {{x}^{2}}}\left( c_{m+1}^{\prime }-2c_{m}^{\prime }+c_{m-1}^{\prime } \right) \\ \bar{\kappa }_{m}^{yy}\frac{{{\partial }^{2}}}{\partial {{y}^{2}}}c_{m}^{\prime } & \approx & \frac{\bar{\kappa }_{m}^{yy}}{\Delta {{y}^{2}}}\left( c_{m+W}^{\prime }-2c_{m}^{\prime }+c_{m-W}^{\prime } \right) \\ 2\bar{\kappa }_{m}^{xy}\frac{{{\partial }^{2}}}{\partial x\partial y}c_{m}^{\prime } & \approx & \frac{2\bar{\kappa }_{m}^{xy}}{\Delta x\Delta y}\left( c_{m+W+1}^{\prime }-2c_{m}^{\prime }+c_{m-W-1}^{\prime } \right) \\ \frac{\partial \bar{\kappa }_{m}^{xx}}{\partial x}\frac{\partial }{\partial x}c_{m}^{\prime } & \approx & \frac{\bar{\kappa }_{m+1}^{xx}-\bar{\kappa }_{m-1}^{xx}}{4\Delta {{x}^{2}}}\left( c_{m+1}^{\prime }-c_{m-1}^{\prime } \right) \\ \frac{\partial \bar{\kappa }_{m}^{yy}}{\partial y}\frac{\partial }{\partial y}c_{m}^{\prime } & \approx & \frac{\bar{\kappa }_{m+W}^{yy}-\bar{\kappa }_{m-W}^{yy}}{4\Delta {{y}^{2}}}\left( c_{m+W}^{\prime }-c_{m-W}^{\prime } \right) \\ \frac{\partial \bar{\kappa }_{m}^{xy}}{\partial x}\frac{\partial }{\partial y}c_{m}^{\prime } & \approx & \frac{\bar{\kappa }_{m+1}^{xy}-\bar{\kappa }_{m-1}^{xy}}{4\Delta x\Delta y}\left( c_{m+W}^{\prime }-c_{m-W}^{\prime } \right) \\ \frac{\partial \bar{\kappa }_{m}^{xy}}{\partial y}\frac{\partial }{\partial x}c_{m}^{\prime } & \approx & \frac{\bar{\kappa }_{m+W}^{xy}-\bar{\kappa }_{m-W}^{xy}}{4\Delta x\Delta y}\left( c_{m+1}^{\prime }-c_{m-1}^{\prime } \right) \\ \bar{u}_{m}^{x}\frac{\partial }{\partial x}c_{m}^{\prime } & \approx & \frac{\bar{u}_{m}^{x}}{2\Delta x}\left( c_{m+1}^{\prime }-c_{m-1}^{\prime } \right) \\ \bar{u}_{m}^{y}\frac{\partial }{\partial y}c_{m}^{\prime } & \approx & \frac{\bar{u}_{m}^{y}}{2\Delta y}\left( c_{m+W}^{\prime }-c_{m-W}^{\prime } \right). \\ \end{array}\end{eqnarray} \tag{ 17 }$$

Cell number m increases from the bottom-left to the top-right in figure 1 such that $m-1,m+1,m-W$ and $m+W$ are the left, right, bottom and top neighbors, respectively. The terms in (17) together with the relaxation rate ${{\bar{r}}_{m}}$ and boundary conditions comprise the discretized transport operator ${\bf B}$ in (5).

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 $\Delta 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 ${\bf a},{\bf b}$ and ${\bf c}$. The anomaly signal is obtained by removing the average 12-month seasonal cycle from the SST fluctuation time series.

The transport between locations ${\bf a},{\bf b}$ and ${\bf c}$ is quantified by elements ${{G}_{ab}}(t)$, ${{G}_{bc}}(t)$, and ${{G}_{ac}}(t)$ of the true response function ${\bf G}(t)$ as shown (by the solid lines) in figure 2. The ${\bf a}\to {\bf b}$ response arrives sooner and has a larger magnitude than the ${\bf a}\to {\bf c}$ response because ${\bf a}$ is nearer to ${\bf b}$ than ${\bf 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.

Zoom In Zoom Out Reset image size

We now estimate the modelʼs SST transport response function and flow parameters using the SST anomaly time series. We estimate the discretized operator ${\bf \hat{B}}$ using the local method (12) where ${\bf c}^{\prime} (t)$ is a finite length time series of monthly SST anomalies. The estimated response ${\bf \hat{G}}(t)$ is calculated by inserting ${\bf \hat{B}}$ into (10). Figure 2 shows elements of the estimated response ${{\hat{G}}_{ab}}(t)$, ${{\hat{G}}_{bc}}(t)$ and ${{\hat{G}}_{ac}}(t)$ using 50 and 100 years of SST data. The order of magnitude and timescale of transport are accurately captured by the 50 year estimate. With 100 years, the average error in the logarithm of ${\bf \hat{G}}(t)$ is approximately 10%.

Parameters of the flow are estimated by separating ${\bf \hat{B}}$ into symmetric, antisymmetric, and divergent operators (13)–(15). Estimates of the velocity vector, diffusivity tensor, and relaxation are obtained by

$$\begin{eqnarray*}\begin{array}{rcl} \widehat{{\bar{\kappa }}}_{m}^{xx} & = & \Delta {{x}^{2}}{{{\hat{S}}}_{m,m+1}} \\ \widehat{{\bar{\kappa }}}_{m}^{yy} & = & \Delta {{y}^{2}}{{{\hat{S}}}_{m,m+W}} \\ 2\widehat{{\bar{\kappa }}}_{m}^{xy} & = & \Delta x\Delta y{{{\hat{S}}}_{m,m+W+1}} \\ \frac{\widehat{\partial \bar{\kappa }}_{m}^{xx}}{\partial x} & = & \frac{1}{2\Delta x}\left( \widehat{{\bar{\kappa }}}_{m+1}^{xx}-\widehat{{\bar{\kappa }}}_{m-1}^{xx} \right) \\ \end{array}\end{eqnarray*}$$

$$\begin{eqnarray}\begin{array}{rcl} \frac{\widehat{\partial \bar{\kappa }}_{m}^{yy}}{\partial y} & = & \frac{1}{2\Delta y}\left( \widehat{{\bar{\kappa }}}_{m+W}^{yy}-\widehat{{\bar{\kappa }}}_{m-W}^{yy} \right) \\ \frac{\widehat{\partial \bar{\kappa }}_{m}^{xy}}{\partial x} & = & \frac{1}{2\Delta x}\left( \widehat{{\bar{\kappa }}}_{m+1}^{xy}-\widehat{{\bar{\kappa }}}_{m-1}^{xy} \right) \\ \frac{\widehat{\partial \bar{\kappa }}_{m}^{xy}}{\partial y} & = & \frac{1}{2\Delta y}\left( \widehat{{\bar{\kappa }}}_{m+W}^{xy}-\widehat{{\bar{\kappa }}}_{m-W}^{xy} \right) \\ \widehat{{\bar{u}}}_{m}^{x} & = & -2\Delta x{{{\hat{A}}}_{m,m+1}}+\frac{1}{2\Delta x}\frac{\widehat{\partial \bar{\kappa }}_{m}^{xx}}{\partial x}+\frac{1}{2\Delta y}\frac{\widehat{\partial \bar{\kappa }}_{m}^{xy}}{\partial y} \\ \widehat{{\bar{u}}}_{m}^{y} & = & -2\Delta y{{{\hat{A}}}_{m,m+W}}+\frac{1}{2\Delta y}\frac{\widehat{\partial \bar{\kappa }}_{m}^{yy}}{\partial y}+\frac{1}{2\Delta x}\frac{\widehat{\partial \bar{\kappa }}_{m}^{xy}}{\partial x} \\ {{\widehat{{\bar{r}}}}_{m}} & = & {{{\hat{J}}}_{mm}}. \\ \end{array}\end{eqnarray} \tag{ 18 }$$

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%.

Zoom In Zoom Out Reset image size

Figure 4 shows the convergence of ${\bf \hat{G}},\widehat{{\bar{u}}},{{\widehat{{\bar{\kappa }}}}^{xx}}$ and $\widehat{{\bar{r}}}$ as the length ℓ of the SST time series increases from 10 to 10,000 years. The error decreases as $1/\sqrt{\ell }$ 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.

Zoom In Zoom Out Reset image size

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 $\overline{{\bf c}^{\prime} {\bf c}{{^{\prime} }^{{\rm 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 $\overline{{\bf c}^{\prime} {\bf c}{{^{\prime} }^{{\rm 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.