Climatology of Borneo Vortices in the HadGEM3-GC3.1 General Circulation Model

: Borneo vortices (BVs) are intense precipitating winter storms that develop over the equatorial South China Sea and strongly affect the weather and climate over the western MaritimeContinent because of their association with deep convectionandheavyrainfall.In thisstudy,theabilityof theHadleyCentreGlobalEnvironmentModel3–GlobalCoupled, version 3.1 (HadGEM3-GC3.1), global climate model to simulate the climatology of BVs at different horizontal resolutions isexaminedusinganobjectivefeature-trackingalgorithm. TheHadGEM3-GC3.1attheN512(25km)horizontalresolution simulates BVs with well-represented characteristics, including their frequency, spatial distribution, and lower-tropospheric structures when compared with BVs identiﬁed in a climate reanalysis, whereas the BVs in the N96 ( ; 135km) and N216 ( ; 65km)simulationsaremuchweakerandlessfrequent.Also,theN512simulationbettercapturesthecontributionof BVs to the winter precipitation in Borneoand the Malay Peninsulawhen compared with precipitation froma reanalysis data and from observations, whereas the N96 and N216 simulations underestimate this contribution because of the overly weak low-levelconvergence of thesimulatedBVs.The N512simulationalso exhibitsan improvedabilityto reproducethemodulation of BV activity by the occurrence of northeasterly cold surges and active phases of the Madden–Julian oscillation in the region,includingincreasedBVtrackdensities,intensities,andlifetimes.Asufﬁcientlyhighmodelresolutionisthusfoundto be important to realistically simulate the present-climate precipitation extremes associated with BVs and to study their possible changes in a warmer climate.


Introduction
During boreal winter, the western ''Maritime Continent,'' including the regions of Borneo, the Malay Peninsula, and Sumatra, experiences pronounced rainfall seasons due to the domination of the northeasterly winter monsoon (Johnson and Houze 1987;Chang et al. 2005a). During the rainfall season, extreme precipitation events associated with deep cumulus convection are frequently observed when the northeasterly winds in the southern South China Sea appear to intensify and lead to strong equatorward surges of cold air (Lim et al. 2017). Such cold surge (CS) outbreaks with strong and concentrated low-level winds are considered to be driven by the sharp pressure gradient induced by the southward propagation of the continental Siberian high in the midlatitudes (Ding 1990;Wu and Chan 1995). One of the main weather systems behind the connection between CSs and precipitation extremes in the western Maritime Continent is the near-equatorial cyclonic disturbances that develop offshore in the downstream region of the CSs near the north of Borneo (Cheang 1977;Chang et al. 1979Chang et al. , 2003Johnson and Houze 1987;Chen et al. 2015a,b;Isnoor et al. 2019). These disturbances, usually referred to as Borneo vortices (BVs) or cold surge vortices (Chen et al. 2002(Chen et al. , 2013(Chen et al. , 2015a(Chen et al. ,b, 2017, are features with a closed cyclonic circulation in the lower troposphere and associated with intense rainfall with deep cumulus convection in the vicinity of the cyclone center (Koseki et al. 2014). Such synoptic systems usually have a shallow vertical structure (Robertson et al. 2011;Trilaksono et al. 2012) and their motion is usually quasi stationary (Chang et al. 2005b).
Previous studies have explored the physical relationship between the formation of BVs and CSs in winter. For instance, some studies have considered BVs to develop as a result of the high vorticity background created by the horizontal cyclonic shear of the CS (Chang et al. 2005b;Ooi et al. 2011). Koseki et al. (2014) performed an analysis of the vorticity tendency of BVs and suggested that the development of BVs is mainly through the vortex stretching effect of the CS. BVs are also facilitated by the low-level convergence of relative vorticity flux and moisture associated with the CS winds blocked and channeled by the adjacent orography of the southern South China Sea (Tangang et al. 2008;Koseki et al. 2014). Additionally, some studies attribute BVs to the downstream Denotes content that is immediately available upon publication as open access. effects of CSs associated with the development of upstream synoptic disturbances in the midlatitudes (Lau and Lau 1984;Chen et al. 2002). Such an effect was further explored by Lim and Chang (1981) using the linearized beta-plane shallow water equations to simulate the dynamical responses to the lowlevel high pressure surge in the midlatitudes. The simulation initialized by the midlatitude pressure surge (which can be considered as the southward extended Siberian high in winter) revealed that the meridional dispersion of the zonal Rossby wave group from the midlatitude forcing can result in a BVlike cyclonic disturbance near the equator after around 6 days (Lim and Chang 1981). Although BVs are considered to be tightly linked to the occurrence of CSs, BVs independent of CSs have also been discussed in previous studies (Koseki et al. 2014;Paulus and Shanas 2017), which are considered to originate from other factors such as westward-propagating waves from the western Pacific (Chang et al. 1979).
BVs are important agents of high-impact weather in the western Maritime Continent as their presence can induce pronounced low-level convergence of moisture, leading to sharp increases in precipitable water and strong mesoscale convective rainfall adjacent to the BV center (Koseki et al. 2014). Although the study of Chen et al. (2013) indicated that BVs contribute only a minority (approximately one-third) of the total precipitation during winter in Borneo and the Malay Peninsula, cases studies have suggested a close relationship between BVs and the occurrence of extreme flood events in the Malay Peninsula (Juneng et al. 2007;Tangang et al. 2008). A typical case of BVs associated with extreme floods was also observed over western Borneo in January 2017, which left heavy precipitation of up to 154.1 mm in one day over the Ketapang region, Indonesia (Isnoor et al. 2019).
Because of the significant socioeconomic impact of BVs across the western Maritime Continent, it is important to improve our understanding of BV climatologies and how they may change under future climate scenarios using numerical models. However, accurate simulation of BVs is a challenge for current climate models because of the need for high model resolution to explicitly resolve the highly localized features of BVs, such as the sharp pressure gradient and convergence of horizontal winds associated with the extreme precipitation near the BV center (Koseki et al. 2014). Also, understanding of the fine-scale dynamics of BVs is still limited by the coarse model resolutions used in climate reanalysis datasets, which are commonly used to investigate the historical climatologies of BVs through objective identification approaches (Chang et al. 2005b;Braesicke et al. 2012;Koseki et al. 2014;Paulus and Shanas 2017). Moreover, to correctly simulate BV climatologies, models should be able to reproduce the interactions of BVs with various large-scale environmental factors, especially the complicated pattern of changes in precipitation related to the interactions among BVs, CSs and the large-scale convection associated with the active phase of the Madden-Julian oscillation (MJO) in the western Maritime Continent (Tangang et al. 2008;Paulus and Shanas 2017;Saragih et al. 2018). So far, research on the ability of climate models to simulate BV climatologies has received little attention, though such studies are important to provide confidence in using climate models to investigate the BV activity in the future under climate change. Understanding the bias in the simulation of BVs also helps improve the simulation of precipitation in Southeast Asia (Loh et al. 2016).
The main aim of this study is to examine how well low-, medium-, and high-resolution general circulation models represent the climatology of BVs, including their spatial distribution, intensity, lifetime, and the associated large-scale environments in the western Maritime Continent. To address this, we first diagnose the characteristics of BVs in the ERA5 reanalysis dataset and a set of atmosphere-only present-climate simulations performed using the Hadley Centre Global Environment Model 3-Global Coupled, version 3.1 (HadGEM3-GC3.1), global climate model (GCM) at horizontal resolutions of N96 (;135 km), N216 (;60 km), and N512 (;25 km) produced for the High Resolution Model Intercomparison Project (HighResMIP; Haarsma et al. 2016;Roberts et al. 2019) of CMIP6. The ability of HadGEM3-GC3.1 to simulate BV climatologies is evaluated with respect to the BVs identified from the ERA5 reanalysis data. The secondary aim is to explore the effect of the model resolution on the simulated BV-associated precipitation and their modulations by CSs and MJO through the intercomparison among the simulations with different horizontal resolutions and precipitation from observations and ERA5.
The article continues with section 2 describing the data used and the methods for the diagnoses of the BV characteristics and their associated large-scale environments. Section 3 presents the results of the comparisons of the BV climatologies between ERA5 and HadGEM3-GC3.1. Section 4 summarizes and discusses the results and describes future work.

a. Verification data
The ECMWF ERA5 reanalysis data (Hersbach and Dee 2016) for the main BV season of October-March (ONDJFM) for the period 1979-2014 (34 BV seasons in total) are used to verify the ability of the HadGEM3-GC3.1 simulations to represent the climatology of BVs and their environments. ERA5 is produced using 4D-Var data assimilation and a spectral model at a spectral horizontal resolution of TL639 (;31 km) with 137 hybrid sigmapressure vertical levels. In this reanalysis, the climatological average precipitation is in good agreement with that observed over Southeast Asia (ECMWF 2019), hence it is suitable to study BVs and the associated precipitation in the region.
The data used for verifying the ability of HadGEM3-GC3.1 to simulate the precipitation associated with BVs include the daily precipitation from ERA5 and two observed daily accumulated precipitation datasets, one based on in situ rain gauge and one based on satellite data. The rain gauge data used come from APHRODITE (Asian Precipitation-Highly-Resolved Observational Data Integration toward Evaluation; Yatagai et al. 2014), which is a gridded daily precipitation dataset over land at 25-km resolution based on rain gauge data developed by the Research Institute for Humanity and Nature and the Meteorological Research Institute of the Japan Meteorological Agency. The satellite-based data are from the Integrated Multi-Satellite Retrievals for Global Precipitation Measurement (IMERG-GPM, referred to herein as GPM) implemented by NASA and the Japan Aerospace and Exploration Agency (Huffman et al. 2015). GPM provides daily precipitation estimates at a spatiotemporal resolution of 0.18 3 0.18 since 2000 over the domain 608N-608S. Comparisons among APHRODITE, GPM, and ERA5 are presented to consider the uncertainties in precipitation products (Kotsuki and Tanaka 2013).

b. Model data
The ability of HadGEM3-GC3.1 (Williams et al. 2018;Menary et al. 2018;Kuhlbrodt et al. 2018) to simulate the climatology of BVs during ONDJFM for the period 1979-2014 is investigated. The HadGEM3-GC3.1 simulations used in this study are the present-climate atmosphere-land-only simulations forced by the observed sea surface temperature, sea ice, and historic forcings, which were performed for the CMIP6 HighResMIP project (Haarsma et al. 2016;Roberts et al. 2019). The atmospheric model of HadGEM3-GC3.1 applies nonhydrostatic Navier-Stokes equations (Davies et al. 2005). The advection scheme applies the Eulerian scheme for density and the semi-Lagrangian discretization (Davies et al. 2005) for prognostic variables (the three-dimensional wind components, potential temperature, Exner pressure, dry density, and the mixing ratios of water vapor). The semi-implicit method is used for the time integration (Davies et al. 2005). The horizontal grids of the simulations are set as regular latitude-longitude grids at the N96 (192 3 144 longitude/latitude, referred to as Had-LM), N216 (432 3 324; Had-MM), and N512 (1024 3 768; Had-HM) resolutions. The vertical resolution of the simulations in the atmosphere is 85 levels with a top at 85 km and four soil levels. All the simulations were run with time steps of 20 min. The model uses the Joint U.K. Land Environment Simulator (JULES) land surface model (Walters et al. 2017) at the atmospheric resolutions. It also includes a full set of physical parameterizations (radiation, convection, boundary layer mixing, and cloud microphysics) and aerosol interactions (Hewitt et al. 2011). The present-climate atmosphere-landonly simulations of the model cover the period 1950-2014. Here we only verify the simulations of BV-associated precipitation for the years 2000-14 given the limitation of the common time span of GPM and the simulations. For verifying the other simulated characteristics of BVs including their frequencies, spatial distributions, structures, and the associated large-scale environments, the common time span of ERA5 and the simulations  is chosen as the validation period.

c. Detection and tracking of Borneo vortices
In this study, BVs are identified and tracked in the ERA5 reanalysis data and the HadGEM3-GC3.1 simulations using the objective tracking algorithm (TRACK) described by Hodges (1994Hodges ( , 1995Hodges ( , 1999. This algorithm has previously been used in studies of high-impact synoptic systems including extratropical cyclones (e.g., Bengtsson et al. 2007;Zappa et al. 2013;Hawcroft et al. 2017;Catto 2018;Priestley et al. 2020), tropical cyclones (e.g., Bengtsson et al. 2007;Manganello et al. 2012Manganello et al. , 2014Manganello et al. , 2016Strachan et al. 2013;Bell et al. 2013Bell et al. , 2014Roberts et al. 2015Roberts et al. , 2020, tropical easterly waves (e.g., Thorncroft and Hodges 2001;Serra et al. 2010;Gomes et al. 2015), polar lows (Zappa et al. 2014;Yanase et al. 2016), and Tibetan Plateau vortices (Curio et al. 2018(Curio et al. , 2019. To identify and track BVs, the TRACK algorithm first spectrally filters the 6-hourly 850-hPa relative vorticity fields, with total wavenumbers less than 5 removed and the field truncated to T63 as well as applying the spectral tapering of the spectral coefficients using the method of Hoskins and Sardeshmukh (1987). This acts to reduce the high-frequency noise (such as spurious numerical oscillations) and removes the low-frequency background in the vorticity fields so that BV tracks can be coherently obtained.
Initially all tropical cyclonic systems are identified and tracked. To do this, all positive vorticity centers in the filtered 850-hPa RV fields are identified as grid point maxima that exceed 5 3 10 26 s 21 , then the off-grid vorticity maxima are obtained by performing a B-spline interpolation of the vorticity field and steepest ascent maximization using the grid point maxima as the starting locations (Hodges 1995). The tracking is then performed by first initializing a set of tracks from the vorticity maxima using a nearest neighbor method, these tracks are then refined by minimizing a cost function for track smoothness subject to suitable adaptive constraints (Hodges 1995(Hodges , 1999. The intensities of the tracked vortices [in terms of the maximum wind speed at 10 m and 925 hPa and the minimum mean sea level pressure (MSLP) near the tracked center] are determined following the tracks. The wind intensities are determined through a direct search of the full-resolution grid point values within a (geodesic) radius of 48 of the storm center. The minimum MSLP is determined through the B-spline interpolation and steepest descent to search for a minimum within 48 of the storm center. The BVs are identified among all of the tracked vortices by requiring that they have a pressure minimum for at least one day within a targeted area (2.58S-7.58N, 102.58-117.58E) near the northwest coast of Borneo; that is, only tracks that are quasi stationary within the Borneo region are considered, although they can come from outside the region and initially be more mobile. The genesis and lysis of BVs are defined respectively as the first and last time steps of the identified tracks wherever these occur.
Following the identification of BVs, the mean composite structures (including MSLP, 850-hPa horizontal winds, and divergence) during the mature phase (i.e., the time step corresponding to the lowest central MSLP during the lifetime) of the 50 most intense BVs within the region defined above are computed. Following Koseki et al. (2014), the composite structures of BVs are presented in a system-centered natural longitude-latitude coordinate. The 50 most intense BVs in terms of their central minimum MSLP are chosen to focus on the most extreme and impactful systems and avoid the large increase in the required computation and data storage (Manganello et al. 2014) for a larger sample.

d. BV-associated precipitation
The contribution of BVs to the precipitation (including the total precipitation and the extreme precipitation amount during ONDJFM) in the region and how well the model is able to represent this is investigated. The precipitation associated with BVs is accumulated on the model grids for all precipitation occurring within a 3.58 geodesic radius of the identified BV centers. Then, the fractional contributions of BVs to the total precipitation and precipitation extremes are calculated. The search radius of 3.58 is approximately 400 km, covering the intense precipitation near the BV center as well as the precipitation extending farther away during the developing and decaying stages (Koseki et al. 2014). To examine the simulations of BV-related extreme precipitation, on each model grid the occurrences of extreme precipitation for daily precipitation . the 95th percentile maximum daily precipitation for rainy days (daily precipitation . 0.1 mm) during ONDJFM is determined.

e. BV-associated environmental variables
The simulation of the large-scale environment associated with BVs is an important aspect in understanding the ability of HadGEM3-GC3.1 to simulate BV climatologies. The study of Lim and Chang (1981) revealed a close relationship between the formation of BVs and the midlatitude high pressure forcing (i.e., the southward extension of the Siberian high pressure in winter). Hence, the MSLP fields indicating the location of the Siberian high pressure during the genesis phase of BVs are investigated. Additionally, the 850-hPa wind field indicating the CS winds and the 850-hPa relative vorticity and divergence fields during the genesis phase of BVs are analyzed because of their close associations with the formation of BVs (Chang et al. 2005a;Tangang et al. 2008;Ooi et al. 2011;Koseki et al. 2014). These environmental fields are examined in the different HadGEM3-GC3.1 versions and compared with those of ERA5 to understand the simulated dynamical factors that affect BVs and their biases in the model.

f. Northeasterly cold surge
The interactions between BVs and the northeasterly CS in winter over the South China Sea can result in strong deep cumulus convection and heavy rains in the west of the Maritime Continent (Koseki et al. 2014;Paulus and Shanas 2017). Therefore, investigating the simulated BVs and the associated precipitation during CS events is helpful in understanding the ability of the model to simulate the contribution of BVs to winter precipitation. Here, the occurrences of BVs and BVassociated precipitation during CS events are analyzed. The CS events are diagnosed based on the method of Lim et al. (2017) and Xavier et al. (2020), which defines CS events as the days when the daily regionally averaged 850-hPa northeasterly wind speed within the region 58-108N, 1078-1158E exceeds 0.75 standard deviations above its long-term mean value.

g. Diagnosis of Madden-Julian oscillation
Near the western Maritime Continent, the eastward-propagating convective envelope of the MJO occurs during phases 3, 4, and 5, whereas suppressed convective activity occurs for the other phases. During winter and early spring, phases 3-5 of the MJO are found to facilitate the low-level convergence and vertical velocity over the western Maritime Continent (Wu and Hsu 2009;Oh et al. 2012) and cause positive anomalies of deep convection and precipitation over the region (Peatman et al. 2014). Enhanced BV activity and associated precipitation are also found during MJO phases 4 and 5 (Saragih et al. 2018). To investigate the ability of HadGEM3-GC3.1 in simulating the interactions between BVs and the active phases of the MJO, the simulated fraction of BV occurrences during MJO phases 3-5 with respect to the total BV occurrences during ONDJFM is computed and compared with that in ERA5. The impact of the MJO on the BV-associated precipitation is also analyzed by calculating the composite difference of the BV-associated daily precipitation rate during MJO phases 3-5 relative to the ONDJFM-mean values.
To determine the days corresponding to the active phases of the MJO in the HadGEM3-GC3.1 simulations, we compute the daily real-time multivariate MJO (RMM) index using the daily outgoing longwave radiation and zonal winds in the upper (200 hPa) and lower (850 hPa) troposphere using the method of Wheeler and Hendon (2004). The index is used to define the eight different MJO phases and amplitude for each day during ONDJFM. The MJO-related BV occurrences and BV precipitation will be analyzed during the days corresponding to MJO phases 3-5.

a. BV frequency
The monthly mean BV frequency for 1979-2014 from ERA5 and the HadGEM3-GC3.1 simulations is shown in Fig. 1. For ERA5, the mean BV frequency during the main BV season (ONDJFM) is 38.8 BVs. The peak in BV frequency (7.8 BVs per month) occurs in January, and the minimum of the monthly BV frequency is in October (with 3.6 BVs per month). For the HadGEM3-GC3.1 simulations, the simulated frequency is 8.5 BVs per season for Had-LM, 27.1 for Had-MM, and 35.4 for Had-HM. Therefore, the HadGEM3-GC3.1 simulations with relatively coarse resolution generally show lower seasonal frequency of BVs. Had-HM simulates reasonably well the mean BV frequency during ONDJFM, with a slight overestimation of the BV frequency (by 3.4 BVs per season). All the HadGEM3-GC3.1 simulations reproduce well the minimum BV frequency in October. However, Had-LM and Had-HM simulate an earlier peak of BV frequency in December, in contrast to the peak in January for ERA5. For Had-MM, the peak of BV frequency in January is not distinct.

b. Spatial distribution of BVs
Evaluation of the ability of HadGEM3-GC3.1 to simulate the spatial distribution of BVs, in terms of the genesis, track, and lysis densities, during ONDJFM is displayed in Fig. 2 . 2b) shows the majority of Borneo Island is affected by at least 9 BV tracks per season per 3.58 spherical area (ONDJFM). The local maximum of the BV track density (around 27-30 tracks per season per 3.58 spherical area) is seen in the coastal regions of northwest Borneo. The track density for the landfalling BVs affecting the southeast of the Malay Peninsula is about 6-9 tracks per season per 3.58 spherical area. Also, about 3-6 tracks per season per 3.58 affect western Sumatra. The HadGEM3-GC3.1 simulations generally reproduce well the location of the maximum track density in northwest Borneo. Relative to ERA5, Had-LM (Fig. 2e) has a pronounced underprediction of the local track density maximum (by more than two-thirds) resulting from the underestimate of the BV genesis. Had-MM (Fig. 2h) exhibits an improved simulation of the local maximum track density and shows well-simulated track density in the Malay Peninsula when compared with ERA5. Had-HM (Fig. 2k) presents the best performance in reproducing the local maximum track density relative to Had-LM and Had-MM, whereas the model overestimates the track density (by around 3 tracks per season per 3.58 spherical area) in the Malay Peninsula and the southern Philippines. Figure 2c shows the BV lysis densities in ERA5. The location where the maximum lysis density is approximately 6.58 to the southwest of the main genesis location (i.e., the local maximum genesis density; Fig. 2a), implying a southwestward propagation for most of BVs. In general, Had-MM and Had-HM simulate the local maximum of lysis density and its displacement from the main genesis area reasonably well, while Had-LM simulates only about one-fifth of the local maximum lysis density relative to ERA5. The lysis densities from Had-HM exhibits the closest magnitude with ERA5 relative to Had-LM and Had-MM, whereas Had-HM still underestimates the local maximum of lysis density by around 2 BVs per 3.58 spherical area. Had-HM also overestimates the lysis density by about 1 BV per 3.58 spherical area in the Malay Peninsula and western Sumatra. The distribution of the minimum central MSLP for all BVs from ERA5 shows peak percentages at 1006-1008 hPa (above 40%; Fig. 3a). For the HadGEM3-GC3.1 simulations, Had-LM and Had-MM generally simulate stronger BVs relative to ERA5 as the peak percentage shifts to lower pressures by about 2 hPa. Had-HM correctly reproduces the maximum percentage .40% at 1006-1008 hPa; however, it also overestimates the BV intensity as it doubles the percentages of BVs at 1002-1006 hPa and reproduces only one-third of the percentage at 1008-1010 hPa relative to ERA5. For the maximum MSLP gradient (Fig. 3b), Had-LM and Had-MM generally present weaker gradients than ERA5. Although Had-HM shows the best distribution when compared with other simulations, the gradient value corresponding to the highest percentage is approximately one-half of that in ERA5. simulations. However, Had-MM and Had-HM generally produce larger and smaller percentages for BVs intensities of greater and less than 16 m s 21 , respectively, as compared with ERA5 and Had-LM. Figure 4 shows the distribution of BV lifetimes calculated by the feature point number of the identified BV tracks. For ERA5, the maximum percentage is shown at the shortest lifetime and the percentage of BVs decreases exponentially with longer lifetimes; over 25% of the BVs have a lifetime shorter than 3 days, and around 20% of the BVs have a lifetime more than a week. For the HadGEM3-GC3.1 simulations, the models are able to reproduce the decreasing percentage with increased lifetime. However, Had-LM overestimates (underestimates) the percentages for lifetime longer (shorter) than 6 days, indicating a general overestimation of BV lifetimes. Had-MM shows improved simulations of BV lifetimes relative to Had-LM, with less overestimation of the percentages at .6 days and less underestimation at 6 days and less. The simulated BV lifetimes for Had-HM shows similar results to Had-MM.

d. BV structure
BVs are generally characterized by closed circulations with shallow structures from the surface to 700 hPa (Trilaksono et al. 2012). Here, we examine how well HadGEM3-GC3.1 is able to simulate the structures of BVs in the lower troposphere by computing the horizontal composite fields near the centers of the selected 50 most intense BVs (in terms of the 50 lowest records of the minimum central MSLP) from ERA5 and the HadGEM3-GC3.1 simulations.
Figures 5a-d show the composite structures of MSLP and its horizontal gradient from ERA5 and the simulations of HadGEM3-GC3.1. The central minimum of MSLP is broadly consistent across ERA5 and the HadGEM3-GC3.1 simulations. However, all of the HadGEM3-GC3.1 simulations produce lower central pressure values by around 0.5-1.0 hPa relative to ERA5, which is similarly shown in Fig. 3a. The horizontal gradient of MSLP of ERA5 shows a steep pressure gradient (.15 3 10 24 Pa m 21 ) within 28 of the BV centers. Had-LM and Had-MM poorly resolve the steep pressure gradient near the BVs centers, although they capture a relatively high pressure gradient (.5 3 10 24 Pa m 21 ) in the northern sector of BVs. In contrast, Had-HM shows the best performance in simulating the steep pressure gradient near the centers when compared with other simulations, which is consistent with result shown by the distribution of BVs with the maximum MSLP gradient (Fig. 3b).
The horizontal structure of wind fields at 850 hPa near the BV centers is displayed in Figs. 5e-h. The wind composites from ERA5 (Fig. 5e) show two local maxima of wind speed (.8 m s 21 ) at 28 and 78 from center, which are partly explained by the steep pressure gradient near the center and at the northwestern sector respectively (Fig. 5a). Relative to ERA5, Had-LM (Fig. 5f) and Had-MM (Fig. 5g) show a poorly resolved near-center maximum wind speed, although they are able to capture the relatively high winds at 48-78 from the center. Had-LM and Had-MM also present underestimated near-center winds in the eastern sector, leading to a semi-enclosed cyclonic flow around the BV centers. Had-HM, however, captures reasonably well the maximum wind speed within 28 and the eye with light winds.
The intense rainfall associated with BVs is found to be closely associated with the low-level convergence of the warm, moist southeasterly flow and the cool northeasterly surge (Tangang et al. 2008;Paulus and Shanas 2017). The low-level convergence also facilitates the maintenance of BVs as meso-a cyclonic systems (Koseki et al. 2014). The ERA5 reanalysis data (Fig. 5i) indicate that the most intense convergence at 850 hPa is located over the northern sectors and a divergence zone is seen in the southern sectors near the cyclone center. We further calculate the decomposed divergence of the zonal (meridional) wind by computing the partial derivative of zonal (meridional) wind speed with respect to zonal (meridional) distance. The divergence of the zonal and meridional winds (Fig. 5m) show a quadrupole structure within 28 from the BV centers and the easterly (southerly) wind components show relatively strong convergence in the northeastern (northwestern) sectors. Such structures are also found in the simulated BVs in the semi-idealized framework (Koseki et al. 2014), which is considered to be the result of the convergence tendency due to the deviatoric stress and self-maintenance via nonlinear process (i.e., the flow deformation tensor between the northeasterly surge and the southwesterly cyclonic wind in the northern sectors of BVs). When compared with ERA5, all of the model simulations (Figs. 5j-l) are able to capture the convergence in the northern sectors and the quadrupole structures of the convergence by the zonal and meridional winds (Figs. 5n-p). However, Had-LM and Had-MM simulate overly weak convergence (,4 3 10 26 s 21 ; Figs. 5j,k) and the maximum centers of the convergence by the different wind components occur farther away from the BV centers (Figs. 5n,o) relative to ERA5. In contrast, Had-HM demonstrates the best performance in representing the convergence zone in the northern sectors (Fig. 5l) and the maximum convergence by the zonal and meridional winds within 28 from the BV centers (Fig. 5p).
For the distribution of precipitation near the center, ERA5 (Fig. 5q) shows that the local maximum precipitation is located in the northern sectors, which is consistent with the strong convergence. Relative to ERA5, Had-LM (Fig. 5r) poorly resolves the near-center maximum precipitation. Had-MM (Fig. 5s) presents an improved representation of the location of the maximum precipitation. Had-HM (Fig. 5t) simulates the strongest precipitation relative to other simulations, which is partly due to the strong low- 3408 level convergence. Also, the simulated precipitation is more concentrated with underestimated precipitation (by about 24 mm day 21 ) across 28-48 west and 48-68 north from the center relative to ERA5.

e. BV-associated large-scale environments
Investigating the simulations of the large-scale environmental fields associated with BVs is an important aspect of understanding the model's ability to reproduce BV climatologies. Here we evaluate the HadGEM3-GC3.1 simulations of the anomaly fields of the large-scale environments associated with BVs in comparison with ERA5. To do so, the anomaly fields (including MSLP, wind fields, and their divergence and relative vorticity at 850 hPa) during the genesis stage of the 200 selected most intense BVs (in terms of the lowest records of the minimum central MSLP) relative to the seasonal mean fields during December-February (DJF) are computed to display the winter circulations associated with well-developed BVs, including the Siberian high and the northeasterly winter monsoon.
During 24 h before the formation of the BVs, ERA5 shows an increase in northeasterly winds of the CS at 850 hPa in the northern South China Sea, with wind speed anomalies up to 1.8 m s 21 to the east of Vietnam (Fig. 6a).
The relatively intense CS winds and their southward deflection due to the channeling of orography (especially to the east of Indochina Peninsula and the Malay Peninsula; Tangang et al. 2008;Chang et al. 2016) lead to increases in low-level convergence (up to 1.2 3 10 26 s 21 ; Fig. 6e) and relatively high vorticity (up to 4 3 10 26 s 21 ; Fig. 6i) to the south of 108N, which can provide a favorable environment for the formation and maintenance of BVs (Chang et al. 2016;Paulus and Shanas 2017). Figure 6b shows that Had-LM correctly simulates the location with intensified northeasterly in the South China Sea. However, the model simulates a strong southwesterly anomaly (up to 2.4 m s 21 ) near the south of the Philippines, which induces an overestimated cyclonic flows and associated convergence (Fig. 6f) and relative vorticity (Fig. 6j) in the area, which partly explains the longer mean BV lifetime relative to ERA5 as indicated by the mean track of BVs and the distribution of BVs with lifetime (Fig. 4). g,k) presents less overestimated southwesterly anomalies and associated convergence and vorticity relative to Had-LM. However, the maximum convergence and vorticity shift to the northeast relative to ERA5, which partially induce a northeastward shift of the mean genesis location of BVs. Similar biases are presented by Had-HM (Figs. 6d,h,l). Also, Had-HM simulates the largest increase in low-level convergence over the southern South China Sea relative to the other simulations, which is consistent with its highest BV frequency relative to other simulations ( Figs. 1 and 2).
Although BVs are closely associated with the interactions between CS and the nearby terrain, the southeastward extension of the Siberian high is also found to trigger equatorial cyclonic flows and facilitate the genesis of BVs (Lim and Chang 1981;Ooi et al. 2011;Chang et al. 2016). This process is associated with the geostrophic adjustment of Rossby wave group responses forced by the strengthened and southeastward extended Siberian high according to the barotropic equatorial beta-plane CS theory (Lim and Chang 1981), which is independent of the orographic effects (Lim and Chang 1981;Chang et al. 2016). Figure 7 show the distributions of DJF-mean MSLP indicating the mean location of the Siberian high and their anomalies during 24 h prior to the genesis phase of BVs. The MSLP contours inside the Tibetan Plateau are masked out as the Siberian high and associated cold air outbreak are usually active outside the plateau (Ding 1990). According to ERA5 (Fig. 7a), the center of the Siberian high indicated by the DJF-mean MSLP is located around 908E and 508N. To the south of the Siberian high center, positive anomalies (.2.5 hPa) of MSLP are found over the north and east of the Tibetan Plateau before the genesis of BVs, implying a strengthening and southeastward extension of the Siberian high. These changes lead to an increase in the MSLP gradient in the northern South China Sea (Fig. 7e), which drives the increased northeasterly as shown in Fig. 6a. All the HadGEM-GC3.1 simulations reproduce the location of the Siberian high reasonably well (Figs. 7b-d). However, Had-LM (Fig. 7b) and Had-MM (Fig. 7c) simulate a generally weakened Siberian high with negative MSLP anomalies near the high pressure center, although the positive MSLP anomalies to the east of the Tibetan Plateau are correctly simulated. These biases are in agreement with the inhibited BV activity in Had-LM and Had-MM relative to ERA5 ( Figs. 1 and 2). In contrast, Had-HM (Fig. 7d) displays the best performance in simulating the positive MSLP anomalies near the high pressure center, although overestimation is seen to the east of the Tibetan Plateau. Besides, in contrast to ERA5, all the simulations exhibit a negative MSLP anomaly in the northeast of the South China Sea (Figs. 7f-h), which partly explains the overestimated cyclonic flows and the northeastward shift of the mean BV genesis locations in the models (Figs. 6i-l). These biases may also lead to a longer distance of BV propagation to the equator, which prolong the BV lifetimes in the models relative to ERA5 (Fig. 4).

f. BV-associated precipitation
To evaluate the ability of HadGEM3-GC3.1 to simulate the BV-associated precipitation, the distributions of mean daily precipitation during the days when BVs occur are presented using ERA5 and all the HadGEM3-GC3.1 simulations (Fig. 8). APHRODITE (Fig. 8a), GPM (Fig. 8b), and ERA5 (Fig. 8c) show that the regions with the highest precipitation rates (.9 mm day 21 in ONDJFM) associated with BVs occur near the northwest coast of the Malay Peninsula and the Iran Mountains region of Indonesia across the northwest of Borneo. In these regions, GPM and ERA5 show a higher precipitation rate when compared with APHRODITE by about 2-4 mm day 21 , which is possibly related to the commonly found underestimates of precipitation in APHRODITE when compared with other precipitation products (e.g., Tong et al. 2014;Tan et al. 2015Tan et al. , 2017Immerzeel et al. 2015;Luo et al. 2020;Ji et al. 2020). In addition, the precipitation rate to the north and south of Borneo in GPM is higher (by around 3-5 mm day 21 ) than that in ERA5. Figure 8d shows that Had-LM reasonably simulates the local maximum precipitation associated with BVs on the east coast of the Malay Peninsula and the northwest of Borneo, while it underestimates the precipitation rate by around 2-7 mm day 21 in these regions relative to ERA5 and GPM. Had-MM (Fig. 8e) simulates well the BV-associated precipitation on the east coast of the Malay Peninsula when compared with APHRODITE, GPM, and ERA5 and improves the simulated precipitation associated with BVs in the northeast of Borneo relative to Had-LM. Had-HM (Fig. 8f) also shows reasonably simulated locations of the maximum precipitation associated with BVs when compared with ERA5. However, relative to APHRODITE, GPM, and ERA5, some underestimates of precipitation are still seen on the east coast of the Malay Peninsula and overestimation is shown in the south of Borneo. As well, for all the simulations, the local maxima of precipitation rate in Borneo and western Sumatra are found to shift inland relative to APHRODITE and GPM. Figure 9 presents the fractional contribution of BVs to the total precipitation amount and the accumulation of extreme precipitation when the daily precipitation is greater than the 95th percentile of maximum daily precipitation. For APHRODITE, GPM, and ERA5, BVs are associated with up to 50% of the total precipitation during ONDJFM across the western coast of Borneo (Figs. 9a-c). Also, up to 15%-25% of the total precipitation in the southeast of the Malay Peninsula is related to BVs. The BV-associated extreme precipitation, in terms of the 95th percentile maximum daily precipitation, accounts for up to 60%-70% of the seasonally accumulated extreme precipitation for ONDJFM on the western coast of Borneo and up to 20%-25% in the southeast of the Malay Peninsula (Figs. 10a-c). Therefore, the fractional contributions of BVs to the accumulation of extreme precipitation are generally greater than the contribution to the total precipitation. Relative to APHRODITE, GPM, and ERA5, Had-LM underestimates the fractional contributions of BVs to the total precipitation by around 25% in the midwest of Borneo (Fig. 9d). Underestimates are also found for the accumulations of extreme precipitation for daily precipitation greater than the 95th percentile (by about 30%; Fig. 10d) in the midwest of Borneo. Had-MM presents well-simulated contributions of BVs to the total precipitation (Fig. 9e), while the contributions to the extreme precipitation are still underestimated by around 10% for the daily precipitation greater than the 95th percentile maximum daily precipitation (Fig. 10e) in the midwest of Borneo. Had-HM generally simulates well the contribution of BVs to total precipitation (Fig. 9f), although it slightly overestimates the contribution of BVs to total precipitation (by 5%-10%) in western Borneo. In the northwest of Borneo, Had-HM also demonstrates some abilities in simulating BVs' contributions to the precipitation extremes (Fig. 10f) relative to APHRODITE, GPM, and ERA5, although underestimates still exist.

g. Impacts of CS and MJO on BVs
The daily mean track densities of BVs in ONDJFM during the presence of CS, MJO phases 3-5, and the combination of both CS and MJO phases 3-5 are shown in Fig. 11. The daily mean track densities are calculated as the track density of BVs during ONDJFM divided by the number of days. All the HadGEM3-GC3.1 simulations reasonably capture the averaged number of days during ONDJFM for CS events, MJO phases 3-5, and both CS and MJO phases 3-5, with biases of ,5 days when compared with ERA5. For BVs associated with CS events, ERA5 shows a track density of up to 32% per 3.58 spherical area per day over the west of the Maritime Continent (Fig. 11a), which is about 1.5 times the daily mean track density during ONDJFM. Such a facilitated occurrence of BVs by the CS is consistent with the result in Fig. 6a. For MJO phases 3-5 (Fig. 11b), an increase in track densities (by about 20%) is also seen but with smaller magnitudes relative to that during CS events. In contrast, the occurrences of CS combined with MJO phases 3-5 approximately double the daily mean track density during ONDJFM (Fig. 11c). The results are consistent with the finding of Lim et al. (2017) suggesting that the increased shear vorticity off Borneo during CS events and the facilitated conditional instability by the MJO can facilitate the development of BVs. Relative to ERA5, all the simulations (Figs. 11d-l) produce reasonable locations of the maximum daily mean track density and the increased densities during CS events. For Had-LM (Figs. 11d-f), the daily mean track densities and their relative increases during CS events are underestimated by about 80%. Also, Had-LM shows no apparent change in the mean track density when MJO phases 3-5 occur. Had-MM (Figs. 11g-i) exhibits an improvement in the simulation of the changes in track densities during CS and MJO phases 3-5, although underestimated changes (by about 20%) are seen during CS events (Fig. 11g) relative to ERA5. In comparison with other simulations, Had-HM FIG. 10. As in Fig. 9, but for daily precipitation greater than the 95th percentile.
(Figs. 11j-l) exhibits the best performance to simulate the changes in track densities associated with CS and MJO.
The effects of CS events and MJO phases 3-5 on the intensity and lifetime distributions of BVs are displayed in Fig. 12. During CS events, ERA5 indicates decreases in the fractions of BVs with relatively low intensities (i.e., the minimum central MSLP greater than 1006 hPa and maximum wind speed lower than 10 m s 21 at 10 m and 14 m s 21 at 925 hPa; Figs. 12a-c) and increased fractions at the higher intensity ranges for the CS-related BVs relative to all the BVs in ONDJFM. This implies an increase in BV intensities on average during CS events relative to all the BVs in ONDJFM. Figure 12d shows increased (decreased) fractions of BVs with lifetime greater (less) than 5 days, indicating a prolonged BV lifetime on average. Similar to the changes during CS events, phases 3-5 of MJO also indicate increased fractions of BV with relatively strong intensities (Figs. 12e-g) and prolonged lifetimes (Fig. 12h) as shown by ERA5. When both the CS events and phases 3-5 of MJO occur, BVs from ERA5 also tend to become more intense (Figs. 12i-k) and last longer (Fig. 12l). The distribution changes for the maximum wind speed at 10 m and 925-hPa (Figs. 12j,k) and lifetimes (Fig. 12l)  categorized by the central MSLP during CS events (Fig. 12a), while for Had-HM the simulated increase in fractions for the lower MSLP compares best to ERA5. In comparing Had-HM with Had-MM and Had-LM, it is seen that the changes in the fractions categorized by the maximum wind speeds at 10 m (Fig. 12b) and 925 hPa (Fig. 12c) shift toward stronger winds, which can be partly explained by the increased BV intensities with increased horizontal model resolution (Figs. 3b,c). For the changes during MJO phases 3-5, the magnitudes of change for the BV fractions categorized by the central MSLP are generally underpredicted by all the simulations (Fig. 12e). For the changes in BV lifetimes shown in Figs. 12d and 12l, Had-LM and Had-HM shift toward longer lifetimes (by around 1-2 days) relative to ERA5, possibly due to the overestimated fractions of BVs with relatively long lifetimes (.7 days; Fig. 4) in these simulations.
Previous studies have revealed a considerable increase in precipitation over Southeast Asia due to the effects of CS events and MJO (Lim et al. 2017;Fauzi and Hidayat 2018;Xavier et al. 2020). Here we discuss the modulations of BVassociated precipitation by CS events and MJO phases in ONDJFM for the period 2000-14 and consider whether the HadGEM3-GC3.1 simulations are able to represent such modulations. APHRODITE is not used for such analyses due to the absence of data over the ocean.
We first analyze the effects of CS events and MJO phases 3-5 on BV-associated precipitation rates in GPM (Figs. 13a,f,k). Increased BV-associated precipitation (by around 5.6 mm day 21 or greater) relative to the average over all BVs is shown in most of the southern South China Sea during CS events (Fig. 13a), especially over the east coast of the Malay Peninsula (.12 mm day 21 ). As shown in Fig. 13f, the presence of MJO phases 3-5 also leads to increased BV-associated precipitation to the north of Borneo (.4.2 mm day 21 ), while it is statistically insignificant and generally weaker than the increase by CS events. The presence of both CS events and MJO phases 3-5 further facilitates the BV-associated precipitation with statistically significant increases (greater than 7 mm day 21 ) in most of the southern South China Sea (Fig. 13k). Significant decreases in the BV-associated precipitation are also seen around the increased precipitation. This pattern of changes can be partly explained by the strengthening of BVs (Fig. 12), which can cause further concentrated water vapor transport and heavy rainfall in the south of the South China Sea (Lim et al. 2017). When compared with GPM, ERA5 presents a similar pattern of changes in BV-associated precipitation related to CS events and MJO phases [3][4][5]g,i), while the magnitude of precipitation increase is generally smaller. In Had-LM, the increased BV-associated precipitation due to the effects of CS events and MJO phases 3-5 are poorly simulated relative to both GPM and ERA5, with very limited areas of increased precipitation surrounded by negative anomalies (Figs. 13c,h,m). Had-MM presents an improved representation of the precipitation increase in the southern South China Sea during CS events (Fig. 13d), while the simulated precipitation increases are statistically insignificant and their magnitudes are underpredicted in northern Borneo by around 4.2 mm day 21 during CS events (Fig. 13d) and 2.8 mm day 21 during MJO phases 3-5 (Fig. 13d) relative to ERA5 (Figs. 13b,g). The limited abilities of Had-LM and Had-MM to simulate the increased BV-associated precipitation may be related to the poorly resolved low-level convergence of horizontal winds (Figs. 5j,k,n,o), which inhibits the concentration of water vapor transport and weakens the precipitation increase associated with CS and MJO (Lim et al. 2017). j,o) also shows apparent underestimation of the precipitation increases in the southern South China Sea. Such biases may be associated with the narrower near-center precipitation of the simulated BVs (Fig. 5t).

Conclusions and discussion
In this paper, the climatologies of BVs in the ERA5 reanalysis datasets and the HadGEM3-GC3.1 general circulation model at three different horizontal resolutions are examined using an objective feature-tracking algorithm (Hodges 1994(Hodges , 1995(Hodges , 1999  BVs. However, other large-scale dynamical features, such as the 850-hPa divergence and relative vorticity, do not explain the improved simulations of BVs. Hence, the main reason for the improvement in the simulation of BVs is possibly related to the finer grid space allowing the fine-scale winds to be further strengthened by the interaction between the constraint of fluid continuity and the emergent scaling properties of winds (Rauscher et al. 2016).
d According to APHRODITE, GPM, and ERA5 precipitation, the occurrence of BVs during ONDJFM can lead to relatively strong precipitation (.13 mm day 21 ) on the east coast of the Malay Peninsula and Borneo. Precipitation associated with BVs accounts for up to 50%-55% of the total precipitation in ONDJFM over the western coast of Borneo and up to 20%-25% in the southeast of Peninsular Malaysia. BVs contribute up to 60%-70% of the seasonally accumulated extreme precipitation for daily precipitation greater than the 95th percentile in Borneo. The HadGEM3-GC3.1 simulation with N512 resolution captures well the precipitation associated with BVs and their contributions to the seasonal total precipitation. The simulations with N96 and N216 underestimate the BV-associated precipitation and their fractional contributions to the seasonal precipitation and extremes, which are possibly due to the overly weak low-level convergence of BVs at these resolutions.
d When compared with both GPM and ERA5, the HadGEM3-GC3.1 simulations show similar modulations of BVs by CSs and active MJO phases (i.e., MJO phases 3-5) near the western Maritime Continent, which include the increase in the daily mean track density of BVs and increased BV intensities and lifetimes. However, over the southern South China Sea, all of the HadGEM3-GC3.1 simulations underestimate the increased BV-associated precipitation due to CSs and active MJO phases relative to GPM and ERA5. On the basis of the HadGEM3-GC3.1 simulations, this study presents improvements of the BV simulations with higher model resolutions. Given the close linkage between BVs and precipitation extremes, the high-resolution simulations are more beneficial for model-based research on extreme precipitation processes in the western Maritime Continent. Although in this study the HadGEM3-GC3.1 simulation at the coarse N96 resolution presents the simulated frequency and dynamical structures of BVs poorly, it still exhibits some ability to capture the large-scale environmental features associated with BVs, including the mean location of the Siberian high and the distribution of the intensified low-level CS winds and associated vorticity fields. Therefore, the N96 coarse-resolution simulations may be useful for studies of the large-scale environments associated with BVs and performing downscaling simulations of BVs based on regional climate models.
The limitations of the study include the lack of analysis of the multimodel uncertainty in simulating BVs due to the use of a single-GCM experiment. This will be addressed in future work using the HighResMIP multimodel ensemble (Roberts et al. 2020). Also, the impacts of different model configurations, such as the vertical resolution, physical parameterizations (especially for the cumulus convection parameterizations that strongly influence the simulation of precipitation over the Maritime Continent; Gianotti et al. 2012), and dynamical core, on the simulation of BVs remain unknown, and require further study.
This study evaluates the BV simulations of HadGEM3-GC3.1 based on comparison with only one reanalysis dataset, which can be affected by the uncertainties in the data assimilation, model formulation, and the quality of the underlying observations used in the assimilation process (Bengtsson et al. 2004). To improve this, a multi-reanalysis climatology of BVs will be investigated in future works to gauge reanalysis uncertainty. Also, datasets based on high-resolution operational analyses can be useful for the evaluation of storm simulations (e.g., Zappa et al. 2014) due to their generally much higher resolutions, which will be studied in the future.
In comparison with the simulations with relatively coarse resolutions, the HadGEM3-GC3.1 simulation with N512 resolution shows improved simulations of the southward extensions of the Siberian high prior to the formation of BVs, which possibly explain the improved BV simulations. However, due to the limited data availability (especially for the 6-hourly pressure-level variables) of the released CMIP6 HighResMIP data, the simulated downstream effects leading to the formation of BVs, such as the Rossby wave packets from the perturbation of the Siberian high (Lim and Chang 1981), are not analyzed here and will be explored in the future model-based studies.
In this study, HadGEM3-GC3.1 exhibits some abilities to simulate the response of BVs to CS and MJO. This includes the increased BV track densities and prolonged BV lifetimes. The facilitated conditional instability of MJO and the low-level convergence on the windward side of the terrain (Lim et al. 2017) can create a favorable environment for the maintenance of BVs, which may explain the prolonged BV lifetime by CSs and active MJO phases. However, the environmental factors related to the lifetime of BVs are not sufficiently understood and are worth studying in the future. In addition, the interannual variabilities of BVs associated with different atmospheric teleconnections, such as El Niño-Southern Oscillation and the equatorial quasi-biennial oscillation (Anip 2012), and how well current GCMs represent such variability is still not well understood. Future work will include analyses of the interannual variability of BVs and the associated physical mechanisms and simulations using GCMs. Besides, the diurnal cycle of precipitation is one of the prominent mechanisms of precipitation in the Maritime Continent (e.g., Oh et al. 2012;Peatman et al. 2014); however, its impact on the BV-associated precipitation remains unclear and needs further research. Also, the linkage between the precipitating storms in Southeast Asia and the equatorial waves including Kelvin waves, equatorial Rossby waves, and westward-moving mixed Rossby-gravity waves (Ferrett et al. 2020) will also be of interest in the future modelbased research.
The study indicates some credibility in the use of HadGEM3-GC3.1 for projecting the possible future changes in BV activity under climate change. Future work will use the HadGEM3-GC3.1 simulations from the HighResMIP project to study the projection of BV climatologies under the scenario of the representative concentration pathway RCP8.5 (Moss et al. 2010). The model evaluation in this article will provide helpful information for interpreting the future projection of BV activity using HadGEM3-GC3.1.