The linear behavior of the joint initial-boundary-value predictability of the climate system

The primary sources of predictability of the Earth climate system come from initial conditions and external radiative forcings, which is a joint initial-boundary-value problem. It is crucial to clearly understand the contributions of initial conditions and external forcings (i.e. initial-value and boundary-value) to the total predictability, for which the traditional predictability study examines the growth of system uncertainties is however difficult to quantitatively identify. This study uses the variations of prediction skills in the space of prediction lead times to systematically examines the linear behaviors of initial-value, boundary-value and joint initial-boundary-value predictability, which are measured by anomaly correlation coefficients between system states in a series of perturbed experiments on initial and boundary conditions. The lead time scale at which the prediction skill induced from external forcing exceeds the skill from the initial condition is a key measurement of system uncertainty, after which the boundary-value predictability efficiently extends the total predictability. These results provide a guideline to understand how the accurate setting of external forcing caused by human activities promotes the prediction skill when a comprehensive Earth climate system model is initialized to make climate predictions.


Introduction
The dynamical and physical processes in the complex Earth system could transform the radiative effects of greenhouse gas emissions and natural aerosols (GHGNAs) into climate signals; for example, modulating the variability of the El Niño and Southern Oscillation (ENSO) (Collins et al. 2010), weakening the Atlantic Meridional Overturning Circulation (AMOC) (Parker et al. 2015), and even influencing the likelihood and intensity of extremes in the atmosphere and ocean (Goddard et al. 2015;Murakami et al. 2020) such as tropical cyclones, heatwaves, winter storms, droughts, floods, and coastal sea-level rise (Solomon et al. 2007). One critical aim of climate impact and vulnerability studies is to detect the influences of radiative boundary conditions (BCs) on an initialized climate signal, addressing the boundary-value (second kind) predictability that is mixed with the initialvalue (first kind) predictability (Lorenz 1975).
Determined by system dynamics, predictability of the climate system refers to the traceable nature of signals when the system evolves over time (Schneider and Griffies 1999), from which one can develop methods to predict the future state of certain signals. This requires an understanding of the associated dynamics and physics that occur during the period. The predictability sources of the Earth climate system reside in both initial and boundary conditions (Boer 1 3 2000) i.e., predicting the future climate state is a joint initial-boundary-value predictability problem (Lorenz 1975;Schneider and Griffies 1999) given the estimated initial and boundary conditions. Here, for the Earth climate system as a whole, "boundary" conditions refer to the effects of external radiative forcings (e.g. Lorenz 1975;Collins and Allen 2002) such as solar radiation and effects of GHGNAs.
Traditional predictability studies, such as weather forecasts and interannual climate forecasts, usually explore the law of error growth and spread, with a focus on initial-value predictability (Thompson 1957;Lorenz 1963;Leith 1973;Madden 1976;Barnett 1981;Hayashi 1983;Trenberth 1985;Li et al. 2018), while climatic phenomena with long timescale impacts are more associated with boundary-value predictability (Reichler and Roads 2003). Previous studies have explored the comprehension of climate system predictability by reducing uncertainties caused by natural fluctuations, models, and scenario changes (Hawkins and Sutton 2009). Although estimating the relative magnitude of external signals and internal noise gives an estimation of the influences of boundary conditions on the predictability of climate signals, how to clearly identify the propagation and transformation of boundary signals and thus understanding boundary-value predictability remains a topic that is worth to be further studied. The Earth climate system is a complex coupled system consisting of various components, such as ocean, sea ice, atmosphere and land processes, and containing motions with multiple time scales from hourly to centennial (Hurrell et al. 2013). To clearly understand the mechanisms of climatic impacts, we choose the slow-varying ocean as the research target. By using three climate models from conceptual models to complex coupled general circulation models (CGCMs), we examine the time scale of the ocean prediction skill and how it evolves in these models. Previous studies employed such as root mean square error (e.g. Lorenz 1965), Lyapunov index (e.g. Toth and Kalnay 1997;Ding et al. 2008) and attractor radius (e.g. Li et al. 2018) to study the predictability of a system, which includes both linear and nonlinear components. However, it is difficult to use such error growth predictability study approaches to quantitively identify contributions of initial-value and boundary-value for predictability of a climate signal. Based on the work of Reichler and Roads (2003), this study explores a systematical method that uses anomaly correlation coefficients (ACCs) (e.g. Wilks 1995) to measure prediction skills and examine the variations of prediction skills in the space of prediction lead times to understand the characteristics of the linear behavior of predictability. The variation of ACCs between system states started from (forced by) a set of control initial (boundary) conditions and a set of perturbed initial (boundary) conditions can be used to understand the characteristics of the linear behavior of initialvalue (boundary-value) predictability, while the combination of both can show the characteristics of the linear behavior of joint initial-boundary-value predictability. It is expected that initial-value predictability would exhibit a decrease in skill due to the growth of the initial error, while boundary-value predictability would show an increase in skill due to the boundary information carried into the internal system over evolution of the states (Reichler and Roads 2003). For the behavior of joint initial-boundary-value predictability, the skill decreases more slowly than initial-value alone predictability due to the contribution of boundary signals.
In this study, we try to answer two questions for advancing the understanding of the problem of joint initial-boundaryvalue predictability: (1) how initial-value predictability and boundary-value predictability superimpose together? (2) For a specific climatic aspect, what is the timescale at which the externally forced boundary signal begins to play a more important role than the decaying initial signal provided by initial conditions, and what is the timescale by which the boundary signal becomes saturated?
The introduction of three coupled models and experimental details are given in Sects. 2 and 3. The results are analyzed and discussed in Sect. 4. Finally, summaries and discussions are given in Sect. 5.

The 5-variable coupled model
Due to the complexity of the physical process and the huge amount of calculation, it is convenient to explore the various kinds of predictability of the climate system starting from a simple model. Using a conceptual model can easily study the nature of two kinds of predictability and lay a theoretical foundation for the predictability examination. Based on the butterfly model (Lorenz 1963) coupled with a slab ocean (Zhang 2011) and the three-term balance mechanism of mean meridional flow (Gnanadesikan 1999), Zhang (2011) developed this 5-variable conceptual pycnocline predictive model (5VCM), taking the form as.
The variables x 1 , x 2 , x 3 originated from the Lorenz63 model and , , b are the original atmospheric parameters that maintain the chaotic regime of the atmosphere, w and represent the upper ocean and deep ocean respectively. (1a) The oceanic and atmospheric components, the deep ocean and the upper ocean are coupled through constant parameters ( c 1 , c 2 , c 3 , c 4 , c 5, c 6 ) . O m and O d denote the heat capacity and damping rate, respectively. Γ represents time scale of the deep ocean changes. The external forcing term of this model is represented by S m + S s cos 2 t∕S pd ; S m and S s stand for annual mean and seasonal cycle of external forcing, and the period of the external forcing is represented by S pd . The meaning of all other parameters can be found in Zhang (Zhang 2011) where shows the default values of all the parameters ( , , b , c 1 , c 2 , O m , O d , S m , S s , S pd , Г, c 3 , c 4 , c 5 , c 6 ) = (9.95, 28, 8/3, 10 -1 , 1, 1, 10, 10, 1, 10, 100, 10 -2 , 10 -2 , 1, 10 -3 ). The 5VCM is starting from ( x 1 , x 2 , x 3 , w, ) = (0, 1, 0, 0, 0) and first spun up for 10 4 time units (TUs). Using a second-order Runge-Kutta differencing scheme with Δt = 0.01 nondimensional time unit (TU) to forward the simple 5-variable coupled model (5VCM), it can simulate the fundamental features of the real climate system. For example, the upper ocean interacts with the transient atmospheric attractor above and the deep ocean below (e.g. Zhang 2011), and the decadal and longer time scale variability of deep ocean has the impact on the transient atmospheric attractor through the feedbacks of the upper ocean to the atmosphere (e.g. Zhao et al. 2019).

The Fast Ocean-Atmosphere Model
To further verify the predictability of the model, we use a coarse-resolution coupled general circulation model (CGCM) called Fast Ocean-Atmosphere Model (FOAM) (Jacob et al. 2001). The FOAM consists of an atmosphere model with Gaussian grids of 4.5 • lat × 7.5 • lon and 18 vertical layers, and an ocean model with the horizontal resolution of 2.8 • latitude × 1.4 • longitude and 24 layers in z-coordinate. The external forcing of the coupled model simply consists of the solar radiation that imposes on the sea surface temperature (SST) as ΔSST solar = β × S 0 /(C p × ρ × h) where S 0 is the sunlight constant and C p is the specific heat of sea water at the constant pressure, ρ is the water density and h is the layer thickness whileβis the air-sea exchange coefficient. It also has a land surface component and a sea-ice component, as well as a river transport model. In FOAM, one model year equal 360 model day.
The FOAM was developed to study the long-term natural variability of the climate system. Compared with 5VCM, it has more realistic air-sea interactions and atmospheric and oceanic circulation structures. Compared to other more comprehensive CGCMs such as the CM2.1 that is used in this study, FOAM is a much coarse resolution model with simpler physics such as fixed vegetation etc. (Jacob et al. 2001). The advantage of this model is that it is very fast and convenient to perform large ensemble experiments to detect initial and boundary predictability as well as joint initial-boundary-value predictability.

The Coupled Model Version 2.1
The Coupled Model Version 2.1 (CM2.1) (Delworth et al. 2006) is a coupled climate model developed by Geophysical Fluid Dynamics Laboratory (GFDL) for seasonal and annual forecasts and multi-centuries of global climate change research. The coupled model includes four modules: ocean, atmosphere, land, and sea ice. The ocean model uses MOM4 (GFDL Modular Ocean Model Version 4), with a horizontal resolution of 1° × 1°, and gradually increases from 30° north to south latitude to the equator to 1/3°, vertically divided into 50 layers, of which the upper 220 m are equally divided into 22 layers. The atmosphere and land models use AM2.1 (GFDL Atmosphere Model Version 2.1) and LM2.1 (GFDL Land Model Version 2.1) respectively, with a horizontal resolution of 2°lat × 2.5°lon, and the atmospheric model has 24 vertical layers. The model has performed well in the previous assessments of the Intergovernmental Panel on Climate Change (IPCC) (Randall 2007) and has been applied to the report of the US Climate Change Science Program (US CCSP). Compared to the FOAM, the CM2.1 is a more comprehensive CGCM which describes the Earth climate system more realistically.

Experimental design
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: where F(x) is the reservoir of the time tendency of the system state, ∂x/∂t, while F B represents the processes associated with BCs (radiative external forcing in coupled systems). For the 5VCM, the state vector x = (x 1 , x 2 , x 3 , 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 S m + S s cos(2 t/S pd ) (Eq. 1d), which is referred to as F B . 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 F B refers to radiative effects of the coupled system as external forcings [solar radiation for the FOAM and Greenhouse gases mainly represented by carbon dioxide (CO 2 ) for the CM2.1]. Once the ICs and BCs are given for the climate (2) ∕ t = ( ) + F B , given | t=0 = I as ICs, and BCs, 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: is the covariance that measures the degree to which both variables covary, and E [⋅] 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 initialboundary-value predictability are described as follows.

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 I i , 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 = I i (I i p ), the "Control" ("Perturbed") system ensemble at time t is defined as , i = 1…N, where I i p represents the perturbed version of I i (red dotted curve in the upper-left panel of Fig. 1a).
Letting A', B' be the anomalous fields of x(I) (t) and x(I p ) (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, I i (i = 1…20) are taken from a 10 6 -TU long integration with an interval of 10 3 TUs. A Gaussian noise is added to the "Control" initial values of variable x 1 -x 3 to form the Perturbed ensemble ICs, I i p , 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, I i (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 I i p are the "same" restart files but only the atmospheric temperature are replaced by the next day counterparts. Then, setting A' = x(I) (t) and B' = x(I p ) (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.

For detecting boundary-value predictability
Different from detecting the behavior of 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).
We first copy two independent ICs (chosen from a long time model integration) into M copies to form the M-member ensemble of ICs denoted as I j (1) and I j (2) , j = 1…M, the superscript (1) and (2) denotes two independent initial states, as shown by the blue dotted and red dotted curves in the upper-left panel of Fig. 1b as I (1) , I (2) . 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. Notice that the BCs were kept constant throughout the experiment.
Letting A', B' be the anomalous fields of x(I (1) ) (t) and x(I (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 as 0 (upper left panel of Fig. 1b).

For detecting joint initial-boundary-value predictability
We also further investigate contributions of initial signals and boundary signals to the joint initial-boundary-value Fig. 1 Schematic illustration of initial-value (boundary-value) predictability determined by the initial (boundary) signals. a The anomaly correlation coefficient (ACC) between the ensembles of "perturbed" (red solid curve connected by the red dots, denoted by x(I i p ) (t) , i = 1…N) and "control" (blue solid curve connected by the blue dots, denoted by x(I i ) (t) , i = 1…N) decreases over the lead time (t). The ensemble members of perturbed simulations are initialized from the perturbed initial conditions (ICs) (red solid curve connected by the red dots at t = 0, denoted by I i p , i = 1…N) from a set of control ICs (blue solid curve connected by the blue dots at t = 0, denoted by I i , i = 1…N) that initialize the ensemble of control simulations. b The ACC between 2 sets of ensemble simulations forced (as shown by green arrows) by the same set of boundary conditions (BCs) but initialized from independent ICs increases over the lead time (t). Here, 2 sets of ensemble simulations are shown as the solid and dotted curves connected by blue and red dots, denoted by x((I j ) (1) ) (t) (BCs-forced simulations-1) and x((I j ) (2) ) (t) (BCs-forced simulations-2), j = 1…M, respectively, each of which is forced by the same set of BCs (dotted-green solid curves connected by the green dots, denoted by B j , j = 1…M) that sample signals from the bound box in the lower part of the upper panel but initialized from one of two independent initial states (denoted by the superscript as (1) or (2)). Here, N (M) is the ensemble size in detecting initial-value (boundary-value) predictability 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 Fig. 2a, b).
First, same as the experiments for detecting initial-value predictability, N independent ICs (I i , i = 1…N) are chosen from a long timeseries of system states, from which M sets of copies are made. Corresponding to each I i , 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 I i (j) (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 this M × N ICs, denoted as I i p(j) (j = 1,…M, i = 1,…N) (see left panel in Fig. 2b).
Letting A′, B′ be the anomalous fields of x(I i (j) ) (t) , x(I i p(j) ) (t) respectively, 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 (Fig. 2a, b). Then, 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.

Results
Focusing on ocean with high heat capacity and low dissipation, we are able to better investigate the impact of initial conditions and external radiative forcing on the linear behavior of predictability. We use a conceptual 5VCM to gain the first glimpse of the fundamental picture of the behavior of joint initial-boundary-value predictability of a climate system. Then based on two CGCMs, the FOAM and CM2.1, we systematically examine the initial-value predictability and boundary-value predictability, as well as the joint initialboundary-value predictability of the coupled climate system. To analyze the predictability of ocean phenomena at different time scales, a high-frequency term w (SST, sea surface temperature) and a low-frequency term (OHC u100 , the sum of 0-100 m ocean heat capacity) in the 5VCM (CGCMs) are chosen.

The behavior of initial-value predictability
The skill of initial-value predictability measures the decaying rate of the initial signals due to the growing of the system error (Lorenz 1975). The linear behavior of such predictability can be quantified by the change in ACCs between numerical simulations with perturbed ICs and control (unperturbed) ICs, as shown in Fig. 3a-c. The prediction skill induced from initial signals decreases due to the growth of the initial error over the lead time as the system subsequently evolves (Reichler and Roads 2003). Figure 3a shows that the characteristic scale of is ten times that of w in the 5VCM, which is a simple conceptual model and the sea-air coupling is not as strong as CGCM (Zhang 2011). It is worth to mention that compared with the other two models, the 5VCM shows a persistent behavior at the beginning period (0-5 Time Units in Fig. 3a) due to the weak sea-air coupling effect. For the CGCMs, the predictability limits of SST and OHC u100 determined by the initial-value problem are seasonal and interannual (Fig. 3b, c). The "high"-frequency variables in both the 5VCM (w) and CGCMs (SST) are more sensitive to the initial signals than low-frequency variables ( in the 5VCM and OHC u100 in the CGCMs) (here, the ACCs of SST and OHC u100 are their global mean). For the ACCs of SST and upper ocean heat content, the slow change is governed by mixing processes and circulation changes in the upper ocean.
In the CGCMs, the predictability exhibits clear regional characteristics (Barnett 1981;Hawkins and Sutton 2009). The left panels of Figs. 4 and 5 further show the spatial distribution of the skill measuring initial-value predictability of SSTs over time in the FOAM and CM2.1, respectively. Due to a large difference in the resolution and dynamical processes between these two CGCMs, the ACC spatial distributions of their SSTs are different. As the forecast time increases, the ACC reflecting initial-value predictability of the FOAM decreases more rapidly, except over the Southern Ocean. In comparison, CM2.1 has longer predictability over most of the global domain, and the predictability is higher in the Northern Hemisphere, especially in the Atlantic region. Therefore, the declining trend of ACCs in CM2.1 is slower than that in the FOAM (Fig. 3b, c). This may reflect the different mixing mechanisms in the two model systems.

The behavior of boundary-value predictability
The behavior of boundary-value predictability can be measured by the carry-in rate of boundary information as the forced signal into the dynamical system as it evolves over the lead time (Lorenz 1975) and therefore shows as the increase of ACCs (Reichler and Roads 2003). These ACCs can be evaluated by two sets of system states starting from , i = 1…N, j = 1…M. Right panel shows the evolution of the joint 2-D state space at t = T t , similarly in the right panel b. b Based on panels a) (at t = 0) but all states are perturbed as I i p(j) (i = 1…N, j = 1…M). Then, the ACC between [(x(I i (j) ) (t) , j = 1…M), i = 1…N] and [(x(I i p(j) ) (t) , j = 1…M), i = 1…N] is a measurement of joint initial-boundary-value predictability. c Illustration of three kinds of predictability. The dotted-red, dotted-green and purple line respectively represents the linear behavior of initial-value predictability, boundary-value predictability and joint initial-boundary-value predictability independent ICs but with the same set of BCs, as shown in Fig. 3d-f. We can see that the skill from BCs increases over the lead time to a modest but stable value. To ensure that the change of ACCs comes from the contribution of boundary signals, we turn off the boundary signal by changing BCs back to the uniform default value after a period of lead time; for example, 80 TUs for w, 800 TUs for (here, a TU is a nondimensional time unit in the 5VCM, the time scale of completing an attractive lob), 5 model years for SST, and 40 model years for OHC u100 in the FOAM. To save computing resources, we turn off the variation signal of BC 5 model years for CM2.1. Then, we monitor the change in the ACCs.
The models take a long time to capture and store the boundary signals. The low-frequency terms in the 5VCM and slow-varying components in the CGCMs have a high heat capacity and low energy dissipation and thus can store information for a long time scale. Therefore, the response of the boundary signals in the slow-varying component is significantly stronger than in fast-varying (high-frequency) variables. Given a specific BC, the upper limit of the Fig. 3 The initial-value (boundary-value) predictability behaving as a decreasing (increasing) ACC curve in the space of forecast lead time, with the variations in the ACCs of a, d w and in the 5VCM, sea surface temperature (SST) and the sum of 0-100 m ocean heat capacity (OHC u100 ) in b, e FOAM and c, f CM2.1 with forecast lead times. Left panels a-c are for the initial-value predictability and right panels d-f are for the boundary-value predictability. The shaded area denotes the uncertainties bounded by a positive/negative standard deviation. In the 5VCM, the standard deviation is evaluated by 20 sets of ensembles differing in initial conditions, while in the CGCMs, it is evaluated based on the variations in ACC values at different latitudes boundary-value predictability can be determined, i.e., the upper limit measured by the ACC of w is 0.22, and is 0.58 (Fig. 3d). The time required to reach this upper limit is also related to the frequency of the variable, low-frequency terms taking longer, i.e., 65 TUs for w and 642 TUs for . In addition, we can see that once the periodic boundary signal is turned off, the skill induced from the boundary signal generally decreases. However, it is clear that the skill decreases much faster in the 5VCM and FOAM (Fig. 3d, e) than in the CM2.1 (Fig. 3f). It is understandable that such a predictability behavior has model dependence (Hawkins and Sutton 2009). For the comprehensive model CM2.1 with complex dynamics and physics, the periodic boundary signal that has been sufficiently mixed with the whole system appears gradually decaying over the subsequent integration even after it is turned off, resulting in a slow pronounced downward trend of skill (5-10 model years in Fig. 3f). For the low-order models such as the 5VCM and FOAM, it quickly disperses. The effect of boundary signal in comprehensive models is worth to be further studied in the future.
External radiative forcings, such as solar radiation and GHG radiative effects that affect the system heat budget, are gradually carried into the system through model dynamics and physics, such as atmospheric dynamical processes and air-sea interactions etc. From the spatial distribution of ACCs measuring boundary-value predictability (middle Fig. 4 The spatial distribution of the ACCs reflecting initial-value (left, denoted as INI), boundary-value (middle, denoted as BOU) and joint initial-boundary-value (right, denoted as JOI) predictability of SST in the FOAM. The red (left), green (middle) and purple (right) boxes represent the initial-value predictability, boundary-value predictability and joint predictability, respectively. Panels a-e represent the first month, sixth month, second year, fifth year and tenth year respectively. Note that after the fifth year, the boundary signal is turned off panels of Figs. 4 and 5), the responses of both CGCMs to the BCs generally start from the tropical area and western boundary area, but pronounced differences exist in some detailed characteristics, which may be associated with model resolution and corresponding dynamics.
For the FOAM with a horizontal oceanic resolution of 2.8 • latitude × 1.4 • longitude, the storage of this boundary information (variational solar radiation) in tropical and subtropical areas is prominent. For the CM2.1 with a horizontal oceanic resolution of 1 • latitude × 1 • longitude, the GHG BC mainly affects the global radiation budget, thus leading to an increase of the ACC value in the areas with strong air-sea interactions, i.e., the tropical Indian and Pacific Oceans, as well as the North Atlantic (especially in the active AMOC area of CM2.1). Panels e in the middle box of Figs. 4 and 5 are the results of the fifth year after the variation signal of BCs is shut off, showing that the boundary signals in the FOAM have been dissipated almost completely, and in CM2.1, the high ACC value reflecting strong boundary signal is still retained, again, due to the model complexity. We also noticed that some areas in CGCMs have a strong negative ACC, especially in CM2.1, indicating that over these areas, the response signal of the comprehensive model to boundary information is out of the phase of the boundary forcings. The corresponding response mechanism should be further examined in follow-up studies.

The behavior of joint initial boundary value predictability
The ACCs between two sets of system states that are initialized from the control and perturbed IC sets, but forced by the same set of variational BCs (see Fig. 2), can be used to show the linear behavior of joint initial-boundary-value predictability, as in Fig. 6. We can see that such ACC values also decrease as the system evolves over the lead time, but different from the initial-value predictability case, the initial and boundary signals play different roles at different lead times. The linear behavior of joint predictability has similar characteristics in different models. Thus, we use the ACC of w of the 5VCM (as shown in Fig. 6a) to look at its general behavior. Since the contribution of the boundary signals takes a longer time, the skill in the early stage mainly reflects the effect of the initial signal, which is similar to the behavior of initial-value predictability in the lead time period of 0-8 TUs. When the model captures and stores the boundary signal more, the skill measuring the joint predictability is affected by the combined effects of the decreased initial signals and increased boundary signals, manifested as a gradual increase in ACCs in the follow-up 10-40 TUs. The contribution of the initial signal decreases over the lead time, while in the period of 40-160 TUs, as the storage of the Fig. 6 Joint initial-boundary-value predictability (purple) combining the processes of initial signal decay (red) and boundary signal carryin (green), thus behaving as decreasing first, then sustaining and finally decreasing to a trivial amount as the boundary signal vanishes, with the variations of the ACCs of a, d w and in the 5VCM, SST and OHC u100 in b, c, e, f) CGCMs with forecast lead times. In each panel, the ACC between any two subsets is plotted in dotted lines as a reference, which evaluates the initial-value (boundary-value) predictability, as shown in Fig. 3. Similar to Fig. 3, the shaded area denotes the uncertainties represented by a positive/negative standard deviation as the upper/lower bound boundary signal increases, the behavior of the joint predictability is close to that of the boundary-value predictability. The spatial distribution of the ACCs reflecting joint predictability (the right panels of Figs. 4 and 5) is used to compare with the initial-value and boundary-value predictability (the left and middle panels of Figs. 4 and 5). In general, due to the different properties of the two CGCMs and their different sensitivity to the given BCs, the spatial distribution of joint predictability is more influenced by the BCs in the FOAM but more influenced by ICs in the CM2.1.
In the first month, the spatial distribution of joint predictability is the same as that of initial-value predictability. With the carry-in of the boundary signal, the joint predictability is affected more by the boundary information, and then it appears to be the superposition of the initial-value and boundary-value predictability as the lead time increases. In the period of a few years of lead times, although the joint predictability reflects the combination effects of the increase in boundary signal and the decrease in initial signal, the joint predictability in most of the regions is consistent with boundary-value predictability, especially in the tropical and subtropical areas in the FOAM and the areas with strong airsea interactions in CM2.1, for example, the Atlantic basin and tropical Pacific and Indian Oceans.
By analyzing the results from the compound experiments above, we attempt to answer the questions raised in the introduction. As the boundary signal begins to exert a certain influence, the changing trend of the ACCs reflecting initial-value predictability and joint initial-boundary-value predictability appears different. When the carry-in of the boundary signal gets saturated, the joint initial-boundaryvalue predictability is dominated by boundary-value predictability. Our results show that the prediction of low-frequency phenomena could be significantly enhanced by considering the influence of the boundary signal from external radiative forcings, helping predict long-term scale phenomena.

Conclusions and discussions
With appropriate experimental design and a series of coupled models, we have examined the linear behavior of joint initial-boundary-value predictability of the coupled climate system, measuring the evolution of anomaly correlation coefficients (ACCs) due to the contributions of initial conditions and external radiative forcings. The decreasing ACC at a certain lead time from a perfect value of 1 at the initial time, is a measure of the growth rate of the initial error, characterizing the linear behavior of initial-value predictability. The higher (lower) the frequency is, the faster (slower) the variable loses its initial-value predictability due to the error growth mechanism of the dynamical system (Lorenz 1975;Zhang 2011). The increasing ACC at a certain lead time from a very trivial value at the initial time is a measure of the carry-in rate of the boundary signal, characterizing the linear behavior of boundary-value predictability. The change of ACCs with the lead time shows an upward trend first due to the carry-in of boundary signals (Reichler and Roads 2003) and then appears flattened as the stored boundary signal gets saturated. The lower the frequency is, the stronger the signal storage, and the longer it takes to get saturated.
In the real world, the predictability sources come from both the initial signals and boundary signals jointly. We have shown that the linear behavior of joint predictability as a combination of linear behaviors of initial-value predictability and boundary-value predictability. Under the circumstance, the decrease of prediction skills due to the growth of system errors can be compensated by the carriedin boundary signal after some time scale by some degree, hence extending the total predictability as the joint value of initial conditions and boundary conditions. This study has certain limitations. First, the present study only focuses on the linear behavior of predictability of the climate system on sea surface temperature and ocean heat content. More characteristic physical variables, such as thermocline and pycnocline depths, subtropical and subpolar gyres, and the Northern Hemisphere subtropical high and North Atlantic oscillation indices etc., also need to be examined. Second, the nonlinear behavior of the joint initialboundary-value predictability of the climate system needs to be addressed. In this regard, exploring the causes of negative ACCs in some areas as the responses of boundary signals is a starting point. Based on comprehensive CGCMs, such as CM2.1 and the community Earth system model, examining the change of strength and frequency on large-scale climate phenomena, such as the El Niño and Southern Oscillation, Atlantic Meridional Overturning Circulation and Antarctic Circumpolar Current as the responses of changes of initial and boundary conditions, can deepen our understanding of the behavior of joint initial-boundary-value predictability.
Finally, it should be noted that the external forcing induced predictability depends on knowledge of the future change of greenhouse gas and aerosols. It also has uncertainties and may involve additional components of the Earth system, such as those involving the carbon cycle, as well as the future behavior of our society. So, this study only serves as a first step, and more further studies shall be followed up.

Conflict of interest
The authors declare that they have no conflict of interest.