How well does MPAS-atmosphere simulate the characteristics of the Botswana High?

The Botswana High is a prominent mid-tropospheric system that modulates rainfall over subtropical southern Africa, but the capability of a global climate model (GCM) to reproduce it remains unknown. This study examines the capability of a GCM with quasi-uniform resolution (Model Prediction Across Scales, hereafter MPAS) in simulating the characteristics of the Botswana High. The MPAS is applied to simulate the global climate at 240 km quasi-uniform resolution over the globe for the period 1980–2010. The model results are validated against gridded observation dataset (Climate Research Unit, CRU), satellite dataset (Global Precipitation Climatology Project, GPCP), and reanalysis datasets (Climate Forecast System Reanalysis, CFSR; the National Oceanic and Atmospheric Administration, NOAA; and the European Centre for Medium-Range Weather Forecasts version 5, ERA5). In general, MPAS replicates all the essential features in the climatology of temperature, rainfall, 500 hPa geopotential height and vertical motion over southern Africa, reproduces the spatial and temporal variation of the Botswana High, and captures the influence of the Botswana High on droughts and deep convections over the sub-continent. In addition, the model reproduces well the anomalies in vertical motion over subtropical southern Africa during +ve and −ve phases of the Botswana High. However, the model struggles to reproduce the precipitation pattern associated with the positive and negative modes of the Botswana High. The results of this study have an application in understanding the characteristics of the Botswana High and in improving MPAS for seasonal forecasting over southern Africa.


Introduction
Subtropical highs are semi-permanent features over the subtropics during summer and are known to influence tropical cyclone tracks (Stowasser et al. 2007;Wu et al. 2005;Kasahara 1959) and play a major role in the formation of the world's subtropical deserts and zones of Mediterranean climate (Miyasaka and Nakamura 2010;Wu et al. 2009;Rodwell and Hoskins 2001). During the austral summer, midlevel (~ 500 hPa) subtropical highs are observed over the Southern Hemisphere subtropical landmasses. These semipermanent anticyclones are typically located over Bolivia in South America, Botswana/Namibia in southern Africa, and over western Australia and are understood to be thermally induced by condensational heating from tropical regions of high rainfall located to their northeast (Reason 2016). The Botswana High, a prominent feature in the mid-levels (~ 500 hPa) over subtropical southern Africa (in austral spring, summer and early autumn), is an example of these semipermanent anticyclones. Reason (2016) showed that the Botswana High forms around August, strengthens throughout spring, and reaches its peak in the late austral summer (January-March). Thereafter, the high becomes weaker in March and finally dissipates around April. Driver and Reason (2017) linked the changes in the Botswana High locations with variation in precipitation over the Congo Basin and the Intertropical Convergence Zone (ITCZ), which lie northeast of the high.
Studies have shown that the Botswana High has a marked interannual variability that correlates with El Niño Southern Oscillation (ENSO; Reason 2016;Driver and Reason 2017). The high tends to be stronger during El Niño (associated summers over southern Africa) and weaker during La Niña events (wet summer). However, the magnitude of the rainfall response is not always proportional to the severity of ENSO (Driver et al. 2018). In addition, the variability of the Botswana High is still high during some neutral phases of ENSO. Furthermore, Driver and Reason (2017) showed that the summer rainfall over southern Africa has a higher correlation with the Botswana High index than with the ENSO index. This suggests that the Botswana High may be useful in predicting southern African rainfall, especially during neutral ENSO years. Hence, a better understanding of the Botswana High variability can give additional insight into the climate variability over southern Africa and foster better seasonal forecasting of summer rainfall over the region.
Several studies have documented the impact of the Botswana High on climate variables in southern Africa (Matarira 1990;Unganai and Mason 2002;Driver and Reason 2017), but there is a dearth of information on the impacts of the Botswana High on droughts over the region. Matarira (1990) and Unganai and Mason (2002) examined the relationship between the Botswana High and rainfall over Zimbabwe and showed that a strong Botswana High was behind the major droughts of 1982-1984 and 1991-1992. These studies used only rainfall to characterise drought (which depends on rainfall and potential evapotranspiration; Abiodun et al. 2018), and its focus was on Zimbabwe. Driver and Reason (2017) showed that the Botswana High variability impacts the frequency of dry spells during summer, the maximum surface air temperature, daily surface air temperature range and days of maximum temperature extremes, all of which are strongly related to drought. However, the direct influence of the Botswana High remained unstudied. Given the devastating impacts of droughts on socio-economic activities in southern Africa (which mostly depend on rain-fed agriculture), it is essential to know how the variability of the Botswana High influences the characteristics of southern African droughts at a regional scale. That is one of the motivations for the present study.
The applicability of Global Climate Models (GCMs) to simulate subtropical highs is well documented in the literature (Paek et al. 2015;Dong et al. 2017;He et al. 2017). For example, Paek et al. (2015) showed that the Central Weather Bureau Atmospheric General Circulation Model (AGCM) gives a realistic simulation of the summer Western Pacific Subtropical High (WNPSH) variability in the high-frequency band (2-3-year) but underestimates it in the low-frequency band (3-5-year). Analysing the simulations of 22 GCMs from the Coupled Model Intercomparison Project 5 (CMIP5), Dong et al. (2017) showed that most CMIP5 models capture the sub-seasonal variability of the WNPSH and reproduce the influence of the high on Sea-level-pressure and 850-hPa wind fields. A similar result was obtained by He et al. (2017) on the ability of 30 CMIP5 models in simulating subtropical anticyclones.
However, no there is no study on the capability of GCMs in simulating the Botswana High.
Most studies that evaluated the capability of GCMs in simulating subtropical highs have focused on GCMs that either use latitude-longitude grid (e.g. He and Zhou 2015) or spectral grids (e.g. Lyu et al. 2017) for horizontal grid discretisation. The limitations of both approaches are well documented in the literature (Washington et al. 2008). For example, with latitude-longitude grid discretisation, the longitude lines are typically clustered at the poles, creating potentially severe Courant-Friedrichs-Lewy (CFL) limitations on the time step and a massive impact on computing time and resources (Zarzycki and Jabłonowski 2014). Although some techniques have been introduced to manage the poles problems, the techniques usually degrade parallel scalability by effectively limiting the model to one-dimensional domain decomposition strategies (Taylor et al. 2008;Prusa 2018). On the other hand, spectral grids discretisation provides a natural spatial filter with a given spectral truncation (e.g. T42); thus the pole problem does not exist; however, the spectral transform method is difficult to implement efficiently on massively parallel computer systems as physical processes (such as radiative forces, convection and precipitation) cannot be included unless the transform method is used (National Center for Atmospheric Research Staff (Eds) 2017; Washington et al. 2008). However, a new GCM called the Model for Prediction Across Scales (hereafter MPAS) overcomes the limitations of both pole clustering and scaling on massively parallel systems through the use of quasi-uniform grids called Spherical Centroidal Voronoi Tessellations (SCVTs) (Du et al. 1999;Ringler et al. 2008). Taylor et al. (2008) illustrated that GCMs with quasi-uniform grids do not require polar filtering scale well on petascale machines and are more computationally economical and efficient than GCMs with latitudelongitude-based grids. Using MPAS, Skamarock et al. (2013) showed the model could produce results similar to the state-of-the-art hydrostatic model at hydrostatic scales (dx > 10 km) as well as nonhydrostatic structures (dx ~ few km) such as convective systems, similar to those produced by nonhydrostatic regional models. While the capability of MPAS to simulate different atmospheric features has been shown (Kramer et al. 2018), no studies have shown how well MPAS reproduces the characteristics of the Botswana High. Hence, the present study aims to examine how well the MPAS model simulates the characteristics of the Botswana High, especially the influence of the high on drought over southern Africa. Section 2 describes the observation data, reanalysis and simulated data used in the study, and the methods employed in the study. Section 3 presents and discusses the results, while Sect. 4 presents the conclusions.

Data
This study analysed observation, reanalysis and model simulation datasets for the period 1980-2010. The observation data include monthly temperature data from the Climatic Research Unit Version 4.03 (Harris et al. 2014) dataset (hereafter CRU) and is accessible from https:// cruda ta. uea. ac. uk/ cru/ data/ hrg/. The CRU dataset covers all land areas (excluding Antarctica) at a 0.5° resolution and is used to assess the skill of the European Centre for Medium-Range Weather Forecasts Reanalysis 5 (ERA5) and MPAS data in simulating the temperature characteristics over southern Africa. Monthly observed rainfall was obtained from the Global Precipitation Climatology Project Version 2.3 dataset (hereafter GPCP; https:// www. esrl. noaa. gov/ psd). The GPCP dataset combines global observations and satellite precipitation data into 2.5° × 2.5° global grids and is used to evaluate the ERA5 and MPAS simulated rainfall over southern Africa.
Observed SSTs data (6 hourly) were retrieved from the Climate Forecast System Reanalysis (CFSR) data products provided by the National Centers for Environmental Prediction (NCEP; https:// rda. ucar. edu/ datas ets/ ds093.1). The SST dataset has a grid resolution of about 0.31° × 0.31° and is used to evaluate how well the model simulates surface temperatures over southern African ocean basins and assesses the correlation between the Botswana High and global SSTs. Monthly SST indices were obtained from www. esrl. noaa. gov/ psd/ gcos_ wgsp/ Times eries. The timeseries data includes the Nino 3.4 and Dipole Mode Index and is used to evaluate the relationship between the Botswana High and SST modes over the Pacific and the Indian Ocean.
Monthly reanalysis temperature, rainfall, geopotential height, outgoing longwave radiation (OLR) and omega data were obtained from the Earth System Research Laboratory 20th Century reanalysis II (hereafter 20C), National Oceanic and Administration (NOAA; https:// www. esrl. noaa. gov/ psd) and the fifth generation European Centre for Medium-Range Weather Forecasts (ECMWF) atmospheric reanalysis of the global climate (ERA5) dataset (https:// clima te. coper nicus. eu/ clima te-reana lysis). Both datasets have regular latitude-longitude grid resolution of 0.25° × 0.25° and were used to evaluate the MPAS model simulation.

Model description and experiment
The model evaluated in this study is MPAS (version 5.2). MPAS is a fully compressible non-hydrostatic global climate model with variable resolution capabilities and uses Advanced Research WRF (ARW) model physics. The MPAS variable resolution uses Voronoi tessellation, which is predominantly composed of hexagons, pentagons and heptagons that are occasionally present (Fig. 1a). The numerical schemes used in MPAS are similar to those in the ARW High area (blue rectangle; 15°-22° E; 20°-25° S). Biases in representation of the topography in the MPAS simulation are shown in contours model with the major differences in that ARW uses rectangular meshes and a hydrostatic vertical coordinate while the fully compressible non-hydrostatic equations in MPAS are derived in terms of a geometric-height vertical coordinate, and the solvers employ a split-explicit time integration scheme which is described in Klemp et al. (2007). The timeintegration scheme uses a 3rd Order Runge-Kutta method with a large time step for meteorologically significant modes and a forward-backwards method with a smaller time step for acoustic modes. For the present study, the MPAS model grid is set to a quasi-uniform resolution of 240 km (10,242 grid cells) over the globe (Fig. 1a). The model simulation uses 41 vertical levels with an integration time step of 960 s and an output rate of 3-hourly data.
The model simulation is initialised with NCEP CFSR data obtainable from https:// rda. ucar. edu/ datas ets/ ds093.1/ and, to give a more realistic simulation, the model was configured to periodically update with 6 hourly CFSR SST and sea-ice data. The model topography was obtained by interpolating the US Geological Survey dataset into the MPAS grid. The model simulation started from the 1st of January 1979 to the 31st December 2011, however, the model results from 1979 were discarded as spin-up and only data from the 1st of January 1980 till the 31st December 2010 was used for the study.

Model validation
To evaluate how well the model simulates the climate over southern Africa, we compared simulated data (e.g. temperature and precipitation) to the observed data. However, the known observation temperature datasets that were accessible have temperature over either land or ocean. Therefore, the CRU data (which has observed surface temperature over the land only) was combined with CFSR-SST (which has observed surface temperature over the ocean) data to have a more complete observed temperature over the study area. In addition to that, we evaluated the model in simulating the climatology of the 500 hPa geopotential height over southern Africa and compared the results to reanalysis data.

Identification of Botswana High
Empirical Orthogonal Functions (EOF) analysis of monthly reanalysis and model data were performed over southern Africa between 1980 and 2010 to extract the Botswana High feature from the 500 hPa geopotential height. EOF analysis was performed for each dataset (i.e. ERA5, 20C and MPAS) over the study period and to verify that the three datasets were showing the same pattern, the reanalysis and model data were combined and another EOF analysis was conducted on the combination of the datasets. The domain for the EOF analysis was chosen for the region lying between 10° E and 40° E and from the equator to 33° S to exclude influence from surrounding oceans and their anticyclones. In addition to that, choosing the southern boundary as 33° S minimises the South Annular Mode impact (Reason 2016).

Performance Measures
We performed a wavelet analysis to evaluate how well the model simulated the Botswana High periodicity compared to reanalysis. We also compared the model and reanalysis periodicities to ENSO periodicities in order to assess the annual variability between the Botswana High and ENSO. The wavelet analysis was performed using monthly SST time series data and Botswana High indices for the ERA5, 20C and MPAS data. The Botswana High indices were computed from the EOF time amplitude functions of the reanalysis and model data. Pearson correlation was used to measure the strength of association as well as the direction of the relationship between the MPAS simulation, reanalysis datasets, and SST data.

Drought Quantification
In this study, the 3-month Standard Precipitation Index (SPI, McKee et al. 1993) and Standard Precipitation-Evapotranspiration indices (SPEI, Vicente-Serrano et al. 2010) (Pearson III type distribution) were used to quantify drought over southern Africa. SPI uses only precipitation to characterise droughts, while SPEI which uses both precipitation and potential evapotranspiration, which may be a better way of characterising drought risks . SPI and SPEI anomalies were calculated from ERA5, 20C reanalysis and MPAS data. For SPEI, the Potential Evapotranspiration (PET) was estimated using the Thornthwaite (1948) method (SPEI_TH), which has the advantage of requiring only the monthly mean temperature and latitude to calculate PET (McEvoy et al. 2012), unlike the Penman-Monteith (SPEI_ PM) which requires more parameters (e.g. solar radiation, wind speed, humidity) for its computations. Despite that, Meque (2015) concluded that both SPEI_TH and SPEI_PM showed a strong correlation and overall good agreement over southern Africa with minor discrepancies.

Climatology over Southern Africa
This section discusses how well MPAS reproduces the climatology over southern Africa with a focus on the spatial distribution of temperature (Fig. 2a, c and e), rainfall (Fig. 2b, d and f) and 500 hPa geopotential height (Fig. 3). Figure 2 shows that MPAS simulates well the spatial distribution of temperature over southern Africa as compared to the observation and ERA5 reanalysis. The MPAS model has a strong spatial correlation (r > 0.9) with the observation, indicating that it reproduces the dominating features in the observed temperature field. It captures the temperature gradient between the tropics and midlatitudes and the relatively cold surface temperatures associated with the Benguela current near the western regions of the continent and the warm surface temperatures associated with the Agulhas Current in the eastern parts of the continent. However, the model simulation shows a cold bias of about 2 °C over much of the Indian Ocean, which may be due to the coarser resolution in the model data (~ 240 km) as compared to the observation (~ 100 km). This bias suggests a weaker response of the MPAS atmosphere to SSTs over the Indian Ocean. The model shows a cold bias over the land, more especially over Botswana, Angola, Zambia and western parts of Tanzania. This bias, which also features in ERA5 results, may be due to how the MPAS resolves the topography over these high elevation areas (see Fig. 1b).
The MPAS simulated rainfall field also has common features with GPCP observation (Fig. 2b, d and f). GPCP shows an area of minimum precipitation over much of the Atlantic Ocean and a band of maximum rainfall (associated with South Indian Convergence Zone (SICZ)) over the Indian Ocean. The SICZ associated rainfall band shows a local maximum over the northern parts of the Indian Ocean (5-15° S; 60-80° E) and Madagascar (about 12 mm/day). MPAS reproduces the tropical rainfall band but with a wet bias of about ± 1 mm/day and fails to capture the local maximum rainfall over Madagascar in comparison to GPCP. The failure of the model to capture the local maximum over Madagascar may be attributed to the overestimation of deep convection over the Mozambique Channel, as the strong convection over the channel can suppress convection over the Malagasy Island, leading to a lack of local maximum rainfall. Additionally, the deep convection over the Mozambique Channel may be due to the dry bias in MPAS simulation over Mozambique and southern Tanzania. Despite that, MPAS shows good agreement with GPCP (r = 0.89) over the southern African continent, where all datasets show more rainfall in the tropics Fig. 2 The climatology of JFM surface temperature (°C) and precipitation over southern Africa  (left and right columns, respectively), as depicted by CRU CFRS and GPCP (top panels), b ERA5 reanalysis (middle panels), and MPAS simulation (bottom panels). Biases in ERA5 data and MPAS simulation (with respect to CRU CFSR and GPCP) are indicated with contours. The corresponding correlation (r) of the spatial patterns are also indicated (6-8 mm/day) and lesser rainfall near the subtropics (about 2-3 mm/day).
MPAS generally agrees with the two reanalysis datasets (ERA5 and 20C) on the spatial pattern of the 500 hPa geopotential height (Fig. 3). All the datasets show areas of high geopotential height over the tropics and low geopotential height poleward of 30°S. However, the MPAS model gives higher values of the 500 hPa geopotential height over much of the study area compared to the ERA5 reanalysis. The MPAS bias is more than 10 m for much of the area, but it is comparable with the difference between the two reanalyses (i.e. reanalysis uncertainty). There are also inconsistencies between the model simulation and the reanalysis over parts of Sub-Saharan Africa and over the Indian Fig. 3 The climatology of 500 hPa geopotential height (m) over southern Africa (1980Africa ( -2010, as depicted by ERA5 a reanalysis, b 20C, and c MPAS simulation. Biases in the 20C data and MPAS simulation (with reference to the ERA5 Reanalysis) are shown with con-tours and the corresponding spatial correlation (r) are also indicated. The blue box shows the Botswana High area (15°-22° E; 20°-25° S) while the cyan box shows the EOF domain (10°-40° E; 0°-33° S) Ocean. For instance, the model shows a high-pressure centre (which resembles the Botswana High) farther east as compared to the reanalysis, which shows the high-pressure centre between Namibia and Botswana's border. Furthermore, MPAS struggles to capture the overall pattern of the 500 hPa geopotential height over the Indian Ocean as compared to the reanalysis. These inconsistencies may be due to the low 240 km resolution in the model simulation. Using a higher resolution may improve the model results. Regardless of that, MPAS features a high spatial correlation (r = 0.98) with the ERA5 reanalysis.

Spatio-temporal distribution of the Botswana High
There is a good agreement between MPAS and reanalysis data on the variability of the pressure gradient over southern Africa as well as the pattern and location of the Botswana High, although with some discrepancies (Fig. 4). In each of the datasets, EOF analysis indicates that mode 1 (EOF1) has a pattern reminiscent of the Botswana High, which is centred between 10° S and 25° S with an area of low pressure south of the ridge. Compared to the actual summer position of the Botswana High, EOF1 is displaced slightly equatorward as in Reason (2018). Despite that, EOF1 in the MPAS model shows an area of high pressure centered between the borders of Botswana, Namibia, Angola and Zambia and accounts for approximately 80% in the variance of the 500 hPa geopotential height in the region, which is similar to ERA5 (85%) and 20C (83%). The 80% variance has also been found by Reason (2018). However, the MPAS high shows a broader centre located farther north in comparison to the reanalysis whilst 20C shows a weaker and narrower centre which is located more eastward compared with the ERA5 and MPAS. Additionally, the combination of the model data and reanalysis EOF analysis (Fig. 4d) shows a similar pattern to Fig. 4a-c,   Fig. 4 The leading factor of EOF analysis (EOF 1) of JFM 500 hPa geopotential height using monthly a ERA5, b 20C reanalysis, c MPAS simulated data, and d the combination of the three datasets for 1980-2010. The numbers in the bottom right denote the percentage of variance explained by the leading factor but more comparable with the ERA5 high and explains 82% variance in the 500 hPa geopotential height over the region.
The time series of EOF1 shows a strong agreement among the three datasets (MPAS, ERA5 and 20C reanalysis) in depicting the interannual variability of the Botswana High. However, the 20C reanalysis exhibits a weaker Botswana High in comparison to ERA5 and the MPAS results (Fig. 5). In addition, even with the combination of the datasets (as in Fig. 4d), the time-series results (figure not shown) are comparable to Fig. 5, indicating that EOF1 shows the same feature in both the model and reanalysis. In agreement with the previous studies by Reason (2016Reason ( , 2018, Fig. 5 shows a strong correlation between ENSO and the Botswana High Index during the study period. All three datasets indicate that the Botswana High was strongest in the years 1983, 1998 and 2010, which is consistent with Reason (2016). In addition, all three cases of the strongest Botswana High occurred during the mature phase of El Nino summers (1983, 1998 and 2010), while the majority of the weakest Botswana Highs occurred during La Nina summers except for 1994, which occurred in a neutral ENSO year. However, MPAS simulation indicates that the weak Botswana Highs occurred in 1985, 1986, 1989, 1999, while ERA5 shows the weak high in 1985, 1986, 1989, 1996and 2000. On the other hand, the 20C showed the highest number of weak highs (1985, 1986, 1989, 1994, 1996, 2000, 2001 and 2008), possibly due to the generally weaker high in the reanalysis. However, the time-series shows that the MPAS results strongly correlated with that of ERA5 (r = 0.83, p < 0.0001), meaning that the MPAS model has a strong skill in simulating the interannual variability of the Botswana High.
A wavelet analysis of EOF1 in the reanalysis and model simulation indicates that the Botswana High is dominated by a strong 4-5-year cycle band before the early 1990s and a weak 2-3-year cycle band afterwards (Fig. 6). All datasets show similarity in the 4-5-year cycle between 1981 and 1989; however, the power spectrum is stronger in the ERA5 than in 20C and MPAS. There are also some discrepancies between the MPAS and reanalysis on the 2-3-year cycle. For example, before the 1990s, the 2-3-year periodicity was weaker in MPAS than the reanalysis but stronger between 1995-2001 and 2007-2010. The dominance of the 4-5-year and 2-3-year periodicities of EOF1 compares well with the ENSO wavelet power spectrum (Fig. 6d) and suggests that the Botswana High variability may be related to forcing from ENSO that also has leading periodicities in these two frequency bands.

The influence of SSTs on the Botswana High
The correlations between the JFM Botswana High index and SST anomalies in different seasons (OND, JFM, and AMJ) indicates a good agreement between MPAS, ERA5 and 20C (Fig. 7). All the datasets feature a high correlation between the Botswana High index and global SST anomalies. During early summer (OND), the most significant correlations occur over the North Atlantic as well as parts of the South-West and North-East Pacific. However, during concurrent summer (JFM), the correlation becomes stronger over the tropical Pacific, tropical Atlantic, South-West Pacific, tropical Indian Ocean, Enderby Plain (60° S, 40° E) as well as over the Kerguelen Plateau (43° S, 72° E). In the Autumn (AMJ), the correlation decreases over   High is highly correlated with the early development phase of ENSO and progressively becomes weaker during the mature phase.  (1983,1998,2010) and −ve Botswana High Years (1989Years ( , 1994Years ( , 2008. The mean and standard deviation used in calculating the anomalies were obtained using the 1980-2010 data To quantify the link between the Botswana High and ENSO, Table 1 shows the coefficient of correlation between the Botswana High and Nino 3.4 seasons. All datasets indicate a strong correlation between the Botswana High and ENSO in JFM; however, the correlation becomes rapidly weaker when negative lags are introduced, indicating a very weak relationship between the Botswana High and decaying phase of ENSO. On the other hand, introducing positive month lags, the correlation becomes progressively weak but remains strong (0.61-0.69) at 95% significance up to a lag of 2, which is related to the ENSO development phase. At lag 3, which is indicative of the mature phase of ENSO, the correlation becomes weaker heading to the decaying phase of ENSO in lag 4 and 5. Despite the strong correlation between the Botswana High and the tropical Indian Ocean SST (as shown Fig. 7), there seems to be no link between the high and the Indian Ocean dipole index (Table 1).
To further assess the relationship between the Botswana High and ENSO, Fig. 8 shows composite anomalies of JFM 200 hPa velocity potential and stream function (contours) during strong Botswana High Years (1983, 1998) and weak Botswana High Years (1989, 1994. The results show that +ve Botswana High years are characterised by anomalous upper tropospheric divergence (negative velocity potential) west of the prime meridian (i.e. 0° longitude) and anomalous upper tropospheric convergence (positive velocity potential) east of the prime meridian and vice-versa during −ve years. This upper-level convergencedivergence pattern during +ve (−ve) Botswana High years is reminiscent of the weakening (strengthening) of the Walker Circulation and typically forms during El Nino (La Nina) years. Furthermore, the increase (decrease) in upper-level convergence over southern Africa may lead to increased (decreased) subsidence over the region and the strengthening (weakening) of the Botswana High. The 200 hPa eddy stream function during +ve (−ve) Botswana High years are characterised with anticyclonic(cyclonic) anomalies in the upper troposphere over the central and eastern pacific. This pattern is characteristic of a Gill-Matsuno type response, which typically occurs during El Nino (La Nina) Years due Fig. 9 Composite of the standardized anomalies of JFM 500 hPa Omega during +ve phase of Botswana High Years (1983,1998,2010) and −ve phase of Botswana High Years (1989,1994,2008).
The mean and standard deviation used in calculating the anomalies were obtained using the 1980-2010 data. The green box shown in the region is used for the vertical-cross section plot in Fig. 10 to warm (cold) SST anomalies in the tropical Pacific Ocean (Gore et al. 2020). This response indicates a weakening (strengthening) of the upper-level cyclonic flow in the eastern Pacific and the weakening (strengthening) of the Walker Circulation. During +ve (−ve) years, the upper-level stream function shows anomalous cyclonic (anticyclonic) flow over the Indian Ocean, Asia and southern Africa, which corresponds to the increase (decrease) in convergence over the area leading to lower-level subsidence (convection). Figures 9 and 10 show the anomalous 500hpa geopotential height and vertical cross-section of the geopotential height during +ve and −ve Botswana High years. In general, the position of the Botswana High varies between the datasets during +ve years (Fig. 9a, c and e). The MPAS model shows a Botswana High located over northwestern South Africa and southeastern Botswana (Fig. 9e), while ERA5 shows the high over much of Namibia (Fig. 9a). However, 20C reanalysis fails to adequately capture the net sinking motion (ω > 0) associated with the Botswana High over Namibia, Botswana or South Africa (Fig. 9c). The reason for this is that the 20C Botswana High is weaker (Fig. 10c) as compared with ERA5 (Fig. 10a) and MPAS (Fig. 10e). This may also explain why the 20C EOF 1 (Fig. 4b) was weaker as compared to ERA5 and MPAS (Fig. 4a and c). Nevertheless, all the datasets show sinking motion associated with the Botswana High, which is characterised by rising air to the west of the High, which sinks and manifests as upper-level convergence in the vicinity of the high, leading to increased subsidence in the mid-levels. Figure 9b, d and f also show variability in the position of the Botswana High during −ve years, however, the position of the Botswana High is more similar as compared to +ve years. Despite that, MPAS agrees well with ERA5 reanalysis on the depth of the −ve Botswana High mode (Figs. 10b, d and f) and shows a −ve Botswana High that is characterised by lower divergence which might lead to increased convection.  (1983, 1998) and −ve Botswana High Years (1989, 1994. The mean and standard deviation used in calculating the anomalies were obtained using the 1980-2010 data. The green box shows the location of the continent while the red box shows the Botswana High

Positive phase of the Botswana High
Figures 11a, c, e and 12a, c, e show the anomalous 500 hPa geopotential height and anomalous surface temperatures, respectively, during +ve Botswana High years. The composite of standardised anomalies for 500 hPa geopotential height during +ve Botswana High years shows a good agreement between the MPAS simulation and reanalysis (20C, ERA5), in a sense that all datasets indicate enhanced geopotential height over the tropics and little to no change poleward of 30° S. Higher geopotential heights are typically associated with warm air masses and increased subsidence and may lead to an increase in temperature over the tropics. When an air mass descends, the pressure on the air mass increases. Because of the increase in pressure, the volume of the air mass decreases, increasing its internal energy, which manifests itself in the increase in temperature of that mass of air. This increase in the geopotential height agrees with Fig. 12, which shows that +ve Botswana High years are generally characterised by warmer temperatures over the land surface as well as over the tropical Atlantic and the tropical Indian Ocean. During +ve Botswana High years, the observational data shows above normal temperatures mostly over the eastern parts of the continent (Zimbabwe, Zambia, and Mozambique) with slightly warmer conditions over the northeast regions of South Africa (Fig. 12a, c, and e). The reanalysis shows a similar pattern; however, the warmer conditions extend throughout the eastern half of South Africa. On the other hand, the model simulation shows above-normal temperatures for much of the eastern parts of the continent as well as over South Africa. In addition to that, the positive anomaly extends along the coastline of Namibia and Angola as compared to the CRU CFSR and ERA5. The above-normal temperatures in the MPAS simulation may be attributed to the anomalous local maximum in the geopotential height over the South African Development Community (Fig. 11e), which is stronger in MPAS as compared to the reanalysis ( Fig. 11a and b).

Fig. 11
Same as Fig. 9 but for 500 hPa geopotential height MPAS and reanalysis datasets indicate that +ve Botswana High years are typically characterised by neutral to positive OLR anomalies (Fig. 13a, c, e), lesser convection (Fig. 9a, b and c) and lesser rainfall (Fig. 14a-c). Positive OLR anomalies indicate suppressed convection, lesser cloud coverage and more radiation emitted back into space, which is consistent with the increased subsidence in the area. However, there are notable discrepancies between 20C, ERA5, and the model simulation of OLR, especially over the ocean. The differences in OLR may be due to the Tiedtke convective parameterisation scheme used in the model, which removes convective instability before resolved-scale motions can respond to it. Another reason could be that the Tiedtke parameterisation is sensitive to entrained air from the free troposphere, and thus convection can be reduced by dry, free troposphere (Ali et al. 2015).
The discrepancies in the position of the Botswana High sinking motion may impact rainfall patterns in the datasets over southern Africa (Fig. 14). In actuality, there is poor agreement between the datasets on the +ve phase rainfall anomalies over land. GPCP shows below-average rainfall for much of the southern African continent and mildly wet conditions over South Africa, northern Angola, and Tanzania. The reanalysis shows an agreement with GPCP over southern Africa, Angola, and Tanzania but shows disagreement with the observation over Madagascar, Zambia and the Democratic Republic of Congo (DRC). The MPAS model also shows disagreements with GPCP and ERA5 over Botswana and northern Namibia, where it simulates above average rainfall, whereas the other datasets show below-average rainfall. The MPAS model's inability to adequately capture the precipitation anomalies may be due to the low resolution (240 km) used in the simulation. A higher resolution simulation could allow for better representation of orography and surface fields vital for the initiation of convection in complex terrains (Hohenegger et al. 2008).
MPAS shows a broad agreement with reanalysis on drought patterns during +ve Botswana High years, although with some discrepancies. Foremost, all datasets show a warm bias in PET over much of the subcontinent (Fig. 15a,  c and e); however, the model shows higher anomalies than reanalysis. Higher PET anomalies indicate enhanced  Fig. 9 but for surface temperature evaporation, leading to more severe drought conditions over the region in the absence of precipitation. Despite the good agreement among the dataset on the PET pattern, they disagree on the meteorological drought patterns (i.e. SPI and SPEI) during +ve Botswana High years (Figs. 16a,c,e and 17a,c,e). For example, both MPAS and 20C show normal to wet conditions over Namibia, while ERA5 shows parched conditions. On the other hand, both reanalysis datasets show a relatively dry bias over Botswana, Zimbabwe and the northeastern regions of South Africa, while the model shows a wet bias. Despite that, some areas show agreement between the datasets, such as Angola, parts of Tanzania, and central parts of South Africa. However, MPAS shows better agreement with 20C in SPEI anomalies than with ERA5. For instance, both MPAS and 20C show +ve SPEI anomalies over central Namibia, Angola and northeastern South Africa, while ERA5 shows severe drought over much of the region.

Negative phase of the Botswana High
MPAS and reanalysis data show that −ve Botswana High years are characterized by a decrease in the geopotential height ( Fig. 11b, d, and f) which is typically associated with cooler air masses and is consistent with the decrease in the surface temperatures over land (Fig. 12b, d, and f). During −ve years, the most obvious disagreement between the reanalysis and MPAS is the local minimum geopotential heights simulated over the Atlantic and the Indian Ocean as compared to the reanalysis. In addition to that, the reanalysis shows a slight increase in the geopotential height over the subtropics whilst the MPAS simulation indicates a decrease in the geopotential height (Fig. 11b, d, and f). In connection to the decrease in geopotential height over the region, both CRU CFSR and ERA5 show below normal temperatures over much of the continent, however, the model simulation shows not much change in temperature except over Zambia, the DRC, Tanzania and northern Mozambique. Despite this, the model shows a good comparison with observation and reanalysis regarding surface temperatures over the ocean. The datasets show two warm tongues over the Atlantic and a local minimum surface temperature south of Madagascar.
There is a better agreement in the OLR, omega, and rainfall patterns between the datasets during −ve Botswana High years as compared with +ve years. However, there are notable discrepancies between the model simulation and ERA5 reanalysis, especially over the ocean. In general, the datasets show negative OLR anomalies during −ve years are indicative of enhanced convection, which may lead to more cloud coverage and rainfall (Fig. 13b, d, and f). More cloud coverage may increase the albedo effect leading to colder temperatures over the study areas as shown by Figs. 12b, d, and f, however, the effect of cloud cover on temperature is not singular as the role of horizontal advection is also important. Additionally, more convection implies higher and colder cloud tops that emit less radiation back into space. 500 hPa Omega anomalies during −ve Botswana High years (Fig. 9b, d and f) indicate areas of convection in the mid-levels (~ 500 hPa) in Namibia, southern Angola, Zambia, Botswana and parts of the interior of South Africa. The anomalous convection in those areas is consistent with the increased rainfall over the region (Fig. 14b, d, and f). Both GPCP and ERA5 show wetter conditions for much of the southern continent except for southern Mozambique, Zimbabwe, and parts of Angola. Similarly, MPAS simulates wetter conditions over much of the southern African continent except over Namibia, however, these are weaker as compared to the observation and reanalysis. This result is consistent with the result of Fig. 13f, which shows weaker negative OLR anomalies for MPAS as compared to Fig. 13b and d).
The agreement between MPAS and reanalysis on drought patterns is better during the −ve Botswana High years than during +ve years. As expected, the datasets generally show a cold bias in PET anomalies over the region, which is indicative of cool (Fig. 12b, d, and f), cloudy (Fig. 13b, d, and f) and rainy conditions (Fig. 12b, d, and f). However, the MPAS model shows near-normal PET anomalies for South Africa, Botswana and Namibia, which disagrees with the reanalysis results. This discrepancy may be due to the near-normal temperatures simulated by the MPAS model in Fig. 12f, which the PET calculation is dependent on. The SPI also shows a good agreement between the datasets especially over South Africa, Botswana, Angola, Zambia, southern DRC, northern Mozambique and much of Tanzania (Fig. 16b, d, and f). However, the model fails to capture the wet (dry) spells near Namibia (Zimbabwe and southern Mozambique) as in the reanalysis. Nonetheless, SPEI also shows a better agreement between the MPAS and reanalysis data during −ve Botswana High Years. Again, the datasets indicate +ve anomalies over much of the region except over Angola, and the southern coastline of South Africa (Fig. 17b, d, and f).

Conclusions
This study has evaluated the capability of MPAS in simulating the characteristics of the summertime Botswana High. The global simulation covers a period of 31-years . Focusing on the southern African domain, EOF analysis was used to extract the Botswana High feature from the 500 hPa geopotential height in the reanalysis and simulation datasets. To evaluate the skill of the model in capturing the spatial and temporal variability of the southern African climate and the Botswana High, the simulated results were compared to the observed and reanalysis data. The results of the study can be summarised as follows: 1. In general, MPAS performs well in simulating all the climate variables over southern Africa. The simulated rainfall and temperature bias are about 1 mm/day (positive) and 2 °C (cold), respectively. However, the spatial variations in the 500 hPa geopotential height are larger in 20C reanalysis (about − 16 to − 18 m) and MPAS simulation (about − 9 to 10 m) as compared to ERA5 reanalysis. 2. All the datasets agree that the Botswana High accounts for about 80% of the variability in the 500 hPa geopotential height over the Botswana area. However, the spatial pattern of the Botswana High simulated by MPAS fea- tures a better agreement with the ERA5 reanalysis than with 20C reanalysis. 3. MPAS reproduced the interannual variability of the Botswana High as in the reanalysis. The variability, which shows a strong connection with ENSO, suggests a link between the Botswana High and ENSO. 4. MPAS captures the 500 hPa vertical velocity over southern Africa with variations in the position of the Botswana High sinking motions. Despite that, the vertical cross-section of omega anomalies shows that the +ve Botswana High mode is characterised by upper level convergence while the −ve mode is characterised up upper level divergence. 5. The MPAS model captures well the drought patterns over southern Africa as compared to reanalysis, however there is a better agreement during −ve Botswana High years than +ve years. In general, +ve Botswana High Years are characterised by increased evaporation while −ve years indicate increased convection.
These results suggest that MPAS can be applied for seasonal forecasting over southern Africa, especially for predicting the influence of the Botswana High variability on the annual variation of summer rainfall. However, more simulations and model improvements are still needed to establish the robustness of the results. For instance, due to the limited computational resources available for the study, we could only perform the MPAS simulation at about 240 km quasi-uniform horizontal grid over the globe. To resolve more regional-scale rainfall-producing systems that may be influenced by the Botswana High, future studies can take advantage of the variable grid resolution of MPAS to perform more simulations with higher resolution over southern Africa. Such studies could improve knowledge on the influence of uniform grid vs. variable-resolution grid in simulating the characteristics of the Botswana High. In addition, our results on the relationship of the Botswana High with ENSO and drought, which is based on correlation and composite analysis, cannot determine causality in the relationships. Therefore, there is a need for a series of sensitivity experiments with MPAS to demonstrate and quantify the cause and effect relationships. However, the present study has demonstrated the capability of MPAS in simulating the characteristics and influence of the Botswana High in southern Africa.