We formulate the joint initial-boundary-value predictability problem as a forwarding equation of vectorized state **x** with given initial conditions (ICs) and boundary conditions (BCs) as:

∂**x**/∂*t* = **F**(**x**) + FB, given **x**|t=0 = I as ICs, and BCs, (2)

where **F**(**x**) is the reservoir of the time tendency of the system state, ∂**x**/∂*t*, while FB represents the processes associated with BCs (radiative external forcing in coupled systems). For the 5VCM, the state vector **x =** (x1, x2, x3, w,η), while for the CGCMs, the FOAM and CM2.1 the state vector consists of atmospheric and oceanic state variables, such as temperature, wind, currents, moisture, and salinity. In the 5VCM, **F**(**x**) represents all terms except for the external forcing Sm+Ss cos(2πt/Spd) **(Eq. 1d**), which is referred to as FB. For the CGCMs, **F**(**x**) consists of dynamical terms, such as advection, pressure gradient forces and inertial force, as well as momentum, energy and mass budgets, caused by physical processes, such as cumulus convection and radiation, as well as interacting processes, such as air-sea interactions, while FB refers to radiative effects of the coupled system as external forcings [solar radiation for the FOAM and Greenhouse gases mainly represented by carbon dioxide (CO2) for the CM2.1]. Once the ICs and BCs are given for the climate system expressed by **Eq. 2**, the joint initial-boundary-value predictability measures the effect of both ICs and BCs.

Based on appropriately designed ensemble experiments, we can simply use the change in ACCs in the space of lead times to examine the linear behavior of predictability of a climate signal. ACC refers to the correlation coefficient between two sets of data after removing climatology (e.g. Wilks 1995). Equation 3 shows the formulation of the correlation coefficient (r) between two variables, *A**’* and *B**’*, as:

$$r\left({A}^{\text{'}},{B}^{\text{'}}\right)= \frac{Cov\left({A}^{\text{'}},{B}^{\text{'}}\right)}{\sqrt{Var\left[{A}^{\text{'}}\right]Var\left[{B}^{\text{'}}\right]}}$$

3

where\(Cov\left(A\text{'},B\text{'}\right)=E\left[\right({A}^{\text{'}}-E\left[{A}^{\text{'}}\right]\left)\right(B\text{'}-E\left[B\text{'}\right]\left)\right]\) is the covariance that measures the degree to which both variables covary, and \(E\left[·\right]\) and \(Var[·]\) denote operations of the expectation (calculated by the arithmetic mean in this case) and variance, respectively.

The detailed strategies of the experimental design for detecting initial-value, boundary-value and joint initial-boundary-value predictability are described as follows.

## 3.1 For detecting initial-value predictability

The principle of experiments for detecting the behavior of initial-value predictability is to measure the system damping rate of initial signals (illustrated in the lower panel of Fig. 1a) due to the existence of initial error and its growth mechanism in the climate system (see the upper three panels of Fig. 1a). For doing that, we first use N independent ICs (chosen from a long time model integration, for instance) to form the N-member ensemble of ICs denoted as Ii, i = 1,…,N where N is the ensemble size (blue dotted curve in the upper-left panel of Fig. 1a). Starting from **x**|t=0 **=** Ii (Iip), the “Control” (“Perturbed”) system ensemble at time t is defined as **x**(Ii)(t) [**x**(Iip)(t)], i=1…N, (red dotted curve in the upper-left panel of Fig. 1a) where Iip represents the perturbed version of Ii. Letting *A’*, *B’* be the anomalous fields of **x**(I)(t) and **x**(Ip)(t) respectively, and substituting into Eq. 3, we can calculate the ACC between them in terms of their variations in ensemble space as a function of time t (upper panels of Fig. 1a). Then, the variation of ACC in the space of lead times measures the change of damping rate of initial signals by the lead time, and thus we can use it to comprehend the behavior of initial-value predictability.

For the 5VCM, 20 Control ICs, Ii (i = 1…20) are taken from a 106-TU long integration with an interval of 103 TUs. A Gaussian noise is added to the “Control” initial values of variable *x*1-*x*3 to form the Perturbed ensemble ICs, Iip, with a standard deviation as 5% of the Control value. For this simple model, we test the sensitivity of results to the ensemble size using different ensemble sizes as 10, 20, 100 and 1000, and it shows that no significant difference is found in predictability detection when the ensemble size is no less than 20. For the CGCMs, Ii (i = 1,…,20) are taken from the restart states of a 100-yr pre-industrial control simulation with 1850 fixed-year forcings of greenhouse gas and natural aerosol (GHGNA) every 5 year, while Iip are the “same” restart files but all the atmospheric variables are replaced by the next day counterparts. Then, setting *A’*= **x**(I)(t) and *B’=* **x**(Ip)(t). For 5VCM (CGCMs), the state variable **x** is selected as *w**a* or η*a* (SST*a* or OHC*a*) (similarly hereinafter), where the subscript “*a*” represents the anomaly, and substituting into Eq. 3 the variation of the ACCs in the space of lead times is computed.

## 3.2 For detecting boundary-value predictability

Different from detecting initial-value predictability that measures the system damping rate of initial signals, the experiments for the behavior of boundary-value predictability shall detect the growth rate of system signals (as illustrated in the lower panel of Fig. 1b) transmitted from boundary signals by the dynamical and physical mechanisms in the system (see upper three panels of Fig. 1b). As illustrated in the middle panels of Fig. 1b, the M BCs (denoted as **B**j, j=1,…M) that sample a boundary signal (a sine curve for instance in this case) force two system ensembles as **x**(**I**j)(1)(t) [**x**(**I**j)(2)(t)], j=1…M, where M is the ensemble size, and the superscript (1) and (2) denotes two independent initial states, as shown by the blue dotted and red dotted curves in the upper panels od Fig. 1b as I(1), I(2).Notice that the BCs were kept constant throughout the experiment.

For the 5VCM, Bj is set as Sm[1 + αsin(2πj/20)] (see **Eq. 1d**), j=1…20, α = 0.15 (0.0 is the value after canceling the boundary signal), and two independent ICs, I(1) and I(2), being any two model states in the 106-TU integration with the 103-TU interval as stated before. For the CGCMs, I(1) and I(2) are any two restart states in the 100-yr pre-industrial control simulation with 5-year interval described before, and Bj is set as a sine variation of solar radiation for the FOAM (greenhouse gases (GHG) mainly represented by carbon dioxide (co2) for the CM2.1) as Bj = ΔSSTsolar[1 + αsin(2πj/20)] (Bj = ΔGHG[1 + αsin(2πj/20)]), j=1,…20, α = 0.15 (0.0). Letting *A’*, *B’* be the anomalous fields of **x(B**(1))(t) and **x(B**(2))(t) respectively and substituting into Eq. 3, we can calculate the ACC between them as a function of time t (upper panels of Fig. 1b), in terms of variations in the ensemble space of sampling boundary signal. Then, the variation of ACC in the space of lead times can measure the growth rate of system signals transmitted from the boundary signal, and thus we can use it to comprehend the behavior of boundary-value predictability. Here, we set the ACC of the two independent initial conditions copies at the initial moment to 0 (upper left panel of Fig. 1b).

## 3.3 For detecting joint initial-boundary-value predictability

We also further investigate contributions of initial signals and boundary signals to the joint initial-boundary-value predictability by combining the aforementioned two experimental design techniques together, applying a small perturbation to the initial states with the given periodic BCs. To incorporate the experiments described before for detecting the initial-value predictability and boundary-value predictability together (see Fig. 1), a joint 2-dimensional “space” of the system states is constructed, each dimension of which corresponds to the variation of states induced by different initial conditions or different boundary conditions (see left panels of **Figs. 2**ab). First, same as the experiments for detecting initial-value predictability, N independent ICs (Ii, i = 1…N) are chosen from a long timeseries of system states, from which M sets of copies are made. Corresponding to each Ii, the M system states are forwarded using the corresponding BCs that sample a boundary signal (a sine curve in this case) to form an M×N-member ensemble of forwarded system states, denoted as (**x**(Ii(j))(t)., j=1,…M), i=1,…N (see left panel in Fig. 2a). Meanwhile, the other M×N-member ensemble of system states is produced following the same fashion but using the perturbed version of N ICs (Iip(j), i = 1,…N), denoted as (**x**(Iip(j))(t), j=1,…M), i=1,…N (see left panel in Fig. 2b).

Same experimental design as the previous two predictability experiments, for the 5VCM (CGCMs), we set **B**j = Sm[1 + αsin(2πj/20)], j=1…20 [ΔSSTsolar (1 + αsin(2πj/20)) for the FOAM and ΔGHG[1 + αsin(2πj/20) for the CM2.1, j=1…20] (See **3.1** for details of the initial perturbation). Then, choosing a state variable **x** and setting *A’* = **x***a*(Ii(j))(t), *B’* = **x***a*(Iip(j))(t) i = 1…20,j=1…20 where the subscript “*a*” represents the corresponding anomaly, and substituting into Eq. 3, we can calculate the ACC between them as a function of time t, in terms of variations in the M×N ensemble space (**Figs. 2**ab).The variation of ACC in the space of lead times can measure both the initial signal damping rate and the boundary signal transmission rate (Fig. 2c), and thus we can use it to comprehend the behavior of joint initial-boundary-value predictability.