Development and evaluation of a new 4DEnVar-based weakly coupled ocean data assimilation system in E3SMv2

The development, implementation, and evaluation of a new weakly coupled ocean data assimilation (WCODA) system for the fully coupled Energy Exascale Earth System Model version 2 (E3SMv2) utilizing the four-dimensional ensemble variational (4DEnVar) method are presented in this study. The 4DEnVar method, based on the dimension-reduced projection four-dimensional variational (DRP-4DVar) approach, replaces the adjoint model with the ensemble technique, thereby reducing computational demands. Monthly mean ocean temperature and salinity data from the EN4.2.1 reanalysis are integrated into the ocean component of E3SMv2 from 1950 to 2021 with the goal of providing realistic initial conditions for decadal predictions and predictability studies. The performance of the WCODA system is assessed using various metrics, including the reduction rate of the cost function, root mean square error (RMSE) differences, correlation differences, and model biases. Results indicate that the WCODA system effectively assimilates the reanalysis data into the climate model, consistently achieving negative reduction rates of the cost function and notable improvements in RMSE and correlation across various ocean layers and regions. Significant enhancements are observed in the upper ocean layers across the majority of global ocean regions, particularly in the north Atlantic, north Pacific, and Indian Ocean. Model biases in sea surface temperature and salinity are also substantially reduced. For sea surface temperature, cold biases in the north Pacific and north Atlantic are diminished by about 1–2 °C, and warm biases in the Southern Ocean are corrected by approximately 1.5–2.5 °C. In terms of salinity, improvements are observed with bias reductions of about 0.5–1 psu in the north Atlantic and north Pacific and up to 1.5 psu in parts of the Southern Ocean. The ultimate goal of the WCODA system is to advance the predictive capabilities of E3SM for subseasonal to decadal climate predictions, thereby supporting research on strategic energy-sector policies and planning.

1 Introduction

Climate predictions are essential for understanding and mitigating the impacts of climate variability and change. The accuracy and reliability of climate predictions depends strongly on the initialization of the climate models, which requires realistic and high-quality initial conditions (ICs) for skillful predictions (Dirmeyer et al., 2018). Data assimilation (DA) techniques are important for providing realistic ICs by integrating observational data into the model, thereby enhancing the predictive capabilities of climate models (Tardif et al., 2014). The efficacy of DA techniques has been demonstrated through enhanced predictability on subseasonal to decadal timescales (Zhou et al., 2024).

Numerous studies have focused on the initialization of climate models for decadal predictions (Branstator and Teng, 2012; Polkova et al., 2019). Climate models integrate multiple components, including the atmosphere, ocean, sea ice, and land. For the initialization of climate models in decadal predictions, DA methods can be categorized into uncoupled data assimilation and coupled data assimilation (CDA). In the uncoupled method, DA is performed independently within the uncoupled atmosphere, land, and ocean models rather than in a coupled model. The optimal analyses from these uncoupled models are then integrated together to establish the ICs for the climate model's predictions (Yao et al., 2021). For example, some studies directly utilize existing reanalysis data to initialize climate models for decadal predictions (Yeager et al., 2012; Tian et al., 2021). Nevertheless, the uncoupled DA method may lead to imbalances between different model components, potentially inducing initial shocks and diminishing the reliability of climate predictions (Smith et al., 2015; Zhang et al., 2020). Therefore, there is a growing interest in exploring and developing CDA methods to enhance the coherence and accuracy of the ICs for climate predictions.

Many research groups and institutions are actively engaged in the development and refinement of CDA methods. In CDA, the assimilation process is conducted directly within a coupled model. Compared to uncoupled DA, CDA provides balanced ICs that are more coordinated across multiple components of coupled models (Zhang et al., 2014). Previous studies have demonstrated that CDA enhances interannual climate predictions more effectively than uncoupled DA (Zhang et al., 2005; Shi et al., 2022). CDA techniques are divided into weakly coupled data assimilation (WCDA) and strongly coupled data assimilation (SCDA). In the WCDA system, reanalysis data are assimilated independently within each component of the coupled model. However, through the coupled model integration, reanalysis information from one component is transmitted to other components through interactions across multiple systems (Browne et al., 2019; He et al., 2020a). Sequential DA is distinctly partitioned into two primary stages: the analysis and forecast steps. During the WCDA analysis step, reanalysis information from one component cannot directly influence other components due to the lack of cross-component background error covariances. Nonetheless, the coupled model is employed during the forecast step to transfer reanalysis information from a single component to others through the integration of the coupled system (Laloyaux et al., 2016; Carrassi et al., 2018). The primary distinction between WCDA and uncoupled DA is the use of the coupled model during the forecast step (Zhang et al., 2020). Recent studies have developed WCDA systems that separately assimilate reanalysis data from the atmosphere (Li et al., 2021), land (Shi et al., 2024), and ocean (He et al., 2017) into coupled models. On the other hand, SCDA employs cross-component background error covariances during the analysis step to directly exert an instantaneous impact of reanalysis information from a single component on the state variables of other components, treating all Earth system components as an integrated whole (Sluka et al., 2016). Moreover, SCDA also allows the reanalysis information from a single component to propagate to other components during the forecast step through the coupled model integration (Yoshida and Kalnay, 2018). Therefore, SCDA offers potential benefits, including reduced model drift and enhanced forecast accuracy (Smith et al., 2015). Nevertheless, the development of SCDA presents considerable obstacles, primarily due to the complexity of accurately establishing cross-component background error covariances (Penny and Hamill, 2017). As a result, most existing CDA systems continue to employ the WCDA systems.

This study presents the development and implementation of the weakly coupled ocean data assimilation (WCODA) system for the fully coupled Energy Exascale Earth System Model version 2 (E3SMv2), utilizing the four-dimensional ensemble variational (4DEnVar) method. The 4DEnVar method is based on the dimension-reduced projection four-dimensional variational (DRP-4DVar) approach, notable for its innovative application of 4DVar by replacing the adjoint model with the ensemble approach (Wang et al., 2010). Previous studies have shown that 4DVar-based methods outperform simpler schemes (e.g., nudging or 3DVar) by maintaining dynamical consistency with the model and minimizing initial shocks in the forecasts (Sugiura et al., 2008; Zhang et al., 2020). In the WCODA system, monthly mean ocean temperature and salinity data from the EN4.2.1 reanalysis are incorporated into the ocean component of E3SMv2 to provide realistic ICs for decadal predictions. Although the assimilation process during the analysis step is conducted independently within the ocean component, the fully coupled E3SMv2 model is employed during the forecast step to transmit reanalysis information from the ocean to other components (e.g., atmosphere and land) through multi-component interactions. Consequently, the reanalysis information assimilated into the ocean ICs affects other model components through the integration of the fully coupled model, emphasizing the operation of this system as a WCDA system. The primary objective of this WCODA system is to advance our understanding of the ocean's role in climate predictability. Shi et al. (2024) implemented a weakly coupled land data assimilation in E3SMv2 for isolating the land's role in climate predictability. By improving the accuracy of ICs for both land and ocean, we aim to advance the predictive capabilities of E3SM for decadal predictions, ultimately supporting research on energy-sector policy and planning.

This study presents and evaluates the 4DEnVar-based WCODA system for E3SMv2. Section 2 provides a detailed description of the E3SMv2 model, the ocean reanalysis data, and the framework for implementing the 4DEnVar-based WCODA system. Section 3 evaluates the assimilation performance of the WCODA system. Finally, Sect. 4 provides the conclusions.

2 Methodology

2.1 E3SM overview

Developed by the US Department of Energy, the Energy Exascale Earth System Model version 2 (E3SMv2) is a state-of-the-art climate model to advance our understanding of climate variability and its future changes (Leung et al., 2020). E3SMv2 integrates multiple components to simulate the complex interactions within the climate system, encompassing the atmospheric, sea ice, ocean, land, and river transport components. The atmospheric component (EAMv2) represents turbulence, clouds, and aerosol processes (Zhang et al., 2023) and features a nonhydrostatic dynamical core (Taylor et al., 2020). It operates on a dynamic grid with a horizontal resolution of approximately 110 km and includes 72 vertical layers that extend to the stratosphere (Golaz et al., 2022). The sea ice component (MPAS-SI) simulates the formation, evolution, and melting of sea ice, with detailed thermodynamics and dynamics processes (Turner et al., 2022). The ocean component (MPAS-O) is responsible for modeling the physical state and biogeochemical processes of the ocean, including detailed simulations of ocean currents, temperature, and salinity (Reckinger et al., 2015). MPAS-O operates at a horizontal resolution of ∼ 60 km in the mid-latitudes and ∼ 30 km at the Equator and poles, differing from the atmospheric model's resolution of 110 km. It is configured with 60 vertical layers, with finer resolution (∼ 10 m) near the surface and coarser resolution (∼ 200 m) at depth. The vertical mixing scheme employed is the K-profile parameterization, as described by Van Roekel et al. (2018). The land component (ELMv2) encompasses various land surface processes, including biophysical processes, soil processes, and surface hydrology (Golaz et al., 2019). These simulations are crucial for understanding land–atmosphere interactions and their impact on climate variability. Additionally, the river transport component (MOSARTv2) simulates the hydrological dynamics of water flow through river basins, providing insights into freshwater resources, flood risks, and sediment transport (Li et al., 2013). The CPL7 coupler dynamically integrates all five components by regulating the exchange of energy, water, and momentum fluxes between different components (Craig et al., 2012). A comprehensive evaluation of the E3SMv2 model is presented by Golaz et al. (2022).

2.2 Ocean reanalysis dataset

The ocean temperature and salinity data in this study are derived from the EN4.2.1 ocean reanalysis dataset. Produced by the Met Office Hadley Centre, the EN4.2.1 product is developed based on quality-controlled ocean temperature and salinity profiles from four input sources: Argo, ASBO (Arctic Synoptic Basin-Wide Oceanography), GTSPP (Global Temperature and Salinity Profile Programme), and WOD09 (World Ocean Database) (Good et al., 2013). The EN4.2.1 dataset includes observations from a wide range of profiling instruments, such as Argo floats, expendable bathythermographs (XBTs), and mechanical bathythermographs (MBTs) (Chen et al., 2020). According to Good et al. (2013), observations in EN4.2.1 are most abundant in the upper 100 m, with vertical resolution refined to ∼ 1 m in the top 100 m. Spatially, data density is high in regions such as the north Atlantic and western Pacific but decreases significantly in high-latitude and deep-ocean regions. This distribution in data availability influences the assimilation results. Areas with denser observational coverage, such as the upper north Atlantic, are expected to show greater improvements through assimilation, while regions with sparse observations may exhibit limited improvements.

To initialize decadal climate predictions, monthly mean ocean temperature and salinity data from the EN4.2.1 reanalysis are assimilated into the fully coupled E3SMv2 model across all 60 ocean layers from 1950 to 2021. The choice to utilize monthly mean reanalysis data is based on two primary reasons: firstly, data with higher temporal resolution (less than 1 month) might produce unwanted noise, potentially compromising the accuracy of decadal predictions. Secondly, the initialization for decadal predictions requires assimilation cycles spanning several decades, and assimilating complex, real-time observations over such extended periods would be computationally prohibitive. Therefore, in line with most existing studies that use reanalysis data for initializing decadal predictions (Pohlmann et al., 2019; Tian et al., 2021), this study assimilates the monthly mean EN4.2.1 reanalysis through the WCODA system for decadal predictions.

2.3 Implementation of the 4DEnVar-based WCODA system

The 4DEnVar method employed by the WCODA system is derived from the DRP-4DVar assimilation approach. The DRP-4DVar technique addresses the high computational demands of traditional 4DVar by employing an ensemble approach rather than utilizing the adjoint model (Wang et al., 2010). Zhu et al. (2022) demonstrated that the DRP-4DVar method significantly reduces computational time by approximately 50 % compared to traditional 4DVar systems. This advanced method enhances computational efficiency by projecting the high-dimensional state space onto a lower-dimensional subspace defined by an ensemble of historical samples. DRP-4DVar achieves an optimal solution within this sample space by aligning observations with model-generated historical time series over a four-dimensional window (Wang et al., 2010). The DRP-4DVar approach has been effectively implemented across multiple numerical models, demonstrating its accuracy and effectiveness (Zhao et al., 2012; Shi et al., 2021; Zhu et al., 2022). A comprehensive explanation of the DRP-4DVar method is provided by Wang et al. (2010). The DRP-4DVar method has also been implemented in a weakly coupled land data assimilation system in E3SMv2 (Shi et al., 2024).

Figure 1 Workflow of the 4DEnVar-based WCODA system utilizing the DRP-4DVar method for the E3SM model (modified from Fig. 1 in Shi et al., 2024).

Figure 2 Schematic diagram of the DRP-4DVar assimilation process within the 4DEnVar-based WCODA system for E3SM. The model background (x_b) includes atmospheric ($x_{\mathrm{b}}^{\mathrm{atm}}$), land ($x_{\mathrm{b}}^{\mathrm{lnd}}$), ice ($x_{\mathrm{b}}^{\mathrm{ice}}$), river ($x_{\mathrm{b}}^{\mathrm{river}}$), and oceanic ($x_{\mathrm{b}}^{\mathrm{ocn}}$) components of the fully coupled E3SMv2. The observational background ($y_{\mathrm{b}}^{\mathrm{ocn}}$) is defined by the model outputs of monthly mean ocean temperature (${\stackrel{\mathrm{‾}}{T}}_{\mathrm{b}}^{\mathrm{m}}$) and salinity (${\stackrel{\mathrm{‾}}{S}}_{\mathrm{b}}^{\mathrm{m}}$) using x_b as the initial state. The ocean observation ($y_{\mathrm{obs}}^{\mathrm{ocn}}$) represents the observed monthly mean ocean temperature (${\stackrel{\mathrm{‾}}{T}}_{\mathrm{obs}}^{\mathrm{m}}$) and salinity (${\stackrel{\mathrm{‾}}{S}}_{\mathrm{obs}}^{\mathrm{m}}$) from the EN4.2.1 reanalysis. The observational innovation ($y_{\mathrm{obs}}^{\prime }$) is calculated as the difference between the observed ocean temperature and salinity ($y_{\mathrm{obs}}^{\mathrm{ocn}}$) and the model's observational background ($y_{\mathrm{b}}^{\mathrm{ocn}}$). $x_{\mathrm{a}}^{\prime }$ denotes the analysis increment. The optimal analysis (x_a) encompasses both the optimal analysis of the ocean component ($x_{\mathrm{a}}^{\mathrm{ocn}}$) and the background states of other components. This optimal analysis (x_a) is used as the initial condition to produce the next month's forecast, transferring ocean reanalysis information to other components.

Figure 1 illustrates the workflow of the 4DEnVar-based WCODA system utilizing the DRP-4DVar approach within the fully coupled E3SMv2 model. The DRP-4DVar algorithm requires three primary inputs: observational innovation (${\stackrel{\mathrm{̃}}{y}}_{\mathrm{obs}}^{\prime }$), model background (x_b), and perturbation samples. Initially, a fully coupled E3SMv2 simulation is conducted for 1 month to generate both the model background (x_b) and the observational background (y_b). Specifically, the model background (x_b) refers to the monthly initial condition prior to the assimilation, while the observational background (y_b) denotes the monthly mean model output. Subsequently, the observational innovation (${\stackrel{\mathrm{̃}}{y}}_{\mathrm{obs}}^{\prime }$) is calculated as the difference in monthly mean ocean salinity and temperature between the EN4.2.1 reanalysis (y_obs) and the monthly mean model output (y_b). From 100 years of balanced pre-industrial control (PI control) simulations, 30 sets of monthly mean forecast samples (${\stackrel{\mathrm{̃}}{y}}^{\prime }$) are selected based on their highest correlations with the observational innovation. More specifically, the monthly mean forecast samples are computed by removing the long-term PI control monthly climatology from the selected PI control monthly mean output, which is then divided by the observational error. Correspondingly, 30 sets of monthly initial condition samples (x^′) for the monthly mean forecast samples are derived. The analysis increment is calculated within the perturbation samples, which consist of 30 monthly initial condition samples and their corresponding monthly mean forecast samples. Due to the limited number of samples and to diminish the influence of spurious correlations between distant grid points, the localization procedure is incorporated into the assimilation process (Wang et al., 2018). Finally, the DRP-4DVar algorithm solves for the analysis increment within the sample space, which is then added to the model background (x_b) to produce the optimal analysis (x_a).

Figure 2 delineates the assimilation process using the DRP-4DVar method within the 4DEnVar-based WCODA system for the fully coupled E3SMv2 model. This assimilation process includes both the analysis and forecast steps through each 1-month assimilation window. In the initial stage, the fully coupled E3SMv2 model employs the model background (x_b) as the monthly initial condition to run for 1 month, producing the monthly mean model outputs for ocean temperature and salinity ($y_{\mathrm{b}}^{\mathrm{ocn}}$). During the analysis step, the observational innovation ($y_{\mathrm{obs}}^{\prime }$) is computed by comparing the discrepancies between the EN4.2.1 reanalysis ($y_{\mathrm{obs}}^{\mathrm{ocn}}$) and the model's monthly mean outputs ($y_{\mathrm{b}}^{\mathrm{ocn}}$) for ocean temperature and salinity. The DRP-4DVar algorithm then utilizes this observational innovation and the PI control samples to compute the optimal analysis of the ocean component ($x_{\mathrm{a}}^{\mathrm{ocn}}$) at the start of the assimilation window. During the subsequent forecast step, the optimal analysis (x_a) includes both the optimal ocean analysis ($x_{\mathrm{a}}^{\mathrm{ocn}}$) and the background states of other components prior to assimilation. This optimal analysis serves as the new initial condition for the fully coupled E3SMv2 model to run for 1 month to generate the next month's forecast. During this fully coupled model integration, reanalysis information from the ocean component is transmitted to the other model components through interactions across multiple systems. Although the assimilation is directly applied to the ocean component, the use of the initial conditions of all components from the optimal analysis and the fully coupled climate model during the forecast step ensures that the reanalysis information from the optimal ocean analysis influences other components through interactions across multiple systems. Therefore, according to the definition of the WCDA system from previous studies (Carrassi et al., 2018; Zhou et al., 2024), this assimilation system is designated as the WCODA system. Using the same DA approach, Shi et al. (2024) documented the implementation of DRP-4DVar as a weakly coupled land data assimilation system in E3SMv2.

2.4 Experiment design

Two distinct numerical experiments are performed in this study to assess the effectiveness of ocean data assimilation within the 4DEnVar-based WCODA system. (1) The control simulation (CTRL) is a free-running fully coupled integration over a 72-year period from 1950 to 2021, driven exclusively by observed external forcings (e.g., solar radiation and greenhouse gas and aerosol concentrations). The observed external forcings, prescribed according to the CMIP6 protocol (Eyring et al., 2016), directly influence the atmospheric component and subsequently affect other components (e.g., land and ocean) through their coupling with the atmosphere. This free-running simulation allows for unrestricted interactions among the various Earth system components, including the atmosphere, land, and ocean. The CTRL simulation serves as a baseline for evaluating the assimilation effectiveness of the WCODA system. (2) The assimilation experiment (ASSIM) incorporates monthly mean ocean temperature and salinity data from the EN4.2.1 reanalysis into the ocean component of the fully coupled E3SMv2 model across all 60 ocean layers spanning the entire ocean depth. This assimilation is conducted using a 1-month assimilation window, covering the same 72-year period from 1950 to 2021. The assimilation run is initialized directly from the historical run in 1950 using the fully coupled state at the start of the simulation. At the beginning of each monthly assimilation window, the EN4.2.1 reanalysis information is incorporated into the ocean state variables, after which the fully coupled model continues with free integration. During this free integration process, the reanalysis information assimilated into the ocean ICs influences other model components through interactions across multiple systems. The historical external forcings for both the ASSIM and CTRL experiments are derived from the CMIP6 protocol (Eyring et al., 2016).

Figure 3 Temporal variation of the reduction rate of the cost function (in %) in the WCODA system based on the 4DEnVar method from 1950 to 2021.

Figure 4 Spatial patterns of root mean square error (RMSE) differences in ocean temperature (in °C) between ASSIM and CTRL across nine ocean layers from 1950 to 2021. The RMSE differences are shown for nine different ocean depths: (a) 5 m, (b) 45 m, (c) 85 m, (d) 135 m, (e) 327 m, (f) 528 m, (g) 708 m, (h) 879 m, and (i) 1106 m. Dotted areas represent statistical significance at the 95 % confidence level.

2.5 Assessment criteria

To comprehensively evaluate the effectiveness of the WCODA system, multiple quantitative metrics are employed, including the root mean square error (RMSE), correlation coefficient, and reduction rate of the cost function. The reduction rate of the cost function serves as a fundamental measure to assess the assimilation system's accuracy and is calculated using the formula:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\frac{\mathrm{1}}{\mathrm{2}}(y_{\mathrm{obs}}-y_{\mathrm{a}}{)}^T{\mathbf{R}}^{-\mathrm{1}}(y_{\mathrm{obs}}-y_{\mathrm{a}})-\frac{\mathrm{1}}{\mathrm{2}}(y_{\mathrm{obs}}-y_{\mathrm{b}}{)}^T{\mathbf{R}}^{-\mathrm{1}}(y_{\mathrm{obs}}-y_{\mathrm{b}})}{\frac{\mathrm{1}}{\mathrm{2}}(y_{\mathrm{obs}}-y_{\mathrm{b}}{)}^T{\mathbf{R}}^{-\mathrm{1}}(y_{\mathrm{obs}}-y_{\mathrm{b}})}}.\end{array}$$

Here, y_obs denotes the EN4.2.1 reanalysis, y_b represents the pre-assimilation observational background, y_a indicates the post-assimilation monthly mean model analyses, and R denotes the observation error covariance matrix. In this study, R is assumed to be diagonal, and its diagonal elements are statistically computed based on the variance of the EN4.2.1 ocean temperature and salinity data. The characteristics of R directly influence the assimilation process, where larger values reduce the relative weight of the EN4.2.1 reanalysis and smaller values increase it. Negative values of the reduction rate of the cost function signify the successful integration of reanalysis data into the model's state variables. To validate the correctness of this assimilation system, the EN4.2.1 reanalysis continues to be utilized as the reference data for evaluation.

Figure 5 Similar to Fig. 4 but for ocean salinity (in psu).

3 Results

3.1 Reduction rate of the cost function

In Fig. 3, the monthly variation in the reduction rate of the cost function for the 4DEnVar-based WCODA system is presented for the 72-year period from 1950 to 2021. As noted earlier, negative values of the reduction rate of the cost function indicate the successful incorporation of reanalysis data into the coupled model. However, the reduction rate is presented here as positive percentages to represent improvements due to the assimilation. The reduction rate of the cost function reaches 12.03 % in the first month. Over the entire 72-year period from 1950 to 2021, the average monthly reduction rate of the cost function is 4.20 % for all months in ASSIM. This average reduction rate of 4.20 % is comparable to the 4.4 % reduction rate reported by He et al. (2020a), who used a similar 4DEnVar-based assimilation system in a different climate model, further supporting the effectiveness of the 4DEnVar approach. The initial sharp reduction rate of the cost function reflects the rapid adjustments made by the model to align with the reanalysis data. As the assimilation progresses, subsequent iterations refine these adjustments, resulting in a slower rate of reduction. More importantly, the reduction rate of the cost function remains below the zero line in each month of assimilation, indicating consistent improvements due to the assimilation. These findings demonstrate the successful implementation of the WCODA system, confirming that the EN4.2.1 reanalysis data have been effectively integrated into the fully coupled model.

3.2 Performance of RMSE differences

Figure 4 illustrates the RMSE differences in monthly ocean temperature between ASSIM and CTRL from 1950 to 2021 across nine ocean layers. Negative values indicate a reduction in RMSE, signifying improvements due to assimilation, while positive values denote an increase in RMSE, indicating degradations. Overall, the assimilation from the WCODA system leads to marked improvements in ocean temperature simulations across most global regions. Both upper and deeper ocean layers exhibit widespread negative RMSE differences, indicating improvements after assimilation, particularly in the tropical and mid-latitude ocean regions. Notable regions of improvement include the north Atlantic, tropical and north Pacific, and parts of the Southern Ocean. However, increased RMSE values are observed near strong ocean currents and upwelling regions, such as the Gulf Stream, Agulhas Current, and California coast. These regions are characterized by strong horizontal gradients and mesoscale variability, which are not well captured by MPAS-O at relatively coarse resolution and hence present challenges for the assimilation system and likely contribute to diminished performance. In the upper ocean layers, RMSE performance is better during winter compared to summer in some regions, such as the tropical Pacific (Figs. A1 and A2). In the deeper layers, the assimilation still shows notable improvements in regions such as the north Pacific and parts of the Southern Ocean, though with more pronounced degradation observed in the equatorial Atlantic and parts of the Indian Ocean. This degradation in the deeper layers may be attributed to larger observational errors in these regions or the inherent complexity of deeper ocean processes that pose challenges for assimilation (Wunsch and Heimbach, 2007; Balmaseda et al., 2013).

Figure 6 Spatial patterns of the differences between ASSIM and CTRL for their correlations of ocean temperature with observations across nine ocean layers. Dotted regions indicate statistical significance at the 95 % confidence level. Panels (a) to (i) represent different ocean depths: (a) 5 m, (b) 45 m, (c) 85 m, (d) 135 m, (e) 327 m, (f) 528 m, (g) 708 m, (h) 879 m, and (i) 1106 m.

The RMSE differences for ocean salinity between ASSIM and CTRL across various ocean layers are presented in Fig. 5. The majority of ocean regions display notable improvements for ocean salinity after assimilation. In the upper ocean layers, significant enhancements are particularly evident in the north Pacific and parts of the north Atlantic. However, certain areas exhibit degradation in RMSE, particularly in parts of the south Pacific. In the deeper layers, the improvements are less extensive but remain evident in regions such as parts of the north Atlantic and north Pacific. However, RMSE degradation becomes notable in the equatorial Atlantic and parts of the Indian Ocean, highlighting the need for further improvements in these regions. The degradation in the deeper ocean layers can be attributed to two main factors: observational data limitations and challenges in representing deep-ocean processes in the model. For the EN4.2.1 reanalysis, the coverage and quality of observations tend to decrease with depth, potentially resulting in greater uncertainties in the deep ocean. This sparse observational coverage limits the constraints that assimilation can impose on the model state. Furthermore, in the E3SMv2 model, the complexity of simulating deep-ocean processes, such as vertical mixing and bottom water formation, may contribute to biases that are difficult to correct through assimilation.

Figure 7 Similar to Fig. 6 but for ocean salinity.

3.3 Performance of correlation differences

Figure 6 illustrates the differences between ASSIM and CTRL in their correlations with observed monthly ocean temperature from 1950 to 2021 across nine ocean layers. The seasonal cycle and linear trend have been removed before computing the correlations. Positive values denote an increase in correlation following assimilation, indicating improvements, whereas negative values suggest a decrease in correlation. In the upper ocean layers, the assimilation has led to improved correlations for ocean temperature across many ocean regions. Notably, the equatorial Pacific Ocean exhibits substantial improvements, indicating potential enhancements in modeling phenomena such as the El Niño–Southern Oscillation (ENSO). Further analysis of the winter Niño 3.4 index (Fig. A3) confirms that the assimilation improves the representation of ENSO variability, with the correlation coefficient increasing from 0.06 in CTRL to 0.79 in ASSIM. Moreover, parts of the north Pacific also exhibit noticeable improvements. In the deeper layers, improvements are observed in the western Pacific and parts of the Southern Ocean. However, certain areas exhibit diminished performance, possibly due to sparse observational coverage introducing higher uncertainty into the assimilation process or imbalances between ocean state variables during the assimilation (Edwards et al., 2015; He et al., 2020b). In summary, ASSIM has enhanced ocean temperature simulations by reducing RMSE (Fig. 4) and improving correlation (Fig. 6) across many ocean regions, with notable improvements in the upper ocean layers, including the equatorial Pacific and north Pacific.

The correlation differences for ocean salinity between ASSIM and CTRL across various ocean layers are depicted in Fig. 7. In the upper ocean layers, the majority of global ocean regions exhibit marked improvements for ocean salinity, with positive correlation differences dominating. Noteworthy improvements are evident in the tropical Pacific, north Pacific, and parts of the north Atlantic. In the deeper layers, the improvements in correlation become more localized, primarily concentrated in the western Pacific and parts of the Southern Ocean. Meanwhile, reductions in correlations are observed in parts of the equatorial Pacific and the south Atlantic, indicating the need for further improvements. Overall, ASSIM has improved simulations of ocean salinity by reducing RMSE (Fig. 5) and improving correlation (Fig. 7) in many ocean regions, with notable enhancements in the upper ocean layers, particularly in parts of the north Pacific and the western Pacific.

3.4 Vertical and temporal analysis of RMSE and bias for ocean temperature and salinity

Figure 8 presents the vertical profiles of the globally averaged RMSE of ocean temperature and salinity comparing ASSIM and CTRL. Negative values in the RMSE difference indicate a reduction in the global mean RMSE due to assimilation. For ocean temperature, the RMSE differences are relatively small but become more negative within the upper 85 m of the ocean. As the depth increases beyond 135 m, the RMSE differences become significantly negative, indicating a marked improvement in ocean temperature after assimilation. Unlike temperature, the salinity RMSE differences show substantial deviations in the upper 155 m of the ocean, indicating notable improvements. The RMSE differences gradually decrease as depth increases from 155 to 305 m, but a slight increase is observed between 305 and 1106 m. This suggests that the assimilation of salinity data has a more pronounced effect in the upper ocean than in deeper layers, possibly due to larger observational errors in these layers (Jacobs et al., 2021; Wang et al., 2015). The extended profiles in Fig. A4 indicate that below 1106 m, the RMSE differences between ASSIM and CTRL gradually decrease for both ocean temperature and salinity, suggesting the limited impact of assimilation in the deeper layers. In summary, these results emphasize the capability of the WCODA system to enhance the simulation accuracy for both ocean temperature and salinity.

Figure 8 Vertical profiles of the globally averaged RMSE differences between ASSIM and CTRL for (a) ocean temperature (in °C) and (b) ocean salinity (in psu) over the period from 1950 to 2021.

The temporal evolutions of the global mean bias and RMSE for vertically averaged ocean temperature and salinity in the top 1000 m are illustrated in Fig. 9. The temperature bias (Fig. 9a) in CTRL is persistently positive, indicating a systematic overestimation of ocean temperature. This overestimation in ocean temperature primarily originates from depths below 300 m (Figs. A5 and A6). In contrast, ASSIM consistently reduces this bias, with values approaching the zero line. Similarly, the temperature RMSE (Fig. 9b) highlights a significant decrease in RMSE for ASSIM compared to CTRL, reflecting a more accurate alignment with observed temperature. For ocean salinity, the salinity bias (Fig. 9c) reveals that CTRL maintains a consistent negative bias, suggesting an underestimation of ocean salinity. This salinity bias in CTRL is already prominent in the upper 300 m (Figs. A5 and A6). However, ASSIM effectively mitigates this bias, bringing the bias values closer to the zero line. Furthermore, the salinity RMSE (Fig. 9d) is notably lower in ASSIM than CTRL, indicating enhanced model performance and a closer match to observed salinity. Notably, it takes approximately 10–15 years for the biases in both temperature and salinity to stabilize near the zero line, reflecting an adjustment period where the assimilation system equilibrates. Overall, ASSIM exhibits superior performance relative to CTRL in reducing bias and RMSE for both ocean temperature and salinity.

Figure 9 Temporal variations in the global mean bias (a, c) and RMSE (b, d) for ocean temperature (in °C) and salinity (in psu) averaged over the upper 1000 m from 1950 to 2021. The red lines represent ASSIM, while the blue lines represent CTRL.

3.5 Climatological mean differences for sea surface temperature and salinity

Figure 10 presents the climatological mean differences between CTRL and observation, as well as between ASSIM and observation, for both sea surface temperature (SST) and salinity (SSS). Pronounced cold biases are evident in the SST difference between CTRL and observation (Fig. 10a), particularly in the tropical and north Pacific, north Atlantic, and parts of the Indian Ocean. Significant warm biases are observed in the Southern Ocean and parts of the south Atlantic. In contrast, these SST biases found in CTRL are substantially reduced by ASSIM (Fig. 10b), with cold biases in the north Pacific and north Atlantic diminished by approximately 1–2 °C and warm biases in the Southern Ocean corrected by about 1.5–2.5 °C. The SSS difference between CTRL and observation highlights a global pattern of salinity biases (Fig. 10c). The CTRL simulation generally underestimates salinity across most global oceans, indicating a widespread lower salinity. This fresh bias is particularly pronounced in the north Atlantic and north Pacific. Notably, in the Mediterranean Sea, CTRL exhibits a large positive salinity bias exceeding 2.5 psu. Compared with CTRL, ASSIM significantly reduces the overall fresh biases in CTRL (Fig. 10d). Notable improvements are observed in the north Atlantic and north Pacific, where salinity biases are reduced by approximately 0.5–1 psu, and in parts of the Southern Ocean, where reductions reach up to 1.5 psu. In summary, ASSIM demonstrates marked improvements in both SST and SSS biases compared to CTRL, emphasizing the importance and effectiveness of the WCODA system in enhancing model accuracy and reliability.

Figure 10 Climatological mean differences in sea surface temperature (a, b; in °C) and salinity (c, d; in psu) from 1950 to 2021. Panels (a) and (c) show the differences between CTRL and observation, while (b) and (d) show the differences between ASSIM and observation. Dotted areas indicate regions where the differences are statistically significant at the 95 % confidence level.

4 Conclusions

This study documents the development and assessment of the new 4DEnVar-based WCODA system in the fully coupled E3SMv2 model, employing the DRP-4DVar method. The DRP-4DVar approach significantly reduces computational demands by replacing the traditional adjoint model with the ensemble technique. As a weakly coupled assimilation system, the WCODA system independently assimilates ocean reanalysis data within the ocean component during the analysis step. However, during the subsequent forecast step, the reanalysis information from the optimal ocean analyses is propagated to other components of the Earth system through interactions across multiple systems, thereby enhancing the coherence of ICs across different components of the climate model.

Monthly mean ocean temperature and salinity data from the EN4.2.1 reanalysis are integrated into the ocean component of E3SMv2 from 1950 to 2021, which can be used to provide realistic ICs for decadal climate predictions. The effectiveness of the WCODA system has been assessed using several metrics, including the reduction rate of the cost function, RMSE differences, correlation differences, and model biases. The reduction rate of the cost function consistently shows negative values in each month over the 72-year period, indicating successful assimilation of the EN4.2.1 reanalysis data into the climate model. Compared to CTRL, ASSIM achieves significant reductions in RMSE and enhancements in correlation in the upper ocean layers, with notable improvements observed in parts of the north Atlantic, north Pacific, and Indian Ocean. ASSIM substantially mitigates model biases for SST and SSS observed in CTRL, particularly reducing cold biases in the north Pacific and north Atlantic by approximately 1–2 °C, correcting warm biases in the Southern Ocean by about 1.5–2.5 °C and significantly increasing salinity estimates to reduce the model fresh biases by approximately 0.5–1 psu in the north Atlantic and north Pacific and up to 1.5 psu in parts of the Southern Ocean.

Despite these advancements, the WCODA system exhibits limitations in certain regions, particularly in the deeper layers of the southern Pacific Ocean and south Atlantic. The reliance on the EN4.2.1 product could pose limitations to the assimilation process due to the sparse salinity observations and potential for static instabilities in data-sparse regions. Reanalysis products such as ORAS5 and GLORYS provide promising alternatives for mitigating these limitations. Future efforts should explore incorporating these reanalysis products into the WCODA system to improve the assimilation performance in challenging areas. Furthermore, expanding the application of the WCODA system to other components of the climate model, such as the atmosphere and sea ice, could enhance overall predictive skill. These developments are essential for providing more accurate and reliable long-term climate predictions, ultimately aiding in the formulation of energy-sector policies and management strategies.

Appendix A

Figure A1 Spatial patterns of RMSE differences between ASSIM and CTRL for ocean temperature (in °C) during summer. Results are presented for nine different ocean layers: (a) 5 m, (b) 45 m, (c) 85 m, (d) 135 m, (e) 327 m, (f) 528 m, (g) 708 m, (h) 879 m, and (i) 1106 m. Dotted areas represent statistical significance at the 95 % confidence level.

Figure A2 Similar to Fig. A1 but during winter.

Figure A3 Time series of the winter Niño 3.4 index from 1950 to 2021 for the observation (black line), ASSIM (red line), and CTRL (blue line). The correlation of the Niño 3.4 index with the observation in ASSIM and CTRL is also shown.

Figure A4 Vertical profiles of the globally averaged RMSE differences between ASSIM and CTRL for (a) ocean temperature (in °C) and (b) ocean salinity (in psu) with depths from 1106 to 5375 m.

Figure A5 Temporal variations in the global mean bias (a, c) and RMSE (b, d) for ocean temperature (in °C) and salinity (in psu) averaged over the upper 300 m. The red lines represent ASSIM, while the blue lines represent CTRL.

Figure A6 Similar to Fig. A5 but averaged over the upper 700 m.