Taiwan Earth System Model Version 1: description and evaluation of mean state

The Taiwan Earth System Model (TaiESM) version 1 is developed based on Community Earth System Model version 1.2.2 of National Center for Atmospheric Research. Several innovative physical and chemical parameterizations, including trigger functions for deep convection, cloud macrophysics, aerosol, and threedimensional radiation–topography interaction, as well as a one-dimensional mixed-layer model optional for the atmosphere component, are incorporated. The precipitation variability, such as diurnal cycle and propagation of convection systems, is improved in TaiESM. TaiESM demonstrates good model stability in the 500-year preindustrial simulation in terms of the net flux at the top of the model, surface temperatures, and sea ice concentration. In the historical simulation, although the warming before 1935 is weak, TaiESM captures the increasing trend of temperature after 1950 well. The current climatology of TaiESM during 1979–2005 is evaluated by observational and reanalysis datasets. Cloud amounts are too large in TaiESM, but their cloud forcing is only slightly weaker than observational data. The mean bias of the sea surface temperature is almost 0, whereas the surface air temperatures over land and sea ice regions exhibit cold biases. The overall performance of TaiESM is above average among models in Coupled Model Intercomparison Project phase 5, particularly in that the bias of precipitation is smallest. However, several common discrepancies shared by most models still exist, such as the double Intertropical Convergence Zone bias in precipitation and warm bias over the Southern Ocean.


Introduction
The Earth system model (ESM) is a state-of-the-art tool that can simulate the long-term evolution of the climate system including the atmosphere, ocean, land, and cryosphere and provide future projections from the scientific aspect to study the impact of global climate change on the natural environment, ecosystem, and human society (IPCC, 2013). Because of the constraint of computing power, the spatial resolution of ESMs participating in the Coupled Model Intercomparison Project Phase 5 (CMIP5; Taylor et al., 2012) is generally on the order of approximately 100 km. However, this coarse resolution is unsuitable for climate studies in the Taiwan area because this island is 400 km long and 150 km wide and occupies only several grid boxes in these ESMs. For the Taiwanese scientific community, developing a global model to provide climate data in various future scenarios with high temporal resolutions -daily or hourly -for dynamical or statistical downscaling is desirable. Taiwan's National Science Council (now Ministry of Science and Technology) has accordingly launched a project to increase climate modeling capability and capacity in Taiwan, the core component of which is Taiwan Earth System Model (TaiESM) development.
In Taiwan, man power and expertise for climate research are limited; thus, we could not create an ESM from scratch. Therefore, TaiESM version 1 is developed on the basis of the Community Earth System Model version 1.2.2 (CESM1.2.2; Hurrell et al., 2013)  Model version 5.3 (CAM5), Community Land Model version 4 (CLM4), Parallel Ocean Program version 2 (POP2), and Community Ice Code version 4 (CICE4). We replace or modify existing parameterizations in CAM5, including new trigger functions for the deep convection scheme (Y.-C. , a new cloud macrophysics scheme for cloud fraction calculation , and a three-moment aerosol scheme (Chen et al., 2013). A novel parameterization for the impact of three-dimensional (3D) radiation-topography interactions  is added to CLM4. In addition, a one-dimensional (1D) mixed-layer ocean model with a high vertical resolution (Tsuang et al., 2009) is used for CAM5 with slab ocean simulation in TaiESM.
An object of TaiESM development is to improve the simulations of climate variability in various spatial and temporal scales for more reliable climate projections in Taiwan. Weather and climate in Taiwan are deeply affected by capricious East Asia and western North Pacific monsoon and typhoons. In addition, because of its small size and steep terrain, predicting the frequencies of severe weather and heavy precipitation in Taiwan is highly difficult (Hsu et al., 2011). Therefore, the parameterizations selected for TaiESM are for enhancing variability simulation. The trigger functions for the deep convection scheme in TaiESM, adopted from the National Centers for Environmental Prediction (NCEP) Global Forecast System (GFS) with the Simplified Arakawa-Schubert scheme (SAS; Pan and Wu, 1995;Han and Pan, 2011), aim to improve the timing of convective precipitation occurrence. As demonstrated by Lee et al. (2008), by using GFS, these trigger functions are key to improved simulations of the diurnal rainfall cycle over the Southern Great Plains (SGP) in the United States. The parameterization for 3D radiation-topography interactions accounts for the effects of shadows and reflections from subgrid topographic variation on the surface solar flux (Lee et al., 2011) designed for the application to general circulation models (GCMs). The high-resolution 1D mixed-layer model can resolve fast change in the skin temperature of the sea surface (Tu and Tsuang, 2005).
The organization of this paper is as follows: Sect. 2 describes TaiESM, particularly the new and modified schemes different from CESM1.2.2. Section 3 presents the design of model experiments. Sections 4 and 5 provide the description of TaiESM performance in preindustrial and historical simulations, respectively. Summary and conclusions are given in Sect. 6.

Model description
The development of TaiESM is based on CESM1.2.2, in which the ocean, sea ice, and river components, as well as the infrastructure of the model, remain unchanged. For the atmosphere, several physical and chemical parameterizations are modified, as two trigger functions are added to the default deep convection scheme, and cloud macrophysics and aerosol schemes are replaced. A parameterization of surface albedo adjustment is added to CLM4 to account for the topographic effect on surface solar radiation. In addition, a 1D mixed-layer ocean model is integrated to TaiESM for simulations of CAM5 coupled with a slab ocean.

Atmosphere
The atmosphere model in TaiESM is based on CAM version 5.3 (Neale et al., 2010). The dynamic core is Finite-Volume (Lin, 2004) in a hybrid sigma-pressure vertical coordinate. The Rapid Radiative Transfer Model for GCMs (RRTMG; Iacono et al., 2008) with two-stream approximation, correlated k-distribution, and Monte Carlo Independence Column Approximation (McICA; Pincus et al., 2003) is employed to calculate radiative fluxes and heating rates in the atmosphere. The shallow convection and moist turbulence schemes are based on those reported by  and Bretherton and Park (2009), respectively. A two-moment cloud microphysics scheme (Morrison and Gettelmen, 2008) is used to predict changes in the mass and number of cloud droplets and to diagnose stratiform precipitation.

Trigger function for deep convection
Convective trigger function is a critical part of the cumulus parameterization scheme to determine the initiation of precipitating convection and thus has a critical role in rainfall variability simulation. With the Zhang-McFarlane scheme framework (Zhang and McFarlane, 1995;Neale et al., 2008), TaiESM has adopted two convection triggers proposed by Y.-C. : unrestricted launching level (ULL) and convective inhibition (CIN). By modifying the deep convection scheme in CESM1.0.3-CAM5.1 Y.-C.  reported significant improvements in the diurnal rainfall peak at the Atmospheric Radiation Measurement (ARM) SGP site, mainly because of the suppression of daytime spurious convection by the CIN trigger and initiation of nighttime mid-level convection by the ULL trigger. ULL may also aid in improving diurnal rainfall phase in many other areas worldwide when implemented in the newly developed Energy Exascale Earth System Model version 1 (E3SMv1) of the U.S. Department of Energy (Xie et al., 2019).
Similar to that in E3SM, improvement in the diurnal rainfall cycle is found in TaiESM. Figure 1 displays local times (LTs) of the diurnal rainfall peak occurrence, referred to as the peak phase from the 11-year (2001-2011) Tropical Rainfall Measuring Mission (TRMM) merged satellite data (Huffman et al., 2007) and the historical model runs during 1979-2005. Areas with an amplitude of diurnal rainfall cycle smaller than 0.5 mm d −1 are masked to emphasize the regions with strong diurnal rainfall signals. Two distinct changes in diurnal rainfall cycle are found in TaiESM com-  2 (1979-2005), and (c) TaiESM . Areas with an amplitude of diurnal precipitation smaller than 0.5 mm d −1 are masked out. pared with those in CESM1.2.2. First, the simulated diurnal rainfall peak over the tropical land areas is improved. For example, the observed peaks in the central Africa (Box A) and the Amazon Basin (Box B) occur around 20:00-22:00 and 18:00-20:00 LT, respectively. These peaks are delayed from 12:00-14:00 LT in CESM1.2.2 to 14:00-18:00 LT in TaiESM. A similar delay is also found in islands such as Borneo and Sumatra. As a result, nocturnal rainfall in TaiESM is increased compared with that in CESM1.2.2.
Second, the propagation of convective organizations is better simulated. Propagating convective systems originating from coastlines or topographical regions could be demonstrated by the gradual phase change in Fig. 1, such as the eastern slope of the Rocky Mountains (Box C) and northern South America (north of Box B). More specifically,  vation (Carbone and Tuttle, 2008). In CESM1.2.2, convection occurs in the early afternoon and peaks before midnight, but it is stationary at the same location. TaiESM successfully captures the eastward propagation of the rainfall and a better occurrence time of convection in the late afternoon, as well as the more realistic rainfall intensity. This result is consistent with the single-column model tests of Y.-C. , indicating that their proposed convective trigger may be the cause of these improvements. Furthermore, Wang and Hsu (2019) demonstrate that the improvement of nocturnal rainfall over SGP is mainly from the superior response of the ULL + CIN convective trigger to the low-level convergence between the branch of the mountain-plain solenoid and low-level jet from the Gulf of Mexico. With the horizontal resolution on the order of 100 km, this result suggests that the convective trigger of TaiESM captures the large-scale preconditioning related to the convective organization there (Dirmeyer et al., 2011) rather than only the convective systems itself.

Cloud fraction
The cloud macrophysics scheme used in TaiESM is the GFS-TaiESM-Sundqvist (GTS) scheme . Its prototype was first developed for the NCEP GFS model and has been further modified for the TaiESM. Similar to that in many numerical weather prediction and global climate models, the GTS scheme is based on the Sundqvist scheme (Sundqvist et al., 1989), which calculates changes in cloud condensates in a grid box on the basis of the budget equation for relative humidity (RH) with large-scale advection. The CAM5 macrophysics (Park et al., 2014) follows this approach and assumes empirical values of critical RH (RH c ) as the threshold of condensation. The key difference of the GTS scheme from the CAM5 macrophysics is the re-derivation of the equation relating the change in the subgrid-scale cloud condensate using the distribution width of the mixing ratio of total water (q t ) to replace RH c , as indicated in Tompkins (2005). The unnecessary use of RH c is consequently removed to allow an improved correlation among cloud fraction, RH, and condensates. Figure 3 illustrates cloud fraction as a function of RH of water vapor (q v /q s ) and RH of condensates (q l /q s ) for the CAM5 macrophysics and the GTS schemes with uniform and triangular probability density functions (PDFs) of q t in a grid box. Given the same RH of water vapor, the PDF-based calculation allows larger cloud fraction if more cloud condensates exist in the grid box than the CAM5 macrophysics. The difference in cloud fraction produced by the two PDFs is small, implying that this scheme might not be very sensitive to the shape of the distribution. The triangular PDF additionally provides rapid changes in cloud fraction when the RH of condensates and water vapor changes, and it is used as the default PDF of the GTS scheme.

Aerosol
The aerosol parameterization used in TaiESM is the Statistical-Numerical Aerosol Parameterization (SNAP; Chen et al., 2013). SNAP is a bulk parameterization, and the modal approach (Seigneur et al., 1986;Whitby and Mc-Murry, 1997) is adopted to describe the particle size distribution. In contrast to conventional aerosol parameterizations in most ESMs, changes in the zeroth moment (number), second moment (surface area), and third moment (mass) due to physical processes are tracked in SNAP. The physical processes included in SNAP are emission, nucleation, coagulation, condensation, mixing, and dry and wet deposition. SNAP has been applied to the US EPA Models-3/Community Multi-scale Air Quality (CMAQ; Byun and Schere, 2006) modeling system and been verified by observations (Chen et al., 2013;Tsai et al., 2015) with the Weather Research and Forecasting Model (WRF; Skamarock et al., 2005).

Land
The land model in TaiESM is CLM4 (Oleson et al., 2010;Lawrence et al., 2011). The surface albedo is primarily a function of vegetation, soil moisture, solar zenith angle, and snow reflectivity calculated by the Snow, Ice, and Aerosol Radiative Model (SNICAR; Flanner and Zender, 2006), which considers the aerosol deposition of black carbon and dust, effective size of snow grains, and vertical profile of heating. As the albedo of a grid box is determined, it is then adjusted to include the topographic effect on surface solar radiation.
The parameterization for 3D radiation-topography interactions is to evaluate the impact of topography on surface solar radiation, including insolation on various slopes and aspects, shadow cast by nearby mountains, and reflections be- tween surfaces . It is developed on the basis of numerous "exact" Monte Carlo calculations that simulate the scattering, reflection, and absorption of photons within the 3D atmosphere and surface (Chen et al., 2006;Liou et al., 2007;Lee et al., 2011). The parameterization adjusts surface albedo so that the solar radiation absorbed by the surface in the land model corresponds to the results of the Monte Carlo calculation. Several geographic parameters are used for input, including the slope, aspect, sky view factor, terrain configuration factor, standard deviation of elevation within a grid box, and solar zenith and azimuth angles. Gu et al. (2012) and Liou et al. (2013) demonstrate that this topographic effect can increase the amount of snowpack in the valley and enhance the snowmelt in mountains in the WRF simulations over the western United States. Lee et al. (2015Lee et al. ( , 2019 also demonstrate that incorporating this parameterization into the Community Climate System Model version 4 (CCSM4) can significantly improve the surface energy budget over the Rocky Mountains and the Tibetan Plateau and thus reduce the systematic cold bias in the CMIP5 models.

Ocean, sea ice, and river
The sea ice and dynamic ocean components of TaiESM are from the CICE4 (Hunke and Lipscomb, 2008) and POP2  of Los Alamos National Laboratory, respectively. The River Transport Model (RTM; Oleson et al., 2010) is designed to route liquid and ice runoff to the ocean as one of the freshwater input to POP2. The configurations of CICE4, POP2, and RTM in the fully coupled TaiESM simulations are identical to those in CESM1.2.2. Note that there is no land ice model in TaiESM. Therefore, the formation of sea ice from the discharge of the ice sheet to the ocean is not simulated.
To save computational resources, a zero-dimensional slab ocean model without dynamical process is commonly used to simulate the thermodynamic interaction between the atmosphere and ocean. In TaiESM, an efficient 1D mixed-layer model is coupled with the atmosphere component to reveal the impact of the fast evolution in upper-ocean layers. The one-column ocean model Snow-Ice-Thermocline (SIT; Tu and Tsuang, 2005;Tsuang et al. 2009) is designed to simulate the sea surface temperature (SST) and upper-ocean temperature variations with a high vertical resolution, including cool skin, diurnal warm layer, and mixed layer of the upper ocean. SIT calculates changes in temperature, momentum, salinity, and turbulent kinetic energy driven by vertical fluxes parameterized using the classical K approach. Cool skin is derived by considering merely molecular transport for vertical diffusion of heat in the skin layer, where the skin layer thickness is calculated as described by Artale et al. (2002). Beneath the skin layer, eddy diffusivity is determined according to a second-order turbulence closure approach (Gaspar et al., 1990), and the 1 m vertical discretization is deployed down to a 10 m depth for resolving diurnal warm layer. Because of the lack of ocean circulation in the one-column ocean model, the calculated ocean temperatures are weakly nudged to clima-tology for ocean below 10 m depth to avoid climate drift. SIT and the atmospheric model exchange SST and fluxes at every time step in the tropics (30 • S-30 • N), whereas climatological SST drives the atmospheric model elsewhere. Note that SIT is not integrated into the dynamic ocean model (POP2); therefore, fully coupled TaiESM simulations do not include SIT.

Model tuning
The preliminary version of TaiESM was very cold compared with CESM1.2.2 using the preindustrial greenhouse gas concentrations and aerosol emissions. The most apparent change was the significant increase in cloud cover, particularly low clouds, which was probably induced by GTS cloud macrophysics and SNAP aerosol schemes. Therefore, several parameters associated with cloud formation were adjusted to reduce shortwave cloud forcing. We first explored that aerosol-cloud interactions were very strong in TaiESM with the SNAP scheme. Therefore, the activation rate of aerosols to cloud condensation nuclei was reduced by 10 % in the microphysics scheme to weaken the aerosol indirect effect. The sizes of detrained liquid particles from shallow convection and solid particles from deep convection were increased from 10 µm in CESM1.2.2 to 14 µm and from 15 to 25 µm, respectively. Larger detrained particles have smaller cloud optical depth and shorter suspension time in the air when the detrained water content is the same. Both effects can reduce cloud albedo. Although RH c is removed from the GTS scheme for grid boxes with the presence of condensates, it is still required for the formation of clouds in a cloud-free grid box. The value of RH c was increased from 0.8 in the free atmosphere in CESM1.2.2 to 0.85 in TaiESM to make cloud formation less efficient. After these adjustments, the global mean surface temperature of TaiESM in the preindustrial simulation is comparable to that of CESM1.2.2 while the radiation imbalance at the top of the atmosphere (TOA) is minimized. Note that this model tuning is made only at the spatial resolution of about 1 • . Additional tuning would be required for stable simulations at higher or lower resolutions.

Experiment design
The horizontal resolution of the atmosphere and land in TaiESM is 0.9 • latitude by 1.25 • longitude, with 30 vertical layers and a model top at 2 hPa in the atmosphere. The ocean and sea ice components use the same horizontal resolution with 320×384 grid points (approximately 1 • ) and 60 vertical layers in the ocean. TaiESM is spun-up using CMIP5 preindustrial conditions, such as greenhouse gas concentrations, surface aerosol emissions, solar constant, and land-use types. Because TaiESM is considerably similar to CESM1.2.2, we use the model restart files of CESM1.2.2 for the 1850 control run as the initial condition to reduce the computational effort, particularly for the ocean component that may need more than 1000 years to reach a steady state. The spin-up integration continues for 500 years, and the climate state at the end of year 500 is used as the initial condition for the 500-year preindustrial control (hereafter piControl) simulation. The historical simulation then starts at the end of pi-Control (i.e., year 1000) with observationally based forcing, including changes in the solar constant, greenhouse gas concentrations, surface aerosol emission, and volcanic eruptions, from 1850 to 2005.

Model stability in the piControl run
In this section, the global means of several climatological variables in the piControl run of TaiESM are evaluated. The climate drift from CESM1.2.2 initial conditions to TaiESM equilibrium during the spin-up is also assessed to represent differences between the two models caused by the new or modified physical processes in TaiESM. Figure 4 illustrates the time series of several global mean variables in TaiESM piControl. The long-term global mean TOA net flux is 0.086 W m −2 , and it decreases by 0.0054 W m −2 in 500 years, but insignificantly. Furthermore, the mean surface net flux is 0.081 W m −2 with an almost identical decreasing trend as TOA net flux. The imbalance at TOA causes heating of the whole model system, and the comparatively smaller imbalance at the surface indicates that a smaller part of excessive energy remains in the atmosphere in piControl. Consequently, the long-term trend of surface air temperature (SAT) is 0.0088 K century −1 in 500 years, which is statistically significant. By contrast, the trend of SST is 0.0047 K century −1 , only about half of the SAT trend and insignificant. By breaking down the surface net flux, we found that the energy exchange between the atmosphere and land is less than 10 −5 W m −2 , whereas the net flux into the ocean is 0.114 W m −2 (figures not shown). The excessive energy enters the deep ocean and leads to a steady increase in global mean ocean temperature of 0.030 K century −1 . Therefore, even after a 1000 years' simulation, the system does not reach the thermodynamic equilibrium. In addition, considering that the heat capacity of the entire ocean is approximately 1000 times larger than the atmosphere, the heating rates of the atmosphere caused by the residual net flux (0.005 W m −2 ) are too small compared with the heating rate of the ocean. It implies that an unknown energy leak may exist in the coupling between the atmosphere and ocean, which requires further investigation in programming to fix this problem.

Time series of climate states
The annual mean time series of sea ice area in the Northern Hemisphere (NH) and Southern Hemisphere (SH) are exhibited in the bottom panels of Fig. 4. The Arctic sea ice has  a small but significant trend of −0.01 × 10 6 km 2 century −1 , corresponding fairly well to the slight warming of the entire model. By contrast, the linear trend of the sea ice area in the Southern Ocean over the 500-year span is almost 0, even though the variation is much larger. The minimal change in the sea ice area indicates that the energy gain of the cryosphere could be negligible compared with other model components. The global mean sea surface salinity (SSS) reduces significantly by −0.0036 g kg −1 century −1 . However, it can be found that SSS is almost constant with a slope of about 10 −4 g kg −1 century −1 after year 700. On the other hand, there is a small but significant decreasing trend of the global mean ocean salinity of 1.3 × 10 −4 g kg −1 century −1 , which is very close to the trend of SSS in the last 300 years.

Comparison with CESM
The long-term means of several variables in piControl runs performed by CESM1.2.2 and TaiESM are listed in  LW cloud forcing (LWCF) are 1.67 and 0.59 W m −2 , respectively; therefore, the warmer surface and atmosphere have greater contribution to additional outgoing longwave radiation (OLR) in TaiESM. However, the amount of high cloud in TaiESM is substantially larger than that in CESM1.2.2. This implies that the high clouds in TaiESM could be optically thinner. The different relation between cloud forcing and cloud cover in SW and LW in TaiESM is probably due to the GTS scheme, which can produce a larger fraction but less dense clouds compared with the cloud macrophysics scheme in CAM5.

Historical simulation
In this section, we evaluate the performance of the TaiESM historical simulation against the observation or reanalysis data. The temporal evolution of global mean temperature from the preindustrial period to the present day is assessed. In the historical simulation, the mean states of the current climate, defined as the period of 1979-2005, are used for comparison. The behavior of the El Niño-Southern Oscillation (ENSO) in TaiESM is also evaluated. Pinatubo (1991), in TaiESM is close to those in the observational data, implying that the radiative forcing due to stratospheric aerosols is in good agreement with the observations. After 1950, the change in SAT of TaiESM follows the observations and captures the trend of global warming very well. The warming rate of TaiESM during 1950-2005 is 1.12 K century −1 , comparable with 1.16 and 1.27 K century −1 of BEST and GISTEMP, respectively. 07. Almost all of the Arctic Ocean is overcast in TaiESM, which is approximately 30 % higher than observational data. Cloud fraction is also severely overestimated over the Antarctic continent and the Southern Ocean. TaiESM produces too much cloud over the southern branch of the Intertropical Convergence Zone (ITCZ) in the central and eastern Pacific, implying the prevalence of double ITCZ, which will be discussed in a subsequent section. An excessive amount of clouds is also noted in the Maritime Continent, western equatorial Indian Ocean, and most of the land areas. By contrast, cloud fraction is remarkably underestimated in the Amazon Basin and the subtropical ocean, particularly the stratocumulus near the western coasts of continents. Compared with the synergic CloudSat and Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) data during (Kay and Gettelman, 2009, low clouds in TaiESM are systematically underestimated over the entire tropical and subtropical regions, as shown in Fig. 6c, whereas they are overestimated in high-latitude areas. The total cloud fraction in the tropics is high because of excessive high cloud in the model (Fig. 6b). Clouds can substantially modulate the radiation field because of its high reflectivity in SW and high absorptivity in LW. Figure 7a illustrates the comparison of SWCF in TaiESM with that in Clouds and the Earth's Radiant Energy System-Energy Balanced and Filled data (CERES-EBAF; Kato et al., 2018Kato et al., ) over 2000Kato et al., -2015 In terms of the global mean, SWCF in TaiESM is very close to that of the observational data with a difference of only 0.19 W m −2 . Although there is excessive cloud over the polar regions in TaiESM, such as the Southern Ocean near the Antarctic continent and almost all of the Arctic Ocean, SWCF is not as strong as that in the observational data. It could be contributed from the op- tically thin polar clouds due to the GTS cloud macrophysics scheme and from the positive bias of sea ice albedo in the Arctic Ocean in TaiESM (not shown). In the subtropical and tropical regions, SWCF generally follows the spatial pattern of total cloud fraction -i.e., a larger cloud fraction produces stronger SWCF, such as the storm track in the North Pacific, southern branch of ITCZ, Maritime Continent, western tropical Indian Ocean, and south of the Sahara Desert. However, SWCF is too strong over the Amazon Basin in TaiESM, even though there is an underestimated amount of clouds. By contrast, because of the underestimated total cloud fraction, SWCF in TaiESM is too weak over the stratocumulus areas off the California and Peru coasts as well as over the subtropical Pacific, Atlantic, and Indian oceans in the SH.

Cloud and radiation
The global mean of LWCF in TaiESM is significantly weaker (by 4.31 W m −2 ) than that in CERES-EBAF. As illustrated in Fig. 7b, TaiESM underestimates LWCF worldwide, and the magnitude of the LWCF bias generally follows the bias of high cloud. Positive LWCF bias only exists in some regions over the tropical ocean with too many high clouds in TaiESM. However, although more high clouds exist along the northern branch of the ITCZ, LWCF is weaker in the model. The remarkable negative LWCF bias seems incompatible with the overestimated high clouds because more high clouds should be able to intercept more LW radiation from the surface. This inconsistency is probably due to the lower altitude of the high clouds or the less dense clouds in TaiESM. Figure 8a illustrates the comparison of SST between TaiESM and the Hadley Centre Sea Ice and Sea Surface Temperature dataset (HadISST; Rayner et al., 2003). The regions with a long-term mean sea ice concentration larger than 15 % are not used for calculations of the mean and RMSD. The global mean bias of SST in TaiESM is 0.01 K with an RMSD of 1.05 K. The overestimated SST over the Southern Ocean and subtropical South Pacific is probably induced by additional downward SW radiation because of the inaccurate microphysical properties of polar clouds (Kay et al., 2016) and the negative bias of cloud fraction as shown in Fig. 6a. The warm bias in the major upwelling regions off the western coasts of the Americas and Africa is a common deficiency in many climate models (Griffies et al., 2009), caused by insufficient spatial resolution of the atmosphere and ocean. Warm bias can also be found in the North Atlantic including the coast of North America, the Labrador Sea, and south of Greenland. Negative biases exist in most of the North Pacific and subtropical North Atlantic, probably because of overestimated wind stress in these regions.

Surface temperature
Although the SST bias in TaiESM is very small, the global mean SAT in TaiESM is substantially colder than the observational data (by 0.49 K with an RMSD of 1.68 K). This result indicates that the temperature over land and sea ice in TaiESM is severely underestimated (Fig. 8b). Cold bias exists over most of the polar regions, the Tibetan Plateau, and tropical land areas (e.g., Amazonia, central Africa, and southeast Asia). It must be due to the excessive cloud that reflects excessive sunlight. SAT bias over the ocean generally follows SST bias, except that the SAT bias in the subtropical South Pacific is very small despite the warm SST bias.  (Lin, 2007;Hirota and Takayabu, 2013) and in CESM1 (C.-c. . The precipitation rates of both the northern and southern ITCZ branches are extremely strong. The overly intense convection strengthens the subsidence and consequently produces too little rainfall along the Equator. Precipitation is also overestimated in the Maritime Continent, while it is severely underestimated in Borneo. In TaiESM, the land-sea contrast in precipitation is not as apparent as in the observation over the warm pool region. The South Pacific convergence zone (SPCZ) is also too strong and too parallel to the ITCZ. The dipole bias in the tropical Indian Ocean, excessive rainfall in the western part and scant rainfall in the eastern part, still exists as in NCAR models (Gent et al., 2011). There is also a double ITCZ bias in the Atlantic Ocean in that the southern branch is too strong and the northern branch is too weak. In South America, precipitation over the Amazon Basin is considerably underestimated, whereas excessive orographic precipitation can be found along the Andes (Cook et al., 2012).  (Peng et al., 2013(Peng et al., ), during 1979(Peng et al., -2005 In the NH, TaiESM severely overestimates sea ice concentration over the North Pacific, particularly in the Sea of Okhotsk. TaiESM also overestimates sea ice in the Barents Sea and near the east coast of Greenland but slightly underestimates sea ice in the Labrador Sea. In the SH, sea ice in TaiESM is generally in agreement with the observation. Excessive sea ice is noted in the area south of New Zealand, but in the Indian Ocean region, sea ice is scant. This deviation follows the SST bias presented previously. Figure 11 illustrates the temporal evolution of the annual sea ice concentration in TaiESM compared with that in the CDR. The change in NH sea ice in TaiESM generally captures the trend in the observation before 2002. However, there  -CDR, 1979-CDR, -2005 is an increase in TaiESM in the last 4 years, in contrast to an accelerated reduction in observational data. This sea ice increase could be a fluctuation in a climate simulation, and it requires longer integration for additional investigation. In the SH, a decreasing trend of the sea ice concentration can be found in TaiESM, whereas it remains almost unchanged in observational data. Because there is no land ice model in TaiESM, the discharge of the ice sheet from the Antarctic continent to the Southern Ocean, the major source of SH sea ice, cannot be simulated accurately. Consequently, the sea ice concentration in the SH could be controlled primarily by temperature in TaiESM, leading to an unrealistic temporal evolution. Figure 11. Time series of annual mean total sea ice area for both NH and SH from TaiESM historical run and observation.

ENSO
To evaluate the ENSO behavior during 1976-2005 in TaiESM, the HadISST sea surface temperature and MRE2 reanalysis data in the same period are used. The observed and simulated spectra of the Nino 3.4 index presented in Fig. 12 reveal the adequate ability of TaiESM in reproducing the periodicity of El Niño. The observed Nino 3.4 index exhibits three statistically significant peaks between 2 and 6 years. TaiESM is able to simulate three spectral peaks with slightly shorter periods, while the amplitudes of all three peaks are larger than observations. The anomalies of surface temperature, sea level pressure, and near-surface wind in December-February when the ENSO is at the mature stage are shown in Fig. 13, which is the composite of five and six El Niño events in observation and TaiESM simulation, respectively. The simulated SST anomaly (SSTA) is evidently larger in both amplitude and spatial coverage than the observed one and with the maximum shifted westward to the central equatorial Pacific compared with the observation, which is the common bias in many climate models. The horseshoe-like negative SSTA in the northwest, southwest, and west of the positive SSTA is stronger and covers much larger areas than the observed one. This over-simulated SSTA structure leads to some marked biases in the simulated atmospheric circulation and temper- ature, such as the cold bias in the western North Pacific and Maritime Continent, warm bias in the western Indian Ocean and Bering Sea, and too strong a convergence in the eastern equatorial Pacific.

Comparison with CMIP5 models
The overall performance of the TaiESM historical simulation during 1979-2005 is evaluated by comparing it with other CMIP5 models following the metrics introduced by Gleckler et al. (2008). Figure 14 shows the normalized space-time root-mean-square error (RMSE) of selected variables from TaiESM, several CMIP5 models, and a multi-model ensemble (MME) against reanalysis and observation datasets. The reference data of air temperatures (TA), zonal and meridional wind velocities (UA and VA), and geopotential height (ZG) at various pressure levels, as well as the surface air temperature (TAS), are from the Collaborative Reanalysis Technical Environment (CREATE) Multi-Reanalysis Ensemble version 2 (MRE2; Potter et al., 2018). The observational precipitation (PR) data are from GPCP. Upward longwave radiation in the total sky (RLUT) and clear sky (RLUTCS) and upward shortwave radiation in the total sky (RSUT) and clear sky (RSUTCS) are from CERES-EBAF. It is expected that the errors of CMIP5 MME are generally the smallest. TaiESM has the smallest bias in PR among all CMIP5 models, and its performance in RSUT and RLUT is also very good. The relatively poor performance in TAS is primarily due to the cold bias over land and sea ice areas. The RMSEs of all variables in TaiESM are smaller than the median CMIP5 error, indicating that the performance of TaiESM is above average among all CMIP5 models. The performance of TaiESM is comparable to that of CESM1-CAM5, and they have similar strengths and weaknesses. Note that three variables with an RMSE larger than the median in CESM1-CAM5 are all improved in TaiESM.

Summary and conclusions
This paper documents the TaiESM version 1, developed on the basis of CESM1.2.2, with revised physical and chemical parameterizations, including (1) trigger functions for deep convection, which can improve the variability simulation in convective rainfall; (2) the GTS cloud macrophysics scheme to avoid an artificial RH threshold for cloud formation; (3) the three-moment SNAP aerosol scheme; (4) 3D radiation-topography interactions to account for the impact of shading and reflection on shortwave radiation in moun- Figure 14. The space-time RMSEs of upward longwave radiation at TOA in total sky and clear sky (RLUT and RLUTCS), upward shortwave radiation at TOA in total sky and clear sky (RSUT and RSUTCS), precipitation (PR), surface air temperature (TAS), geopotential height (ZG), meridional wind (VA), zonal wind (UA), and air temperature (TA) from TaiESM, CMIP5 models, and CMIP5 MME. The values of shading represent the magnitude of the normalized error with respect to the median CMIP5 error. For example, a value of −0.2 indicates that the RMSE of a model is 20 % smaller than the median error. tains. A 1D mixed-layer ocean model is incorporated into the atmosphere component to simulate the thermodynamic air-sea interaction, but it is not used for fully coupled simulations.
TaiESM stability is assessed using the 500-year piControl. Although constant imbalance in the net flux at the TOA exists, the drifts of global mean SAT and SST are very small, with long-term trends of 0.0088 and 0.0047 K century −1 , respectively. The excessive energy enters the deep ocean and leads to continuous warming by 0.030 K century −1 . The drifts in the sea ice concentration in both NH and SH are both small because of the nearly zero net energy flux from the atmosphere to sea ice. However, the global mean SSS and total ocean salinity both demonstrate significantly decreasing trends.
For the historical evolution of SAT, the warming of TaiESM from 1850 to 1935 is too weak compared with the observation. After 1950, TaiESM satisfactorily captures the trend of global warming with a heating rate of 1.12 K century −1 comparable to the observation of 1.16 K century −1 .
The current climatology of TaiESM during 1979-2005 is generally in agreement with the observations. The overall performance of TaiESM is better than the median of CMIP5 models, particularly in that the RMSE of precipitation is smallest. There are too many clouds in TaiESM, whereas the SWCF and LWCF are mostly similar to and weaker than the observation, respectively. This result implies that the new cloud macrophysics scheme produces a larger amount of but optically thinner clouds. SST in TaiESM is very close to the observation, whereas SAT is significantly colder, implying remarkably underestimated SAT over land and sea ice surfaces. TaiESM produces excessive precipitation, and the biases of the double ITCZ and dipole in the tropical Indian Ocean exist, whereas there is a severe dry bias in the Amazon Basin. The trend of the NH sea ice concentration in TaiESM follows the observation well, whereas it might not capture the accelerating reduction in the 21st century. For the ENSO simulation, TaiESM is able to reproduce three spectral peaks similar to observation with periods between 2 and 6 years, while the variability of SST, including the magnitude of the anomaly and spatial coverage, is too strong.
This paper focuses on the evaluation of the long-term climatological state and evolution of global mean quantities in TaiESM in preindustrial and historical simulations. The other part of the characteristics of an ESM, climate variability, is also very critical to the performance of a model, and it requires additional in-depth research. Further investigation of climate variability in TaiESM, including the intraseasonal oscillation, monsoon, and extreme precipitation, will be documented in follow-up papers.
Author contributions. HHH is the initiator and the primary investigator of the TaiESM project. WLL is the main model developer and wrote the majority of the paper. YCW is the developer and writer of trigger functions for deep convection. CJS and YCW are the developer and writers of cloud macrophysics. ICT and JPC are the developers and writers of SNAP aerosol scheme. CYT and YYL are developers of 1D mixed-layer model, and CYT is the writer of this section. HLP helped develop the theoretical basis of trigger functions for deep convection and cloud macrophysics. Review statement. This paper was edited by Qiang Wang and reviewed by Ingo Bethke and one anonymous referee.