Changyu LI , Jianping HUANG , Lei DING , Yu REN , Linli AN , Xiaoyue LIU , and Jiping HUANG

1College of Atmospheric Sciences,Lanzhou University,Lanzhou 730000,China

2Collaborative Innovation Center for Western Ecological Safety,Lanzhou University,Lanzhou 730000,China

3Enlightening Bioscience Research Center, Mississauga, L4X 2X7, Canada

Text S1.The ensemble empirical mode decomposition (EEMD) method

The ensemble empirical mode decomposition (EEMD) was developed based on the empirical mode decomposition(EMD). EMD is a one-dimensional data analysis method that is adaptive, has high locality, and can thereby handle the nonlinear and nonstationary nature of data. EEMD adds robustness to the EMD decomposition when data is perturbed by noise, guaranteeing that the physical interpretation of the decomposition result is not sensitive to the noise inevitably contained in real data. In this article, the steps of EEMD can be described as follows:

(1) Add a white noise series with an amplitude fraction (0.2 in this study) of that of the standard deviation of the raw data seriesx(t);

(2.1) Setx1(t)=x(t) and find the maxima and minima ofx1(t). Then obtain the upper envelopeeu(t) and the lower envelopeel(t) using cubic splines to connect the maxima and minima, respectively;

(2.2) Find the local meanm(t)=[eu(t)+el(t)]/2 and then determine whetherm(t) is close to zero at any location based on the given criterion;

(2.3) If yes, stop the sifting process; otherwise, setx1(t)=x(t)-m(t) and repeat steps 2.1 to 2.2;

(2.4) In this manner, we obtain the first intrinsic mode function (IMF), and by subtracting it fromx(t), we obtain a remainder series. If the remainder still contains oscillatory components, we again repeat steps 2.1 to 2.2, but with the newx1(t) as the remainder.

So, each time series is decomposed into different IMFs, which can be expressed as:

whereCj(t)represents the jth IMF, which is an amplitude-frequency-modulated oscillatory component, andRn(t) is the residual of datax(t), which is either monotonic or contains only one extreme.

(3) Repeat steps 1 and 2 again and again but with different white noise series added each time and obtain the (ensemble)means for corresponding IMFs of the decompositions as the final result.

The noise added to data has an amplitude that is 0.2 times the standard deviation of the raw data, and the ensemble number is set to 400.

Fig.S1. The time series of combustion ratios (O:C) from 1990 to 2014.

Fig.S2. The annual mean air-sea O2 flux (mol m-2 yr-1) averaged from 1985 to 2014 in each CMIP6 model.

Fig.S3. Map of the ocean regions. For comparison with other studies, the global ocean is divided into 13 regions. The boundaries of these regions are defined at 58°S,36°S, 13°S, 13°N and 36°N in the Pacific and at 53°N in the Atlantic.

Fig.S4. Standard deviation of the natural variability based on the EEMD method.The natural variability is the sum of IMFs 2-5 from EEMD.

Fig.S5. The spatial differences of mixed layer depth (m) between 1985-2014 and 1950-79 in CMIP6 historical simulations.