Roberto De Almeida: The SACZ/SST feedback

Another open notebook science post. Part I, part II.

The central idea of this work is that the interaction between sea surface temperatures (SST) and the South Atlantic Convergence Zone (SACZ) is characterized by a negative feedback. In this feedback, positive (warm) SSTs enhance the SACZ, while the associated increase in cloudiness acts to damp the initial SST anomaly. This relationship can be modeled as:

$x_{i + 1} = a_{1} x_{i} - {by}_{i} + γ$
$y_{i + 1} = {bx}_{i} + a_{2} y_{i} + η$

Where $x$ is the SST, $y$ is a variable related to the SACZ intensity (like precipitation); $a_{1}$ and $a_{2}$ correspond to the “memory” of each variable between timesteps, and $γ$ and $η$ represent the external forcing which excites the system. The equations are know as the stochastic oscillator (De Almeida et al. 2008).

It’s possible to quantify the negative SACZ/SST feedback by estimating the $b$ coupling term from data. If we plot the spatial distribution of $b$ in the South Atlantic ocean, a pattern of high values appears beneath the SACZ climatological position, with a northward shift due to migration of the SACZ during its lifetime.

In order to estimate $b$ we first calculate $a_{1}$ and $a_{2}$ from the autocorrelation of the variables. After that, we can simply estimate $b$ through a linear regression that minimizes the stochastic noise $γ + η$ :

$x_{i + 1} - a_{1} x_{i} + y_{i + 1} - a_{2} y_{i} = b (x_{i} - y_{i}) + γ + η$

The code is pretty simple:

mapb = zeros((ENSEMBLE, SPACE), 'f')
for member in range(ENSEMBLE):
    # remove seasonal cycle and normalize by interannual variability
    ppta = ppt[member] - mean(ppt[member], axis=0)
    ppta.shape = (TIME, Z, SPACE)
    ppta = ppta/std(ppta, axis=0)
    tsa = ts[member] - mean(ts[member], axis=0)
    tsa.shape = (TIME, Z, SPACE)
    tsa = tsa/std(tsa, axis=0)
 
    # calculate b at each point in space
    for i in range(SPACE):
        # slice for JF (2:4) and FM (3:5)
        s1 = (slice(None), slice(2,4), i)
        s2 = (slice(None), slice(3,5), i)
 
        # estimate memory a1 and a2
        a1 = corrcoef(tsa[s1].flat, tsa[s2].flat)[0,1]    # JFM-FMA
        a2 = corrcoef(ppta[s1].flat, ppta[s2].flat)[0,1]
 
        # linear regression minimizing error
        LHS = tsa[s2] + ppta[s2] - a1*tsa[s1] - a2*ppta[s1]
        LHS.shape = (-1, 1)
        RHS = tsa[s1] - ppta[s1]
        RHS.shape = (-1, 1)
        b, dd, rank, sv = linalg.lstsq(RHS, LHS)
        mapb[member,i] = b[0,0]
mapb.shape = (ENSEMBLE, LATITUDE, LONGITUDE)