Performance of the Taiwan Earth System Model in Simulating Climate Variability Compared With Observations and CMIP6 Model Simulations

Abstract

This study evaluates the performance of the Taiwan Earth System Model version 1 (TaiESM1) in simulating the observed climate variability in the historical simulation of the Coupled Model Intercomparison Phase 6 (CMIP6). TaiESM1 is developed on the basis of the Community Earth System Model version 1.2.2, with the inclusion of several new physical schemes and improvements in the atmosphere model. The new additions include an improved triggering function in the cumulus convection scheme, a revised distribution‐based formula in the cloud fraction scheme, a new aerosol scheme, and a unique scheme for three‐dimensional surface absorption of shortwave radiation that accounts for the influence of complex terrains. In contrast to the majority of model evaluation processes, which focus mainly on the climatological mean, this evaluation focuses on climate variability parameters, including the diurnal rainfall cycle, precipitation extremes, synoptic eddy activity, intraseasonal fluctuation, monsoon evolution, and interannual and multidecadal atmospheric and oceanic teleconnection patterns. A series of intercomparisons between the simulations of TaiESM1 and CMIP6 models and observations indicate that TaiESM1, a participating model in CMIP6, can realistically simulate the observed climate variability at various time scales and are among the leading CMIP6 models in terms of many key climate features.

Introduction

The Taiwan Earth System Model version 1 (TaiESM1) is developed on the basis of the Community Earth System Model version 1.2.2 (CESM1.2.2; Hurrell et al., 2013) by implementing several improvements in the parameterization schemes in the atmospheric component of CESM1.2.2. The modifications include the following: (1) replacing the three-mode version of the Modal Aerosol Module (Liu et al., 2012) aerosol scheme with the Statistical-Numerical Aerosol Parameterization scheme (Chen et al., 2013); (2) replacing the trigger function in the Zhang-McFarlane convection scheme (Zhang & Mcfarlane, 1995) with one that considers convection inhibition and the initiation of elevated instability (Wang et al., 2015); (3) improvement in the cloud fraction scheme (Park et al., 2014; Zhang, 2003) to allow cloud fraction determination based on the distribution of the total water content instead of the relative humidity threshold (Shiu et al., 2021); and (4) implementing a surface radiation scheme that considers the effect of three-dimensional topography on the absorption of shortwave solar radiation (Lee et al., 2013) into the Rapid Radiative Transfer Method for GCMs scheme of CAM5 (Iacono et al., 2008). Detailed descriptions of the developments and tuning of TaiESM1 and the evaluation of its performance based on the piControl run and historical runs conducted with the Coupled Model Intercomparison Phase 5 (CMIP5) setup are provided in an accompanying report (Lee, Kim, et al., 2020).

A basic requirement of a climate model is satisfactory performance in the simulation of mean climatology. In addition, a model used for future climate projections should be able to realistically simulate the observed climate variability at various time scales, which is modulated not only through the long-term mean climate state but also through feedback to the mean state. For this purpose, the TaiESM1 is designed to enhance the ability of simulating variability from diurnal to interdecadal time scales. The basic approach is to improve or implement parameterization schemes such that the modules could more realistically represent the observed temporal and spatial variations. In an accompanying report, Lee, Kim, et al. (2020) demonstrated that the TaiESM1, when driven by the forcing designed for the CMIP5 historical experiments, can simulate long-term climatological mean fields with a score higher than those of most other CMIP5 models. In this study, we further demonstrated the ability of TaiESM1 to simulate the seasonal cycle, monsoon evolution, synoptic and intraseasonal variability, characteristics of precipitation extremes, the diurnal cycle, the El Niño-Southern Oscillation (ENSO), interannual teleconnection variability, and oceanic interdecadal oscillations in the historical experiments of the Coupled Model Intercomparison Project Phase 6 (CMIP6).

The remainder of the paper is organized as follows. The methodology for analyzing climate variability in various time scales is described in Section 2. An evaluation of the mean state and warming of the historical simulation is presented in Section 3. Section 4 presents the evaluation of seasonal evolution and major monsoon systems. Section 5 details intraseasonal and synoptic variability, extremes, and the diurnal rainfall cycle. Interannual-interdecadal variability is reported in Section 6, discussion is in Section 7, and summary is provided in Section 8.

Model, Experimental Setup, and Data

The historical experiment is conducted using TaiESM1, driven by the forcing provided by the CMIP6 (Eyring et al., 2016) for the 1850–2014 period, following the procedure described by Lee, Wang, et al. (2020). The historical run is initiated from the year 671 in the piControl run of TaiESM1 with a horizontal resolution of 0.9° latitude × 1.25° longitude and 30 vertical layers. The performance of the model is evaluated for two data periods: 1915–2014 and 1980–2014. The longer period is used in the evaluation of interdecadal variability, such as the Atlantic Multidecadal Oscillation (AMO) and the Pacific Decadal Oscillation (PDO), whereas the shorter period that covers the satellite observation era is used for evaluating phenomena of shorter time scales, such as seasonal, intraseasonal, synoptic, extreme weather, and interannual scales.

Table 1 presents all the CMIP6 historical model runs used in this study for evaluating the performance of TaiESM1. Models are numbered alphabetically, and the numbers are used in model analysis presented below. The data are downloaded from the CMIP6 archive (https://esgf-node.llnl.gov/search/cmip6/). However, selection of the models used in each analysis is based on the availability of the required variables for analysis. A complete list of the models used in individual sessions of the study can be found in Table S1. Table S2 lists the observational data sets used in this study for validation. Most of the free tropospheric variables are evaluated against the Collaborative Reanalysis Technical Environment Multi-Reanalysis Ensemble version 2 (MRE2; Potter et al., 2018). The MRE2 is a product of ensemble average of 7 reanalysis products, including CFSR (Saha et al., 2010), ERA-Interim (Dee et al., 2011), MERRA (Rienecker et al., 2011), MERRA-2 (Gelaro et al., 2017), JRA-25 (Onogi et al., 2007), JRA-55 (Kobayashi et al., 2015), and 20CRv2c (Compo et al., 2011). For selected atmospheric variables, it is found that the errors of ensemble average are reduced compared with individual reanalysis (Potter et al., 2018). In addition, precipitation data obtained from the Global Precipitation Climatology Project (GPCP V2.3, 1980–2014; Adler et al., 2003; Huffman et al., 2009), radiation fluxes, including outgoing longwave radiation (OLR, 1980–2014) from the Clouds and the Earth's Radiant Energy System-Energy Balanced and Filled (CERES-EBAF; Loeb et al., 2018; Wielicki et al., 1996), and sea-surface temperature (SST; HadSST V1.1, Rayner et al., 2003, 1980–2014 for the ENSO and 1915–2014 for the AMO and PDO) are used in the model evaluation. The historical warming trend is compared with two sets of observations: the Hadley Centre—Climate Research Unit Temperature Anomalies (HadCRUT; Jones et al., 2012) and the Berkeley Earth Surface Temperature (BEST; Rohde et al., 2013).

Table 1 The 45 CMIP6 Coupled Atmosphere-Ocean Climate Models Used in the Historical Warming Analysis

#	Model	Description	References	Spatial resolution (# of lon × # of lat)
1	ACCESS-CM2	Commonwealth Scientific and Industrial Research Organisation (Australia)	Bi et al. (2020)	192 × 144
2	ACCESS-ESM1-5		Ziehn et al. (2020)	192 × 145
3	AWI-CM-1-1-MR	Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research (Germany)	Semmler et al. (2020)	384 × 192
4	AWI-ESM-1-1-LR		Sidorenko et al. (2015)	192 × 96
5	BCC-CSM2-MR	Beijing Climate Center, Beijing (China)	Wu et al. (2019)	320 × 160
6	BCC-ESM1		Wu et al. (2020)	128 × 64
7	CAMS-CSM1-0	Chinese Academy of Meteorological Sciences (China)	Chen, Li, et al. (2019)	320 × 160
8	CanESM5	Canadian Centre for Climate Modelling and Analysis, Environment and Climate Change Canada (Canada)	Swart et al. (2019)	128 × 64
9	CESM2-FV2	National Center for Atmospheric Research, Climate and Global Dynamics Laboratory, (USA)	Danabasoglu et al. (2020)	144 × 96
10	CESM2			288 × 192
11	CESM2-WACCM-FV2			144 × 96
12	CESM2-WACCM			288 × 192
13	CIESM	Department of Earth System Science, Tsinghua University (China)	Lin et al. (2020)	288 × 192
14	CMCC-CM2-SR5	Fondazione Centro Euro-Mediterraneo sui Cambiamenti Climatici (Italy)	Cherchi et al. (2018)	288 × 192
15	CNRM-CM6-1	Centre National de Recherches Météorologiques (France)	Voldoire et al. (2019)	256 × 128
16	E3SM-1-0	Lawrence Livermore National Laboratory, Department of Energy (USA)	Golaz et al. (2019)	360 × 180
17	E3SM-1-1-ECA		Burrows et al. (2020)	360 × 180
18	E3SM-1-1			360 × 180
19	EC-Earth3	EC-Earth consortium, Rossby Center, Swedish Meteorological and Hydrological Institute/SMHI (Sweden)	Döscher et al. (2021)	512 × 256
20	EC-Earth3-Veg			512 × 256
21	EC-Earth3-Veg-LR			320 × 160
22	FGOALS-g3	Chinese Academy of Sciences (China)	Li et al. (2020)	180 × 80
23	FIO-ESM-2-0	First Institute of Oceanography, State Oceanic Administration (China)	Bao et al. (2020)	288 × 192
24	GFDL-CM4	National Oceanic and Atmospheric Administration, Geophysical Fluid Dynamics Laboratory (USA)	Held et al. (2019)	288 × 180
25	GFDL-ESM4		Dunne et al. (2020)	288 × 180
26	GISS-E2-1-G	Goddard Institute for Space Studies, National Aeronautics and Space Administration (USA)	Kelley et al. (2020)	144 × 90
27	GISS-E2-1-H			144 × 90
28	HadGEM3-GC31-LL	Met Office, Hadley Centre (UK)	Williams et al. (2018)	192 × 144
29	INM-CM4-8	Institute for Numerical Mathematics, Russian Academy of Science, Moscow (Russia)	Volodin et al. (2013)	180 × 120
30	INM-CM5-0		Volodin et al. (2017)	180 × 120
31	IPSL-CM6A-LR	Institut Pierre Simon Laplace (France)	Boucher et al. (2020)	144 × 143
32	KACE-1-0-G	National Institute of Meteorological Sciences/Korea Meteorological Administration, Climate Research Division (Republic of Korea)	Lee, Kim, et al. (2020)	192 × 144
33	MCM-UA-1-0	Department of Geosciences, University of Arizona (USA)	Delworth et al. (2002)	96 × 80
34	MIROC6	Japan Agency for Marine-Earth Science and Technology (Japan)	Tatebe et al. (2019)	256 × 128
35	MICROC-ES2L		Hajima et al. (2020)	128 × 64
36	MPI-ESM-1-2-HAM	ETH Zurich, Switzerland; Max Planck Institut fur Meteorologie (Germany)	Tegen et al. (2019)	192 × 96
37	MPI-ESM1-2-HR	Max Planck Institute for Meteorology (Germany)	Müller et al. (2018)	384 × 192
38	MPI-ESM1-2-LR		Mauritsen et al. (2019)	192 × 96
39	MRI-ESM2-0	Meteorological Research Institute (Japan)	Yukimoto et al. (2019)	320 × 160
40	NESM3	Nanjing University of Information Science and Technology (China)	Cao et al. (2018)	192 × 96
41	NorCPM1	Climate Modeling Consortium consisting of Center for International Climate and Environmental Research (Norway)	Counillon et al. (2016) and Wang et al. (2017)	144 × 96
42	NorESM2-LM		Seland et al. (2020)	144 × 96
43	NorESM2-MM			288 × 192
44	SAM0-UNICON	Seoul National University (Republic of Korea)	Park et al. (2019)	288 × 192
x	TaiESM1	Research Center for Environmental Changes, Academia Sinica (Taiwan)	Lee, Wang, et al. (2020)	288 × 192
45	UKESM1-0-LL	Met Office, Hadley Centre (UK)	Sellar et al. (2019)	192 × 144

For evaluation of the climate modes, the empirical orthogonal function (EOF) method is commonly used to extract geographical patterns with maximum variability. However, using the EOF modes derived from models in model performance evaluation presents several challenges in comparison with the observed leading climate modes. For example, Lee et al. (2019) explored the interannual and decadal modes and found that the order of model-derived EOFs may need to be swapped before comparisons are made with the observed EOFs. This problem is more severe in the evaluation of climate modes, such as Pacific-Japan (PJ) and Pacific-North America (PNA) patterns, which do not explain variance as efficiently as other leading modes. To avoid the ordering of climate modes based on EOFs, Lee et al. (2019) proposed the common basis function (CBF) method, wherein model anomalies are projected onto the geographical patterns of the observed EOFs for comparison. We found that the derived model modes based on the CBF method are generally more consistent with the observed modes and that the CBF method provides a more consistent framework for model evaluation.

To quantify the model performance against observations, we use the Pearson correlation coefficients and the skill scores for pattern evaluations. The Pearson correlation R is defined as [Image Omitted. See PDF]where $x_i$ and $y_i$ are the simulation and observation data, respectively (Wilks, 2011). The skill score S (Taylor, 2001) is calculated as follows: [Image Omitted. See PDF]where R represents the spatial pattern correlation coefficient between the observation and model simulation and $\sigma$ is the ratio of the spatial standard deviation of the model simulation relative to that of the observation. $R_0$ is the maximum correlation attainable and is assumed to be 1 here. We then ranked model performances based on these metrics for model intercomparisons in the following analysis. All data used to plot model intercomparison figures can be found in the supplementary materials.

Climatological State and the Evolution of Historical Simulations

Figure 1 presents the model global mean of near-surface (2 m) air temperature (SAT) anomalies of the historical simulations of TaiESM1 with the mean temperatures in 1951–1980 as a reference. Two sets of observations, HadCRUT and BEST, are plotted for comparison. The gray lines in Figure 1 denote the temperature time series of other CMIP6 models. TaiESM1 responds with similar magnitudes of cooling to CMIP6 forcing during major volcanic eruptions, such as those of Krakatoa (1883), Agung (1963), and Pinatubo (1991), compared with observations, implying that the model sensitivity of radiative forcing to stratospheric aerosols is reasonable. From 1850 to 1950, the SAT anomaly in TaiESM1 is ∼0.3°C warmer than the observed anomaly; moreover, decadal fluctuations rather than a warming trend are noted before 1950. During this period, TaiESM1 is at the warmer end in the CMIP6 model spectrum of SAT anomalies. After 1960, the change in SAT simulated by TaiESM1 is close to the observed value in 2014. This feature is highly similar to the SAT evaluation when TaiESM1 is driven by CMIP5 historical forcing (Lee, Wang, et al., 2020). The causes of the warm bias in the beginning of the period are unknown and require further investigation.

Figure 1. Historical warming trend of TaiESM1 with 1951–1980 as a reference (black line). Two sets of observed global temperature data, HadCRUT (blue) and BEST (red), are also shown for comparison. Temperature anomalies calculated from other CMIP6 models are presented as gray lines.

Figure 2 presents the global patterns of the observed annual mean precipitation and surface temperature and the corresponding values obtained using TaiESM1 (Figures 2a–2f) and the performance of TaiESM1 compared with that of other CMIP6 models (Figure 2g). The multimodel ensemble (MME) means of rainfall and surface temperature denote the averages for all CMIP6 models listed in Table 1. The CMIP6 MME overestimates rainfall over the tropical oceans (Figure 2b), especially on the southern side of the equator, but underestimates rainfall over tropical lands such as the Amazon basin and Indian subcontinent, indicating the long-standing rainfall bias in previous CMIP simulations (Stephens et al., 2010). TaiESM1 demonstrates rainfall biases similar to those of other CMIP6 models in the overestimation of warm pool rainfall and the Intertropical Convergence Zone (ITCZ) over the eastern Indian and eastern Pacific oceans (Figure 2c). For near-surface temperature, TaiESM1 demonstrates a warm bias over the southern oceans and the west Eurasian continent but a colder bias in both the Arctic and Antarctica (Figures 2e and 2f).

Figure 2. (a–f) Annual mean rainfall and surface temperature in observational analysis, and corresponding model biases of TaiESM1 and CMIP6 multiple model means with respect to observations. Observations used here are listed in Table 1. (g) Model ranking of basic mean variables compared with observations, following evaluation of the IPCC AR5 report. TaiESM1 is denoted by red crosses, and the CMIP6 multiple model mean is denoted by orange pluses. The air temperature is abbreviated as TA, zonal winds as UA, meridional winds as VA, geopotential height as ZG, longwave outgoing radiation as RLUT, shortwave upward radiation as RSUT, surface clear-sky shortwave upward radiation as RSUTCS, and surface clear-sky longwave upward radiation as RLUTCS.

The model performance rankings, shown in Figure 2g, are evaluated using the metrics introduced by Gleckler et al. (2008) through comparison of the relative performance with the median-performance model member among the CMIP6 models. The normalized space-time root-mean-square-error (RMSE) of selected variables, including air temperature, zonal and meridional wind velocity, and geopotential height at various pressure levels, and the SAT are evaluated against MRE2. Precipitation is evaluated against the GPCP, and radiation fluxes such as in total OLR, clear-sky upward longwave radiation, upward shortwave radiation in the total sky, and clear-sky shortwave radiation are evaluated on the basis of CERES-EBAF. The same evaluation is conducted for all CMIP6 models and the MME. Because the RMSEs of all models are compared with the median-performance model, high-performing models are shown in the lower part of the figure. Overall, TaiESM1 is among the top 50% of all of the CMIP6 models for performance in terms of the evaluated mean variables. Especially, it is in the top group for simulating tropospheric winds and temperature, except for the 200-hPa temperature, among the models.

Seasonal Evolution and Monsoons Seasonal Means

We evaluated the historical simulation (1976–2014) of tropical precipitation and circulation, and focuses on the global monsoon, which refers to the dominant mode of annual variation over the tropical region, coupled with large-scale overturning (Trenberth et al., 2000; Wang & Ding, 2008). The seasonal mean precipitation and circulation in the upper and lower troposphere during June-August (JJA) and December-February (DJF) are presented in Figure 3. TaiESM1 realistically simulates the major characteristics of seasonal mean fields. In the boreal summer (JJA, Figures 3a–3c), major precipitation occurs over the ITCZ and the monsoonal regions of West Africa, South Asia, East Asia-Western North Pacific (EAWNP), and tropical America. As shown in Figure 3b, the monsoonal precipitation in South Asia and EAWNP, the associated planetary divergent flow, and regional monsoon circulation are well simulated. In the boreal winter (DJF, Figures 3d–3f), major precipitation occurs in the Southern Hemisphere, which corresponds to the monsoons in Africa, the Maritime Continent, Australia, and South America. TaiESM1 also well presents the monsoonal winds in the lower troposphere (vectors in Figure 3) and velocity potential in the upper troposphere (contours in Figure 3). It successfully characterizes the southwesterly flow in West Africa, the Arabian Sea, and the South China Sea and the southeasterly flow in North America in summer, and the northeasterly flow in West Africa, the Arabian Sea, the South China Sea, and South America in winter. The upper tropospheric divergence center in winter is located over the Maritime Continent, whereas that in summer is located over South Asia and EAWNP.

Figure 3. Precipitation (mm day−1; shading), 200-hPa velocity potential (106 m2 s−1; contours), and 850-hPa wind (m s−1; vectors) from the observation results (GPCP in 1997–2014; MRE2 in 1979–2005), TaiESM1 outputs (1979–2005), and TaiESM1 minus the observation results in (a–c) JJA and (d–f) DJF. Vectors in (a, b, d, and e) denote wind speeds higher than 3 m s−1, and those in (c, f) denote a wind speed difference greater than 1 m s−1. In (c, f), precipitation is shown when the differences between model outputs and observations have a confidence level of 99%.

Model biases in seasonal climatology can be quantitatively demonstrated through direct comparisons of model outputs with satellite observation and reanalysis data, as shown in Figures 3c and 3f. Differences at a confidence level of 99% are shown. TaiESM1 overestimates the summer precipitation over the equatorial Indian Ocean and central Pacific but underestimates the precipitation south of the Tibetan Plateau and the northern Bay of Bengal (Figure 3c). By contrast, the precipitation in extratropical East Asia is reasonably simulated. The biases are associated with excessively strong divergence in the central tropical and southeastern Pacific and a convergence bias in the tropical Atlantic and Western Africa. This planetary-scale divergence feature indicates the potential global impact of regional bias. As a result, a double ITCZ feature in the western/central tropical Pacific and the easterly winds are stronger than that observed in the tropical eastern Pacific. In winter (Figure 3f), the model overestimates the off-equatorial tropical precipitation with the double ITCZ feature in the eastern Pacific, whereas the precipitation in storm tracks is well simulated.

The ranking of the performance of the CMIP6 models based on the model-observation pattern correlations of surface temperature (TS), OLR (RLUT), and precipitation (PR) over North America, Africa, Asia, Australia, and the global domain is shown in Figure 4. Range of these regions is denoted in Figure S1. In general, MME means have the optimal performance, while the individual models have more difficulty in representing the warm-season rainbands than representing the annual rainfall means (Figures 4a–4c). Among the three variables, rainfall is the variable that the models exhibit the least ability to simulate, in all regions and both cold and warm seasons. Notably, almost all models have good capability to simulate surface temperature in all monsoon regions, especially during the cold season. The reasons underlying the biases in particular regions remain unclear and warrant further investigation. Overall, the ability of TaiESM1 to represent all three variables in all four regions is usually above the median member of CMIP6 and sometimes pretty close to the MME.

Figure 4. (a–o) Model ranking of pattern correlation annually (top) and in JJA (middle) and DJF (bottom) for (a–c) across the globe and four monsoon regions, including (d–f) Asia, (g–i) Africa, (j–l) Australia, and (m–o) North America, among surface temperature (TS), outgoing longwave radiation (RLUT), and precipitation (PR). The multiple model ensemble (MME) means of CMIP6 models are presented in orange, those of individual models are presented in gray, and that of TaiESM1 is presented in red. The four monsoon regions are annotated in Figure S1.

Seasonal Evolution of Monsoonal Precipitation

The seasonal evolution of precipitation in the monsoon regions is presented in Figure 5: West Africa (20°W–30°E), South Asia (50°E−100°E), EAWNP (110°E−140°E), and Central America (110°W–50°W). As shown in Figures 5a–5c, TaiESM1 efficiently captures the seasonal evolution of the African monsoon (e.g., the northward advancement in spring, southward retreat in autumn, and peak in August), but with excessive tropical precipitation in April-June and November-December and slightly weaker precipitation in the peak period. However, the simulation of seasonal evolution in South Asia is less well simulated (Figures 5d–5f). Instead of simulating the northward advancement in spring and the southward retreat in autumn, the model simulates a relatively stationary precipitation pattern between 0°N and 20°N during June-August and a strong stationary rainband over the southern slope of the Himalayas and the Tibetan Plateau (25°N–30°N). These biases are consistent with the excessive precipitation south of the Indian subcontinent and the Tibetan Plateau displayed in Figure 3c; however, the reason for the biases is not well understood. One possibility is that the excessive precipitation south of the Tibetan Plateau, likely due to the oversimulated topographic effect in the region, induces anomalous subsidence over South Asia and prevents the development and northward advancement of precipitation. Another reason can be the undersimulated forcing of the Arakan Mountains in western Myanmar. Wu et al. (2014) reported that the deficiency in resolving the narrow north-south-elongated Arakan Mountains could lead to poor simulation of monsoon onset in the Bay of Bengal and underestimation of precipitation in the northeastern corner of the Bay of Bengal. The degree to which the topographic factors contribute to model biases requires further investigation. For the EAWNP (110°E−140°E; Figures 5g–5i), achieving a realistic simulation of the asymmetric seasonal variation (e.g., strong/fast northward advancement and weak/slow southward retreat) by using climate models is often challenging (Chen, Hsu, et al., 2019; Kusunoki & Arakawa, 2015). Each year, beginning in March, the EAWNP undergoes a series of transitions in precipitation (10°N–40°N in Figure 5g). The East Asian spring rain over subtropical East Asia (20°N–30°N) is followed by Mei-yu/Baiu and its northward migration in May-June. The termination of Mei-yu/Baiu in late July coincides with the onset of the WNP summer monsoon and typhoon season when the monsoon trough is established over the Philippine Sea and the subtropical anticyclonic ridge shifts suddenly northward. In September, the WNP monsoon begins a southward retreat (Chou et al., 2011; LinHo et al., 2008; Murakami & Matsumoto, 1994; Suzuki & Hoskins, 2009; Wu et al., 2009, 2018). The northward advancements in the stage-wise development of precipitation (Figure 5g), which is often a challenge for climate models to simulate, is well simulated in TaiESM1 (Figures 5g–5i). However, TaiESM1 still unrealistically simulates the split rainbands off of the equator during autumn and winter, which are associated with the double ITCZ model bias. The model-simulated American monsoon (110°W–50°W) features presented in Figures 5j–5l indicate that the seasonal evolution is reasonably simulated by TaiESM1, with excessive precipitation during the evolution and dryness outside the precipitation region.

Figure 5. Latitude-time cross section of precipitation (mm day−1, shading) obtained from observations from the GPCP and TaiESM1 and the differences between averages over regions (a–c) 20°W–30°E (Africa), (d–f) 50°E−100°E (South Asia), (g–i) 110°E−140°E (East Asia), and (j–l) 110°W–50°W (Central America) during 1998–2014. (m) Model ranking based on skill scores averaged over four monsoon regions with data simulated by TaiESM1 (red cross), CMIP6 models (gray numbers), and the multimodel ensemble mean (orange plus).

Figure 5m presents the skills of the CMIP6 models to simulate the annual cycles (as shown in Figures 5a, 5d, 5g and 5j) of four major monsoon regions during 1998–2014. The CMIP6 models (denoted by gray numbers) generally have good skill scores (i.e., >0.7) in the simulation of the seasonal evolution of precipitation in Africa, South Asia, East Asia, and Central America monsoonal regions. The CMIP6 ensemble mean (denoted by an orange plus) scores much higher (i.e., >0.9) than most models in all regions. TaiESM1 (denoted by a red cross) performs better than most CMIP6 models, especially over Africa and East Asia, with scores between 0.85 and 0.95.

Intraseasonal Variability, Synoptic Variability, and Extremes Intraseasonal Variability

The Madden-Julian Oscillation (MJO) is the dominant pattern of atmospheric intraseasonal (e.g., 20–100 days) variability in the tropics (Lau & Waliser, 2005; Madden & Julian, 1972; Zhang, 2005). MJO events are characterized by large-scale tropical circulation anomalies that develop over the Indian Ocean and propagate eastward into the western Pacific in 2–3 weeks. The summertime intraseasonal oscillation (ISO) is an important component of the Asian summer monsoon, which involves the movement of convection centers both northward and eastward in the equatorial and northern Indian Ocean and north-northwestward in the WNP (Hsu, 2005; Hsu & Weng, 2001; Lau & Chan, 1988). MJO events have major global impacts on monsoons, tropical storms, extratropical weather, and the ENSO. However, realistically representing the MJO by using the current climate models remains difficult (Hung et al., 2013; Kim et al., 2009, 2020). We evaluate the ability of TaiESM1 to simulate the MJO. The CLIVAR MJO Working Group diagnostics package is used to isolate and analyze intraseasonal variability (CLIVAR Madden-Julian Oscillation Working Group, 2009). Here, two seasons are defined: boreal winter (November to April) and boreal summer (May to October).

The wavenumber-frequency spectra of 850-hPa zonal wind averaged over 10°S–10°N simulated using TaiESM1 are compared with the observation spectra in Figure 6a. TaiESM1 simulates the observed wavenumber-1 structure with a much broader periodicity band and a maximum in the longer period (∼80 days compared with the observed 30–80 days) during boreal winter (Figure 6a, upper). This low-frequency tendency is reflected by the weaker and slower eastward propagation in the time-longitude Hovmöller diagrams (Figure 6b). Meridional propagation is one of the major characteristics of ISO. With very weak coupling between winds and convection, TaiESM1 only simulates the observed northward propagation of the MJO over the northern part of the Indian Ocean without rooting the ISO from the deep tropics. It also undersimulates the southward propagation tendency south of the equator (Figure 6c). The overall ISO performance of CMIP6 models, evaluated on the basis of the indices for the boreal winter and summer, is summarized in Figure 6d. In general, TaiESM1 tends to better simulate the overall amplitude of the intraseasonal variability but fairly simulates the propagation tendency (e.g., eastward/westward component ratio, W1, and northward propagation over the near equatorial region, NEQ S2) and periodicity of MJO.

Figure 6. (a) Zonal wavenumber-frequency spectra for equatorial 850-hPa zonal wind between 10°S and 10°N. Lag-longitude diagrams of intraseasonal rainfall (color) and 10-m zonal wind (contours) averaged over (b) 10°S–10°N correlated with Indian Ocean (10°S–5°N, 75°E−100°E) precipitation and (c) 80°E−100°E correlated with near equatorial region (5°N–10°N, 85°E−90°E) precipitation. (d) Summary diagram. Skill scores and the ratio of the 14 GCMs indicating their fidelity in representing the characteristics of MJO simulations. Each item is described in detail in the Text S1. The left Y-axis shows units of W1 and W2, while the right Y-axis shows units of W3 and S1–S3. The cross denotes TaiESM1, and number marks denote CMIP6 models based on Table 1.

The efficiency of TaiESM1 in simulating the characteristics of convectively coupled equatorial waves over the tropical belt (30°S–30°N) are assessed (Kiladis et al., 2009; Kim et al., 2009; Takayabu, 1994; Wheeler & Kiladis, 1999). Figures 7a and 7b show the space-time spectra of the symmetric component of equatorial precipitation following Wheeler and Kiladis (1999). The observation is characterized by strong variance associated with the MJO, equatorial Kelvin and Rossby waves, and slightly weaker mixed Rossby-gravity waves. The simulation by TaiESM1 realistically represents these three main features but with weaker amplitudes. The Rossby waves and the high-frequency/high-wavenumber Kelvin waves are particularly weak. By contrast, the mixed Rossby-gravity waves are not realistically simulated. Following Dias and Kiladis (2014), we evaluated the model performance by season, and the results for the Kelvin and Rossby waves are presented in Figure 7c. A comparison of the annual eastward and westward wave spectra reveals that TaiESM1 exhibits higher skills than the other models, especially for westward propagation. In general, all models exhibit lower simulation skills in JJA and September-November (SON) and higher skills in March-May (MAM) and DJF. Notably, TaiESM1 exhibits higher skills for representing equatorial Rossby waves in DJF and MAM (i.e., index 07 and 10) compared with other models. This result is consistent with the high efficiency in simulating westward propagation (e.g., index 02).

Figure 7. Symmetric coherence-squared in wavenumber-frequency spectrum between 15°S–15°N of precipitation in (a) observations and (b) TaiESM1 simulations. (c) Summary diagram. Skill scores and the ratio of the 14 GCMs indicating their fidelity in representing the characteristics of convective coupled equatorial waves in simulations. A higher value indicates better simulation ability. Each item is explained in detail in the Text S2. The cross indicates TaiESM1 and numbers denote other CMIP6 models.

Synoptic Eddy Variability

Synoptic eddy activity causes daily temperature and precipitation fluctuation in the extratropics. The two-way energy conversion between the mean state and synoptic perturbations (i.e., the eddy-mean flow interaction) is a key physical process that keeps the balance between mean flows and synoptic eddies and helps maintain the atmospheric general circulation in the extratropics. For example, a synoptic eddy grows at the expense of the mean potential energy in the early stages of the lifecycle and feeds back kinetic energy to the mean flow in the later stages. A climate model that reasonably captures the aforementioned dynamic process and synoptic eddy activity is likely to more realistically simulate a mean state. Therefore, the model ability to simulate synoptic activity must be evaluated. In this study, synoptic perturbations are defined as 1–10-days band-pass-filtered fields such as wind and temperature.

Synoptic meridional momentum flux at 250 hPa and meridional heat flux at 850 hPa are important variables in kinetic and potential energy conversion, respectively, between the mean state and synoptic eddies and are often adopted as proxies to represent eddy activity. The observed 250-hPa meridional momentum flux in the Northern Hemisphere (Figure 8a) is maximized in the central North Pacific and North Atlantic around 30°N with an eastward extension to North America and Europe, respectively. In the Southern Hemisphere, the maxima appears in the Southern Atlantic and Southern Indian Ocean around 40°S. The observed meridional heat flux at 850 hPa (not shown), which is often spatially and temporarily associated with momentum flux, is found mostly located to the westward and poleward side of the meridional momentum flux. The magnitudes and spatial distribution of these major features are realistically simulated by TaiESM1 (Figure 8c). A similar comparison for June-August also reveals the high performance of TaiESM1 in simulating the overall spatial distribution of eddy fluxes (not shown).

Figure 8. Performance of TaiESM1 in the simulation of synoptic eddy variability. (a, b) Observed and (c, d) simulated synoptic eddy momentum flux at 250 hPa during 1980/1981–2013/2014. DJF (a, c) climatology (m s−1) and (b, d) interannual variance (m2 s−2). Pattern correlation (e) and normalized root-mean-square-error (RMSE) (f) for synoptic eddy fluxes of 12 CMIP6 models (see Table S2). The seasonal mean climatology and interannual variance of the eddy momentum flux at 250 hPa (uv250) and eddy heat flux at 850 hPa (vt850) in the Northern Hemisphere storm track (NH, 15°N–75°N) and Southern Hemisphere storm track (SH, 30°S–60°S) during DJF and JJA are evaluated. Red crosses and gray numbers represent TaiESM1 and other CMIP6 models, respectively. The medians of RMSEs are subtracted from the RMSEs, and the difference is then divided by the median. Smaller normalized RMSE values indicate higher performance.

Despite fluctuating in a time scale of <10 days, the overall activity of synoptic eddies varies with a large-scale background environment such as the location and strength of jet streams and temperature gradient, which are strongly affected by known interannual fluctuations such as the ENSO. The interannual variance in meridional momentum for the observation and simulation is presented in Figures 8b and 8d. The interannual variance is defined as the variance of seasonal mean fluxes during 1975–2005. In the observation results, two major active regions of meridional momentum flux in the Northern Hemisphere are the eastern North Pacific and central North Atlantic (Figure 8b). The maxima of the meridional heat flux are found in the west of the maxima in the meridional momentum flux and in the east of Greenland with an eastward extension into the Eurasian Arctic. TaiESM1 reasonably simulates these features with weaker variance in both the North Pacific and Atlantic (Figure 8d). The variance of both fluxes in the Southern Hemisphere is also reasonably simulated but less skillful so compared with the simulation of the Northern Hemisphere. The spatial distribution during JJA is also realistically simulated (not shown). Overall, TaiESM1 can realistically simulate the spatial distribution and temporal (seasonal and interannual) fluctuation in synoptic eddy activity, thus demonstrating above average performance among the CMIP6 models (Figures 8e and 8f).

An interesting contrast in model performance is the overall lower pattern correlation (0.6–0.07) for interannual variance in the Southern Hemisphere compared with the counterpart in the Northern Hemisphere (pattern correlation, ∼0.8). The contrast in model performance reveals that even under the same forcing as the ENSO, the synoptic eddy activity in the Northern Hemisphere is easier to simulate than that in the Southern Hemisphere. The strong control of the significant land-sea contrast and topography in the North Hemisphere, which is absent in the Southern Hemisphere, is likely one of the major reasons. Large-scale orography (e.g., the Himalayas and the Rocky Mountains) and land-sea contrast have dominant importance in forcing stationary waves (Brayshaw et al., 2009; Egger, 1976; Held et al., 2002; Hoskins & Karoly, 1981), which influence the location of the jet stream and synoptic eddies.

Extreme Precipitation Events

We evaluate the performance of TaiESM1 by examining indices associated with extreme precipitation and compare the results with those of CMIP6 models. The skills of TaiESM1 and other models are evaluated against the precipitation indices derived from the 1-degree grid of the GPCP. The study region is the 40°S–40°N tropical belt during 1998–2014, considering the common data period of the GPCP and model outputs.

The indices of simple daily intensity (SDII), maximum 1-day precipitation (RX1day), maximum 5-days precipitation (RX5day), extreme precipitation intensity (PR99), total rainfall occurrence (Totfq, defined as daily precipitation exceeding 1 mm), and consecutive dry days (CDD) are analyzed to examine the representation of precipitation characteristics associated with extreme events.

To demonstrate the performance of TaiESM1 in simulating extreme indices, Figures 9a–9d display the spatial pattern of RX1day in the GPCP and TaiESM1 as an example. Compared with the GPCP, TaiESM1 generally captures the seasonal main feature of RX1day in JJA (Figures 9a and 9c) and DJF (Figures 9b and 9d), despite some degree of overestimation or underestimation over the tropical region (e.g., overestimation in the EAWNP, eastern tropical Pacific, and southern Pacific regions and underestimation over Central Africa, Central America, and South America in JJA). Compared with the other CMIP6 models, TaiESM1 exhibits higher skill scores in the simulation of SDII, RX1day, RX5day, and PR99 in JJA and DJF (red cross in Figures 9e and 9f), similar to the performance in the CMIP6 ensemble. However, the scores associated with rainfall occurrence such as CDD and Totfq are relatively low. The models also tend to obtain lower scores when a lower threshold is selected for defining a wet day is selected (0.1 mm day⁻¹), indicating that the biases associated with too-frequent precipitation in previous model simulations (e.g., in CMIP Phase 5) still exist in the sixth-generation CMIP models.

Figure 9. Spatial distribution of RX1day in (a) JJA and (b) DJF in the GPCP. (c and d) Same as (a) and (b), but for TaiESM1. The purple (gray) dotted region denotes the overestimated (underestimated) precipitation against the GPCP. (e and f) Model ranking based on the skill scores in simulating the extreme indices in JJA and DJF. CDD(0.1) and Totfq(0.1) represent the indices estimated using the wet-day definition of 0.1 mm day−1.

We further examine the performance of TaiESM1 in simulating the extreme precipitation over the EAWNP region (115°E−135°E, 20°N–50°N). In the EAWNP region, the main wet season between the 28th and 54th pentad (i.e., May 16 to September 27) is divided into two rainy periods: the first wet season (1st_wet) is associated with the Mei-yu frontal system, and the second wet season (2nd_wet) is related to typhoon rainfall (Chen & Chen, 2003; Chen, Hsu, et al., 2019; Chou et al., 2009; Hsu et al., 2014; LinHo & Wang, 2002). The current General Circulation Models (GCM) are less skillful in simulating precipitation intensity during the Mei-yu season (Chen, Hsu, et al., 2019; Endo & Kitoh, 2016; Kusunoki, 2018; Kusunoki & Arakawa, 2015) and tropical cyclone activities (Flato et al., 2013; Murakami et al., 2012a, 2012b), during which extreme precipitation is observed. (Chen et al., 2021) compared the performance of CMIP5 and CMIP6 models in precipitation simulation in seasonal evolution and extreme indices in the EAWNP region. CMIP6 models generally have higher skill scores in the EAWNP region, which indicate their improvement over the CMIP5 models. As shown in Figures 10a–10d, despite some degree of inconsistency, TaiESM1 can capture the spatial patterns associated with the Mei-yu rainband in the 1st_wet season and typhoon-related precipitation in the 2nd_wet season, which are similar to those obtained in the GPCP. TaiESM1 performs well in simulating the extreme indices in the EAWNP wet seasons (i.e., skill score >0.6; Figures 10e−10f) and exhibits higher scores than most other CMIP6 models, especially in SDII, RX1day, RX5day, and PR99. These results indicate that the improvement in TaiESM1 might be mainly associated with the improved precipitation intensity rather than rainfall frequency. Model simulations in the EAWNP region present errors associated with rainfall occurrence [i.e., CDD(0.1) and Totfq(0.1)], with more errors in the 1st_wet season than in the 2nd_wet season.

Figure 10. Spatial distribution of RX1day in the WNP-EA region (90°E−180°E, 0°N–40°N) in the (e) 1st_wet and (f) 2nd_wet seasons.

Diurnal Rainfall Phase and Amplitude

The diurnal cycle denotes the prominent oscillation of a climate system forced by the diurnal variation of solar radiation. We conduct evaluations of the diurnal peak phase, which is the local time when diurnal rainfall peak occurs, and the diurnal amplitude of the diurnal rainfall cycle by using TaiESM1 against the 3-hourly Tropical Rainfall Measurement Mission Multisatellite Precipitation Analysis (TMPA; Huffman et al., 2007) 3B42 observations. The evaluation is based on the first harmonic of the climatological diurnal rainfall cycle retrieved from the total data. We use data from 1998 to 2010 for TMPA and 30-years data for the TaiESM1 historical run. The phase and amplitude of diurnal rainfall in TaiESM1 share similar biases as seen in other Atmosphere-Ocean GCMs (Figures 11a and 11b), including the underestimation of diurnal amplitude over the tropical lands and the early peaking time over most of the land region (Covey et al., 2016; Dai, 2006). However, unlike other models, TaiESM1 exhibits improved representations of the propagation behavior in the diurnal peak phase over many topographical regions and coastal regions, including the Southern Great Plains and the coastal regions of the maritime continent. This improvement result is also reported by Lee, Wang, et al. (2020) and can be attributed to the improvements in the convective trigger function designs in TaiESM1 (Wang & Hsu, 2019). Figure 11c presents the performance of the CMIP6 model compared with TMPA observations based on the pattern correlation for the diurnal rainfall phase and amplitude. To compute the pattern correlation for the diurnal phase, the local phase is weighed with the local amplitude to emphasize on the regions with stronger diurnal signals. The pattern correlation coefficients for amplitude and phase between TaiESM1 and TMPA observations are 0.68 and 0.69, respectively (Figure 11c). The performance of TaiESM1 in the amplitude simulation is above average compared with the other CMIP6 models (correlation, 0.5–0.8); however, TaiESM1 is demonstrated to be one of the best models for the simulation of diurnal phase distribution. The higher performance in the phase simulation can be attributed to the more efficient simulation in the rainfall propagation regions.

Figure 11. (a) Diurnal amplitude and (b) peak phase of TRMM 3B42 multisatellite products and historical run of TaiESM1. TRMM is regridded onto the 0.9° × 1.25° grid of TaiESM1 for comparison. (c) Model ranking based on pattern correlations of phase and amplitude. The pattern correlation of phase is only considered for regions with diurnal amplitudes greater than 1.2 mm day−1, which approximately equates the global mean of diurnal amplitude. The results of TaiESM1 are denoted by red crosses, medians are denoted by orange pluses, and the results of the CMIP6 models are numbered based on Table 1.

Interannual-Interdecadal Variability El Niño-Southern Oscillation

The ENSO is one of the most prominent phenomena that contribute significantly to interannual variability. Like the MJO, realistic simulation of the ENSO by using climate models is challenging (Bellenger et al., 2014; Beobide-Arsuaga et al., 2021; Brown et al., 2020). To construct model-simulated ENSO composite, we identify the strong ENSO events simulated with the Nino3.4 index larger than 1.5 standard deviation of the entire period. Five ENSO events are then identified from the TaiESM1 simulation and used for the composite. Figure 12 presents the spatial distributions of surface temperature, mean sea-level pressure, and 1,000-hPa wind anomalies in December-February when the ENSO is in the mature stage; they are based on the composites of five El Niño events by TaiESM1 (Figure 12a) and their differences from the composites of six events in the MRE2 data set (Figure 12b). The TaiESM1-simulated SST anomaly (SSTA) is evidently larger in terms of both amplitude and covered area than the observed values, and the maximum shifts westward to the central equatorial Pacific compared with observations (Figure 12a), which is a common bias seen in many climate models (Bellenger et al., 2014). The horseshoe-like negative SSTA in the northwest/southwest and west of the positive SSTA correspondingly shift westward compared with observations and exhibit a cold bias in the far western tropical Pacific (Figure 12b). This overestimated SSTA structure leads to marked biases in the simulated atmospheric circulation and temperature. The biased strong SSTA induces a stronger-than-observation near-surface convergence toward the central equatorial Pacific. In response to the westward shift of positive SSTA in the equator, the circulation and temperature anomaly patterns shift westward. For example, the westward-shifted western Pacific anticyclonic anomaly and the observed warm-cool SSTA dipole in the subtropical WNP during El Niño induced by the local atmosphere-ocean interaction (Wang et al., 2000) is not efficiently simulated. The region is dominated by negative SSTA, with the positive SSTA restricted over the coastal East Asia. The cool-warm structure in extratropical Asia seen in observations shifts westward toward the interior of the Asian continent in the simulation. The warming in the Indian Ocean, which often occurs following the onset of El Niño, also occurs to the west of the observed location together with the westward shift of the anticyclonic circulation anomaly in the Southern Indian Ocean. The observed northwesterly anomaly to the west of Australia is seen as a southwesterly anomaly in the simulation, which evidently induces upwelling and cooler temperatures along the Australian west coast. By contrast, the temperature and circulation over North America and the North Atlantic are more realistically simulated; however, the Bering Sea is warmer and the northwestern North America is cooler than the observed temperatures.

Figure 12. Spatial composites of surface temperature (shading), mean sea-level pressure (contours), and 1,000-hPa wind anomalies (green arrows) of ENSO for (a) TaiESM1 and (b) differences between TaiESM1 and MRE2 reanalysis ensemble based on the December, January, February after ENSO events (DJF+1; year + 1). (c) Spectrum of Nino 3.4 from TaiESM1 (red), observations (black), and a range of CMIP6 models (25%: green, median: blue lines, 75%: brown dashed lines, full range: gray shading). (d) Normalized skill scores of geological patterns obtained from TaiESM1, multimodel ensemble, and CMIP6 models (see Table 1 for corresponding model numbers) compared with the MRE2 ensemble in the winter season. They are normalized on the basis of the median member of all CMIP6 models. The x-axis shows from June, July, August in year 0 (JJA0) to March, April, May in year + 1 (MAM+1) of ENSO events.

Figure 12c presents the observed and simulated spectra of the Niño 3.4 index. The observed Niño 3.4 index exhibits three statistically significant peaks between 2 and 8 years. TaiESM1 simulates a strong spectral peak at ∼4–5 years and another strong peak at ∼8 years. El Niño simulated by TaiESM1 has a larger amplitude than that simulated by most of the CMIP6 models (e.g., exceeding the 75th percentile of CMIP6 models, dashed line in Figure 12c), which likely leads to the large extratropical responses displayed in Figure 12a. The spectra of all CMIP6 models indicates a wide range of ENSO amplitudes among CMIP6 models (gray shading in Figure 12c). Figure 12d presents the skill scores of the models in simulating the observed El Niño SST in the tropical Pacific for four seasons preceding and following the peak of El Niño. The performance of TaiESM1 is among the best, except in the winter of year + 1 (DJF⁺¹) when El Niño is in its peak. The CMIP6 models tend to perform poorly in the spring of year + 1 (MAM⁺¹) following the peak phase of El Niño.

Atlantic Multidecadal Oscillation and Pacific Decadal Oscillation

Oceanic interdecadal-multidecadal fluctuations are considered the major reason for climate variability beyond interannual time scales. Two well-known oscillations beyond the decadal time scale, the AMO and PDO, are identified for further evaluation of TaiESM1.

A comparison between observed PDO and PDO simulated by TaiESM1 is presented in Figures 13a and 13b. The overall simulated pattern is similar to that of the observation: a negative SSTA in the extratropical North Pacific, a horseshoe-like positive SSTA in the extratropical eastern Pacific, and a positive SSTA in the equatorial central-eastern Pacific. The weaker SSTA structure in the other oceans is also simulated. However, consistent with the bias in the El Niño simulation, the positive and negative SSTA in the equatorial Pacific shift to the west of the observed pattern. The Taylor diagram presented in Figure 13c reveals that all CMIP6 models exhibited a pattern correlation between 0.8 and 0.9 and large disagreements in terms of mode variability. TaiESM1 exhibits a correlation of 0.9 and a variability ratio of ∼1.25, which is at the upper end of that of the CMIP6 models. This bias is consistent with the strong ENSO signal simulated by TaiESM1.

Figure 13. (a–c) Geospatial patterns of Pacific Decadal Oscillation (PDO) derived from HadISST and TaiESM1 and Taylor diagrams between the two modes. (d–f) Similar comparison as (a–c) for the Atlantic Multidecadal Oscillation (AMO) between HadISST and TaiESM1. Model numbers on Taylor diagrams are based on Table 1 and TaiESM is marked as red cross.

Figures 13d–13f display the observed and simulated AMO, which significantly affects the regional and global climate. The major centers of SSTA in the extratropical North Atlantic, tropical Atlantic, central/eastern equatorial Pacific, and extratropical Pacific are realistically simulated by TaiESM1 (Figures 13d and 13e). Figure 13f shows that the performance of TaiESM1 in AMO simulation is average, with a correlation of 0.7, compared with a correlation range of 0.6–0.85 among the CMIP6 models.

Atmospheric Teleconnection

Atmospheric teleconnections are important phenomena that link climate variation in separate remote regions, thereby influencing regional climates. Therefore, climate models should be able to accurately simulate the main characteristics of teleconnection patterns to produce a reasonable global climate variability distribution. The performance of TaiESM1 in simulating well-known teleconnection patterns is evaluated. We adopt the CBF method to extract climate modes of models to avoid problems such as mode swapping. This is important for regional modes, such as the Pacific-Japan (PJ) pattern, which appears as the leading EOF in MRE2 precipitation (110°E−180°E, 5°N–55°N; similar structure in the GPCP, not shown), with a tripolar structure, but as the second EOF in TaiESM1.

Figure 14 presents the spatial structures of the Arctic Oscillation (AO; Thompson & Wallace, 1998) in winter and the PJ pattern (Hsu & Lin, 2007; Kosaka & Nakamura, 2006; Nitta, 1987) in summer as examples. Observations and simulations of the AO, which is defined as the first EOF of seasonal mean sea-level pressure north of 20°N, are shown in Figures 14a and 14b, respectively. The out-of-phase relationship between the polar region and middle latitudes and the location of major centers is reasonably simulated by TaiESM1. However, the Atlantic component of the meridional dipole is weakly simulated, whereas the simulation of the Pacific component is stronger than the observation. The explained variance in simulation is ∼33.4%, which is very close to the observed 33.5%. For the PJ pattern, TaiESM1 realistically simulates the out-of-phase pattern between north and south parts of the East Asia region seen in the observations (Figures 14c and 14d). Although the locations of the maximum are captured, TaiESM1 has a stronger simulated amplitude and a higher explained variance of 25.6% compared with the 24.6% of variance explained by observations.

Figure 14. (a, b) Geospatial patterns of the Arctic Oscillation (AO) based on the first empirical orthogonal function (EOF) of sea-level pressure of regions north of 20°N. (c, d) Geospatial patterns of Pacific-Japan (PJ) teleconnection between MRE2 and TaiESM1 based on vorticity over the East Asia.

Figure 15 presents the Taylor diagrams of four important teleconnections, namely the AO, PJ, and PNA (Wallace & Gutzler, 1981) patterns, and the North Atlantic Oscillation (NAO; Rogers & vanLoon, 1979; vanLoon & Rogers, 1978; Walker & Bliss, 1932) obtained by all CMIP6 models and TaiESM1, on the basis of the CBF method. The Taylor diagram for AO shows that TaiESM1 realistically simulates the AO (red cross) like most CMIP6 models do, with a pattern correlation of 0.9 (Figure 15a). As displayed in Figure 15b, the simulated PJ pattern in TaiESM1 has a pattern correlation of 0.85 with the observed pattern, whereas the CMIP6 models have correlations of 0.8–0.9. The PNA and NAO are two dominant atmospheric intrinsic modes of seasonal and interannual variability and influence interdecadal and longer time scales through interactions with oceans (Battisti et al., 2019). For simulating the PNA pattern (Figure 15c), the TaiESM1 is the second best among the CMIP6 models, with a pattern correlation of 0.96 and an amplitude very close to the observation. Most models simulate the NAO pattern well, with a pattern correlation ranging between 0.85 and 0.9 (Figure 15d). TaiESM1 reasonably simulates the NAO, with a pattern correlation of 0.85. These results indicate the TaiESM1 can reproduce the overall pattern structure of observed teleconnection patterns, but still have regional biases such as those shown in the AO pattern. Notably, the CMIP6 models exhibit a similar ability to simulate the patterns of the major atmospheric modes but have a large spread of variability, which is shown by the model spread over the radial distance in Taylor diagrams. TaiESM1 is also one of the models with the best performance in terms of normalized RMSE and the ratio of variability relative to the observations.

Figure 15. Taylor diagrams for evaluating important geographical features of teleconnection obtained using CMIP6 models. (a) The Arctic Oscillation (AO) in DJF, (b) Pacific-Japan pattern in JJA, (c) Pacific-North American Oscillation (PNA), and (d) North Atlantic Oscillation (NAO) in DJF among TaiESM1 (red cross), CMIP6 models (numbers based on Table 1), and the MRE2 ensemble. Evaluations are based on regressed spatial patterns of sea-level pressure with the AO, PNA, and NAO and vorticity with PJ between TaiESM1 and MRE2.

Equilibrium Climate Sensitivity

We estimated the equilibrium climate sensitivity (ECS) to characterize the surface warming of TaiESM1 in response to increasing greenhouse gas emissions. Compared with CMIP5, CMIP6 models are reported to have an even larger range of ECS between 1.8°C and 5.6°C (Meehl et al., 2020). To better interpret future climate changes projected under various emission scenarios, the range of model sensitivity should be understood. We used the method proposed by Gregory et al. (2004), which contrasts the residual energy and surface temperature obtained from the 4 × CO₂ (years 1–150) and the preindustrial (years 501–650) experiments to evaluate the ECS of TaiESM1 (Figure 16a). Compared with the ECS calculated by Zelinka et al. (2020; their Table S1), the ECS of TaiESM1 is 4.32°C, which is higher than the average of 3.8°C among the CMIP6 models (Figure 16b). The CMIP6 models are considerably more sensitive than the CMIP5 models (Zelinka et al., 2020). The sensitivity of TaiESM1 is at the 59th percentile among the CMIP6 models.

Figure 16. (a) Scatter plot of changes in the top-of-atmosphere energy residual and surface temperatures in TaiESM1 between the 4 × CO2 (years 1–150) and preindustrial (years 501–650) experiments for estimating equilibrium climate sensitivity (ECS). ECS = −intercept/2 × slope of the regression line in (a), following Gregory et al. (2004). (b) Model spread for ECS among CMIP6 models. TaiESM1 is denoted by a red cross; CMIP6 models are denoted by the number marks (see Table 1), and their median members are denoted by an orange plus. The ECSs of the CMIP6 models in Zelinka et al. (2020) are plotted for comparison.

Discussion

To understand the performance differences of TaiESM1 and its base model CESM1, we evaluated the model performances with the historical run of CESM1 from the CMIP5 data archive. It is referred to as “CESM1/CMIP5” in the following discussion. Due to data availability, we can only conduct analysis for mean states, monsoon progression, extreme rainfall, ENSO, and teleconnections. Figure S2 shows the differences of mean states of rainfall and surface temperature between TaiESM1 and CESM1/CMIP5. The RMSE of CESM1/CMIP5 is 1.31 for annual rainfall and 2.9 for annual surface temperature, while those of TaiESM1 is 1.27 for rainfall and 2.7 for temperature. It is notable that TaiESM1 reduces double ITCZ bias and cold surface biases in CESM1/CMIP5, but the two biases still remain significant in TaiESM1 compared to the observations. Similar to the biases found in Figure 5, Figure S3 shows that the seasonal mean states and monsoon progression in CESM1/CMIP5 are very similar to TaiESM1, which indicates spurious equatorial rainfall and weaker northward monsoon progression. From Arica, South Asia, East Asia, and Central America, the skill scores of each region are 0.94,0.88, 0.87, 0.92 for CESM1/CMIP5 and 0.94, 0.87,0.87, 0.92 for TaiESM1.

Figure S4 shows the spatial distribution of RX1day rainfall in the two wet seasons over East Asia in CESM1/CMIP5. In the first wet season, CESM1/CMIP5 simulates very weak RX1day rainfall over the Indo-China peninsula, South China Sea (SCS), and the west Pacific to the east of Taiwan (Figure S4a). TaiESM1 has improved the rainfall intensity over the east coast of Myanmar and SCS, although the intensity is still not strong enough (Figure S4c). In the second wet season, CESM1/CMIP5 also has a weak bias of rainfall intensity over the mainland China and the east part of Indo-China peninsula (Figure S4b). Again, TaiESM1 has improved the extreme rainfall intensity in the East Asia region, especially over the southern China and Taiwan, as well as the east part of Indo-China peninsula (Figure S4d). Compared with GPCP the pattern correlations are 0.72 for CESM1/CMIP5 and 0.77 for TaiESM1 for the first wet season, and 0.69 and 0.77 for the second wet season, respectively.

Figure S5 shows the differences of ENSO composite between TaiESM1 and CESM1/CMIP5 and spectrum of Nino 3.4 of CESM1/CMIP5 similar to Figures 12b and 12c. For ENSO composite in DJF, TaiESM1 has improved surface temperature responses of ENSO seen in CESM1/CMIP5, especially for cold bias over North America and warm bias over Antarctica (Figure S5a). The skill scores of the model-simulated ENSO composite compared to MRE2 are 0.55 for CESM1/CMIP5 and 0.56 for TaiESM1 over the 1980 to 2005 period. Compared with the two strong peaks of TaiESM1 (Figure 12c), CESM1/CMIP5 has variability amplitude similar to the observed data set (Figure S5b). Examining the phase locking behavior of ENSO for both CMIP5 and CMIP6 models, Chen and Jin (2021) also shows both CESM1/CMIP5 and TaiESM1 have similar biased dynamics for ENSO phase transition, but TaiESM1 has stronger intensity. Sensitivity studies are being planned to pinpoint the causes of differences between the two models.

For teleconnection, TaiESM1 improves the tropical responses to both the PDO and AMO over the Pacific of CESM1/CMIP5 (Figure S6). The deviation ratio and pattern correlations of PDO are 0.93 and 0.83 for CESM1/CMIP5, and 1.16 and 0.88 for TaiESM1. Stronger tropical warm anomaly is also shown in TaiESM1 compared with CESM1/CMIP5 (Figures S6b and S6d). For AMO, they are 1.12 and 0.74 for CESM1/CMIP5, and 1.36 and 0.72 for TaiESM1. Figures S6a and S6c suggest TaiESM1 has too strong cooling over the tropical east Pacific compared with CESM1/CMIP5. For AO mode, the ratio and correlation are 1.05 and 0.89 for CESM1/CMIP5, and 0.90 and 0.88 for TaiESM1. The weaker anomaly over the northern Atlantic is the cause of the deteriorated performance in TaiESM1 (Figures S7a and S7c). For PJ mode, TaiESM1 reduces the amplitude biases of the oversimulated latitudinal gradient of vorticity in CESM1/CMIP5 over the East Asia (Figures S7b and S7d). The variability ratio and pattern correlation are 1.26 and 0.87 for CESM1/CMIP5, and 1.08 and 0.81 for TaiESM1.

While detailed sensitivity experiments are needed to understand the cause, we think the improvements in rainfall are very likely to be from better rainfall characteristics due to convective trigger scheme and the better surface radiative heating due to the inclusion of radiative effects of 3D topography over the East Asia region. Previous studies show these two schemes can improve intensity and distribution of land convection (Lee et al., 2013; Wang & Hsu, 2019), which is consistent with improvements of land rainfall in TaiESM1 (e.g., Figure S4). In terms of the climate variability of ENSO, decadal variability, and teleconnection, we think the improvements are likely due to changes of moisture-cloud relationship and cloud radiative properties from GTS macrophysics scheme and convection behaviors of convective trigger scheme, through atmosphere-ocean interactions.

It is also notable that, while we only use one member of the historical simulation of TaiESM1 in our analysis, the analysis presented here does not reflect the internal variability of climate systems in the climate models. To fully address the impacts of internal variability, large number of ensemble members with carefully designed initialization are needed. Projects like the CESM Large Ensemble Project (Kay et al., 2015) are needed for more complete evaluations of interannual and interdecadal variability of TaiESM1. For other evaluations with shorter time scales, including diurnal cycle, MJO, and extreme rainfall, the results are less sensitive to the internal variability of TaiESM1.

Summary

This study evaluated the performance of TaiESM1 in simulating the observed climate variability in the historical simulation driven by CMIP6 forcing. TaiESM1 is developed on the basis of CESM1.2.2, with modifications to the cumulus convection scheme and cloud fraction scheme, replacement of the aerosol scheme with a new one, and implementation of a unique scheme for three-dimensional surface absorption of shortwave radiation that resolves the effects of complex terrains on the surface radiation budget. Most model evaluations focus on the climatological mean, whereas TaiESM1 focuses on climate variability, including precipitation extremes, synoptic eddy activity, intraseasonal fluctuation, monsoon evolution, and interannual and multidecadal atmospheric and oceanic teleconnection patterns. The series of intercomparisons between the simulations of TaiESM1 and CMIP6 models and observations indicate that TaiESM1 is capable of realistically simulating the observed climate variability at various time scales and is favorable to the CMIP6 models in terms of many key climate features. Biases are also identified and discussed.

TaiESM1 is a participating climate model in the CMIP6 and has provided independent simulations and outputs for the Diagnostic, Evaluation, and Characterization of Klima experiments and model intercomparison projects (MIPs) such as the Scenario MIP, Global Monsoon MIP, Aerosol/Chemistry MIP, and Cloud Feedback MIP, which are available at the data portal (https://doi.org/10.22033/ESGF/CMIP6.9684). Through CMIP6 and preparation for the upcoming Sixth Climate Change Assessment report of the Intergovernmental Panel for Climate Change, TaiESM1 will join the efforts to expand our understandings of variability and changes of the earth system with the global scientific community.

Acknowledgments

We thank two anonymous reviewers for their constructive suggestions and comments to make this study more complete. This study is conducted under the Consortium for Climate Change Study and is supported by grants from the Ministry of Science and Technology, Taiwan (MOST 109-2123-M-001-004, MOST 105-2119-M-001-018, and MOST 110-2123-M-001-003). We extend special thanks to the National Center for High-performance Computing for their computational support for TaiESM1 development and participation in CMIP6 activity.

Data Availability Statement

All observational data sets used in this study are available online, and a complete list of download sites and literature is included in Table S2. All computed correlation coefficients and skills scores used in our figures are documented in supplementary materials, too. The data used to produce Figure 16b are based on the data provided in Supplementary Materials of Zelinka et al. (2020). All CMIP6 model data presented in this research, including for TaiESM1, can be downloaded from the data portal for CMIP6 hosted by Lawrence Livermore National Laboratory, Department of Energy (http://esgf-node.llnl.gov/search/cmip6/). A complete list of the CMIP6 model data used for analysis in this study can be found in Table S1. Historical simulations of CESM1/CMIP5 are downloaded from the CMIP5 archive maintained by the Centre for Environmental Data Analysis at UK (https://data.ceda.ac.uk/badc/cmip5/data/cmip5/output1/NSF-DOE-NCAR/CESM1-CAM5/historical/).