Description and basic evaluation of Beijing Normal University Earth System Model ( BNU-ESM ) version 1

An earth system model has been developed at Beijing Normal University (Beijing Normal University Earth System Model, BNU-ESM); the model is based on several widely evaluated climate model components and is used to study mechanisms of ocean-atmosphere interactions, natural climate variability and carbon-climate feedbacks at interannual to interdecadal time scales. In this paper, the model structure and individual components are described briefly. Further, results for the CMIP5 (Coupled Model Intercomparison Project phase 5) pre-industrial control and historical simulations are presented to demonstrate the model’s performance in terms of the mean model state and the internal variability. It is illustrated that BNU-ESM can simulate many observed features of the earth climate system, such as the climatological annual cycle of surface-air temperature and precipitation, annual cycle of tropical Pacific sea surface temperature (SST), the overall patterns and positions of cells in global ocean meridional overturning circulation. For example, the El Niño-Southern Oscillation (ENSO) simulated in BNU-ESM exhibits an irregular oscillation between 2 and 5 years with the seasonal phase locking feature of ENSO. Important biases with regard to observations are presented and discussed, including warm SST discrepancies in the major upwelling regions, an equatorward drift of midlatitude westerly wind bands, and tropical precipitation bias over the ocean that is related to the double Intertropical Convergence Zone (ITCZ).

have been published via an Earth System Grid Data Node located at Beijing Normal University (BNU) and can be accessed at http://esg.bnu.edu.cn, as a part of internationally federated, distributed data archival and retrieval system, referred to as the Earth System Grid Federation (ESGF).
Many studies have utilized CMIP5 results from BNU-ESM, and the model has received comprehensive evaluations. For example, Wu et al. (2013) evaluated the precipitation-surface temperature (P -T ) relationship of BNU-ESM among 17 models in CMIP5 and found BNU-ESM has better ability in simulating P -T pattern correlation than other models, especially over ocean and tropics. Bellenger et al. (2013) used the metrics developed within the Climate Variability and Predictability (CLIVAR) Pacific Panel and additional metrics to evaluate the basic El Niño-Southern Oscillation (ENSO) properties and associated feedbacks of BNU-ESM and other CMIP5 models. BNU-ESM performs well on simulating precipitation anomalies over the Niño-4 region; the ratio between the ENSO spectral energy in the 1-3 year band and in 3-8 year band is well consistent with observational result, but the model has stronger sea surface temperature (SST) anomalies than observational estimates over Niño-3 and Niño-4 regions. Fettweis et al. (2013) reported BNU-ESM can simulate the 1961-1990 variability of the June-August (JJA) North Atlantic Oscillation (NAO) well and the sharp decrease of the NAO index over the last 10 years as observed, and the model projects similar negative NAO values into the future under RCP 8.5 scenario. Gillett and Fyfe (2013) reported no significant Northern Annular Mode (NAM) decrease in any season between 1861 and 2099 in historical and rcp45 simulations of BNU-ESM as with the other 36 models from CMIP5. Bracegirdle et al. (2013) assessed the model's simulation of near-surface westerly winds over the Southern Ocean and found an equatorward bias in the present-day zonal mean surface jet position in common with many of the CMIP5 models. Among other studies, Chen et al. (2013) evaluated the cloud and water vapor feedbacks to El Niño warming in BNU-ESM. Vial et al. (2013) diagnosed the climate sensitivity, radiative forcing and climate feedback of BNU-ESM. Roehrig et al. (2013) assessed the performance of BNU-ESM on simulating the West African Monsoon. Sillmann et al. (2013) evaluated the model performance on simulating climate extreme indices defined by the Expert Team on Climate Change Detection and Indices (ETCCDI). Wei et al. (2012) utilized BNU-ESM in assessment of developed and developing world responsibilities for historical climate change and CO 2 mitigation.
Although the simulation results from BNU-ESM are widely used in many climate studies, a general description of the model itself and its control climate is still not available. Documenting the main features of the model structure and its underlying parameterization schemes will help the climate community to further understand the results from BNU-ESM.

D. Ji et al.: Description and basic evaluation of BNU-ESM 2041
This paper provides a general description and basic evaluation of the historical climate simulated by BNU-ESM. Particular focus is put on the model structure, the simulated climatology, internal climate variability and terrestrial carbon cycle deduced from the piControl and historical simulations submitted for CMIP5. The climate response and scenario projections in BNU-ESM will be covered elsewhere. The paper is organized as follows. In Sect. 2, a general overview of BNU-ESM is provided, elaborating on similarities and differences between the original and revised model components in BNU-ESM. In Sect. 3, the design of the piControl and historical model experiments is briefly presented, as well as the spin-up strategy. In Sect. 4, the general model performance is evaluated by using the Taylor diagram (Taylor, 2001). The following two sections focus on the model performance on simulating physical climatology and climate variability. Several key modes of internal variability on different timescales ranging from interseasonal to interdecadal are evaluated. The terrestrial carbon cycle is evaluated in Sect. 7, and particular focus is put on terrestrial primary productions and soil organic carbon stocks. Finally, the paper is summarized and discussed in Sect. 8.

Atmospheric model
The atmospheric component in BNU-ESM is based on Community Atmospheric Model version 3.5 (CAM3.5), which is an interim version of the Community Atmospheric Model version 4 (CAM4) (Neale et al., 2010(Neale et al., , 2013. Here, the main difference of the atmospheric component in BNU-ESM relative to the original CAM3.5 model is the process of deep convection. The BNU-ESM uses a modified Zhang-McFarlane scheme in which a revised closure scheme couples convection to the large-scale forcing in the free troposphere instead of to the convective available potential energy in the atmosphere (Zhang, 2002;Zhang and Mu, 2005a). On the other hand CAM3.5 adopts a Zhang-McFarlane scheme (Zhang and McFarlane, 1995) modified with the addition of convective momentum transports (Richter and Rasch, 2008), and a modified dilute plume calculation (Neale et al., 2008) following Blyth (1986, 1992). BNU-ESM uses the Eulerian dynamical core in CAM3.5 for transport calculations with a T42 horizontal spectral resolution (approximately 2.81 • × 2.81 • transform grid), with 26 levels in the vertical of a hybrid sigma-pressure coordinates and model top at 2.917 hPa. Atmospheric chemical processes utilize the tropospheric MOZART (TROP-MOZART) framework in CAM3.5 , which has prognostic greenhouse gases and prescribed aerosols. Note that the aerosols do not directly interact with the cloud scheme so that any indirect effects are omitted in CAM3.5, as well as in BNU-ESM.

Ocean model
The ocean component in BNU-ESM is based on the GFDL Modular Ocean Model version 4p1 (MOM4p1) released in 2009 (Griffies, 2010). The oceanic physics is unchanged from the standard MOM4p1 model, and the main modifications are in the general geometry and geography of the ocean component. MOM4p1 uses a tripolar grid to avoid the polar singularity over the Arctic, in which the two northern poles of the curvilinear grid are shifted to land areas over North America and Eurasia (Murray, 1996). In BNU-ESM, MOM4p1 uses a nominal latitude-longitude resolution of 1 • (down to 1/3 • within 10 • of the equatorial tropics) with 360 longitudinal grids and 200 latitudinal grids, and there are 50 vertical levels with the uppermost 23 layers each being 10.143 m thick. The mixed layer is represented by the K profile parameterization (KPP) of vertical mixing (Large et al., 1994). The idealized ocean biogeochemistry (iBGC) module is used in BNU-ESM, which carries a single prognostic macronutrient tracer (phosphate, PO 4 ), and simulates two main representative biogeochemical processes, i.e., the net biological uptake in the euphotic zone due to phytoplankton activity as a function of temperature, light and phosphate availability, and regeneration of phosphate as an exponential function below the euphotic zone.

Sea ice model
The BNU-ESM sea ice component is the Los Alamos sea ice model (CICE) version 4.1 (Hunke and Lipscomb, 2010). The CICE was originally developed to be compatible with the Parallel Ocean Program (POP), but has been greatly enhanced in its technical and physical compatibility with different models in recent years. In particular, supporting tripolar grids makes it easier to couple with MOM4p1 code. In BNU-ESM, CICE uses its default shortwave scheme, in which the penetrating solar radiation is equal to zero for snow-covered ice, that is, most of the incoming sunlight is absorbed near the top surface. The visible and near infrared albedos for thick ice and cold snow are set to 0.77, 0.35, 0.96 and 0.69, respectively, slightly smaller than the standard CICE configuration, as they are used as tuning parameters during model control integration. The surface temperature of ice or snow is calculated in CICE without exploiting its "zero-layer" thermodynamic scheme, and the "bubbly brine" model based parameterization of ice thermal conductivity is used.

Land model
The land component in BNU-ESM is the Common Land Model (CoLM), which was initially developed by incorporating the best features of three earlier land models: the biosphere-atmosphere transfer scheme (BATS) (Dickinson et al., 1993), the 1994 version of the Chinese Academy of Sciences Institute of Atmospheric Physics LSM (IAP94)  (Dai and Zeng, 1997) and the NCAR Land Surface Model (LSM) (Bonan, 1996(Bonan, , 1998. The CoLM was documented by Dai et al. (2001) and introduced to the modeling community in Dai et al. (2003). The initial version of CoLM was adopted as the Community Land Model (CLM) for use with the Community Climate System Model (CCSM). The land model was then developed separately at NCAR and BNU. Currently, the CoLM is radically different from its initial version and the CLM (Dai et al., 2004;Bonan et al., 2011); including the following: (i) improved two stream approximation model of radiation transfer of the canopy, with attention to singularities in its solution and with separate integrations of radiation absorption by sunlit and shaded fractions of canopy. (ii) A photosynthesis-stomatal conductance model for sunlit and shaded leaves separately, and for the simultaneous transfers of CO 2 and water vapor into and out of the leaf. (iii) Lund-Potsdam-Jena (LPJ) model (Sitch et al., 2003) based dynamical global vegetation model and terrestrial carbon cycle, and LPJ-DyN (Xu and Prentice, 2008) based scheme on carbon-nitrogen cycle interactions. Note that in all BNU-ESM's CMIP5 and GeoMIP simulations, carbon-nitrogen cycle interactions are turned off as the nitrogen cycle has not yet been fully evaluated.

Component coupling
The coupling framework of BNU-ESM is largely based on the coupler in NCAR CCSM3.5 (an interim version of NCAR CCSM4), with changes on grid mapping interpolation to allow for the identical tripolar grids used in both ocean and sea ice components. The time evolution of the whole model and communication between various component models are all synchronized and controlled by the coupler in the BNU-ESM. Since MOM4p1 and CICE4.1 are both Arakawa B-grid models, the coupling between them is efficient, and the exchanged fields need no transformation or additional treatment (e.g., vector rotation, grid remapping, grid-point shifting, etc.). The different model components are run simultaneously from their initial conditions. The atmospheric component uses a 1 h time step for atmospheric radiation and 20 min time step for other atmospheric physics. The ocean, sea ice and land components have a 2 h, 1 h and 30 min time step, respectively, while direct coupling occurs hourly among atmospheric, sea ice and land components, and daily with the ocean component without any flux adjustment.
All biogeochemical components are driven by the physical climate with the biogeochemical feedback loops combined. The terrestrial carbon cycle module determines the exchange of CO 2 between the land and the atmosphere. It is coupled to the physical climate through the vegetation distribution and leaf area index, which affects the surface albedo, the evapotranspiration flux and so on. As with the terrestrial carbon cycle module, the ocean biogeochemistry module calculates the ocean-atmosphere exchange of CO 2 , and both are coupled with the TROP-MOZART framework in the atmospheric component to form a closed carbon cycle.

Experiments
Following CMIP5 specifications (Taylor et al., 2009), BNU-ESM has performed all CMIP5 long-term core experiments and part of the tier-1 experiments. The CMIP5 specification requires each model to reach its equilibrium states before kicking off formal simulations, especially for long-term control experiments. BNU-ESM adopted a two-step spin-up strategy to achieve model equilibrium. Firstly, the land component including vegetation dynamics and terrestrial carbon cycle, and the ocean component including biogeochemical module were separately spun-up to yield an initial estimate of equilibrium states. In these off-line integrations of the first step spin-up, surface physical quantities such as winds, temperature, precipitation, moisture, and radiation flux are taken as the climatology of a pre-industrial run of the fully coupled BNU-ESM with carbon cycles turned off. Then, the resultant equilibrated physical and carbon cycle states were fed into the coupled model as initial conditions to do on-line spin-up to achieve final equilibrium states. During the second stage, the coupled model was forced with constant external conditions as specified for CMIP5 pre-industrial control simulation as stated below.   (Rienecker et al., 2011); d ERBE (Barkstrom, 1984); e CERES-EBAF (Loeb et al., 2009); f GPCP (Adler et al., 2003); g CMAP (Xie and Arkin, 1997); h ISCCP-D2 (Rossow and Schiffer, 1999;Rossow and Dueñas, 2004); i CLOUDSAT (L'Ecuyer et al., 2008); j RSS (Wentz, 2000(Wentz, , 2013; k NVAP ; l HadISST (Rayner et al., 2003); m OISST (Reynolds et al., 2002); n NOCS (Josey et al., 1999); o FLUXNET-MTE ; p LDEO (Takahashi et al., 2009). In this paper, we focus on the 559 year (from model year 1450 to 2008) pre-industrial control simulation (piControl) and 156 year historical simulation representing the historical period from year 1850 to 2005. The piControl simulation is integrated with constant external forcing prescribed at 1850 conditions (the solar constant is 1365.885 W m −2 , the concentrations of CO 2 , CH 4 , N 2 O are 284.725 ppmv, 790.979 ppbv, and 275.425 ppbv respectively, CFC-11, CFC-12 and volcanic aerosols are assumed to be zero). In terms of energy balance and model stability, the global mean topof-atmosphere (TOA) net radiation flux over piControl period is 0.88 W m −2 , while the global mean surface net radiation flux is 0.86 W m −2 . The global mean sea surface temperature over piControl period is 17.69 • C with a warming drift of 0.02 • C per century (Fig. 1). The historical simulation is initialized with the model states of 1850 year from pi-Control simulation, and forced with natural variation of solar radiation Wang et al., 2005), anthropogenic changes in greenhouse gases concentrations, stratospheric sulphate aerosol concentrations from explosive volcanoes (Ammann et al., 2003), and aerosol concentrations of sulfate, black and organic carbon, dust and sea salt according to Lamarque et al. (2010). Note that there is no land-cover change related to (anthropogenic) land use because the vegetation distributions evolve according to the model-simulated climate, and the areal fraction of non-vegetated regions (lake, wetland, glacier and urban) are fixed according to the Global Land Cover Characterization (GLCC) Database. Therefore, changes in physical and biogeochemical properties of the vegetation due to actual land-cover changes are excluded by design.

General model performance
To systematically evaluate the general performance of BNU-ESM, we use the Taylor diagram (Taylor, 2001;Gleckler et al., 2008), which relates the "centered" root-mean square (RMS) error, the pattern correlation and the standard deviation of particular climate fields. We selected 24 fields (Table 1) and compared model simulations with two different reference data sets (only one data set was available for gross primary production over land and surface CO 2 flux over ocean). The selection rationale for the fields and reference data sets follows Gleckler et al. (2008), where most of reference data sets are briefly described. One notable difference is that we use ERA-Interim (Dee et al., 2011) and JRA-55 (Ebita et al., 2011) reanalysis data instead of ERA40 and NCEP to reflect recent advances in reanalysis systems. We use estimates of specific humidity from National Aeronautics and Space Administration (NASA) Modern Era Retrospective analysis for Research and Applications (MERRA, Rienecker et al., 2011) instead of the Atmospheric Infrared Sounder (AIRS) experiment, as Tian et al. (2013) indicated MERRA specific humidity probably has a smaller uncertainty than the AIRS data set. The International Satellite Cloud Climatology Project (ISCCP, Rossow and Schiffer, 1999;Dueñas, 2004) D2 andCLOUDSAT (L'Ecuyer et al., 2008) data sets are used to examine the total cloud cover. The Clouds and the Earth's Radiant Energy System -Energy Balanced and Filled (CERES-EBAF) data set (Loeb et al., 2009) is used instead of the CERES observations, because the energy balanced characteristics of CERES-EBAF made it more suitable for the near balanced energetics of the earth system. Two carbon cycle fields (gpp and fgco2) were added to fill the gap between climate system model and earth system model. The reference data used to examine gross primary production (gpp) over land is FLUXNET Model Tree Ensembles (FLUXNET-MTE) estimates , which are restricted to vegetated land surface. The reference data used to examine surface CO 2 flux over ocean (fgco2) is from Lamont-Doherty Earth Observatory (LDEO, Takahashi et al., 2009), this climatology data set was created from about 3 million direct observations of seawater pCO 2 around the world between 1970 and 2007. Figure 2 shows six climatological annual-cycle space-time Taylor diagrams for the 24 selected fields in Table 1 for the tropical (20 • S-20 • N) and the northern extra-tropical (20-90 • N) zones. It is clear from Fig. 2 that the accuracy of the model varies between fields and domains. Some simulated fields over the northern extra-tropics have correlations with the reference data of greater than 0.95 (e.g., zg-500hPa, ta-850hPa, rlut, rsnt, tos), and most of fields have correlations with the reference data of greater than 0.8, whereas one field has much lower correlation of 0.38 (fgco2 over the northern extra-tropics). The amplitude of spatial and temporal variability simulated by the model is reasonably close to that of observationally based reference data. The normalized standard deviations between the simulation and the reference data of most fields have a bias of less than 0.25, and several fields have a bias of less than 0.1 (e.g., ta-850hPa, hus-850hPa, rlut, rsnt, psl, tos). One outlier in Fig. 2 (NHEX G3 and TROP G3) is the sensible heat flux over ocean (hfss) examined with National Oceanography Centre, Southampton (NOCS) reference data (Josey et al., 1999). The model shows better skills when compared to ERA-Interim reanalysis, although the pattern correlations against two reference data sets are both of about 0.6. Previous studies suggest that there are large uncertainties in NOCS data set, and their pattern has better agreement with reanalysis products than the magnitude Each field is normalized by the corresponding standard deviation of the reference data, which allows multiple fields to be shown in each sub-figure. Red/Blue markers represent the simulation field evaluated against the Reference1/Reference2 data defined in Table 1. of their fluxes (e.g., Taylor, 2000). In general, most of fields over the tropics are closer to reference data than those over the northern extra-tropics in Taylor diagrams, but some fields with relatively high correlations in the northern extra-tropics have a lower skill in the tropics. These features are consistent with Gleckler et al. (2008). Figure 3 shows the zonally averaged mean atmospheric temperature, zonal wind and specific humidity for the historical simulation of the BNU-ESM and its deviations from the ERA-Interim reanalysis (Dee et al., 2011). The air temperature in the troposphere is in general cold for both boreal summer and winter, especially during the boreal summer ( Fig. 3a). Near the polar tropopause (about 250 hPa) there is a relatively large cold bias up to 8 K over the Arctic during JJA, and up to 10 K over the Antarctica during December-February (DJF). This tropospheric cold bias is one common problem in many CMIP5 models (Charlton-Perez et al., 2013;Tian et al., 2013). In the lower polar troposphere and specific humidity (c) climatology from BNU-ESM historical simulation (black contours) and bias relative to the ERA-Interim climatology (color filled, color bar is of same units except as % for specific humidity) for 1986-2005. during JJA, there is a notable cold bias over the Antarctic. In the stratosphere, the very low winter temperature at 50 hPa in the Southern Hemisphere associated with the polar night jet is overestimated in the model.

Atmospheric mean state
With respect to zonally averaged winds (Fig. 3b), the seasonal mitigation of the northern tropospheric jet is well captured in the simulation, but the westerlies at 200 hPa in this jet are too strong by up to 4 m s −1 during DJF and 8 m s −1 during JJA compared with ERA-Interim reanalysis. The southern tropospheric jet during DJF is also too strong by up to 12 m s −1 , while the westerlies from the surface to about 100 hPa at 60 • S during DJF are weak relative to the reanalysis. The westerly wind maximum in the Southern Hemisphere during JJA extends upward into the stratosphere at higher latitudes as is observed. In the stratosphere, the polar-night jets in both hemispheres are shifted slightly polewards relative to the reanalysis. Over the equator in the upper tropopause the model overestimates the easterly velocities, the largest biases occur at roughly 50 hPa. Figure 3c shows the modeled zonally averaged specific humidity and their differences relative to the ERA-Interim reanalysis shown as percentages because the relative error provides a better measure of the water vapor's impact on the radiative transfer than does the absolute errors (Soden et al., 2005). The model can simulate the strong meridional and vertical gradients in tropospheric specific humidity that decrease with both latitude and altitude. For example, the specific humidity decreases from around 14 g kg −1 at 1000 hPa near the equator to around 1 g kg −1 at 1000 hPa near the poles and around 0.5 g kg −1 at 300 hPa over the equator. In comparison with ERA-Interim reanalysis, the model has a moist tendency in the southern tropical upper troposphere (above 700 hPa) and a slightly dry tendency in the tropical lower troposphere. In terms of relative difference, the model's dry bias in the tropical lower troposphere approaches 15 %, and the wet bias in the tropical upper troposphere approaches 50 %. This humidity bias pattern is also presented in many CMIP5 models (Tian et al., 2013).
Clouds are always a major source of uncertainty in climate models. In BNU-ESM the total cloud fraction is generally underestimated (Fig. 4a), the global mean value for the years 1976-2005 of the historical simulation gives a bias of −14 % with a root-mean-square error (RMSE) of 18 % compared with the ISCCP observational data set. A notable exception is Antarctica where there are too many clouds. The tropical central eastern Pacific and southern Africa also have more clouds than observations. The latitudinal averaged cloud fraction bias within the tropics and subtropics is much lower than at higher latitudes (Fig. 4b), and is similar to results from the original CAM3.5 and CAM4 at 2 • × 2 • horizontal resolution (Neale et al., 2013). At the same time, the liquid water in clouds over ocean is generally exaggerated in the simulation (Fig. 4c), and is particularly pronounced in the extratropical storm track regions.
Clouds have a significant impact on the global radiative balance that is often assessed using TOA shortwave cloud forcing (SWCF) and long-wave cloud forcing (LWCF) (Ramanathan et al., 1989). In BNU-ESM, the simulated shortwave cooling effect of clouds is too strong in the tropics and too weak in the mid-latitudes (Fig. 5b), especially over oceans, these biases are common in climate models (Trenberth and Fasullo, 2010). BNU-ESM also overestimates LWCF in the tropics due to the presence of a double Intertropical Convergence Zone (ITCZ) (Fig. 5d), and it largely offsets the bias of SWCF in the tropics. In AMIP simulation with sea surface temperature and sea ice boundary conditions specified, the SWCF biases in BNU-ESM (not shown) resemble that in CAM4, except for Eurasian continent (Kay et al., 2012). Over Eurasia, BNU-ESM simulates moderate shortwave cooling effects, while CAM4 simulates opposite warming effects. In South Africa and Amazon regions, both models exhibit strong shortwave cloud cooling effects.  (Rossow and Schiffer, 1999;Rossow and Dueñas, 2004). (b) Zonally averaged total cloud fraction compared with ISCCP D2 retrievals and CLOUDSAT retrievals (L'Ecuyer et al., 2008). (c) Zonally averaged total liquid water path (LWP) compared with Special Sensor Microwave/Imager (SSM/I) retrievals (Wentz, 2000(Wentz, , 2013 over oceans.

Surface temperature and precipitation
The mean observed and modeled climatological annual cycles of surface-air temperature and precipitation for nine representative land regions are shown in Figs. 6 and 7. The most prominent differences from observations in modeled surface-air temperature are a positive bias in Europe of up to 4 • C and negative bias in Eastern Siberia up to nearly 7 • C. In Central Canada, China, and India, the biases are relatively small. In addition to Europe, eight of nine regions exhibit cold biases in annual mean surface-air temperature, and the model generally underestimates the annual temperature over the global land area (excluding Antarctica) by −0.47 • C (−0.28 • C) with an RMSE of 2.25 • C (2.40 • C) compared with CRU TS3.1 (Matsuura and Willmott, MW) data. Compared with two observational precipitation data sets, BNU-ESM has a wet bias at high latitudes. Excessive rainfall during winter seasons in Europe results from too strong mid-latitude westerlies, in particular over the North Atlantic, which carry moist maritime air to the continent. The wet season precipitation in the Amazon exhibits a dry bias, and this tendency extends to August. In Southeastern Asia, the monsoon rainfall in India is more realistic than in China; this is consistent with Sabeerali et al. (2013), who found that the BNU-ESM can simulate a climatologically realistic spatial pattern of June to September precipitation over the Asian summer monsoon region. Globally, BNU-ESM overestimates the annual precipitation over the land (excluding Antarctica) by 0.47 mm day −1 (0.44 mm day −1 ) with a RMSE of 1.42 mm day −1 (1.33 mm day −1 ) compared with CMAP (MW) data. These regional biases may cause dynamic vegetation models in BNU-ESM to produce unrealistic vegetation in affected regions.
In Fig. 8, global surface temperature for the period 1976-2005 of historical simulation is compared with observations. The globally averaged bias is −0.17 • C with a RMSE of 1.83 • C. Over ocean, positive sea surface temperature (SST) biases are seen in the major eastern coastal upwelling regions; probably due to coastal winds that are not favorable for upwelling or underestimation of stratocumulus cloud cover, which is also an issue with other models (e.g., Washington et al., 2000;Roberts et al., 2004;Lin, 2007;Gent et al., 2011). Negative SST biases are mainly found in South Atlantic, South Indian, and subpolar North Pacific Oceans. Another notable negative SST bias is seen in a narrow region associated with East Greenland and Labrador cold currents. In South Atlantic and South Indian Oceans, a tendency for negative SST biases along the northern flank of the Antarctic Circumpolar Current (ACC) are mostly due to insufficient southward transport of heat out of the tropics and a positioning error of the ACC caused by equatorward shift of the westerlies; although there is a small positive bias of the shortwave cloud radiation effect at the cold band between 40 and 50 • S (Fig. 5b). Gupta et al. (2009) noted that relatively small errors in the position of the ACC lead to more obvious biases in the SST. Over continents, the temperature biases are likely consistent with cloud fraction and TOA shortwave cloud forcing (SWCF) biases (Figs. 8b and 5b). Such as the negative temperature bias over South Africa is likely linked to the negative SWCF bias and excessive cloud fraction, and the positive temperature bias over central USA is probably linked to less cloud fraction (Ma et al., 2014). The global average precipitation in BNU-ESM is 0.18 mm day −1 larger over the period of 1979-2005 year ( Fig. 9) than the Global Precipitation Climatology Project (GPCP) data set, which combines surface observations and satellite precipitation data (Adler et al., 2003). While the GPCP data has been claimed to be an underestimate over ocean by Trenberth et al. (2007), the magnitude of tropical precipitation is clearly overestimated by BNU-ESM. In common with many climate models (e.g., Li and Xie, 2014;Lin, 2007), we note a bias in precipitation, characterized by a double Intertropical Convergence Zone (ITCZ) structure over much of the Tropics. This produces excess precipitation over the Northern Hemisphere's ITCZ, Southern Hemisphere's South Pacific convergence zone (SPCZ), the Maritime Continent and the tropical Indian Ocean, together with insufficient precipitation over the equatorial Pacific. BNU-ESM displays the characteristic pattern of the double ITCZ problem with too much precipitation in the central Pacific near 5 • S and too little precipitation in the west and central Pacific between 15 and 30 • S which is similar to CCSM4 (Gent et al., 2011). BNU-ESM underestimates precipitation at 5 • N latitude but overestimates it along the 5 • S parallel in the tropical Atlantic. Compared with observations, the BNU-ESM develops too weak a latitudinal asymmetry in tropical precipitation and SST over the eastern Pacific and Atlantic Oceans. The negative precipitation bias in the South and Northwest Atlantic is closely associated with local negative SST biases (Fig. 8). The band of excessive precipitation over the Southern Ocean between the southernmost of Southern Africa (about at 35 • S, 30 • E) to southwest of Australian is consistent with the spatial pattern of warm SST biases and is along the northern flank of a cold SST bias, which probably produces more convective precipitation. Over continents, there is excessive precipitation in India, northern China, western USA, South Africa and west coast of South America, and less precipitation in southern China and Amazon.
The frequency and intensity of precipitation in the model is highly dependent on the formulation of the convection parameterization (Wilcox and Donner, 2007). Figure 10 shows frequency versus daily precipitation rate over land in the tropics between 20 • N and 20 • S, and compared with the observational estimates from the GPCP 1-degree daily data set  and the Tropical Rainfall Measuring Mission (TRMM) satellite observations (Kummerow et al., 2000). It is clear that BNU-ESM produces a realistic number of precipitation events at a wide range of precipitation rates, although the model has a tendency to underestimate extreme precipitation events (over 50 mm day −1 ). We note that CCSM4 also produces similar precipitation characteristics at 1 and 2 • resolutions (Gent et al., 2011).

Tropical Pacific SST
The tropical Pacific SST is closely associated with the El Niño-Southern Oscillation (ENSO), and exerts a strong influence on the East Asian monsoon Li et al., 2010). Figure 11 shows the 20th century mean and annual cycle of SSTs along the equator averaged between 2 • S and 2 • N in the Pacific Oceans from HadISST observations and the BNU-ESM historical run. The modeled mean SST is colder by about 0.4 • C than the observations over most of the western Pacific and by nearly 1.3 • C over the eastern basin, while warmer than reality at both the western and eastern boundaries of the Pacific (Fig. 11a). These biases are caused by the strong easterly winds in the central and western Pacific and weaker zonal wind at the equatorial boundaries of the Pacific, which result in cold and warm SST biases through enhanced or weakened Ekman pumping in these regions. The different cold SST biases in the central eastern Pacific along the equator result in a stronger equatorial westward SST gradient than observed. In terms of seasonal variation, the observations show a dominant annual cycle in SST in the eastern Pacific Ocean, with anomaly patterns propagating westward across the central Pacific (Fig. 11b). BNU-ESM reasonably reproduces features of the annual cycle structure in the eastern Pacific (Fig. 11c); such as its transition phases and the amplitude and the position of the cold tongue, but the warm season peak is 1 month later in the model than in observations. The westward propagation of positive SST anomaly patterns in BNU-ESM is at about the correct speed between April and November, with 0.5 • C seasonal warming extending to a little west of 160 • W while the observed anomaly remains east of 160 • W. On the other hand, the observed 0.5 • C seasonal cooling near the dateline in March is not seen in the model. The semiannual cycle in SST that dominates in the western Pacific in the HadISST observations is also reasonably simulated in BNU-ESM.

Sea ice extent
Sea ice has long been recognized as a critical aspect of the global heat balance. Unrealistic simulation of sea ice usually exposes deficiencies in both atmospheric and oceanic forcing (e.g., Losch et al., 2010). The observational data used to evaluate the BNU-ESM is monthly climatological sea ice concentrations from the Special Sensor Microwave Imager (SSM/I) data set (Comiso, 1999), obtained from the National Snow and Ice Data Center (NSIDC). We also use the NSIDC's Sea Ice Index (Fetterer et al., 2002), which contains monthly values of sea ice extent and sea ice area. Figure 12 shows the climatological sea ice concentration in the Arctic and Antarctica for the period 1979-2005 of BNU-ESM historical simulation, and the solid black lines are the 15 % mean concentration values from SSM/I satellite observations. The sea ice extent is overestimated in March (Fig. 12a) and slightly underestimated in September (Fig. 12b) following the summer in the Northern Hemisphere (the average mean sea ice extents of March and September are 18.46 and 5.87 million km 2 , while the NSIDC sea ice extents for the same periods are 15.48 and 6.67 million km 2 .). In the Southern Hemisphere both March (Fig. 12c) and September (Fig. 12d) extents are overestimated (the average mean sea ice extents of March and September are 4.96 and 25.94 million km 2 , while the NSIDC sea ice extents are   4.02 and 18.45 million km 2 ). The excessive sea ice extent following the winter in the Northern Hemisphere is mostly due to too much sea ice in the Labrador Sea, Bering Sea, Sea of Okhotsk and adjacent North Pacific. The modeled geographic distribution of ice in the Northern Hemisphere is close to observations in summer. In the Southern Hemisphere, the main overestimation in summer is in Weddell Sea. The too extensive sea ice simulated in both hemispheres is consistent with the cold SST bias found in corresponding areas (Fig. 8). The simulated atmospheric fields are at least partly responsible for the Southern Hemisphere sea ice bias. One notable bias is that the annual average zonal wind stress from about 35 to 55 • S latitudes over ocean is 23.2 % stronger compared with ERA-Interim reanalysis and 42.8 % stronger compared with NCEP reanalysis, which likely inhibits sufficient southward transport of heat, and contributes to cold surface temperatures that are directly linked to a biased ice extent.
In terms of seasonal cycle of sea ice extent, the simulated Arctic sea ice extent for the period 1980-1999 is within the range of 42 CMIP5 models reported by Flato et al. (2013). In Antarctica, BNU-ESM estimates reasonable sea ice extents for February, but overestimates them in September (26 million km 2 ), which is somewhat above the range of 42 CMIP5 models. BNU-ESM and CCSM/CESM adopt similar sea ice schemes, and both models can simulate both the  (Comiso, 1999).
September Arctic sea ice extent and the rate of Arctic sea ice decline over recent decades better than many other CMIP5 models (Liu et al., 2013). While for Antarctica BNU-ESM and CCSM both have a tendency to overestimate sea ice extent.

Ocean meridional overturning circulation
The meridional overturning circulation (MOC) of the global ocean is a system of surface and deep currents encompassing all ocean basins. It transports large amounts of water, heat, salt, carbon, nutrients and other substances around the globe, and is quite important for the chemical and biological properties of the ocean. The Atlantic MOC (AMOC) is an important part of the system and is responsible for a considerable part of northward oceanic heat transport. Figure 13 shows 30 year means of the global MOC and the AMOC over the 1976-2005 period of the BNU-ESM historical run; the overall patterns and positions of cells, water masses, and overturning are similar to observed patterns (Lumpkin and Speer, 2007). above the estimate of 18.7 ± 4.8 Sv for the AMOC strength at the same latitude found by the RAPID/MOCHA monitoring array for the years 2004-2011 (Rayner et al., 2011). Over the historical simulation period , the maximum annual mean AMOC strength at 26.5 • N decreases 12.6 % from 26.9 to 23.5 Sv. The BNU-ESM global MOC possesses a strong Deacon cell of about 40 Sv between 60 and 45 • S, which penetrates to 4 km depth and is a result of increased zonal wind stress driving the ocean. The mean transport of the Antarctic Circumpolar Current (ACC) through Drake Passage is about 101.7 Sv. This is less than the measured value of 134 ± 11 Sv (Cunningham et al., 2003) and at the low end of the range of 90-264 Sv from 23 CMIP5 models (Meijers et al., 2012). One reason for weaker ACC transport through the Drake Passage is that the model-simulated westerly wind stress maximum is shifted equatorward. The mean zonal wind stress over ocean is 26 % lower than ERA-Interim reanalysis products at the latitude of the Drake Passage. Antarctic Bottom Water (AABW) is located north of 50 • S at depths greater than 3.5 km, and the deep MOC in the Southern Hemisphere is about 4 Sv and weak compared with estimates of 8-9.5 Sv from observations (Orsi et al., 1999). 6 Climate variability

Tropical intraseasonal oscillation
The dominant component of the tropical intraseasonal oscillation (ISO) is the Madden-Julian Oscillation (MJO) Julian, 1971, 1972) which affects tropical deep convection and rainfall patterns. During the boreal winter an eastward propagating component affects rainfall over the tropics, while during the boreal summer a northward propagating ISO affects much of southern Asia (e.g., Krishnamurti and Subrahmanyam, 1982;Lau and Chan, 1986;Annamalai and Sperber, 2005;Yang et al., 2008). The MJO plays the prominent role in tropical climate variability, but is still poorly represented in climate models (Lin et al., 2006;Kim et al., 2009;Xavier et al., 2010;Lau and Waliser, 2012;Sperber and Kim, 2012). Here, we adopt the set of community diagnostics developed by the CLIVAR MJO Working Group to examine simulated MJO characteristics. In BNU-ESM, the winter eastward propagation is well detectable in zonal winds at 850 hPa (U850) over a region from the maritime continent to the western Pacific, but is absent over the Indian Ocean and not evident in precipitation ( Fig. 14a and b). Meanwhile, the northward propagation in summer can be realistically simulated; particularly in the off-equatorial region from 5 to 20 • N ( Fig. 14c and d). The quadrature relationship between precipitation and U850 is also well reproduced in northward propagation signals, consistent with observations. The observed MJO (Fig. 15a) exhibits peak power at zonal wavenumber 1 at a period of 30-80 days in both boreal winter and summer (e.g., Weickmann et al., 1985;Kiladis and Weickmann, 1992;Zhang et al., 2006). The power spectrum of BNU-ESM shows that the zonal wave number power distribution is well captured during boreal winter (Fig. 15b); but the eastward propagating power tends to be concentrated at lower than observed frequencies (periods > 80 days). The power density for westward propagation is overestimated, and consequently the east-west ratio of MJO spectral power is smaller than observed. As with BNU-ESM, the power spectra maximum produced by CCSM3.5 using its default convection parameterization is also greater than 80 days , while spectra computed by Zhang and Mu (2005b) for CCM3 adopting the same convection parameterization scheme as BNU-ESM peaks at approximately 40 days. These studies suggest that the ability of a climate model to simulate realistic MJO depends not only on its convective parameterization, but also on interactions between convection and other physical processes in the model. BNU-ESM simulation shows a northward propagating mode of precipitation during boreal summer at wavenumber 1 with a maximum variance between 30 and 50 days (Fig. 15d), but the northward propagating band is weaker than observed (Fig. 15c). Sabeerali et al. (2013) analyzed the boreal summer ISO of BNU-ESM along with 32 CMIP5 models. They found that BNU-ESM is one of six models which captures the three peak centers of boreal summer ISO variance over the Indian summer monsoon region adequately.
We also compared space-time spectra of daily tropical precipitation from BNU-ESM with observed precipitation estimates from GPCP 1-degree daily data set from 1997 to 2005 using the methodology of Wheeler and Kiladis (1999). Figure 16 shows the results of dividing the symmetric raw spectra by estimates of their background spectra. Kelvin, equatorial Rossby (ER), westward inertia-gravity (WIG) waves and the MJO are readily identified in the observational GPCP symmetric spectra. Signals of convectively coupled Kelvin and ER waves appear in the model, and the spectral signature of the MJO is also represented. In observations there is a clear distinction between eastward power in the MJO range (20 day-80 day) and westward power associated with ER waves. The BNU-ESM model exhibits this distinction to some extent, with the eastward power lying at a constant frequency across all wavenumbers and the westward power lying more along the ER dispersion curves. BNU-ESM represents signals of convectively coupled equatorial waves (CCEWs) similarly as CCSM4 (Hung et al., 2013), such as the equivalent depth of the waves and the low power of WIG waves (Fig. 4 in Hung et al., 2013). The powers of eastward propagating components near the MJO spatial and temporal scale in BNU-ESM are more distinctive than that of their westward propagating counterparts compared with CCSM4 (Hung et al., 2013).

El Niño-Southern Oscillation
The El Niño-Southern Oscillation (ENSO) phenomenon is the dominant mode of climate variability on seasonal to interannual time scales (Zhang and Levitus, 1997;Wang and Picaut, 2004;Zhang et al., 2013). Bellenger et al. (2013)  in December and reaches a minimum (1.21 K) in May. BNU-ESM exhibits realistic timing of the seasonal cycle with one peak and one minimum, but the amplitude is much stronger than in observations. Figure 18 shows the power spectra of the normalized time series of Fig. 17 (the detrended SST anomalies normalized by their long-term standard deviation). The observation based Niño-3.4 index has most power between 3 and 7 years, while both BNU-ESM indices have the most prominent variability between 2 and 5 years with a narrow peak at 3.5 years. On timescales longer than 10 year, the piControl and historical simulations have similar power spectra but less power compared with HadISST observations. The presence of variability in the external forcing during the historical simulation does not induce significant changes in decadal and longer period variability. Another aspect of the BNU-ESM ENSO historical simulation, shown in Fig. 19, is the correlation of monthly mean Niño-3.4 SST anomalies with global SST anomalies compared with that from HadISST observations. The figure shows a realistic but narrower meridional width of the positive correlations in the central and eastern tropical Pacific. A horseshoe pattern of negative correlations in the western tropical Pacific is seen in HadISST but is less pronounced in the model. The positive correlation in the western part of the Indian Ocean is well simulated in BNU-ESM, but the extension of this positive pattern into the Bay of Bengal, Gulf of Thailand and South China Sea is missing from the model. The correlation patterns in the Atlantic Ocean are similar between HadISST and BNU-ESM, but more pronounced in the model.
The Southern Oscillation is the atmospheric component of El Niño. Figure 20 shows the Southern Oscillation Index (SOI) from BNU-ESM compared to observation. The observed SOI is calculated using station data from Darwin and Tahiti. For the model, areal averages of mean sea-level Figure 18. Power spectra of the Niño-3.4 index (the SST anomalies of Fig. 17 normalized with the standard deviation) using the multitaper method (Ghil et al., 2002) with resolution p = 4 and number of tapers t = 7. pressure over 125-135 • E, 17-7 • S and 155-145 • W, 22-12 • S (10 • × 10 • areas centered close to the Darwin and Tahiti stations) are used. The interannual variability in the modeled SOI due to ENSO events is well reproduced and shows the expected negative correlation with Niño-3.4 SST anomalies (Fig. 17). The modeled regression coefficient between monthly deseasonalized SOI and Niño3.4 SST anomalies is −0.52 hPa K −1 while the observed is −1.52 hPa K −1 . Hence, the model underestimates the strength of the atmospheric response to ENSO.

Pacific Decadal Oscillation
Another prominent structure of low-frequency climate variability in the North Pacific, with extensions to the tropical Indo-Pacific, is the Pacific Decadal Oscillation (PDO) (Mantua et al., 1997). PDO and ENSO exhibit similar spatial patterns of SST variability but with different regional emphasis Deser et al., 2007). During the positive (negative) phase of PDO, waters in the east tropical Pacific and along the North American west coast are anomalously warm (cool) while waters in the northern, western, and southern Pacific are colder (warmer) than normal. Coupled climate models can simulate some aspects of PDO, although linkages between the tropical and North Pacific are usually weaker than observed (Stoner et al., 2009;Furtado et al., 2011). Figure 21 shows the regression maps of monthly SST anomalies upon the normalized leading principal component time series of monthly SST anomalies over the North Pacific domain (20-40 • N). The first empirical orthogonal function (EOF) mode of BNU-ESM and HadISST observations explains 22.4 and 25.8 % variance respectively. BNU-ESM exhibits generally realistic PDO spatial patterns and its connections to the tropical Pacific are of comparative strength with respect to HadISST observations, but with a narrower meridional extent in the tropical Pacific region. The maximum amplitude of the negative SST anomalies in the North Pacific shifts a little too far west, to the east of Japan, rather than in the central basin. Figure 22 shows time series of the normalized first EOF mode of SST anomalies of BNU-ESM and HadISST observations over the North Pacific domain. It is evident that both patterns show prominent decadal variability. Carbon flux components are hard to measure directly, presenting a challenge in evaluating the model performance. Global products for land gross primary production (GPP) and net primary production (NPP) exist but are model-based and have large uncertainties (Anav et al., 2013;Ito, 2011). Figure 23 shows regional averages of monthly land gross primary production (GPP) for BNU-ESM compared with FLUXNET-MTE estimates . BNU-ESM replicates the annual cycle of GPP in arctic, mid-latitudes, and tropical regions, but the model has a tendency for underestimation during boreal summer, especially over Alaska, the eastern USA, and Europe. Differences between the estimates from our model and those from FLUXNET-MTE may be caused both by differences in the near surface climatology and land cover characteristics, as BNU-ESM dynamically simulates vegetation characteristics as a function of climate and atmospheric CO 2 concentration. In Alaska, the model simulates more C 3 arctic grass and less boreal shrub compared with the observed International Geosphere-Biosphere Programme (IGBP) vegetation distribution (not shown). While in Europe, although the model simulates more broadleaf deciduous temperate tree cover and less grassland, the biased high temperature and low precipitation during boreal summer suppress GPP significantly. In the Amazon, the model simulates a reasonable vegetation distribution of broadleaf and evergreen tropical trees, but the wet season precipitation suffers a dry bias until August (Fig. 7), and the model systematically underestimates GPP. The interannual variability of the GPP estimated by the model is larger than the observational estimates from FLUXNET-MTE and this may be connected with the stronger interannual variability of the physical fields.
The global terrestrial GPP simulated in the BNU-ESM is 106.3 Pg C yr −1 over the period 1986-2005. Various studies estimated the global terrestrial GPP to be about 120 ± 6 Pg C yr −1 over similar periods (Sabine et al., 2004;Beer et al., 2010;Jung et al., 2011). However, these are well below the range of 150-175 Pg C yr −1 from recent observational estimates (Welp et al., 2011). The global simulated NPP over the period 1986-2005 is 49 Pg C yr −1 , which is consistent with the range of 42-70 Pg C yr −1 from earlier studies (Schimel et al., 2001;Gruber et al., 2004;Zhao et al., 2005;Ito, 2011). Net biosphere production (NBP) simulated in the model for the 1990s and 2000-2005 are 1.6 and 1.4 Pg C yr −1 , which is also consistent with estimates of 1.5 ± 0.8 and 1.1 ± 0.8 Pg C yr −1 , respectively reported by Ciais et al. (2013).

Soil organic carbon
Soil organic carbon is a large component of the carbon cycle that can participate in climate change feedbacks, particularly on decadal and centennial timescales (Todd-Brown et al., 2013). The amount of soil organic carbon simulated by models is strongly dependent on their design, especially the number of soil-carbon pools, turnover rate of decomposition and their response to soil moisture and temperature change. Figure 24a, b show the distribution of global soil organic carbon content, including litter, from BNU-ESM compared with the most recent high-resolution observation-based Harmonized World Soil Database (HWSD; FAO/IIASA/ISRIC/ISSCAS/JRC, 2012). The HWSD data provides soil-carbon estimates for topsoil (0-30 cm) and subsoil (30-100 cm) at 30 arc-second resolution. Overall, the ecosystem carbon content follows the precipitation and temperature distribution (Figs. 8 and 9). The BNU-ESM model can capture the large store of soil organic carbon in the boreal and tundra regions of Eurasia and North America, and the small storage in tropical and extra-tropical regions (Fig. 24b). The model underestimates soil-carbon density in the upper 1 m globally compared with the HWSD (Fig. 24a), especially in boreal regions. Soil carbon is overestimated in the model on the Tibetan plateau, because the coarse horizontal resolution does not correctly represent the rugged terrain and overestimates vegetation cover.
The total simulated soil organic carbon, including litter, is 700 Pg C for the period 1986-2005, is well below the 1260 Pg C (with a 95 % confidence interval of 890-1660 Pg C) estimated from HWSD data (Todd-Brown et al., 2013), and 1502 Pg C estimated by Jobbágy and Jackson (2000) for the upper 1 m of soil. However, there is still considerable uncertainty for those observation-based estimates because of limited numbers of soil profiles with organic carbon analyses (Tarnocai et al., 2009). In addition, the soil-carbon sub-model of BNU-ESM is not yet designed to simulate the large carbon accumulations in organic peat soils, or the stocks and dynamics of organic matter in permafrost, a common failure of many CMIP5 models. It is thus to be expected that simulations without these processes underestimate the global soil organic carbon stock. Especially, the temperature sensitivity of soil-carbon decomposition is described by the Q 10 equation (Lloyd and Taylor, 1994) in BNU-ESM, and the environmental controls of moisture and temperature are diagnosed at 0.25 m depth. In Fig. 24c, the zonally averaged soil-carbon density from BNU-ESM is compared with those from HWSD and IGBP-DIS for upper 0.3 m and upper 1.0 m depth ranges. The model simulates substantially less soil carbon than those from the HWSD and IGBP-DIS for the upper 1.0 m, but agrees much better with upper 0.3 m soil-carbon density estimates on magnitude and latitudinal gradients.

Summary and discussion
In this study, the BNU-ESM is described and results for the CMIP5 pre-industrial and historical simulations are evaluated in terms of climatology and climate variability. The climatological annual cycles of surface-air temperature and precipitation generally agree with observations, but with the annual temperature underestimated and the annual precipitation overestimated over global land areas (excluding Antarctica). The sea ice extent of both polar regions agrees better   with the observations in summer seasons than in winter seasons, and the model has a tendency to have excessive ice extent during winter seasons. The global and Atlantic ocean meridional overturning circulation patterns are similar to those observed. With respect to climate variability, BNU-ESM captures some features of tropical intraseasonal oscillation such as the quadrature relationship between precipitation and zonal wind in the northward propagation direction. The MJO signal in large-scale circulation (U850) is not as well simulated as it is in convection (precipitation), but the northward and eastward propagating motions are both weaker than observed. The annual cycle patterns of tropical equatorial Pacific SST, the periods of ENSO, and the leading EOF mode of PDO in the historical simulation are reasonably well simulated. As BNU-ESM has similarities and some heritage in common with CCSM4, in particular for the atmosphere, land and sea ice components, many characteristics in BNU-ESM are probably shared by CCSM4, such as some notable surface climate biases over land (Lawrence et al., 2012) and the dipole precipitation bias in the Indian Ocean.
BNU-ESM has significant biases that need to be improved, such as the tropical precipitation bias over ocean related to the double ITCZ that has long been a problem among many climate models (Lin, 2007). Note that BNU-ESM uses the revised Zhang-McFarlane scheme on deep convection (Zhang, 2002;Zhang and Mu, 2005a), and CCSM4 also uses a re-  (Richter and Rasch, 2008;Neale et al., 2008). It turns out that neither of them eliminates the double ITCZ problem (Gent et al., 2011), so further parameterization improvements are certainly required. Land surface-air temperature simulated for the last few decades of the 20th century exhibit a mean bias greater than 2 • C over significant regions compared with observations, which also shows room for further improvements. Another related discrepancy is that modeled temperatures increase significantly during the last few years of the historical simulation relative to observations (not shown). This is very likely related to the lack of indirect aerosol effects in the atmospheric component (e.g., Gent et al., 2011), and we note that NorESM, which is also based on CCSM4, but which includes indirect of aerosol effects, does not exhibit similar problems (Bentsen et al., 2013). The positive SST biases prevailing at major coastal upwelling regions are clearly related with the relatively coarse horizontal resolution used by the atmospheric component. According to Gent et al. (2010), the most important factor for SST improvements in CCSM3.5 is the finer resolution and better representation of topography, which produces stronger upwelling and favorable winds right along the model coasts rather than being located somewhat offshore. The cold biases in mean SST along the equator in the Pacific Ocean have several causes. One is the stronger easterly winds on the equator which result in stronger equatorial upwelling; another may be weaker activity of tropical instability waves in the ocean. The ocean component MOM4p1 uses the horizontal anisotropic friction scheme from Large et al. (2001), which induces more frictional dissipation and prohibits vigorous tropical instability wave activity (Wittenberg et al., 2006). Stronger activity of tropical instability waves could prevent the cold tongue water from cooling down by mixing with the warm off-equatorial water Menkes et al., 2006;Seo et al., 2006;Zhang and Busalacchi, 2008). The negative SST bias in the southern ocean and excessive sea ice extent in the Antarctic suggest a need to correct the wind stress field to ensure sufficient southern ocean heat transport and proper ocean gyre boundaries.
The strength and frequency of ESNO variability in BNU-ESM highlights potential improvements. The model has a robust ENSO with an irregular oscillation between 2 and 5 years and a peak at about 3.5 years, whereas the HadISST observations show an oscillation between 3 and 7 years. The seasonal phase locking feature of ENSO is well captured in the model, although the standard deviation of Niño-3.4 SST anomalies from the historical simulation is significantly large than in the observations. The causes of biases in ENSO occurrence and amplitude in BNU-ESM may involve many different physical processes and feedbacks. Because of the dominant role of the atmospheric component in setting ENSO characteristics (Schneider, 2002;Guilyardi et al., 2004;Kim et al., 2008;Neale et al., 2008;Wu and Kirtman, 2007;Sun et al., 2009), previous studies have diagnosed the dynamical Bjerknes feedback (Bjerknes, 1969;Neelin and Djikstra, 1995) and the heat flux feedback (Waliser et al., 1994;Jin et al., 2006) during ENSO. Bellenger et al. (2013) found that BNU-ESM underestimates both the positive Bjerknes and the negative heat flux feedbacks by about 45 and 50 % respectively, which could be the major causes of the ENSO biases in the model. This also raises the importance of further improvements on the deep convection parameterization scheme, as the representation of deep convection is central in defining both the dynamical and the heat flux atmospheric feedbacks (Guilyardi et al., 2009). Another possible cause for the excessive ENSO amplitude is the lack of a sufficient surface heat flux damping of SST anomalies in the model, as weaker heat flux damping tends to destabilize and amplify ENSO (Wittenberg, 2002;Wittenberg et al., 2006). Further studies on these topics are warranted.
Despite the drawbacks of the model in simulating some details of the climate system, BNU-ESM has proven to be a useful modelling tool, and is being actively used by many researchers in prognostic simulations for both anthropogenic and geoengineering forcing scenarios. The BNU-ESM represents an addition to the diversity of earth system simulators,

D. Ji et al.: Description and basic evaluation of BNU-ESM
and currently is evolving in many respects. As global biogeochemical cycles are recognized as being evermore significant in mediating global climate change, improvements of BNU-ESM are underway in the terrestrial and marine biogeochemistry schemes. On terrestrial biogeochemistry, the LPJ-DyN based carbon-nitrogen interaction scheme (Xu and Prentice, 2008) will be evaluated and activated in the future. The soilcarbon scheme will be further improved to simulate the large carbon accumulations in organic peat soils, the stocks and dynamics of organic matter in permafrost. A dynamic marine ecosystem scheme will replace the current iBGC module, the new marine ecosystem scheme has improved parameterizations of dissolved organic materials and detritus , a phytoplankton dynamic module that produces a variable of carbon to chlorophyll ratio (Wang et al., 2009a), and refined nitrogen regeneration pathways (Wang et al., 2009b). Additionally, a three-dimensional canopy radiative transfer model (Yuan et al., 2014) will be adopted to replace the traditional one-dimensional two-stream approximation scheme in the land component to calculate terrestrial canopy radiation more realistically. The spatial resolution of the BNU-ESM will be increased to better the simulation of surface physical climate, especially for the atmospheric and land components. Currently a 0.9 • × 1.25 • resolution land and atmosphere components adapted from the finite-volume dynamic core in CAM is being tested. We also note that CAM5 has made significant progress, such as correcting well-known cloud biases from CAM3.5 (Kay et al., 2012). Further discussions of how to incorporate these developments from CAM5 into BNU-ESM are underway.

Code availability
Please contact Duoying Ji (E-mail: duoyingji@bnu.edu.cn) to obtain the source code of BNU-ESM.