Numerical Weather Prediction at 200 m Local Resolution Based on an Unstructured Grid Global Model

It is a great challenge to reach sub-km resolution in numerical weather prediction models without using domain nesting techniques. The Model for Prediction Across Scales – Atmosphere (MPAS-A) uses an unstructured grid to achieve smooth transition of resolution in a global model. However, a number of model physical parameterization schemes and methods to enhance computational efficiency assume resolution down to about 3 km only. We extended the MPAS-A model to support resolution down to 200 m and carried out idealized and hindcast numerical experiments over a domain covering Hong Kong. The Shin-Hong Planetary Boundary Layer scheme was integrated to parameterize partially resolved turbulence in a

Resolution is one of the important factors to enhance the model forecasting skill but applying finer resolution incurs substantially higher computational costs. The higher the resolution desired, the smaller the area such resolution can be applied under practical computational cost constraints, and the transition to a coarser resolution for less interesting parts of the globe needs to be strategically designed. On the other hand, scale-awareness of physics parameterization schemes is also a concern in applying variable-resolution models. MPAS-A provides a 60-to-3 km variable-resolution mesh (largest resolution variation ratio and finest resolution among other downloadable meshes) and a scale-aware cumulus parameterization schemes for convection permitting resolution usage. For planetary boundary layer (PBL) parameterization, the YSU and MYNN schemes provided do not embed scale-awareness formulation. At around 1 km grid spacing, turbulence is only partially resolved; the use of non-scale-aware turbulence schemes thus becomes questionable, if not invalid.
The grid spacing used in real-world applications can be down to around 1 km -200 m, so-called "gray zone" resolution. We intend to extend the resolution supported by the MPAS-A framework down to 200 m. Shin and Hong (Shin & Hong, 2015) modified YSU PBL scheme to support gray zone region. They handled nonlocal transport via strong updrafts and local transport via the remaining small-scale eddies separately. Moreover, an entirely new formulation for nonlocal heat transport term that fits the LES output is constructed and this nonlocal term multiplies a grid-size dependency function. Additionally, the local transport which is still in eddy-diffusivity form multiplies a different grid-size dependency function. They used a Large Eddy Simulation (LES) with 25 m grid spacing as a benchmark. Based on this LES result, a list of reference data with different grid sizes is derived from 25 m to 4 km. For each reference data, the resolved and subgrid parts are separated and the ratio of subgrids over the total was calculated. This Shin-Hong model was implemented as a boundary layer parameterization module in the Weather Forecast and Research (WRF) model. Previous studies of the Shin-Hong PBL scheme using WRF or WRF-LES simulations models showed that the scale-aware Shin-Hong scheme can be used to reproduce consecutive rolls in a desert region (Xu et al., 2018) and an intensification scenario of a tropical cyclone .
In this paper, we elaborate the integration of the Shin-Hong scheme to our modified version of MPAS-A, as well as other practical issues for supporting 200 m resolution. We also verify such global model with idealized tests and real case tests with Hong Kong chosen as the city to apply 200 m resolution, analyze experiment results and highlight interesting findings. To the best of our knowledge, this is the first time 200 m resolution can be used in a global multi-resolution NWP model with unstructured mesh.
The paper is organized as follows. In Section 2, we describe our modified MPAS-A model and the components added or integrated to support 200 m resolution simulation. In Section 3, we present the validation of the dynamical core on the global mesh with the extreme resolution variation to reach 200 m resolution. In Section 4, results from simulations using different resolutions will be compared. In Sections 5, further discussion and conclusion are presented.

Model and Static Input Data
The MPAS-A was developed by the US National Center for Atmospheric Research (Skamarock et al., 2012). Its dynamical core solves a set of fully compressible non-hydrostatic equations using the finite-volume method, discretized horizontally on a C-grid staggered unstructured spherical centroidal Voronoi tessellation (SCVT) (Ringler et al., 2010;Thuburn et al., 2009) and vertically in a geometric-height hybrid terrain-following coordinate (Klemp, 2011). MPAS-A has been used for applications such as medium-range convection-permitting ensemble forecasts (Schwartz, 2019), study of spatial and temporal wind speed variability during open cellular convection (Imberger et al., 2021), diurnal variation of Mei-yu rainfall over East China (Xu et al., 2021), investigation of western North Pacific tropical cyclone tracks and intensities (Lui et al., 2021), etc. Modifying the source code of MPAS-A, we extend the model to cover usages where one would conventionally use WRF with domain nesting for high-resolution simulations around a region of interest.
Inspired by the use of different time steps in nested domains in WRF and remarked as a feature worth exploring for MPAS-A (Skamarock et al., 2012), the MPAS-A dynamical core was modified to use a hierarchy of time-steps which uses the smallest timestep, Δt, for the region with smallest grid spacing. As the grid spacing undergoes transition to coarser and coarser resolution outside the region of interest, the modified dynamical core uses timesteps of 2Δt, ⋯ , 2 i Δt, ⋯ , 2 (n−1) Δt, where i is called the HTS level and n is the number of HTS 10.1029/2022EA002342 3 of 34 levels. The mesh is spatially partitioned into zones with the minimal grid spacing occurring in each zone matching a 2 i Δt timestep that satisfies the Courant-Friedrichs-Lewy condition for numerical stability. The determination of zones and their boundaries is performed with a minimization of computational cost in parallel computing. This hierarchical time-stepping (HTS) technique essentially avoids cells with large grid spacing from using small timesteps unnecessarily and enables MPAS-A to be used on meshes of very large resolution variation ranges with affordable computational resources Ng et al., 2019). For the 200 m-resolution mesh (details in Section 2.1), there are 9 HTS levels, with timestep 0.781 s for the finest resolution zone (with grid spacing range 0.171-0.513 km), timestep 1.563 s for the second highest HTS level zone (with grid spacing range 0.349-1.16 km), so on and so forth to 200.0 s for the coarsest resolution zone (with grid spacing 42.9-59.7 km). Note that the grid spacing range of each HTS level may actually overlap.
Moreover, a technique based on the OLAM Walko & Avissar, 2011) for Customizable Unstructured Mesh Generation (CUMG) enables local mesh refinement in arbitrary shape using user-defined target horizontal resolution in any desired locations. It generates SCVTs that satisfy the original formulation of MPAS-A dynamic core. The generation of the meshes used in this paper is described in Section 2.1.
Both the techniques of HTS and CUMG relax the computational cost, making local high resolution in a global NWP model feasible. The new model is referred to as the Clustertech Platform for Atmospheric Simulations (CPAS) which is implemented on a HPC platform (https://cpas.earth/). The use of HTS and CUMG in CPAS was separately tested to produce very similar results to MPAS-A, which uses uniform timesteps, in a number of meteorological episodes (Ng et al., 2019;Tse et al., 2020). CPAS version 0.17.0, which is based on MPAS-A v7.1, is used in all the experiments in this paper.

Variable-Resolution Meshes
We generated several SCVT meshes with different resolution variation settings. The mesh for idealized test has a wide range of resolution variation 100 km-200 m while two CUMG-generated meshes were generated with configuration shown in Table 1. A web-based tool with a map panel was used to specify the target resolutions of the regions listed in the table. The CUMG algorithm begins from an icosahedron and performs a global multi-section and regional bisections according to the target resolutions. Finally, it performs another multi-section to widen the transition zone and vertex location adjustment to optimize the shapes of the cells in the transition zone. However, this version of CUMG algorithm has the effect of increasing the resolution in the areas around pentagons, including those at the 12 vertices of the original icosahedron (North pole, South pole, 5 on a parallel in the North Hemisphere and another 5 in the South Hemisphere), while it maintains the same vertex connectivities under Delaunay triangulation in vertex location adjustment. This undesirable effect may be eliminated in a future version of the CUMG algorithm. Still it is noteworthy that CUMG-generated meshes used in this paper have no obtuse Delaunay triangle in the triangulation. Each edge of the resulting Voronoi Tessellation intersects with exactly one edge of the Delaunay triangulation dual mesh, and hence the mesh is well-centered (Engwirda, 2017).
In real case testing, we use three meshes, simply called the 200 m, 1 and 3 km meshes. The 3 km mesh is the 60-3 km mesh downloaded from MPAS site and rotated to have the fine resolution domain centered at Hong Kong. They are shown in Figure 1.
With resolution transition in the mesh, undesirable waves and associated reflections are to be eliminated along with the simulation. The Smagorinsky horizontal diffusion scheme with coefficients listed in Table 2, which follows the default values in MPAS-A, was used for filtering of waves. The stability and validity of such settings are to be tested in the experiments. For MPAS-A, a height-based terrain-following vertical coordinate with smoothed coordinate surfaces (Klemp, 2011) is used. A set of default 55 zeta levels provided by MPAS-A is used in simulations and shown in Table A3 in the appendix.

The PBL Parameterization
Turbulence is a chaotic, unsteady, and multi-scale phenomenon and often occurs within the atmospheric boundary layer. A numerical simulation can only resolve the turbulence that is larger than several times of the size of the grid cells and the unresolved part is handled by the PBL parameterization.
For turbulence modeling, denoting w as vertical wind velocity, c as the physical quantity being studied, dash superscript as the unresolved component (fluctuation), and angle brackets or overbar as averaging, the term of turbulent flux ⟨ ′ ′ ⟩ is the extra term when applying Reynolds averaging and an extra condition has to be introduced to close the equation. Then, we need to specify the vertical profile of the turbulent flux, explicitly or via higher-order closures. One method is the eddy-diffusivity model (often called K-theory). Based on this idea, the YSU scheme (Hong & Pan, 1996; uses the K-profile to formulate the functional form of eddy diffusivity coefficient K c within the boundary layer and is reproduced as below for convenience, (1)  where c is variables such as u, v, θ, q and superscripts L and NL denotes local and non-local transport term respectively, K c is the vertical diffusivity coefficient, κ is Von Karman constant (usually set as 0.4), h is the height of the planetary boundary layer and w s is the velocity scale at the surface. The γ c is a correction term to the local gradient, and the ⟨( ′ ′ ) ℎ ⟩ ( ℎ ) 3 term is an asymptotic entrainment flux term at the inversion layer. Moreover, the YSU scheme considers the region above the mixed layer (z > h) to be a free atmosphere.
When the horizontal resolution is up to a few hundred meters, turbulence is partially resolved by the dynamic core, and the YSU scheme is modified to be a scale-aware scheme, namely the Shin-Hong scheme (or SH scheme) (Shin & Hong, 2015). In this study, the scale-aware SH scheme is applied.

Shin-Hong PBL Scheme
The SH scheme introduces a grid size dependency function P that suppresses the parameterized transport as the grid size decreases. The local term (Equation 7 in Shin and Hong (2015)) is shown below, for potential temperature, where Δ is the horizontal resolution, K H is the vertical diffusivity coefficient of heat, P L is a grid-size dependency for the local term and the superscripts S(Δ * ) and L denote subgrid-scale (SGS) and local respectively, shown by the red line in Figure 2. On the other hand, a new modified non-local term ⟨ ′ ′ ⟩ is reformulated by considering three linear profiles as surface layer cooling, mixed-layer heating and entertainment to fits LES results and here this new modified non-local term is labeled as ⟨ ′ ′ ⟩ . Meanwhile, it is also multiplied a grid-size dependency (blue line in Figure 2), where C CS is a stability function and is an empirical function (Shin & Hong, 2013) fit into the reference data for a free convection case and NL denotes nonlocal. Local and non-local momentum terms are both formulated by the same grid-size dependency function P u (green line in Figure 2) as The SH scheme was verified with Large-Eddy Simulation output in the convective boundary layer (CBL) at the gray-zone resolution as a workable scale-aware PBL scheme (Shin & Hong, 2015).

Static Input Data
In meteorological models, land surface processes are one of the crucial issues to conduct an accurate NWP, particularly within the planetary boundary layer. Meteorological variables of the planetary boundary layer are strongly affected by land surface states and surface energy fluxes. MPAS-A uses the same geographical static data from the WPS of WRF. The highest resolution available is 15 s (∼450 m). It is necessary to use higher resolution Figure 2. The grid-size dependency functions P for subgrid-scale local (red) and non-local (blue, C CS = 1) vertical heat transports. P u is for local and non-local momentum terms.
static data for simulations with sub-kilometer resolution. In this study, land use, soil type (both embed land/ water mask), and elevation are prepared from 3 s (∼90 m)-resolution static data sources for the 200 m resolution experiments. The source of dominant vegetation category data is MODIS 2019, the soil type data is from Beijing Normal University and terrain data is from Shuttle Radar Topography Mission. Nevertheless, raster soiltype and land cover/land use data often have wrong pixels along the coast. We implemented a heuristic algorithm which takes the latest coastline information from Open Street Map and extends nearby inland data to the corrected coastline. The region using gray-zone resolution in Hong Kong and the updated static data are shown in Figure 3.
For terrain static data, take the highest mountain, Tai Mo Shan which has an elevation of 957 m, as an example. After populating terrain data into the 200 m resolution mesh, it is represented to have a terrain height of 915.7 m; populated to a 1 km resolution mesh, it is only represented to have a terrain height of 760.9 m. Terrain effects are expected to be better simulated using the 200 m resolution mesh compared to the lower resolution meshes.
For soil-type data which also embeds water/land contrast information, the new data set carries more accurate data for the narrow part of the Victoria Harbor (which was represented as connected lands in the original data set), the airport with a new runway, rivers, and reservoirs.

Weather Stations Near Coastline
Many cities in the world over which high-resolution models have potential applications, like Hong Kong, are at the coastline or riverside. Indeed, many weather stations are located within one or two km from the water front.
There is a practical need to provide numerical forecasts at the exact locations of weather stations, as a basis of district-level forecast dissemination, accuracy evaluation as well as building statistical post-processing models. In this study, barycentric interpolation involving the nearest three cells to map simulated scalars from grid to station locations is used. We define "landmask fraction" as the average landmask value (1.0 for land, 0.0 for water) weighted by the weights used in barycentric interpolation. All weather stations in Hong Kong studied in this paper are built on land, though some are on tiny islands or at the waterfront. Figure 4d shows the landmask fractions for stations having landmask fraction less than 1.0 in Hong Kong, which are more susceptible to landmask misrepresentation problems, leading to erratic micro-scale temperature and humidity diurnal behavior in simulations. It can be seen that for a 200 m mesh very often the landmask fraction is higher than that for the 1 and 3 km mesh, having cells near the stations treated as land.

Results of Idealized Tests
The Jablonowski and Williamson baroclinic wave test with initial conditions that leads to steady flow was conducted (Jablonowski & Williamson, 2006). A dry dynamical core was initialized with a zonal symmetric and a meridional formulation of wind. The mesh used in the test here has a region of 200 m resolution being a circle with 10 km radius centered at (22.30°N, 114.2°E) (i.e., Hong Kong) and the resolution changes gradually (1 km resolution coarsening per 12 km distance outward) to coarsest ∼100 km resolution, denoted as the HK200m-100 km mesh.
Along with time integration, we evaluate the numerical errors evolved due to HTS or artifacts in the computational grids if any. The model zonal wind fields were first interpolated onto a regular 0.225° × 0.225° latitude-longitude grid for plotting and error metric evaluation. Figure 5 shows the evolution of the zonal and meridional winds at the lowest level for the HK200 m-100 km mesh. The flows are zonal initially and meridional wind is zero. However, note that the MPAS-A framework and hence also CPAS uses C-staggered mesh that discretizes wind state at the edges of the unstructured grid, and regrids them as uReconstructZonal and uReconstructMeridional at cell centers back and forth as needed. Even before time integration, the first row of Figure 5 (which further undergoes regridding to the latitude-longitude grid) shows that there is tiny (<0.01 m/s) numerical artifact, which has highest magnitude around the tiny high-resolution area centered at Hong Kong. Nevertheless, tiny deviation from the steady state solution present in the grid formed disturbance in the simulation, which gradually evolved along day 0 to day 5.5 but its magnitude is small (<0.2 m/s). Baroclinic waves are formed near day 9.5 and the meridional components become as large as 1.5 m/s. The baroclinic wave has grown very large when approaching day 15.
The two root-mean-square error (RMSE) norms defined in (Jablonowski & Williamson, 2006) were evaluated: where the overbar denotes the zonal average and λ i , φ j , ζ k denote the longitude, latitude, vertical levels, respectively, given any indices (i, j, k) in the global regular grid. t denotes time. w j and Δζ k denote the summation weights, in which w j = | sin(φ j+1/2 ) − sin(φ j−1/2 )| and Δζ k = (ζ k+1/2 − ζ k−1/2 ). Here, we also substitute u by v for evaluating the norms for meridional wind in addition to zonal wind.
We compare the results of the above mesh with large resolution variation and HTS (HK200 m-100 km with 11 HTS levels), three cases with mesh generated by Lloyd-based algorithm and time-integrated without HTS ) (160-1 km-no-HTS, 60-3 km-no-HTS and 92-25 km-no-HTS), three global near-uniform meshes with HTS (2 HTS levels) (48 , 64 and 96 km) and a mesh NH24-SH48, having 24 km on Northern Hemisphere and 48 km on Southern Hemisphere, with HTS (3 HTS levels).
In the steady flow test, the baroclinic waves are produced when disturbance from numerical artifacts evolves. Figure 6 shows the result of pressure at the lowest level and relative vorticity at 850 hPa. As expected, two small-amplitude waves in the Northern and Southern Hemisphere start to appear on day 9.5. The amplitude of pressure and vorticity wave grows with time. The high-pressure and low-pressure regions on the sides of the waves with vorticity changes indicated the vortex formations along the baroclinic waves. The overall pattern Figure 5. Evolution of re-gridded zonal and meridional winds at the lowest levels from initial condition to day 15 for the HK200 m-100 km mesh. Each sub-plot may use a different color scales and the values marking the ticks in the color bar must be carefully noted. From Figures 7a and 7b, it can be seen that the Norm-1 curve of zonal wind for the HK200 m-100 km mesh (the black dashed one) has values smaller than meshes generated by Lloyd-based methods but slightly larger than those of simple near-uniform or low-variation resolution meshes, and reaches 1 m/s at around the 15th day. The order of magnitude is less than 0.01 m/s for the first 7 days. The Norm-1 of meridional wind is similar (Figures 7d  and 7e).
For Norm-2 in Figure 7c, the RMSE of zonal wind to the initial state remains less than 0.1 m/s for the first 13 days, after which its value increases rapidly. Most of the meshes are in a similar sequence, except that the 160-to-1 km mesh is a bit higher. The Norm-2 of meridional wind is one order of magnitude smaller than that of zonal wind (Figure 7f). Norm-2 measures the deviation of the resulting wind speed from analytical solution along time. With the zontal and meridional Norm-2 being small from start to the 13th day, it can be inferred that the horizontal components of momentum and kinetic energy are nearly conserved in this steady flow case and the NWP model is good for practical use.

Model Configurations and Experimental Setup
In this part, the model physics scheme options include the SH PBL scheme mentioned above, Monin-Obukhov surface layer scheme (Jiménez et al., 2012), the WSM6 microphysics scheme (Hong & Lim, 2006), the new Tiedtke (Bechtold et al., 2014;Zhang & Wang, 2017) cumulus parameterization scheme, the orographic gravity wave drag parameterization (Choi & Hong, 2015) which is turned off over locations where the grid resolution is finer than 10 km, the Xu and Randall cloud fraction for radiation scheme (Xu & Randall, 1996), the Noah land surface model (Chen & Dudhia, 2001a, 2001b, and Rapid Radiative Transfer Model longwave and shortwave (Mlawer et al., 1997).

Observational Data and Evaluation Methodology
For all cases, observational data from the surface weather stations of the Hong Kong Observatory are utilized for result verification. There are 36 (surface temperature), 14 (relative humidity) and 28 (wind) stations distributed over Hong Kong (in Figure A1), and as listed in Tables A1 and A2. Barycentric interpolation is used to obtain the forecast data at station locations. After that, the CPAS output and the HKO data are compared for calculating statistics.

Numerical Stability Test With Long Forecast
A high resolution forecast using meshes with large resolution variation is a challenging task and we push the limit up locally to the ultra-high resolution of 200 m. In this high resolution variation, the first issue we have to check is the overall numerical stability. The idealized test in the above section was carried out without terrain and other static data, nor physics parameterization schemes. In this section, real simulations are carried out for long forecasts to confirm its validity and the software program is repeatedly run in an HPC environment for behavior checking.
The Global Forecast System (GFS) 0.25° data produced by the National Centers for Environmental Prediction at date-time 2021-05-08 00UTC was chosen as the initial condition. The weather conditions during 8 May 2021 to 17 May 2021 were mainly sunny, accompanied with a subtropical high belt covering the South China Sea that brought fine conditions and hot weather to Hong Kong. Under the influence of an anticyclone aloft in the subtropical high, 8 to 10 May were mainly sunny. On 11 May, the anticyclone weakened a little. A southerlies airstream affected the local area, so it was partly sunny on 11 to 14 May with some showers in regional areas and even a few isolated thunderstorms on 13 and 14 May. As the subtropical ridge re-strengthened, it became mainly sunny on 15 to 17 May.
The first verification concerns reproducibility: 3 simulations with 3.5 days forecasts are executed with the same initial condition and static data using the same number of CPUs, all simulations are successfully done and the root-mean-square difference of surface temperature (2 m above ground, t2m) among simulation results are all 0 degree Celsius up to 16 decimal places. In addition 5.5 and 9.5 days forecasts were carried out for testing general model behavior. Checking surface temperature, the first 3.5 days of the 5.5 and 9.5 days forecast are the same as the 3.5 days forecast, and the first 5.5 days of the 9.5 days forecast is the same as that of the 5.5 days forecast as well (up to 16 decimal places).

Hourly Surface Temperature
Next, we compare the 200 m-resolution with 1 km-and 3 km-resolution long forecasts. In the simulations using the 3 meshes introduced in 2.1, the first 12 hr were treated as spin-up and were excluded from analysis. Hourly time-series of surface temperature at each station are produced and compared against observational data. Figure 8 shows two representative stations, (a) Hong Kong Observatory (Tsim Sha Tsui) which is located in the urban area but less than 1 km far from water front, and (b) Tai Mo Shan which is the station with highest altitude.

Table 3 (a) The Average Correlation Coefficient of Surface Temperature of 200 m-, 1 km-, and 3 km-Resolution, and (b) Percentage Difference of Correlation Coefficient Between 200 m-Resolution and Coarser Resolutions (1 and 3 km) Over 36 Stations
shows that the correlation coefficients are improved along refining resolution from 3 km, to 1 km, and 200 m. The correlation coefficients from 200 m mesh is 1%-3% higher than those from 1 km mesh, and are 3%-11% higher than those from the MPAS-A 60-to-3 km mesh.

Daily Minimum, Maximum and Mean Surface Temperature
Daily min, max and mean surface temperature of district stations are the common form of summary local citizens pay attention to, and hence they are evaluated here. Scatter plots of the 36 stations times 9 days predicted versus observational data are shown in Figure 9a. It can be seen that the 200 m resolution simulation results (red dots) are closer to observed data than those based on the 1 km (blue dots) and 3 km (green dots) resolutions. Figure 9b shows that 200 m simulation has smaller RMSE than the 1 and 3 km counterparts. In particular, the 3 km resolution has the largest RMSE of 3.82°C for daily maximum temperature, which is reduced to 2.95°C in 200 m resolution. The two representative stations in Figure 8 can be checked to indicate the source of error. Typically, 200 m resolution resolves land-water contrast around stations much more satisfactorily, and the rise of temperature in the afternoon is more realistically simulated near the vicinity of land-water boundary. 200 m resolution also resolves terrain altitude around stations more accurately and avoid positive temperature bias for stations on mountains in the 3 km resolution.

High-Resolution Statics Data
We now investigate the performance of 3 s static data mentioned in Section 2.3 on 200 m resolution simulation as compared with 15 s static data. Eight episodes without prominent mesoscale weather feature (otherwise the  All experiments were initialized at 00 UTC of the previous day by the GFS 0.25 × 0.25 data and were integrated for 40 hr. The first 16 hr were treated as spin-up and excluded from our analysis. The output interval is 10 min. The surface temperature (t2m) is most concern meteorological variable daily and its bias and root mean square error are shown in Figure 10. It can be revealed that using 3 s static data can improve the accuracy of 200 m resolution on surface temperature. This 3 s static data is used in all other experiments of this paper.

Comparison Among 200 m, 1, and 3 km Resolution
The experiments for the 8 sunny day episodes are repeated with the 200 m mesh replaced by the 1 and 3 km mesh mentioned in Section 2.1. The sensitivity to resolution is studied for surface temperature and relative humidity.

Surface Temperature
The bias and RMSE for each station is shown in Figure 11, sorted according to the 3 km-resolution result so that the improvement by 1 km and 200 m resolution can be easily seen. Note the 4 stations near the right hand side,  The correlation coefficient of the 10-min interval, 24-hr time series of each station for each episode are calculated. Averaging over stations and episodes, the correlation coefficient from 200 m-mesh is 0.894, 1 km-mesh is 0.874 and 3 km-mesh is 0.824. It can be said that the correlation coefficient resulting from temperature prediction with 200 m resolution is 2.36% higher than that resulting from 1 km resolution and also is 8.53% higher than that resulting from an MPAS 60-3 km mesh.

Relative Humidity
In this part, we investigate the relative humidity at 2 m above ground that is calculated from 2 m specific humidity, surface temperature and surface pressure in MPAS-A/CPAS. Time-averaged surface relative humidity, temperature and specific humidity are shown in Figure 14. The spatial variations of relative humidity are contributed from both temperature and moisture concentration variations.
Sorting simulation bias and RMSE and analyzing the locations (details in Appendix B), it is found that the observation stations can be divided into two groups shown in Figure A2, near-water (14 red dots) and inland (6 blue dots) stations. Figure 15 shows the diurnal pattern of two groups. The near-water group shows smaller temporal variability than the inland group, especially for the 3 km model, which could not simulate the drop in relative humidity in the afternoon well. This can be understood as an artifact that the barycentric interpolation in the near-water stations involves more cells regarded as sea in the land surface model (Figure 14 3 km-resolution column). Since the relative humidity is consistently higher above the sea, hence there is an overestimation with reduced diurnal variation. As shown in Table 4, the 3 km model has significantly smaller correlation than the 1 km and 200 m counterparts.
On the other hand, we find that the relative humidity by the 3 km simulations is usually higher than the remaining two higher-resolution simulations. First, the moister conditions by 3 km mesh can be explained by the advection from the sea. For the cases chosen, the local region is either affected by mild easterlies or south-westerlies in  However, it can be seen in Figure 16 that the moisture brought by the breeze in 3 km mesh can travel into deeper inland regions than that in 200 m mesh, especially for the mountain ranges around the Tai Mo Shan, while their moisture contents in the open sea do not differ significantly. Such a phenomenon is likely due to 3 km's coarser features of resolving the topography. Second, it can be observed that coarser-resolution tends to give a moister condition to the lower levels. In Figure 17, the values of the water vapor mixing ratio by the 3 km mesh are usually higher in the lower levels, but sometimes lower in the upper levels, such as in the cases 2021-03-29, 2021-04-13, and 2021-05-09. Under 400 m from the surface, 200 m, 1, and 3 km are the usual order of the water mixing ratios.

Wind Prediction
The effect of resolving topographic features with fine resolution to winds is also studied. Complex terrain features, such as mountain ranges, can be an important factor to the local wind. Current global numerical weather models seldom use a fine resolution mesh to simulate the local or regional meteorology, however, such an effect under a global model will be attempted here. The purpose of this experiment is to study the performance of wind prediction by the 200 m-resolution, by comparing it to the 1 km-resolution.
In Hong Kong, persisting windy scenarios over a day can be due to tropical cyclones or the monsoon. We selected six scenarios with windy conditions shown in Table 5. Four of the cases are tropical cyclones, while the remaining two are monsoon cases. Different wind directions happened in the cases and thus the effect of topographic features on winds can be studied.
The purposes of the verification are to investigate the wind prediction enhancement due to resolved terrain altitude at the station and the sheltering effect to the winds by the surrounding mountains. Observation data from a maximum of 28 stations in each case will be used for verification (See Table A2). Both the output data from the models and the observation data are smoothed to hourly averaged values, instead of instantaneous values.

Wind Speed
The RMSEs of wind speed prediction are shown in Figure 18a. It can be seen that the accuracy of 200 m has improved from the 1 km simulation in all cases.
In Figure 18b, we can see that the 200 m-resolution and 1 km-resolution simulations show a little difference in ranges of wind speed, with the 200 m one having lower values in general. It indicates that the overall intensities of simulated winds are similar to each other, meaning the accuracy in predicting winds in local feature has improved without changing simulated mesoscale meteorological features.

Mountain Top
In fact, for some stations, the correctness of the resolved terrain altitude of the station directly improves the accuracy of predicted wind speed. Higher mountain tops usually have higher wind speed conditions that are underestimated in simulations. Resolving the altitude of the mountain top closer to its reality in static data helps predict the wind more accurately. Figure 19 illustrates the relationships between wind speed biases and biases of the anemometers' altitude. Take the Tate's Cairn Meteorological Station as an example. It is located at nearly the Figure 20 also shows the spatial plots of simulated 10-m wind speeds and resolved terrain altitudes by the two models near the Tai Mo Shan, the highest peak in Hong Kong with an elevation of 957 m. It can be seen that the 200 m-resolution simulation results in very high wind speed along several ridges spanning from the central mountain peak, and lower wind speed in the interlacing valleys, an apparent realism that 1 km-resolution simulation cannot provide.

Sheltering Effect
Apart from better resolving the altitude of stations near mountain tops, a local resolution model to 200 m can improve the wind prediction associated with the mountain blocking effect. To study the effect of mountain blocking, first, for each station, we divide the horizontal space into four quadrants, NW, NE, SE and SW respectively. We define a neighborhood distance, 1,500 m, which is actually picked arbitrarily, and look up the terrain data (ter) for the highest cells in each quadrant within the neighborhood. We define "mountain displacement" as that    Figure 21a illustrates such definition.
We pick several stations with obvious blocking topography, and show the plot of wind speed bias versus "mountain displacement" for 1 km-resolution and 200 m-resolution simulations in Figures 21b and 21c respectively. The data in each of the quadrants only present the period of the corresponding prevailing wind directions. The 2018-09-16 case out of the six cases is not included due to the high proportion of unavailable data for these stations. Comparing each corresponding quadrant in subfigure b and c, from a resolution of 1 km to 200 m, we find that the 200 m model is more capable of resolving the surrounding blocking mountains as expected and reflected by larger "mountain displacement," which often results in smaller wind speed biases.
The station Wong Chuk Hang located at the southern part of Hong Kong Island is a typical example that the complex terrain resolved by 200 m resolution makes improvement in prediction. As seen in Figure 22, the wind speed biases are reduced by 2-7 m/s. In fact, the station is surrounded by up to 400 m hills, however, the 1 km-resolution mesh is seriously misrepresenting them (compare the color contour of terrain in Figures 22a  and 22b). The time-series in Figure 22c shows that the positive bias of wind speed by 1 km-resolution simulation is reduced in 200 m-resolution simulation.

Wind Direction
More than half of the stations show improvement in wind direction forecast from 1 km to 200 m simulation. Take 2017-08-27 as an example. Figure 23 shows that the 200 m simulation has reduced the RMSE in mountain-top stations and the stations having mountain-blockage effects, though the northern part of Hong Kong shows a little degradation of prediction in wind direction. The overall RMSE of the wind direction for the day by 200 m simulation is 34.3°, while that by the 1 km simulation is 35.6°.

Discussion and Conclusion
Ultra-high resolution is a challenging problem in NWP, and it is related to all aspects of modeling including numerical stability, parameterization schemes and static data. In this study, we verify the application of locally 200 m resolution on the CPAS model, which is based on MPAS-A and uses unstructured variable resolution global meshes. The advantage of CPAS on locally ultra-high resolution is OLAM-based customized mesh generation that allows steep yet smooth resolution transition with no ill-staggering, and CPAS' hierarchical timestepping that saves computational cost and makes such simulation possible under practical computational cost constraint. Moreover, scale-aware Shin-Hong PBL scheme is integrated to support gray-zone resolution.
In this study, first, a steady flow dry baroclinic wave test is carried out to verify the numerical behavior of the HTS-enabled dynamical core. Then, the full set of parameterization schemes and static data are incorporated together for test-running 9.5 days long forecast with real world initial condition, after which verification against observation data is done. Meanwhile, we have examined elements of the static data such as land use, soil category and terrain height, and verified that the updated, higher-resolution data is beneficial to the accuracy of 200 m-resolution simulation. The heuristic algorithm for handling land-water boundary in the manipulation of upstream geographical data sources and populating them to unstructured mesh is verified to be effective for Hong Kong, and is applied to the global data in the CPAS platform. We have also compared 200 m resolution simulation results with 1 and 3 km resolution on surface temperature, relative humidity and wind. Overall, the performance of 200 m resolution is better. In addition, mountain top and sheltering effect on wind speed and wind direction are also analyzed simply to demonstrate special wind simulation benefits due to high resolution.
While most previous ultra-high-resolution NWP is performed by regional model with domain nesting, CPAS conducts it in a single global model with a unstructured grid with resolution spanning 2 orders of magnitude (∼52 km -200 m) which yet has smooth resolution transition (no ill-staggered mesh cell). This opens up theoretical interests as well as real-world applications. Being an initial value problem, such locally ultra-high resolution model does not depend on an external driving model (to supply boundary condition data for regional model), hence the model can be considered self-sustainable if a data assimilation cycle can be implemented with observation data supply. It may have long validity of locally high-resolution forecast like 9 days, that an externally driven regional model cannot attain. Without long validity and high resolution together, sub-city scale forecast providing individual numerical temperature, humidity and wind forecast for each small district heavily relies on statistical method like Model Output Statistics or Kalman filter to generate individual adjustment on top of a moderate-resolution NWP model. With a locally ultra-high resolution global model, the NWP may provide reasonable variability across the districts itself, reducing the needs or magnitude of adjustment by statistical methods. For less frequently occurring weather conditions, which are often severe weather events the community demands accurate forecast, statistical methods do not have enough training samples and often cannot output the variability that a physically based model may simulate. On the other hand, the apparent realism of wind simulation can have great potential application to wind farm prediction, especially for those in complex terrain. The vertical profile of simulated wind (10 -∼400 m above ground) should be further studied. On the other hand, rain prediction was not covered in this paper. Further research on ways ultra-high resolution simulation may improve the forecast skill of rain should be pursued.

Table A3
Heights of Defined Zeta Levels Used in the Model Figure A2. The stations for relative humidity. Near-water stations are red dots and inland stations are blue dots.

Appendix B: Details of Relative Humidity Results
For inland group, the bias of each station, sorted according to the 3 km model, is shown in Figure B1a. The diurnal pattern of the bias is shown in Figure B1b. It is seen than 200 m and 1 km underestimate the relative humidity for these 6 inland stations and over the whole diurnal cycle. The inland group has RMSE 11.31%, 9.36%, and 7.12% for 200 m, 1 and 3 km resolution respectively, with the high resolution ones contributed mostly by the bias.
For the near-water group, the bias and RMSE for each station is shown in Figure B2, sorted according to the error in 3 km model. The diurnal pattern and overall values of bias and RMSE of relative humidity are shown in Figure B3. The 3 km model has positive bias. The bias of the 200 m model is closest to zero, and the RMSE is the smallest.

Data Availability Statement
Model experiments for this research are conducted at the high-performance computing cluster provided by Clus-terTech Limited. All raw data for reproducing this study is publicly accessible in the following ways. The raw data for all simulations are available on Zenodo (https://doi.org/10.5281/zenodo.6837936) . In addition, the raw data for clustertech platform for atmospheric simulations (CPAS) simulations are stored on the CPAS cloud platform (https://cpas.earth/) and can be accessed with its Jupyter visualization tools and user interface (see How-to-visualize-CPAS-model-data.pdf on Zenodo).