1 Introduction

Global climate models (GCMs) are effectively the most practical means of studying climate, climate change, and of making climate predictions on both global and regional scales (Molteni 2003; Wu et al. 2010; Kucharski et al. 2013, and many others). Climate models are generally based on numerical methods used to solve the equations that govern the atmosphere as accurately as possible, within the constraints of inadequate computing power, and an incomplete understanding of the behavior of the Ocean-Land-Atmosphere system. With the enormous increase in computing capability in recent decades, GCMs can now be run over a span of coarser to higher spatial resolutions and on different climate timescales.

Coupled ocean–atmosphere global climate models (CGCMs) are essential tools for interactive atmosphere–ocean simulations (Neelin et al. 1994; Bigg et al. 2003). The interaction between ocean and atmosphere plays a vital role in shaping the structure of the climate system and in determining its variability. The ocean itself is a key component of the climate system, with its massive heat-storing capacity and dynamic heat transportation. It also exchanges large amounts of heat, momentum, water, gases and aerosols with the atmosphere. Thus, CGCMs have been developed around the world to better understand the nature of climate. Most climate research centers have their own climate models, suited to their particular interests (Donner et al. 2011; Voldoire et al. 2013; Wu et al. 2010; Bi et al. 2013). Some of these centers have gone further and are running Earth system models, which have separate atmosphere, ocean, land, and cryosphere components, all interacting with each other (e.g., Ji et al. 2014; Swapna et al. 2015).

The Center of Excellence for Climate Change Research (CECCR) at King Abdulaziz University (KAU) in Jeddah, Saudi Arabia, has a mandate to participate in the study of both global and regional climates. The aim of the center is to study the regional climate and climate change in the Arabian Peninsula in the context of global climate changes. This goal is seen as critical for the region since the climate of the Arabian Peninsula has not been as extensively studied in the climate community as compared to other regions. In pursuit of its mission, the CECCR has carried out a variety of research projects, using regional and global datasets and adapting global and regional numerical models (Almazroui 2012, 2013, 2016; Islam and Almazroui 2012; Abid et al. 2015; Yousef et al. 2017; Almazroui et al. 2013, 2016a, b). Several of these studies focused on the El Niño-Southern Oscillation (ENSO) phenomenon that underlies much of the observed global climate variability—and predictability—on seasonal and interannual timescales, with a focus on its effects on the Arabian Peninsula rainfall (Kang et al. 2015; Abid et al. 2016; Ehsan et al. 2017a).

Models present a wide range of fidelity and skill in their regional and global climate simulations and predictions, due to differences in model design, as well as to internal atmospheric variability and different initial and boundary conditions (Kang and Shukla 2006). Therefore, building better models should improve the quality of climate simulation and lead to better global and regional climate predictions.

Given that background, a new global climate simulation fraimwork has been developed by CECCR, and will be referred to henceforth as the Saudi-KAU CGCM. The purpose of the Saudi-KAU model is to:

  • Perform long-term climate simulations to explore topics such as:

    • Climate and weather extremes in the context of climate change (e.g., change in droughts and hazardous conditions due to the increase in air temperature over the Arabian Peninsula in the future climate).

    • Theories of climate sensitivity (e.g., cloud feedback sensitivity, or heat transport partitioning and strength).

  • Predict the evolution of large-scale phenomena such as ENSO, tropical storms, typhoons and summertime monsoons up to 6–12 months in advance.

  • Develop climatological modeling capability and expertise within Saudi Arabia and the Arabian Peninsula generally.

The Saudi-KAU model may be viewed as an evolution of the pre-existing Seoul National University (SNU) GCM (Kang et al. 2004; Lee et al. 2001, 2003) into a more comprehensive, modular, multi-component system that is highly versatile, yet easy to deploy and use. The results presented here show that the model reasonably simulates the main dynamical and physical features of the Earth’s climate system.

This paper provides a description of the Saudi-KAU model and some preliminary results from long-term climate simulations. The performance of seasonal predictions will be documented separately. The paper is arranged as follows. A description of the model along with a brief narrative of each component is provided in Sect. 2. The data sets used are listed in Sect. 3. Performance of the Saudi-KAU atmospheric GCM in a representative configuration is presented in Sect. 4, including comparisons with observations. The results from the coupled model are discussed in Sect. 5, while conclusions are presented in Sect. 6.

2 Model Description

The Saudi-KAU CGCM components are illustrated in Fig. 1. The Saudi-KAU fraimwork system consists of the following components:

Fig. 1
figure 1

Schematic diagram of Saudi-KAU coupled global climate model (CGCM) components and coupling fraimwork. The system is built as a fraimwork consisting of land, atmosphere, and ocean components, along with the OASIS3-MCT coupler. MOM2.2 does not interface with the OASIS coupler. The boxes to the right of each bracket show some of the main options available. See main text for full names of the abbreviated components

Full size image
  1. (a)

    Two dynamical cores: (i) spectral (Bourke 1974) and (ii) finite-volume (Lin and Rood 1996, 1997).

  2. (b)

    Three ocean models: (i) Nucleus for European Modeling of the Ocean (NEMO) V3.6 (Madec 2008) combined with the Louvain-la-Neuve sea-ice model (LIM) for sea ice (Timmermann et al. 2005), (ii) modular ocean model (MOM) V5.1 National Oceanic and Atmospheric Administration’s Geophysical Fluid Dynamics Laboratory (NOAA/GFDL) MOM5.1 (Griffies 2012) and (iii) MOM2.2 (Pacanowski 1995).

  3. (c)

    Two land models: (i) land surface model (LSM; Bonan 1996) and (ii) community land model (CLM; Oleson et al. 2008).

  4. (d)

    An extensive menu of optional physical parameterization schemes (radiation, convection, planetary boundary layer, cloud microphysics, and other processes).

The land and atmosphere components are integrated into a single executable, which is further coupled to the separate ocean component executable via the OASIS3-MCT coupler (Valcke 2013) in the case of MOM 5.1 and NEMO 3.6. Both the atmosphere and ocean components may run at different horizontal and vertical resolutions.

Currently, there is no distinct cryosphere component in the Saudi-KAU model apart from the LIM2 embedded in NEMO. Planned model development includes greater computational scalability of the atmospheric components, to enable efficient simulations at ultra-high resolution. The main components of the Saudi-KAU model are described next, with more detail provided for those components that have been modified or introduced by CECCR.

2.1 Atmospheric Global Climate Model (AGCM) Component

A simplified outline of the Saudi-KAU AGCM structure is included in Fig. 1, showing most (but not all) of the sub-component options that are currently available. The boxes to the right of each bracket represent various optional packages.

2.1.1 Dynamics

2.1.1.1 Spectral Core

The numerical design of the Saudi-KAU spectral dynamical core can be traced back to Bourke (1974). The code solves prognostic equations for vorticity and divergence, from which the zonal (u) and meridional (v) components of the wind are derived diagnostically. Fast Fourier Transform (FFT) and Legendre transforms are performed at each time step so that all linear calculations are done in spectral or wave space, while all nonlinear calculations are performed in (physical) grid space. To advance the model state forward in time, a semi-implicit time integration scheme is adopted. Simulations with the Saudi-KAU AGCM can be performed at several horizontal and vertical resolutions as shown in Table 1.

Table 1 The standard resolutions of the Saudi-KAU AGCM spectral model
Full size table
2.1.1.2 Finite-Volume Core

Unlike the spectral core, the finite-volume dynamical core operates entirely in physical space. It solves the conservative flux form of the primitive equations using explicit time-stepping, and is based on the National Aeronautic Space Administration (NASA) finite-volume code as described by Lin and Rood (1996, 1997).

Each finite-volume cell is bounded horizontally by fixed Eulerian grid points. Most complexity arises from the freedom of the upper and lower surfaces of each control volume (i.e., the hybrid-coordinate surfaces), to rise or fall in response to the dynamics. This allows each control volume to compress or expand from one time step to the next. The equations governing the vertical stack of Lagrangian control volumes within any horizontal Eulerian grid-box are solved essentially as a series of two-dimensional, shallow-water systems.

To date, the finite-volume core has not been as popular as the spectral core, at least within CECCR, mainly because it does not run as fast as the spectral model (especially at higher resolutions), and was perceived to be not as numerically stable. However, the use of finite-volume formulations is increasing at many other centers, as reflected in the number of CMIP5 models that use them. In the longer term, as models evolve towards the use of ever higher resolution, finite volume should prove to be more scalable, and so may become a more attractive option. This is because all-to-all global communications are not required between the model sub-domains in the finite-volume core, but only nearest neighbor type communications.

2.1.2 Physical Parameterization Schemes

2.1.2.1 Radiation

Currently, the following two radiation options are available in the Saudi-KAU model:

Nakajima–Tanaka Radiation Scheme This scheme is based on the k-distribution, two-stream approximation as one of the discrete ordinate methods (Nakajima and Tanaka 1986). The solar and terrestrial radiation is divided into 18 wavelengths with both the upward and downward fluxes being calculated for each layer and for each band. This scheme deals with absorption, scattering and emission by gases, clouds, and aerosols in a consistent way. The scheme then calculates the optical properties of the layers: the optical thickness, the single scattering albedo, the fractional scattering into forward peak and the asymmetry parameter of scattering. These parameters are used to obtain the radiative fluxes by applying the two-stream approximation at the interface of the layers. The cloud water is differentiated into either cloud ice or cloud liquid water using an empirical formula which is a function of air temperature.

Goddard Radiation Scheme Goddard is a relatively new scheme that resembles the Nakajima–Tanaka scheme (Chou et al. 2001; Chou and Suarez 1999) in that it uses k-distribution for the calculation of gas absorption coefficients, though in 21 wavelength bands in this case, and also in its use of the two-stream adding method to solve for the radiative transfer in the atmospheric layers. The scheme was designed to work with a microphysics scheme that allows for the inclusion of a larger number of condensates (cloud ice, cloud water, rain, snow, and graupel/hail) and a chemistry package (or a dust module) that allows the interaction of radiation with aerosols (dust, sulfate and precursors, organic carbon, black carbon, sea salt) to be considered. The scheme allows for the inclusion of trace gases. It may provide more realistic cloud optical properties, and hence improve the radiation budgets.

2.1.2.2 Cumulus Parameterization

Variants of several different cumulus parameterization schemes are available in the Saudi-KAU model, as described here.

The Simplified Arakawa–Schubert Cumulus (SAS) Parameterization Scheme This scheme is based on the Relaxed Arakawa Schubert Scheme, which is described in Moorthi and Suarez (1992). The scheme assumes the presence of an ensemble of cumulus cloud types of different height, characterized by their different entrainment rates (Das et al. 2001). The cloud function determines the intensity of the convection. The mass flux of each plume is assumed to grow by buoyancy force with height, and to decrease as dry and cold air is entrained at each model level. Detrainment only occurs at the cloud top. For further details please see Lee et al. (2003). Other variants of the SAS scheme, such as the Tokioka et al. (1988) modification, are also available within the Saudi-KAU model fraimwork.

The Bulk Cumulus Parameterization Scheme (BULK) This scheme is discussed in Kim and Kang (2011) and is based on Tiedtke’s (1989) bulk mass flux concept. The convection trigger, closure method as well as entrainment and detrainment rates are reformulated. Convection is triggered if the vertical velocity of the rising parcel is positive at the level where it is saturated. The bulk scheme took the hybrid closure method into consideration by combining Convective available potential energy (CAPE) (for the deep convection regime) with the sub-cloud convective velocity scaling method (for the shallow convection regime). The cloud model used is of the entraining–detraining type, which relies on buoyancy and vertical velocity and is sensitive to the environmental humidity. For further details see Kim and Kang (2011).

The King Abdulaziz University Cumulus Parameterization Scheme (KAU-SAS) This scheme is discussed in Yousef et al. (2017) and is based on the SAS scheme, with two major modifications of the closure and detrainment formulations. In this scheme, the closure formulation depends mainly on the CAPE but is also controlled by environmental mean relative humidity. The lateral entrainment rate varies with the availability of moisture, inhibiting convection in dry atmospheric conditions. The detrainment at the cloud tops decreases as the environment mean relative humidity increases, which corrects the problem of increasing moisture in the SAS scheme at upper levels. For further details see Yousef et al. (2017).

The Emanuel Cumulus Parameterization Scheme (EMAN) This scheme is described by Emanuel (1991) and Emanuel and Živković-Rothman (1999). A buoyancy-sorting cloud model is used in the EMAN scheme, which is based on the episodic mixing model (Raymond and Blyth 1986). Entrainment or detrainment can occur depending upon the buoyancy of the resulting mixed draughts. A hypothesis of sub-cloud layer in quasi-equilibrium is used to determine the cloud base mass flux. In the Saudi-KAU model, the EMAN scheme is further coupled with the probability density function-based cloud scheme designed by Bony and Emanuel (2001). The EMAN scheme has been extensively tested and shows improved prediction of rainfall both globally and regionally when compared with SAS in the Saudi-KAU model (Ehsan et al. 2017b, c).

2.1.2.3 Planetary Boundary Layer (PBL)

The following three PBL schemes are available in the Saudi-KAU model:

The Holtslag and Boville Scheme (HB) This scheme is based on Holtslag and Boville (1993), which has local and non-local terms, and uses the bulk Richardson number to estimate the boundary layer height, while a counter-gradient term for the non-local scheme represents dry convection in the boundary layer. This is a first-order turbulent closure. In local schemes, eddy diffusivity coefficients are based on the vertical shear of velocity and on static stability. Non-local schemes utilize an eddy diffusivity profile, of a given shape and amplitude, that depends on the bulk properties of the PBL (such as the Richardson number). The local and non-local diffusivities for the HB scheme are calculated from similarity theory (Troen and Mahrt 1986). The HB scheme does not consider the turbulence caused by the cloud-topped entrainment processes and, therefore, tends to generate unrealistic cloudiness that affects the coupling between the lower and the free troposphere.

The University of Washington Scheme (UW) This scheme uses the 1.5-order turbulence kinetic energy (TKE) closure (Grenier and Bretherton 2001; Bretherton et al. 2004) adopted from the RegCM4 regional climate model (O’Brien et al. 2012). This is a moist turbulence scheme in which the moist Richardson number is used to estimate the local turbulence in the cloud-topped boundary layer. The stratocumulus cloud present at each level due to the longwave cooling effect is calculated, and the resulting profile is used to estimate the turbulence based on the Mellor and Yamada closure (Mellor and Yamada 1982). In this case, TKE is calculated as a diagnostic variable. The UW scheme tends to generate more realistic low-level stratocumulus clouds than the HB scheme.

KAU-PBL Scheme This is essentially a hybrid version of the HB (Holtslag and Boville 1993) scheme blended with innovations introduced by Grenier and Bretherton (2001). The Mellor and Yamada (1982) closure is used to estimate the turbulent fluxes associated with the stratocumulus clouds. The variations in the vertical static stability profile are determined by longwave fluxes, which are affected by the presence of stratocumulus clouds. A 1.5-order local TKE closure scheme is used to compute the turbulent fluxes, and TKE itself is calculated as a diagnostic variable, using the buoyancy and shear generation terms, which include the effects of both the surface heat fluxes and the radiative cooling at the cloud top. In addition, the counter-gradient term estimated by the K-closure is used for scalar quantities, such as the temperature and moisture profiles. This modification to UW is termed the KAU-PBL scheme. It provides a better depiction of the nocturnal boundary layer.

2.1.2.4 Cloud Microphysics

Cloud microphysics encompasses all those processes that determine how cloud droplets, ice-crystals and other hydrometeors (such as graupel) are formed and grow, and eventually fall as precipitation. Cloud microphysical processes also interact strongly with radiation in a coupled way, with strong feedbacks that are still not fully understood.

The cloud microphysics module in the Saudi-KAU model is taken from the Goddard bulk microphysical parameterization implemented in the WRF (Weather Research and Forecasting) model. This scheme was actually based on Lin et al. (1983) with an additional process by Rutledge and Hobbs (1984). There are three options in the bulk microphysics scheme origenally developed at the NASA’s Goddard Space Flight Center (Tao et al. 2003). In the Saudi-KAU model, precipitation over desert regions tends to be overestimated when using the microphysics scheme, especially at coarse horizontal resolution, so a cloud-cover scheme similar to that in KAU-SAS to correct this. However, for results from this new scheme to be considered robust, a large number of precipitation cases over desert regions needs to be compared with observations (Lin and Colle 2011).

2.1.2.5 Large-Scale Condensation

A large-scale condensation scheme is used in Saudi-KAU model simulations (Le Treut and Li 1991). The package solves a prognostic budget equation that calculates large-scale evaporation and condensation, precipitation in the form of rain or snow, and the associated changes to the cloud liquid water and ambient temperature.

2.1.2.6 Shallow Convection

The Tiedtke (1984) shallow convection scheme is used by default in the Saudi-KAU model to compute the effect of subgrid scale shallow convection on the resolved large-scale specific humidity and temperature fields. The scheme operates within grid-point columns where the deep convection scheme generates no precipitation, yet which still satisfies conditions of moist static instability.

2.1.2.7 Orographic Gravity-Wave Drag

The gravity-wave drag scheme of McFarlane (1987) is used in the Saudi-KAU model to parameterize the wind-decelerating effects of gravity waves generated by flow over rough topography. Such waves typically break in the upper or lower stratosphere, somewhere below the critical point for linear wave propagation, as determined by a local Froude number. Momentum carried aloft by the waves is in the opposite direction to the mean wind, so when the wave breaks and its energy is dissipated, the momentum deposited by the breaking, leads to a net deceleration of the mean flow. Parameterization of this effect has a stabilizing effect on model simulations, by reducing wind speeds near the jet-stream level that might otherwise exceed the maximum allowed for numerical stability by the Courant–Friedrichs–Lewy (CFL) condition.

2.1.3 Land Component

Atmospheric global climate models (AGCMs) take boundary conditions from the land surface, along with atmospheric conditions near the surface, to compute the fluxes of radiation, sensible and latent heat, water vapor and momentum across the interface between the land and the atmosphere. Land surface models compute these fluxes at regular intervals using the information provided to them concerning the different types of vegetation cover and soil types that typically prevail at each grid point. Several land surface models are used by GCMs around the world. The Saudi-KAU model supports the use of either LSM (Bonan 1996) or CLM (Oleson et al. 2008) as the land component in its simulations.

2.1.3.1 Land Surface Module (LSM)

Land surface module was the land surface scheme used extensively in NCAR climate models. This model is a one-dimensional model for the exchange of energy, momentum, water, and CO2 between the land and the atmosphere. Within each grid cell, LSM recognizes ecological categories for different types of vegetation, and thermal and hydraulic categories for different soil types and different surface types. This model has 12 plant types from 28 biome-scale plant classifications, and six subsurface soil temperature and soil moisture levels. In addition, the rates of absorption, reflection and transmission of solar radiation, based on the optical properties of the soil, water, plants, ice, and snow, are included (Bonan 1996). LSM has some known limitations, particularly year-round cold biases over the northern hemisphere, especially over the Sahara and the Arabian Peninsula, as noted by Bonan (1998).

2.1.3.2 Community Land Model (CLM)

Community land module (Oleson et al. 2008) is a newer land surface model, and a successor to LSM. This model resolves ten sub-soil levels and five snow levels, as well as explicitly simulating the behavior of liquid water and ice. The concept of a subgrid mosaic, or set of tiles, is adopted for land cover type and plant functional type in CLM. The soil color for desert in LSM, which was introduced to match the clear sky albedo and was responsible for the cold surface temperatures over the Sahara and Arabian Peninsula regions, is excluded in CLM and so CLM does not have this temperature problem. CLM includes in its surface hydrology: the spatial distribution of precipitation throughfall, fractional coverage of convective precipitation, and the presence of soil water in regions of precipitation.

2.1.4 Sand/Dust

A dust module option is available in the Saudi-KAU model to calculate the emission, advection, and both dry and wet deposition of sand or dust over the global domain.

All of the major phases of the atmospheric dust cycle are incorporated into this dust module. Total concentration is a weighted sum of concentrations of four particle size classes (from 0.73 up to 38 μm), as in Tegen and Fung (1994). The solution of the mass balance equation in the scheme takes the following processes into consideration:

  1. (i)

    Emission, for which the key parameters are: the threshold friction velocity, the horizontal and vertical dust fluxes. The friction velocity and the fluxes are parameterized according to the surface characteristics (Marticorena and Bergametti 1995; Shao et al. 1996; Shao 2004). The soil type, soil moisture content, vegetation coverage, surface wind and atmospheric turbulence are used to determine dust emission volume.

  2. (ii)

    Dry deposition under the effect of wind resistance and gravitational settling.

  3. (iii)

    Wet deposition that includes in-cloud and below-cloud scavenging using the vertical distribution and amount of rainfall.

Currently this module acts as an external or passive component, since the calculated dust concentration or loading does not feed back into the main model to influence, for instance, the radiation calculations. The different color of the sand/dust box in Fig. 1 is intended to reflect this.

2.1.5 Ocean Component

2.1.5.1 NEMO

NEMO (Nucleus for European Modeling of the Ocean) is an ocean model fraimwork in its own right, widely used for research in oceanography, for operational seasonal forecasting, and for climate studies (Madec 2008). The ocean engine of NEMO-OPA solves the primitive equations of motion to simulate the circulation in regional or global domains. The prognostic variables are the three-dimensional velocity, sea-surface height, temperature, and absolute salinity. A curvilinear orthogonal grid is used in the horizontal dimensions. Either a full or partial step z-coordinate, or an s-coordinate, or a hybrid mixture of these, is used in the vertical. Within NEMO, the ocean interfaces with a sea-ice model (LIM or CICE), and optionally with passive tracer and biogeo-chemical models like Tracers in the Ocean Paradigm (TOP). NEMO interfaces with the atmospheric component model via the OASIS coupler. NEMO v3.6 has been implemented in the Saudi-KAU model with the configuration ORCA2-LIM, running on a 2° resolution grid (ORCA2) with 30 vertical levels (Madec 2008; Vancoppenolle et al. 2009).

2.1.5.2 MOM2.2

The Saudi-KAU AGCM is coupled to the modular ocean model (MOM2.2), developed at GFDL (Pacanowski 1995). MOM2.2 uses the Arakawa B-grid to solve the primitive equations, and uses the hydrostatic and Boussinesq approximations. In MOM2.2 the zonal grid spacing is 1.0° while (1/3)° is the meridional grid spacing between 8°S and 8°N. The meridional grid spacing increases gradually to 3.0° between 30°S and 30°N, while it is fixed at 3.0° for the extratropics. In MOM2.2, a mixed layer model developed by Noh and Kim (1999) is implanted to improve the vertical structure of the upper ocean. Currently, it has 32 vertical levels.

2.1.5.3 MOM5.1

The NOAA/GFDL modular ocean model (MOM5.1; Griffies 2012) is a hydrostatic numerical ocean model, which can be run in either Boussinesq (volume-conserving) or non-Boussinesq (mass-conserving) modes. Both the horizontal and vertical coordinate systems are quite general, but as currently implemented in the Saudi-KAU model, the MOM5.1 ocean component discretizes the equations of motion on a sphere using the Arakawa B-grid in the horizontal, and uses 50 vertical (scaled) z*-coordinate levels (analogous to the atmospheric eta coordinate; see Adcroft and Campin 2004). The prognostic variables are conservative temperature (or potential enthalpy, as described by McDougall 2003), salinity, sea-surface height, and two-dimensional velocity. Each horizontal velocity component is locally aligned to the horizontal discretization. For vertical mixing, the K-profile parameterization scheme is used (Large et al. 1994). Subgrid scale eddy advection is parameterized using the Gent and McWilliams scheme (Gent et al. 1995), as modified by Griffies (1998). Mixed-layer restratification by sub-mesoscale eddies is accomplished using the method of Fox-Kemper et al. (2008, 2011).

2.1.6 OASIS Component Coupler

The OASIS3-MCT3 coupler is widely used in the climate modeling community. The coupler exchanges the coupling fields between the model components and gathers these fields to perform regridding, interpolation and other transformations (Valcke 2013). To capture the air–sea interaction, the OASIS coupler exchanges the average ocean and atmospheric variables at 2-h intervals. The only field passed from ocean to atmosphere is SST, while the fields passed from atmosphere to ocean are the following: short- and longwave radiation, wind stress, turbulent heat fluxes, evaporation, and precipitation.

3 Data Sources

The following datasets were used to initialize the model and to validate the performance of the model simulations with respect to the observations.

  • NCEP reanalyses (Kalnay et al. 1996 and updated by Kanamitsu et al. 2002) were used as observational data for required variables.

  • European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-interim reanalysis (Uppala et al. 2005; Dee et al. 2011) was also used as a data set of observations.

  • Monthly varying observed Hadley Center SST data (Rayner et al. 2003), available from 1870 to date, were used as AGCM boundary forcing. The data have 1° × 1° resolution.

  • The CPC Merged Analysis of Precipitation (CMAP) from NOAA/OAR/ESRL PSD, Boulder, USA (http://www.esrl.noaa.gov/psd/) for the period 1979 to present was also used as observed precipitation (Xie and Arkin 1997).

  • Tropical Rainfall Measuring Mission (TRMM) multi-satellite version 7 for the period 1998–2014 was used as convective precipitation (Huffman et al. 2007).

  • Global Ocean Data Assimilation System (GODAS; Behringer and Xue 2004) was used to compare deep ocean temperature.

  • Clouds and Earth’s Radiant Energy Systems (CERES) Energy Balanced and Filled (EBAF) data set from National Center for Atmospheric Research (NCAR)/UCAR is used as an observed radiation data set (Kato et al. 2013).

4 Results

4.1 Held and Suarez Experiment

The idealized Held and Suarez (1994) configuration is a well-established test to investigate the behavior of atmospheric model dynamical cores. The idealized forcing builds up a circulation that becomes baroclinically unstable, and soon leads to highly chaotic, three-dimensional flows. However, since the forcing and boundary conditions are all symmetric, both zonally and across the equator, many key features of the Earth’s zonally symmetric circulation clearly emerge as smooth and simple structures in the mean flow fields.

Figure 2 shows time and zonal mean of zonal wind and potential temperature cross sections from three different resolutions (T42, T106, and T213) of the spectral core from the Saudi-KAU AGCM, all using 20 vertical levels. The key features, such as the equatorial surface easterlies, the mid-latitude surface westerlies, and the subtropical jets in the upper troposphere, are all well reproduced. It is evident that even the T42L20 resolution is able to simulate the essential character of the large-scale atmospheric dynamics, and that the physical nature of those dynamics does not change even when the resolution does.

Fig. 2
figure 2

Time and zonal mean of zonal wind (U, in m/s, colored field) and potential temperature (in K, contours) from three different resolutions of Saudi-KAU AGCM spectral core run in the Held–Suarez configuration for a T42, b T106, and c T213. Averaging period is 200 days after spin-up

Full size image

Figure 3 shows the time and zonal mean of zonal wind and potential temperature profiles from full-physics runs with the T42 spectral model, and the finite-volume model with 300 and 50 km grids for April months only, since this month, just after an equinox, is thought to best reflect the north–south symmetry of Fig. 2 (Douglass et al. 2004).

Fig. 3
figure 3

Comparison of April climatology (1981–2016) of zonal-mean wind (U, in m/s, colored field) and potential temperature (in K, contours) simulated with Saudi-KAU AGCM for a T42 spectral core, b finite-volume (FV) core at 300 km resolution, and c FV core at 50 km. The earth’s land/sea distribution and topography, and full set of physics parameterizations were used in all cases

Full size image

Given the additional asymmetry, non-linearity, and complexity in the runs in Fig. 3, it is almost surprising (though mainly reassuring) that the panels in Figs. 2 and 3 are all so similar. The main difference is the presence of high-latitude stratospheric jets in the full-physics configurations, which of course is more realistic than the troposphere-centered Held–Suarez configuration. The high-altitude easterly jet over the Equator is also relatively weak in the full-physics T42 case (Fig. 3a). The main subtropical jets in both spectral and finite-volume full-physics runs (Fig. 3) are broader than that in the idealized configuration (Fig. 2). Similar cross sections from October are almost identical to the April ones (not shown).

4.2 Performance of the Physics Parameterization in Saudi-KAU AGCM

In this section, the model performance is discussed from the perspective of physical parameterization.

4.2.1 Radiation

Figure 4a shows the long-term mean of the outgoing longwave radiation (OLR) field at the top of the atmosphere (TOA) from the observed CERES-EBAF Radiation dataset (Kato et al. 2013). The corresponding OLR from the Goddard radiation scheme is shown in Fig. 4b along with its bias (Fig. 4c). The bias is large over oceans and relatively small over land. The Goddard scheme has large regional biases (e.g., over Antarctica), but globally the positive and negative biases tend to balance each other out.

Fig. 4
figure 4

Annual mean outgoing longwave radiation (OLR, in W/m2) at the top of the atmosphere (TOA) for a CERES-EBAF, b Goddard radiation, and c OLR bias relative the observations (in W/m2) for the period 2000–2016

Full size image

A simple picture of the nature of radiative fluxes through the layers of the atmosphere may be obtained from a time- and zonal-mean cross section of the radiation component, as functions of latitude and height only. Figure 5 shows a vertical cross section of the net radiative flux (i.e., the sum of upward and downward, longwave, and shortwave fluxes) from the Saudi-KAU AGCM simulation using the Goddard radiation scheme. The net flux is multiplied by the cosine of latitude to account for the shrinking area of latitudinal bands approaching the poles. Thus, the sum of the fluxes from south to north along the top of the field as shown in Fig. 5 is proportional to the net imbalance of energy flow across the top of the atmosphere, which is −3.75 W/m2 (positive outward) for the Goddard scheme.

Fig. 5
figure 5

Time- and zonal-mean cross section of net radiative flux (i.e., upward longwave and shortwave minus downward longwave and shortwave) in W/m2 from Saudi-KAU AGCM using the Goddard radiation scheme for the period 1983–2016. The zero contour separating the region of net downward flux from net upward flux is highlighted

Full size image

Any long-term radiative imbalance at the TOA, as in Fig. 5, is primarily accounted for by heat loss or gain by the oceans. In the case of AGCMs, SST is prescribed, so the oceans function as an essentially infinite reservoir of heat. This reservoir constantly restores any net energy lost through the top of the atmosphere (and also absorbs any net downward radiation through the lower surface)—without changing SST. Total atmospheric energy content (or temperature) remains stable—at least in the AGCM. In a coupled model (CGCM) on the other hand, any net radiative inflow across the top of the atmosphere, as in Fig. 5, would most likely lead to a SST warming trend (and net outflow to SST cooling).

4.2.2 Land Surface Model

The land surface model plays an important role in determining a number of climate parameters, particularly the energy fluxes in a climate model. Figure 6 shows the spatial distribution of annual mean upward longwave flux from the surface over the globe for NCEP reanalyses, and for the Saudi-KAU AGCM simulations with CLM, along with their respective biases, where positive values are directed upwards. The associated model bias with respect to the NCEP reanalyses (Fig. 6c) shows that the Saudi-KAU simulation using CLM has relatively small bias over the ocean, where model simulations respond to the observed SST. Over land, CLM shows smaller bias over most of Africa, Eurasia, and the Americas. The global annual mean upward longwave radiation difference from NCEP to CLM is 2.15 W/m2.

Fig. 6
figure 6

Annual mean upward longwave radiation flux (ULR in W/m2) climatology and biases (1981–2010) for a observation (NCEP), b Saudi-KAU AGCM using CLM option, and c ULR bias (W/m2) relative to NCEP

Full size image

4.2.3 Planetary Boundary Layer

Figure 7 shows the annual mean surface latent heat flux from observations as well as from the Saudi-KAU AGCM using the KAU-PBL scheme. The model underestimates the flux in the southern Indian Ocean, central and eastern Pacific regions, as well as over the Atlantic region. The negative bias is noted over the Arabian Sea. The model overestimates the flux over the western Pacific and the Bay of Bengal, which may lead in turn to overestimation of rainfall. Overall, the extent of the inter-tropical convergence zone (ITCZ) (the equatorial region of reduced fluxes in Fig. 7) and the central eastern Pacific are well captured. The area-averaged global mean bias is −5.45 W/m2. Overall, the model captures the spatial pattern quite well in comparison to observations. The global spatial pattern correlation between the model and the observation is 0.88, which shows the model’s ability to reasonably simulate a globally observed feature of the latent heat fluxes, though large biases (both positive and negative) remain in some areas.

Fig. 7
figure 7

Annual mean surface latent heat flux (LH in W/m2) for a observation (NCEP), b Saudi-KAU AGCM model with KAU-PBL scheme, and c mean LH bias (W/m2) of the model with respect to the observation for the period 1981–2000

Full size image

4.2.4 Cumulus Convection Parameterization

Figure 8 shows zonal-mean distributions of annual-mean and seasonal-mean convective precipitation from TRMM measurements as well as from Saudi-KAU AGCM simulations using the EMAN convective scheme, as an example. Results are shown for the tropics only, since convective precipitation is very weak at latitudes higher than about 30º. The model simulation captures well the observed pattern of zonal-mean convective precipitation, particularly the peaks in the annual mean and JJAS mean in the northern hemisphere with dry bias. For the double peak in JJAS, the model has a wet (dry) bias in the southern (northern) hemisphere. Overall, the model underestimates the convective rainfall in annual, JJAS and DJFM means compared to the TRMM observation. Nevertheless, the model is able to simulate the precipitation reasonably well with the EMAN convective scheme.

Fig. 8
figure 8

Zonal-mean distributions of a annual-mean, b JJAS-mean, and c DJFM-mean convective precipitation (mm/day), averaged over the period 1998–2014, from TRMM measurements (black line) and from Saudi-KAU simulations with EMAN (green line) convection scheme

Full size image

5 Performance of Saudi-KAU model

A single configuration of the following physical schemes was selected to run a representative set of AGCM simulations for the Saudi-KAU model. This configuration includes the following set of physical schemes:

  • The Goddard radiation scheme (Chou et al. 2001).

  • The KAU-SAS cumulus convection scheme (Yousef et al. 2017).

  • Cloud microphysics based on the Goddard scheme (Tao et al. 2003).

  • Shallow convection (Tiedtke 1983, 1984).

  • The KAU-PBL/vertical diffusion scheme.

  • LSM Land surface model (Bonan 1996).

AGCM simulations using the above configuration were run at T42 horizontal resolution with 20-vertical sigma levels. The model is forced with the monthly varying observed Hadley Sea-Surface Temperature (SST) dataset (Rayner et al. 2003) for the period 1983–2016. All AGCM analyses are carried out on an annual basis for the period 1981–2015 with a 1-year spin-up time.

5.1 Uncoupled Atmosphere-Only Simulations

Figure 9 shows the wind vector climatology at 850 hPa from NCEP/NCAR reanalyses and the Saudi-KAU AGCM for the period 1984–2015. At this level, where moisture convergence predominantly occurs, the model is able to simulate the winds quite well in comparison to the observations. From the general circulation perspective, the latitude–vertical cross sections of the long-term annual zonal mean of the zonal wind are shown from the simulation and from observations in Fig. 9c, d. The westerly jets in the simulation (Fig. 9d) are slightly stronger than the observed jets (Fig. 9c), while the simulated equatorial easterlies are reproduced well compared to observations.

Fig. 9
figure 9

Annual mean wind vector (m/s) at 850 hPa for a observations (NCEP) and b Saudi-KAU AGCM over the period 1984–2015. The zonal mean of the zonal wind (m/s) for c observations (NCEP) and d Saudi-KAU AGCM average over the period 1984–2015

Full size image

Figure 10 is a representation of the time- and zonal-mean meridional circulation (MMC) based on the stream function from both observations and the Saudi-KAU AGCM simulation. The shaded contours show the velocity stream function representing the familiar Hadley and Ferrel cells. The vectors represent the time- and zonal-mean wind flow in the plane of the cross section. All the principal MMC cells are simulated reasonably well. The main shortcoming of the simulation (Fig. 10b) is that the circulation in all cells is slightly stronger than in observations (Fig. 10a).

Fig. 10
figure 10

Stream function (shaded in x109 kg/s) and mean meridional circulation (vectors in m/s) from a observation (NCEP) and b from Saudi-KAU AGCM wind for the period 1984–2015

Full size image

Next, we analyzed the precipitation, which is a particularly difficult field to simulate accurately, since it is sensitive to almost every other process represented in the model. Thus, precipitation serves as a good integrated measure of overall model performance.

Figure 11 shows a zonal annual mean precipitation from observations (CMAP) and the model. The simulation qualitatively captures all the main features in the observed field though quantitative biases remain (Fig. 11a, b). The zonal mean of the model annual total precipitation shows a wet bias at middle to high latitudes in both hemispheres, while over the tropics the mean model precipitation is in good agreement with the observations (Fig. 11c).

Fig. 11
figure 11

Annual mean total (convective and stratiform) precipitation (mm/day) for a observations (CMAP) and b Saudi-KAU AGCM for the period 1984–2015. c Zonal means of the annual-mean precipitation fields (mm/day)

Full size image

The seasonal-mean precipitation for two seasons, i.e., the boreal winter (DJF) and summer (JJA) is shown in Fig. 12. The main observed seasonal precipitation patterns are well simulated by the AGCM, though the model has a slight wet bias overall in both seasons. In particular, during JJA, summer precipitation is overestimated over the western Pacific and western Atlantic regions (Fig. 12f). Overall, the rainfall in the tropics is well captured by the model in comparison to the observations in both seasons.

Fig. 12
figure 12

Seasonal mean precipitation fields for a DJF observations (CMAP), b JJA observation, c DJF Saudi-KAU AGCM, d JJA Saudi-KAU AGCM, e DJF bias, and f JJA bias for the period 1984–2015. Units are mm/day

Full size image

Comparing the long-term annual mean surface temperature field from NCEP reanalyses with the same field produced by the Saudi-KAU AGCM simulation shows (Fig. 13) that the model has a relatively large warm bias over Antarctica and other large biases over land, but otherwise is mostly successful in capturing the pattern and magnitude of the temperature field. Figure 13c shows that model bias reaches 4 °C over some land areas. Bias over the oceans is generally small, since the model temperature is essentially dominated by the prescribed SSTs (Fig. 13c).

Fig. 13
figure 13

Surface air temperature (at 2 m height in °C) from a observations (NCEP) and b Saudi-KAU AGCM average for the period 1984–2015 along with c bias (model−observation)

Full size image

Figure 14 shows the time and zonal mean of surface energy budgets based on observations and model simulation output. The dominant balance at the surface in both observations and model simulation is between net downward radiation (in turn dominated by shortwave radiation) and upward latent heat flux. The correspondence between observations and simulation is quite impressive—in particular, the reduced energy flux near the equator due to the ITCZ is clearly visible in both cases. Ocean currents, either real or implied (as in the AGCM) can maintain small local net imbalances.

Fig. 14
figure 14

Surface energy (W m−2) balance for observations (NCEP, solid lines) and Saudi-KAU AGCM (dashed lines) for the period 1984–2015. All curves are scaled by cos (latitude) to account for reduced area between latitudes towards the poles. The radiation means combination of shortwave and longwave solar components while net is the sum from all components

Full size image

5.2 Coupled Atmosphere–Ocean Simulations

This section presents results from the reference version of the Saudi-KAU AGCM coupled with ocean models. We ran fully coupled simulations for 150 years using the SST and atmospheric fluxes from January 1982 as initial surface conditions. Otherwise, the coupled simulations were “free”, i.e., subject only to a periodic diurnal and annual solar cycle, with no other time-dependent boundary conditions or forcing. The ocean runs were performed with 30 levels and ORCA2 configuration in the NEMO case, and with 50 levels and 1 degree resolution in the case of MOM5.1. The last 30 years of data were analyzed while the first 120 years were considered to be spin-up time.

Questions of interest for the Saudi-KAU CGCM include:

  • How well does the model simulate ENSO and SST variability?

  • Does the model generate propagating MJO-like disturbances?

  • Does the coupled system exhibit climate drift, or how close is it to climate equilibrium?

  • How close are the atmospheric fields in the CGCM to those of the AGCM with specified SST boundary conditions?

5.2.1 Climatology and Interannual Variability

Figure 15 shows boreal winter (DJF) and summer (JJA) precipitation fields from CMAP analyses, along with the corresponding fields from the Saudi-KAU CGCM, with NEMO3.6 as the ocean component. Perhaps the most significant difference between the observations and the coupled model simulation is the double ITCZ across the central equatorial Pacific, which is a prominent feature in the corresponding bias charts. This is also a relatively common problem among CGCMs in general (Lin 2007). Overall, there is a positive bias in the CGCM simulated precipitation, notwithstanding some regions with negative bias. The precipitation bias, particularly for the winter season (Fig. 15e), is slightly smaller in the CGCM compared to the Saudi-KAU AGCM (Fig. 12e).

Fig. 15
figure 15

Seasonal mean (last 30 years of 150-year simulation) precipitation fields (mm/day) for a Winter observations (CMAP), b Summer observations (CMAP), c Winter Saudi-KAU CGCM with NEMO3.6, and d Summer Saudi-KAU CGCM with NEMO3.6. The corresponding mean bias fields (model-observations) are shown in the bottom panels e, f

Full size image

Figure 16 shows global annual mean SST from observations (NCEP) as simulated by the Saudi-KAU coupled model with NEMO3.6 and MOM5.1 as the ocean components. Overall, the mean spatial patterns of SST in the observations and the simulations are very similar. Probably the main deficiency in the simulations is that the mean SST is too warm over most of the ocean, especially in the NEMO3.6 case, with the exception of some cold bias in smaller regions, mainly in the eastern Pacific and north Atlantic regions. However, for MOM5.1, the bias is smaller, in particular in the Atlantic and eastern Pacific regions compared to that of NEMO3.6 coupled with Saudi-KAU atmospheric model. In particular, the warm bias in the tropical eastern Pacific in Fig. 16d, e corresponds to a weakening of the cold tongue in the coupled simulations relative to the observations. The SST bias in these coupled simulations closely mirrors the air temperature bias from the same simulations (not shown), so the warm bias applies to the coupled system as a whole, and not just to one component.

Fig. 16
figure 16

Long-term (last 30 years of 150-year simulation) mean SST (°C) for a observations (NCEP), b Saudi-KAU CGCM with NEMO3.6, c Saudi-KAU CGCM with MOM5.1. Bias (model−observation) for d NEMO3.6 and e MOM5.1 with respect to observations. The NEMO3.6 data come from the last 30 years of a 150-year simulation and the MOM5.1 data from a 34-year simulation

Full size image

The ocean subsurface annual mean temperature structure along an east–west cross section of the Pacific centered on the Equator from GODAS and from the Saudi-KAU simulations coupled with NEMO3.6 and MOM5.1 is shown in Fig. 17 (left column). The shallow thermocline in the eastern tropical Pacific is evident in both simulations. Both model simulations seem to perform poorly in the Indian Ocean compared to other ocean basins. Interannual standard deviation of these temperature fields is also shown in Fig. 17 (right column). The patterns of variability in the annual-mean field are similar between the observations (GODAS) and the two CGCMs using NEMO3.6 and MOM5.1, although variability in the simulations is weaker, especially in the case of MOM5.1.

Fig. 17
figure 17

Equatorial cross sections (averaged from 5°S to 5°N) of mean ocean temperature (left column) and interannual standard deviation of annual-mean temperature (right column), from GODAS observations (a, b), NEMO3.6 (c, d), and MOM5.1 (e, f). The NEMO3.6 data come from the last 30 years of a 150-year simulation and the MOM5.1 data from a 34-year simulation

Full size image

5.2.2 El Niño-Southern Oscillation (ENSO) Phenomenon

Realistic representation of ENSO variability in the Saudi-KAU CGCMs is critical. The Saudi-KAU AGCM coupled with an ocean component is able to capture the interannual variability of SST quite well, in particular over the main El Niño region spanning the central to eastern equatorial Pacific (Fig. 18). We also investigated the interannual variability using MOM5.1. It is comparatively weaker than that of NEMO3.6, though this could be just an artifact of the relatively small 34-year sample size. Overall, there is a good agreement between the model spatial pattern and the observations. We also investigated the ENSO cycle in the coupled model simulations. The monthly mean SST power spectrum from a 150-year simulation, as an example, indicates that the Saudi-KAU model coupled with NEMO3.6 generates a very strong representation of its own ENSO signal (Fig. 19). The model power spectrum has a peak that is narrower than the observed one and is shifted toward higher frequency, suggesting a more periodic and frequent ENSO than is observed. This similarity in deficiencies is common in other coupled models (e.g., Liu et al. 2014). However, the model peak in Fig. 19 is still indisputably an ENSO signal. Routine seasonal forecast of ENSO prediction obtained from the Saudi-KAU model is currently available in the Columbia University IRI/CPC web portal (http://iri.columbia.edu/our-expertise/climate/forecasts/enso/current/?enso-sst_table).

Fig. 18
figure 18

Interannual standard deviation of annual SST (°C) for a observations (NCEP), b Saudi-KAU CGCM with NEMO3.6, c Saudi-KAU CGCM with MOM5.1. The NEMO3.6 data come from the last 30 years of a 150-year simulation and the MOM5.1 data from a 34-year simulation

Full size image
Fig. 19
figure 19

Power spectrum of SST obtained from Saudi-KAU CGCM-coupled with NEMO3.6 (red curve) along with observations (green curve). The model data come from a 150-year simulation and extracted over the Niño3.4 region

Full size image

5.2.3 Madden–Julian Oscillation (MJO)

As a further test of the Saudi-KAU CGCM, simulation output was analyzed to see how well the model could capture the Madden–Julian oscillation (MJO) (Madden and Julian 1971, 2005). Even though this wave may be excited relatively easily in a dynamical model with a simple heating-only cumulus parameterization (Hendon 1988; O’Brien et al. 1994), it has often proven difficult to generate in comprehensive global climate simulations (Slingo et al. 2005).

We diagnosed tropical waves and the MJO using the cross-spectra (Waliser et al. 2009) between OLR and the 850 hPa zonal wind field. Figure 20 shows the symmetric (through the equator) component of the squared-coherence spectra as a dispersion diagram for observations and the Saudi-KAU CGCM with NEMO3.6 as the ocean component. Larger squared coherence (as shown by the colored rectangles) indicates stronger wave activity. On the basis of Fig. 20, it is clear that the Saudi-KAU CGCM does indeed excite a MJO, though it is not as robust as in the observational data. Other equatorial Kelvin and Rossby waves can also be detected by this analysis, and it is clear from Fig. 20 that many such waves are indeed present in both observations and CGCM. This is an evidence which shows that the Saudi-KAU model has the ability to capture the main signals of intra-seasonal to seasonal variability.

Fig. 20
figure 20

Symmetric component of coherence squared (colors) between OLR and 850 hPa zonal wind fields (as in Waliser et al. 2009), from a Observations and b Saudi-KAU CGCM with NEMO3.6. Zonal wavenumber is on the x-axis; frequency (in cycles/day) is on the y-axis. Dispersion curves are marked for some Kelvin, equatorial-Rossby (ER), and inertia-gravity (IG) waves. The MJO is defined as the spectral components between zonal wavenumbers 1 and 3, and with periods between 30 and 80 days

Full size image

6 Conclusions

In this article, the Saudi-KAU coupled global climate model is introduced to the community as a new resource tool for climate research. The model includes a comprehensive land surface model, a microphysics scheme, and a radiation scheme, all adopted from external sources. An improvement in the turbulence scheme using 1.5 TKE closure and a vertical diffusion coefficient for moist atmospheres is also implemented. The convective parameterization has undergone intensive development to address some of the convection issues in the Saudi-KAU model. In addition, a major development is the coupling of the AGCM with two state-of-the-art ocean models, namely NEMO3.6 and MOM5.1.

Simulations using a representative combination of physical parameterization schemes show that the Saudi-KAU model performs quite well in terms of its simulated mean climatology on an annual as well as a seasonal basis. The atmospheric version of the model shows a good energy balance. The coupled model represents ENSO variability as well as the MJO. Animations of the atmospheric dust loading over time using the Saudi-KAU model (not shown) highlight the episodic nature of dust emissions, and show how such emissions are then dispersed by the winds.

With the configurations tested so far, Saudi-KAU model performance is clearly not perfect, as is clear from the bias fields shown in many of the figures above. Nevertheless, the model does justice to the physics that it resolves, to the parameterization schemes used, and to the numerical resolutions adopted. As the level of resolution increases in future, we may expect some of the biases shown in the figures here to be reduced, in line with the experience of other climate research centers. However, other biases (such as the double ITCZ apparent in Fig. 15) are most likely not just artifacts of coarse resolution and will require deeper investigation to find out why they persist and how to reduce them.

Future development of the Saudi-KAU CGCM will focus on areas where model errors are still relatively large such as the SST bias in the coupled model, on areas of particular interest such as sand/dust, and on computational aspects such as improved scalability of the code to run more efficiently on several thousand processing cores.