An improved ENSO simulation by representing chlorophyll-induced climate feedback in the NCAR Community Earth System Model

The El Niño-Southern oscillation (ENSO) simulated in the Community Earth System Model of the National Center for Atmospheric Research (NCAR CESM) is much stronger than in reality. Here, satellite data are used to derive a statistical relationship between interannual variations in oceanic chlorophyll (CHL) and sea surface temperature (SST), which is then incorporated into the CESM to represent oceanic chlorophyll -induced climate feedback in the tropical Pacific. Numerical runs with and without the feedback (referred to as feedback and non-feedback runs) are performed and compared with each other. The ENSO amplitude simulated in the feedback run is more accurate than that in the non-feedback run; quantitatively, the Niño3 SST index is reduced by 35% when the feedback is included. The underlying processes are analyzed and the results show that interannual CHL anomalies exert a systematic modulating effect on the solar radiation penetrating into the subsurface layers, which induces differential heating in the upper ocean that affects vertical mixing and thus SST. The statistical modeling approach proposed in this work offers an effective and economical way for improving climate simulations.

. A schematic diagram illustrating the inclusion of the CHL-induced bio-feedback using the NCAR CESM, consisting of the original atmosphere model (CAM5) and the OGCM (POP2), and an embedded statistical representation for interannual CHL variability developed in this work. The total CHL field is separated into its climatological component (CHL clim ) and interannual anomaly component (CHL inter ), written as CHL = CHL clim + CHL inter , which directly affect the vertical distribution of the incoming solar radiation in the upper ocean. The statistical relationships between the interannual variations in SST and CHL are derived from satellite observations and used to explicitly calculate CHLAs in the CESM to represent CHL-induced heating effects.  Table 1. Experiments performed to demonstrate the CHL-induced heating effects using the CESM. CHL clim is a climatological run that uses CHL as its seasonally varying climatology prescribed from satellite observations. CHL inter is an interannual run that considers the interannually varying CHL effects (interannual CHL anomalies calculated from the SST anomalies by using the statistical CHL model). In addition, an El Niño case in the model years 153 and 154 is selected for two more experiments that use the CESM designed to illustrate direct effects of CHLAs on SST differences and the underlying processes responsible for the differences in two runs; CHL clim -EN is a run that uses CHL as its seasonally varying climatology, and CHL inter -EN is a run that considers the interannually varying CHL effects. These two runs are restarted from the same initial state at model year 153 and are time integrated for 2 years. The results from these two runs are compared with each other to identify the effects on the detailed evolution of some related fields.
averaged SST index is shown in Fig. 3. The amplitude of ENSO in the CHL inter run is substantially smaller than that in the CHL clim run, especially during the first several decades of the 50-year integration period. The damping effect of CHLAs on ENSO occurs not only during El Niño events but also during La Niña events. The commonly used Niño3.4 SST series are used to quantify the dominant time scales of ENSO (Fig. 4). The interannual variability has a sharp peak at 4 years in the CHL inter run compared with 5 years in the reference CHL clim run, with a difference of approximately 12 months in the oscillation periods.
Quantitatively, Fig. 5 shows the standard deviation of the SSTA for the observations and the CHL clim and CHL inter experiments. The standard deviation can be used as a metric for measuring ENSO amplitude. The strongest interannual variability of the observed SSTA appears in the eastern equatorial Pacific, which has a standard deviation exceeding 1.0 °C (Fig. 5a). In the CHL clim run, the model well captures the spatial distribution of the SST variability; however, the amplitude of the variability is significantly overestimated, with the maximum standard  deviation being over 1.5°C (Fig. 5b). This bias is a long-standing problem for the CESM, and it can be attributed to many factors 1,2 . When considering the interannually varying effect of CHLAs in the CHL inter run, the SST amplitude bias can be effectively reduced, with the maximum standard deviation being reduced to 1.25°C (Fig. 5c). The area averaged standard deviation of SSTAs in the Niño3 and Niño4 regions is summarized in Table 2. In the Niño4 region, the standard deviation is 1.03 for the CHL inter run and 1.22 for the CHL clim run, whereas it is 0.71 for the corresponding observation. Thus, the inclusion of the interannual varying CHL effect leads to a damped ENSO amplitude. A similar reduction in the SST variability also occurs in the Niño3 region, where the standard deviation of SSTAs is 1.17 for the CHL inter run and 1.46 for the CHL clim run, and the standard deviation is 0.93 for Figure 4. The power spectra for the Niño3 SST anomalies estimated for the observations (black), the CHL inter run (red), and the CHL clim (red dash) run. The dot-dashed line is the 95% significance level for these runs, assuming a white noise process. Since the interannual variability in the tropical Pacific has a strong interaction with the mean states in the region, the effects on mean states and seasonal cycles are presented in the Supplementary section. It is evident that there is little change in the mean thermocline and seasonal cycle of SST in the equatorial Pacific due to the feedback of CHL and SST. For example, the seasonal cycle of SST in the eastern Pacific has been changed about 0.2 K when the CHL-induced feedback on SST is taken into account, with slight amplification of seasonal cycle of SST in the feedback (CHL inter ) run (a warming in the spring and a cooling in the fall) compared with the non-feedback (CHL clim ). Previous studies have suggested that an enhanced seasonal cycle tends to reduce interannual variability in the tropical Pacific 12 . It is unlikely that such a slight increase in the seasonal cycle in the feedback (CHL inter ) run can be responsible for the reduction in the interannual variation of SST as seen in the two runs.
The processes responsible for the reduced ENSO amplitude as seen in the CHL inter run are examined (the details are presented in the Supplementary section, including a calculation of the SST budget). CHL in the ocean affects the penetration of the incoming downward shortwave flux (SWDN) at the sea surface, and two heating terms are directly related with the bio-effects in the upper ocean. One is the SWDN component that is absorbed within the first layer (Q abs ) and another is the SWDN component that penetrates into the subsurface layers through the bottom of the first layer (Q pen ). It is found that the Q abs and Q pen fields are two terms that are directly modulated by interannual anomalies of CHL 19 . Therefore, the relationships among these related fields can be analyzed to understand the effects of the CHL-related feedback and the involved processes. As detailed in the Supplementary section, two pathways are possible for the variations in CHL that can affect SST in the equatorial Pacific, depending on the relative dominance of which term is modulated most pronouncedly, systematically, and coherently. On one hand, if a change in CHL induces a coherent and systematic change to Q abs , a direct heating effect on SST can be a dominant process that causes a corresponding change to SST; this is referred to as a direct effect. Then, it follows that the differences in CHL and Q abs between the two runs should be well matched with each other. On the other hand, if a change in CHL causes a systematic and significant modulation of Q pen , differential heating can be then induced vertically between the first layer (Q abs ) and subsurface layers (Q pen ). These processes can further modify the stratification and vertical mixing in the upper ocean and thus affect SST; this is referred to as an indirect effect on SST. Therefore, the extent to which these two heating terms are affected by interannual CHL anomalies implies that different dominant processes are involved in the bio-effects. How Q abs and/or Q pen are affected by interannual CHL anomalies can be a good indicator of which influence pathways are taken and the underlying processes operating in association with the CHL-induced heating feedback in the NCAR CESM.
To identify the effect of the reduced CHL concentration on the total SWDN reaching the ocean surface, Q abs and Q pen are analyzed during the El Niño evolution in the feedback and non-feedback runs (the Supplementary section). Interannual CHL anomalies can have a direct effect on the total SWDN reaching the ocean surface and the vertical penetration in the upper ocean. It turns out that the total SWDN is not systematically and coherently affected by CHLAs. Similarly, the difference fields in Q abs do not show coherent and systematic relationships with CHLAs. Indeed, the differences in Q abs and SST between the two runs do not well match each other. As a result, the direct heating effects within the first layer induced by CHLAs (Q abs ) cannot be a dominant process that causes the systematic differences in the SST simulations between these runs. In fact, the differences in Q abs are negatively correlated with those in SST. This is consistent with the understanding that on interannual time scales associated with ENSO, heat flux and related direct solar radiation heating/cooling play a damping role on SST warming during El Niño evolution. As the SST difference between the two runs is spatially well defined, with consistent cooling (warming) in the eastern (western-central) equatorial Pacific (the Supplementary section), oceanic dynamical processes must play a dominant role in regulating SST variability in the tropical Pacific.
Indeed, it is found that although the difference in the total SWDN between the feedback and non-feedback runs is not systematic, the difference in Q pen exhibits a well-defined and coherent pattern that is well matched with CHLAs (the Supplementary section). For example, a consistent increase in Q pen is seen during El Niño events in the feedback run, which is accompanied with a low CHL concentration. The negative CHLAs increase penetration into the subsurface later throughout the bottom of the first layer (Q pen ) and reduce the component that is absorbed within the first layer (Q abs ). So, the increased Q pen field warms the temperature of the subsurface layers. Such vertical redistribution of the absorbed SWDN (i.e., warming in subsurface layers but cooling in the first layer) acts to enhance the vertical mixing and further decrease the SST. As a result, the direct effect of SWDN absorbed in the first layer is unlikely to be a major contributor to the SST differences; instead, the indirect effect of the differential heating induced by Q pen and Q abs resulting from interannual anomalies of CHL can be a major contributor. As seen, relative to the non-feedback run, the warm SST anomalies in the feedback run are weakened due to the inclusion of CHLAs in the ocean, leading to a damping effect on the El Niño event.
In short, Q pen is the term that is systematically affected by interannual CHL anomalies, and it can be responsible for the differences in the ENSO simulations seen in the feedback and non-feedback runs. Through the  dynamic response induced by modulation of Q pen in the upper ocean, the SSTs are modulated. Therefore, the SST difference between the two runs can be attributed to the effects of interannual CHL anomalies on Q pen . This is consistent with previous modeling studies that demonstrate that the CHL anomalies act to modulate ENSO mainly through modulation of solar radiation penetration throughout the bottom of the oceanic mixed layer 17,19 .

Discussion
Recent studies indicate interactive coupling between ocean ecosystem and climate system on various space-time scales [24][25][26] . As a main focus in the tropical Pacific, ENSO has been extensively investigated; various observations-based analyses have demonstrated significant feedbacks of the ocean biology-induced heating on interannual variability [27][28][29][30][31] . Because interactions between ocean biology and physics are very complicated, numerical models provide a powerful tool to study the related feedback and coupling [32][33][34] . However, our current understanding of the influences of oceanic biology-induced heating feedback on ENSO and of the related physical processes remains elusive in ocean and climate models. In particular, the relationships between the existing ENSO biases in the NCAR CESM and how bio-feedback is represented are not well known. Previously, great efforts have focused on the physical aspects, but with less attention to the oceanic biological attributions. In this study, an SVD-based statistical model for interannual CHL anomalies is derived to identify its relationship with SST 23 , which is then incorporated into the NCAR CESM. We show that such an SVD-based statistical model can well reproduce the main characteristics of interannual CHL variability compared with observations. Numerical runs with and without the interannually varying CHLA effects are performed using the CESM and compared to each other. These results clearly demonstrate that climate biases in the current CESM can be greatly reduced by considering the bio-feedback. Specifically, the CHLA has obvious damping effects on ENSO by shortening the oscillation periods. To explain this, two more experiments are performed with a focus on one specific El Niño event. The difference fields between the two runs are analyzed to reveal dominant terms that are directly modulated by interannual CHL anomalies and thus the underlying processes responsible for the differences in ENSO simulations. The effects and related processes during the El Niño evolution are clearly illustrated in association with interannual CHL anomalies. Q abs (the absorbed component of the incoming solar radiation in the first layer) and Q pen (the component penetrating into the subsurface layers through the bottom of the first layer) are two terms that are directly affected by CHLAs. The differences in Q abs between the two runs are not well matched with CHLAs; Q abs is not systematically modulated by CHLAs, so it is not a major factor directly responsible for the SST difference. Instead, the Q pen field is systematically modulated by CHLAs. The differences in CHL and Q pen between the two runs are well matched with each other. As shown, the Q pen field causes a dynamic response in the upper ocean, leading to modulations of SSTs. During an El Niño event, for instance, the negative CHLA allows solar radiation to penetrate more into the subsurface layers throughout the bottom of the first layer and to be absorbed less within the first layer. This directly leads to warming in the subsurface layers and cooling in the first layer. Furthermore, the temperature decrease in the first ocean model layer and temperature increase in the subsurface layers induce a vertical differential heating, which enhances vertical mixing and indirectly leads to SST decrease and so, the El Niño event is damped. As Q pen is responsible for the damping effect on SST, the effects indirectly induced by CHLAs play a more dominant role in the damping effect on ENSO. At present, obvious large biases exist in many of the state-of-the-art fully coupled models that explicitly include a biogeochemical component. For example, climatological CHL simulations still have large discrepancies in the tropical Pacific. In this study, an anomaly approach is taken in which the CHL climatological field is prescribed from observation, and its anomaly field is derived from SST anomalies. It is evident that the SVD-based statistical model for CHLAs is very efficient computationally and can well reproduce the main characteristics of the CHLAs compared to those of the observations. In this study, the feasibility and effectiveness of this statistical modeling approach are clearly justified by its success in improving ENSO simulations when this simple model is implemented in the NCAR CESM. It is expected that this simple model for interannual CHL variability can be used in other climate models to consider bio-feedback adequately, with climate simulations improved significantly.
This paper is intended to illustrate ENSO modulations that are induced by ocean biology-induced feedbacks within the tropical Pacific climate system. Numerical experiments are performed using the empirical model with the parameter representing feedback intensity that is tunable. Thus, the results presented in this paper are preliminary and need to be validated using more realistic and comprehensive models. For example, there is uncertainty/ sensitivity in the represented amplitude of the bio-feedback to SST changes as indicated by the feedback strength factor α. In the current modeling study to reasonably represent the intensity of the OBH effects, α is doubled (Note that when taking α = 2 in the feedback simulation, the standard deviations of interannual H p variability are well comparable to those from satellite-based estimates and thus, the OBH effect can be reasonably represented in the model simulation. Sensitivity experiments testing the effects of varying α on ENSO simulations are performed (Supplementary section). Indeed, within the statistical modeling context of H p , the choice of this parameter is rather arbitrary, and the modeling results are sensitive to the related parameters used. Additionally, a long term integration is also performed and the results are presented in Supplementary section.

Methods
The physical model used in this study is the Community Earth System Model version 1.0.5 (CESM1.0.5) developed at the National Center for Atmospheric Research (NCAR); see the Supplementary section for details 21,22,35 . It is notable that the oceanic model in the NCAR CESM already includes a biogeochemical component that can be used to describe the biogeochemical processes and their interactions with the ocean physical system. Its performance, however, is still not accurate enough for practical uses. Additionally, the CESM with an interactive and explicit biogeochemical component included is extremely time consuming and impractical to run extensively. Instead, CHL is commonly prescribed to represent seasonally varying climatology derived from satellite data; we use this CHL setting as our reference for the CESM simulation (CHL clim ), in which the oceanic biogeochemistry model is not activated, so there is no interannual bio-coupling between the oceanic physical model (POP2) and biogeochemical model. As such, the CHL clim simulation does not include interannually varying CHL-induced climate feedback and the interactions between the physical and biogeochemical systems. As shown in this work, there are large biases in ENSO simulations when using this prescribed CHL clim setting for the CESM.
Based on our previous work, we adopted a statistical modeling approach to represent the interannually varying CHL-induced climate feedback. This is based on the fact that on interannual time scales, ENSO is a major driver in the tropical Pacific, and interannual anomalies of CHL are represented as a response, which in turn can have feedback effects on ENSO. In particular, observed changes in ocean biology (e.g., CHL) closely follow physical changes (e.g., SST) during ENSO cycles. Based on the analysis conducted by Zhang et al. 23 , in this study, we adopted a singular value decomposition (SVD) method to capture the co-variations between the historical CHLA and SSTA during Sep. 1997-May 2016; a corresponding statistical model is then derived to relate the SSTA to the CHLA. Then, given an SSTA, the statistical model is used to calculate CHLAs with the first 5 SVD modes being retained, accounting for approximately 65% of the variance. This statistical model for CHLAs covers only the tropical Pacific (25°S-25°N, 124°E-76°W).
The derived statistical model for CHLAs is then incorporated into the CESM to represent CHL-induced feedback (Fig. 1). The total CHL field in the CESM can be described as where α is a scaling factor that is introduced to represent the amplitude of the derived CHLA, and f is a relationship between the interannual variations in SST and CHL that is determined using SVD analysis. In this way, the CHL in the CESM is composed of the seasonally varying climatological component (chl clim ) and the interannually varying component (chla). chl clim is prescribed from the historical satellite-based observations, and chla is calculated at each ocean time step based on the instant SSTA through the SVD-based statistical model (Fig. 1). The SSTA is calculated at each ocean time step by subtracting the SST climatology field (SST clim ) from the modeled SST (the SST clim is derived from a long-term simulation with prescribed CHL climatology). As such, an anomaly coupling is adopted for the physical (SST) and biological (CHL) systems with the climatological CHL field being prescribed. When only the first five SVD modes are maintained in the statistical model, some variance is inevitably lost because the SVD modes contain only approximately 65% of the variance. In this study, α = 2.0 is taken so that the amplitude of the CHLA simulated from the statistical model can be comparable to the observations (Fig. S1 in the supplementary section). Note that such an SVD-based statistical model is only applied to the tropical Pacific. Using this statistical method, interannually varying CHL-induced climate feedback can be considered in the NCAR CESM even though the ocean biogeochemical component is disabled. Thus, the net downward shortwave flux at the ocean surface (SWDN) from the atmospheric model is attenuated in the upper ocean due to the CHL field, which includes the interannual variation in the tropical Pacific. The vertically redistributed SWDN in the upper ocean acts to affect the SST. Interactively, the changed SST can also influence the atmosphere and the SWDN. Thus, the interactions between the CHLA and the physical system can be represented in the NCAR CESM (Fig. 1).
To assess the bio-effects, numerical experiments are conducted using the NCAR CESM (Table 1). First, the CESM with prescribed chl clim is conducted for 150 years as a spin-up; then, two runs continue for an additional 50-year run (years 151-200). The CHL clim run only includes the chl clim field, and the CHL inter run includes both the prescribed chl clim and chla fields derived from the SVD-based model ( Table 1). The results from the CHL clim and CHL inter runs are analyzed and compared with each other. As the only differences in the two runs are interannual CHL variability included or not in the tropical Pacific, the effects of CHLAs on ENSO can be isolated in a clear way.
Additionally, two case studies are conducted to illustrate the detailed evolution of the related anomaly fields and the physical processes by focusing on a specific El Niño event that appears in the time integration during the model years 156-157 (Table 1). CHL clim -EN is a run in which CHL is used as its seasonally varying climatology; CHL inter -EN is a run in which interannually varying CHL effects are considered. These two runs are restarted from the same initial state at model year 153 and are time integrated for 2 years only. Detailed analyses are presented in the Supplementary section.