Stripe patterns in wind forecasts induced by physics-dynamics coupling on a staggered grid in CMA-GFS 3.0

An unphysical stripe pattern is identified in low-level wind field in China Meteorological Administration Global Forecast System (CMA-GFS), characterized by meridional stripes in u-component and zonal stripes in v-component. This stripe noise is primarily confined to the planetary boundary layer over land. The structural mismatch between static field variations and the observed 2Δx noise amplitude suggests that locally forced mechanisms from surface inhomogeneity alone cannot explain the wind stripe patterns. Meanwhile, pure dynamical core simulations exhibit no such noise, confirming that the dynamical core itself does not generate these patterns. These results suggest that staggered-grid mismatch in physics-dynamics coupling is likely the primary mechanism. Idealized two-dimensional experiments demonstrate that combining one-dimensional dynamic-core advection and physics-based vertical diffusion on a staggered grid generates 2Δx-wavelength spurious waves when surface friction is non-uniform. One-dimensional linear wave analysis further confirms that staggered-grid coupling between dynamic advection and inhomogeneous damping forcing induces dispersion errors in wave solutions. Sensitivity tests validate that eliminating grid mismatch in physics-dynamics coupling removes this stripe noise. These findings collectively indicate that while staggered grids benefit the dynamic core's numerical stability and accuracy, their inherent grid mismatch with physics parameterizations requires specialized coupling strategies to avoid spurious noise. Potential solutions to remedy this issue are discussed.

1 Introduction

The wind field, as a fundamental state variable in atmospheric dynamics, governs not only the accuracy of circulation forecasts but also controls the spatiotemporal distributions of temperature, humidity, aerosols, and cloud microphysical properties through advection processes. Numerical weather prediction (NWP) systems simulate three-dimensional wind field evolution through numerical integration of atmospheric governing equations with multi-scale physical parameterizations. The momentum equations describe wind field evolution through three key forces: the pressure gradient force, Coriolis force, and nonlinear advection terms. The numerical discretization of advection terms is particularly important as it directly affects solution stability (Durran, 2010). Parameterized physical processes at subgrid scale significantly influence wind field predictions. Surface stress not only exerts substantial control over the wind intensity in both the near-surface and planetary boundary layers (PBL) through turbulent momentum transport (Blackadar, 1957; Stull, 1988), but also significantly influences the location of the surface westerlies and upper-level jets in the midlatitudes (Robinson, 1997; Chen et al., 2007). When flow encounters subgrid-scale orography, it excites gravity waves that propagate vertically, altering large-scale winds through momentum flux divergence (McFarlane, 1987; Kim et al., 2003; Kim and Doyle, 2005). Mesoscale topographic blocking and turbulent form drag modify low-level wind fields, which in turn affect global circulation patterns (Beljaars et al., 2004; Sandu et al., 2016). In addition to dynamics and physical processes, the consistency of dynamic-physics coupling is equally critical for overall model performance (Bauer et al., 2015; Gross et al., 2016, 2018).

In NWP models, numerical noise typically stems from non-physical high-frequency oscillations induced by dynamical process discretization (Wurtele, 1961; Arakawa, 1966; Mesinger and Arakawa, 1976). Additionally, spatial grid mismatch in physics-dynamic coupling can generate computational noise (Chen et al., 2020). While consistent spatial grids between dynamical and physical processes are theoretically preferable for seamless coupling, practical implementations often employ differing horizontal and vertical grid configurations. For vertical discretization, two conventional grid arrangements are employed: (1) the Lorenz grid (L-grid; Lorenz, 1960), which collocates thermodynamic variables (e.g., potential temperature and moisture) and horizontal wind components at the half-levels (mass levels), with the vertical velocity and geopotential typically located at the full-levels (interface levels); and (2) the Charney–Phillips grid (CP-grid; Charney and Phillips, 1953), which staggers the thermodynamic variables at the full-levels above and below the corresponding half-levels of the horizontal wind, while the vertical velocity is also defined at the full-levels. Chen et al. (2017) demonstrated that Lorenz-grid physics coupled with CP-grid dynamics generates non-physical computational modes in PBL scheme solutions. Chen et al. (2020) confirmed that unified grid coupling significantly improves the stratocumulus representation and overall forecast skill while eliminating this spurious noise. These findings underscore that physics-dynamics coupling inconsistencies not only introduce numerical noise but also substantially degrade model performance.

In global models, physical processes are typically computed in single columns, independent of horizontal grid configuration. Thus, horizontally distributed noise cannot originate from physical parameterizations if the underlying static physical data is noise-free. In two-dimensional rectangular grids, horizontal grid configurations are categorized into non-staggered (A-grid) and staggered (B grid to E grid) types, distinguished by their arrangement of wind components (u, v) and mass variables (height z or pressure p) (Arakawa and Lamb, 1977). Among these, the Arakawa C-grid – where mass variables are collocated at grid centers while velocity components are staggered at cell interfaces – exhibits optimal accuracy when the grid spacing resolves scales smaller than the Rossby radius of deformation (Batteen and Han, 1981; Xu and Lin, 1993; Randall, 1994). This superiority stems from its inherent conservation properties and reduced numerical dispersion in simulating geostrophic adjustment processes (Arakawa and Lamb, 1977; Randall, 1994). The Arakawa-C grid has been widely implemented as a preferred discretization framework in structured-grid NWP systems, including the Weather Research and Forecasting (WRF) Model's Advanced Research core (Skamarock and Klemp, 2008) and operational systems such as the UK Met Office Unified Model (Walters et al., 2017), the NOAA Global Forecast System (GFS) (Sela, 2010), the Environment and Climate Change Canada (ECCC) Global Environmental Multiscale (GEM) model (Côté et al., 1998), and the China Meteorological Administration's (CMA) operational GRAPES system (Chen et al., 2008). Figure 1 illustrates the wind field distribution on an Arakawa C-grid, where the physics-parameterized winds are located at the central mass points while the dynamically computed winds are staggered at the cell interfaces. This inherent wind field grid mismatch is typically handled via two-step interpolation during physics-dynamics coupling. However, as demonstrated by Chen et al. (2020), the grid inconsistency during physics-dynamics coupling introduces spurious computational noise.

Figure 1 Wind configuration in Arakawa C grid; the subscript “p” denotes the physics grid points.

High-order horizontal diffusion (Xue, 2000) and filters (Shapiro, 1970, 1975; Raymond and Garder, 1988) are standard techniques for suppressing high-frequency computational noise in NWP models. These approaches were successfully implemented in the GRAPES model and effectively mitigated small-scale noise (Wang et al., 2008). Compared to CMA-GFS 2.4 (originally GRAPES_GFS V2.4), the new CMA-GFS 3.0 has removed fourth-order horizontal diffusion to improve small-scale system simulation capabilities (Shen et al., 2023). Following the upgrade to CMA-GFS 3.0, stripe patterns became evident in the low-level wind fields – an issue absent in previous versions.

This study systematically investigates the underlying causes of these stripe patterns in CMA-GFS 3.0, particularly examining physics-dynamics coupling mismatches as the probable source. Section 2 details physics-dynamics coupling scheme on a staggered grid for the wind field prediction in CMA-GFS 3.0 and characterizes the observed horizontal wind noise patterns. Section 3 identifies noise sources via idealized experiments and analytical wave solutions, revealing fundamental discretization constraints in physics-dynamics coupling. Section 4 designs sensitivity experiments to verify the origin of horizontal wind noise in CMA-GFS 3.0, conclusively attributing it to physics-dynamics coupling mismatches. Conclusions and discussions are presented in Sect. 5.

2 The wind prediction in the PBL in CMA_GFS

2.1 Model description and physics-dynamic coupling in wind prediction

The CMA-GFS (formerly known as GRAPES_GFS) was developed based on a two-time-level semi-implicit semi-Lagrangian (SISL) non-hydrostatic dynamic core. Its development began in 2007, and the system became operational in late 2015 (Shen et al., 2020). The model employs a horizontal Arakawa C-grid on a regular latitude-longitude grid, combined with a vertical CP grid using a height-based terrain-following coordinate system (Chen et al., 2008).

Physics package includes the Rapid Radiative Transfer Model for GCMs (RRTMG) for longwave and shortwave radiation (Pincus et al., 2003; Morcrette et al., 2008), the double-moment microphysics and an explicit prognostic cloud cover scheme (Ma et al., 2018), the common land model (CoLM) (Dai et al., 2003), the subgrid-scale orographic (SSO) scheme including gravity wave drag (GWD) (Chen et al., 2016) and turbulence orographic form drag (TOFD) (Beljaars et al., 2004), the new Simplified Arakawa Schubert (NSAS) convection, and the new medium-range forecast (NMRF) for vertical diffusion (Han and Pan, 2011). The implementation of vertical diffusion and cloud physics schemes on the CP grid eliminated computational noise induced by physics-dynamics coupling mismatches, leading to significant improvements in stratocumulus cloud simulations and overall forecast skill scores (Chen et al., 2020).

In CMA-GFS's dynamical core, horizontal winds $u_{i\pm \frac{\mathrm{1}}{\mathrm{2}},j}$ and $v_{i,j\pm \frac{\mathrm{1}}{\mathrm{2}}}$ are staggered at cell interfaces, while physics-parameterized winds u_p and v_p are collocated at mass center (Fig. 1). In our coupling scheme, linear interpolation bridges the grid mismatch for physics-dynamics coupling:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\displaystyle }P\left(u_{i,j}\right)& {\displaystyle }=P\left({\displaystyle \frac{u_{i-\frac{\mathrm{1}}{\mathrm{2}},j}+u_{i+\frac{\mathrm{1}}{\mathrm{2}},j}}{\mathrm{2}}}\right)\\ {\displaystyle }P\left(v_{i,j}\right)& {\displaystyle }=P\left({\displaystyle \frac{v_{i,j-\frac{\mathrm{1}}{\mathrm{2}}}+v_{i,j+\frac{\mathrm{1}}{\mathrm{2}}}}{\mathrm{2}}}\right)\end{array}\end{array}$$

where P denotes the physics parameterization operator. When feeding physics tendencies back to the dynamical core, the coupling follows:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\displaystyle }\mathrm{\Delta }{u}_{i+\frac{\mathrm{1}}{\mathrm{2}},j}^{\mathrm{P}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[P\left(u_{i,j}\right)+P\left(u_{i+\mathrm{1},j}\right)\right]\mathrm{\Delta }t\\ {\displaystyle }\mathrm{\Delta }{v}_{i,j+\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{P}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[P\left(v_{i,j}\right)+P\left(v_{i,j+\mathrm{1}}\right)\right]\mathrm{\Delta }t\end{array}\end{array}$$

where ΔV^P denotes the physics-induced wind increment to be fed back to the dynamic core. However, such interpolation-based coupling schemes often fail to mitigate computational noise induced by physics-dynamics grid mismatches (Chen et al., 2020).

As theoretically anticipated, lower-level winds of CMA-GFS 3.0 exhibit distinct stripe patterns. We therefore systematically investigate these anomalies in detail, first characterizing their spatiotemporal structures, and then evaluating potential links to physics-dynamics coupling mechanisms.

2.2 Spatiotemporal distribution of wind forecast noise: a case study

In this case study – hereafter referred to as Exp_Ctrl – we utilize CMA-GFS 3.0 with a horizontal resolution of 0.25° (approximately 25 km at the equator) and 87 vertical levels extending up to 1 hPa. The simulation spans from 12:00 UTC 1 July 2021, initialized with European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5 reanalysis data (Hersbach et al., 2020; data available at https://cds.climate.copernicus.eu/datasets/, last access: 3 November 2025). Hourly forecast fields are saved directly on the model's Arakawa C-grid at model levels throughout the 5 d simulation.

Global distribution of the wind field noise reveals pronounced land-sea contrasts, with predominant stripe patterns over land while relatively smooth over oceans. Figure 2 displays the 18 h forecast (valid at 06:00 UTC) of near-surface winds at the lowest model level over South Asia (5–34° N, 88–113° E), a region encompassing both land and ocean. Distinct stripe patterns dominate the wind field over land areas, where the u-component displays meridional stripes and the v-component exhibits zonal stripes. With the model's 0.25° resolution, these structures are clearly resolved with a wavelength of 0.5° (2Δx), demonstrating characteristic 2Δx wave properties.

Figure 2 18 h forecast of (a) u- and (b) v-components at the lowest model level from CMA-GFS 3.0.

Figure 3 displays vertical cross-sections of u- and v-components along 30° N and 100° E, respectively, for the region shown in Fig. 2. Pronounced 2Δx wave patterns are evident throughout the PBL in both components, whereas the wind field becomes markedly smoother in the free atmosphere above the PBL. The strongest 2Δx signatures are concentrated in the bottom of the PBL.

Figure 3 Vertical distribution of 18 h forecasted winds: (a) u-component along 30° N; (b) v-component along 100° E. Black solid lines denote the PBL height.

Figure 4 presents the 120 h (5 d) evolution of surface-layer u- and v-components along the zonal and meridional transects shown in Fig. 3. The 2Δx oscillations exhibit continuous presence throughout the period. Their amplitude shows strong diurnal variations but no systematic growth with forecast days, suggesting that this noise is unlikely to directly induce model instability. This persistent noise pattern shows daytime intensification and nighttime weakening, suggests a potential linkage to the inherent diurnal cycle of the development of PBL.

Figure 4 Temporal evolution of winds at the lowest model level (0–120 h): (a) u-component along 30° N, (b) v-component along 100° E.

The vertical structure of low-level winds and their diurnal variations in the surface layer imply a connection between the stripe patterns and momentum-related physical processes, particularly turbulent transport within the PBL. However, the single-column design of physical schemes permits PBL diffusion to affect such horizontal wind noise solely via surface friction forcing. In the following section, we systematically examine the mechanisms for these 2Δx stripe patterns.

3 Tracing the cause of wind stripes predicted in CMA-GFS 3.0

3.1 Examination of associated surface static fields

The surface static fields directly governing surface friction in the wind field are surface inhomogeneity descriptors, principally characterized by subgrid orographic variability and aerodynamic roughness length. Figure 5 displays the spatial distributions of subgrid orography standard deviation and momentum roughness length across the domain shown in Fig. 2.

Figure 5 Horizontal distributions of (a) Standard deviation of sub-grid orography; (b) Roughness length over the study region (domain same as Fig. 2).

The distributions of the two static fields in Fig. 5 differ markedly from the wind stripe patterns in Fig. 2, indicating that surface parameter inhomogeneity is not the direct origin of the stripe noise. However, visual comparison with Fig. 2 demonstrates a distinct spatial correspondence between the stripe noise locations and the surface roughness distribution: the 2Δx noise patterns predominantly occur over areas with higher roughness, while remaining absent over smooth surfaces. The Andaman Islands (11–13° N, 93° E), located between the Bay of Bengal and the Andaman Sea, generate localized surface inhomogeneity in this predominantly oceanic region. Consistent with this inhomogeneity, Fig. 2 exhibits obvious meridional stripes in the u-component precisely over this archipelago. These observations imply that although the surface inhomogeneity parameters themselves exhibit no stripe patterns – and thus cannot directly produce the wind stripes through physical processes – the spatial organization of 2Δx noise systematically correlates with surface inhomogeneity.

3.2 Dynamical core test: Noise assessment in physics-off simulation

Section 3.1 demonstrated that the surface-friction-related static fields lack horizontal noise patterns (Fig. 5), definitively excluding surface-mediated physical processes – which are computed in single columns – as possible sources of the observed 2Δx wind stripes. We now evaluate the dynamical core's potential contribution through a physics-free simulation using identical configuration described in Sect. 2.2. Figure 6 displays 18 h forecasts of near-surface u- and v-components over the same domain as in Fig. 2. The absence of surface friction yields stronger winds compared to the default physics-on simulation (Fig. 2), and more critically, the wind fields exhibit smooth distributions completely free of 2Δx stripes. This definitive evidence excludes the dynamical core as the source of the observed low-level wind noise.

Figure 6 Same as Fig. 2, but for the physics-off experiment.

Having definitively excluded both the dynamical core and physical parameterizations as direct sources of these horizontal noise patterns, the residual evidence points to the physics-dynamics coupling interface as the sole remaining candidate. As schematized in Fig. 1, CMA-GFS employs staggered grids between physics and dynamics components, requiring linear interpolation for wind field coupling. This inherent grid mismatch is now the most probable mechanism generating the observed 2Δx stripe noise. We proceed with idealized experiments and theoretical analysis to examine whether the grid mismatch between physics and dynamics could induce the computational noise.

3.3 Two-dimensional idealized test

The governing equation for the idealized simulation combining one-dimensional horizontal advection with vertical diffusion is

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial u}{\partial t}}+u{\displaystyle \frac{\partial u}{\partial x}}=-{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}}{\partial z}}\end{array}$$

where u is flow speed and $\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}$ the turbulent vertical momentum flux. In Eq. (3), the nonlinear term $u\frac{\partial u}{\partial x}$ represents the dynamical advection process, while the right-hand side $-\frac{\partial \stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}}{\partial z}$ represents the physical forcing process, specifically parameterized turbulent vertical transport. The momentum flux can be expressed as:

$$\begin{array}{}\text{(4)}& \stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}=-K{\displaystyle \frac{\partial u}{\partial z}}\end{array}$$

where K(z) is given by:

$$\begin{array}{}\text{(5)}& K\left(z\right)=\left\{\begin{array}{ll}\mathrm{10},& z\le H\\ \mathrm{0},& z>H\end{array}\right.\end{array}$$

and H=500 m denotes the prescribed boundary layer height.

The governing Eq. (3) is numerically solved using the discretization strategies and physical configurations summarized in Table 1.

Table 1 Configuration for idealized simulation.

The surface friction configuration in Table 1 creates deliberate localized forcing, enabling controlled investigation of non-uniform surface effects on wind noise generation mechanisms in the simulated boundary layer flow.

Based on the configurations in Table 1, we conducted two contrasting idealized numerical experiments:

Ideal_Ctrl (Control Experiment):

• Solves both horizontal advection and vertical diffusion terms at collocated grid points (x_i)

Ideal_Test (Test Experiment):

• Solves horizontal advection at x_i

• Computes vertical diffusion at staggered grid points ($x_i+\mathrm{\Delta }x/\mathrm{2}$),with surface friction velocity (u_∗) prescribed at corresponding staggered positions

• Couples solutions through linear interpolation between grids

The two experiments are designed to isolate the impact of grid staggering on noise generation under inhomogeneous surface forcing. Both experiments are simulated for 24 h (288 steps). Figure 7 displays the spatio-temporal distribution of the simulated u-wind from both cases. The Ideal_Ctrl simulation (Fig. 7a) exhibits localized wind reduction at the surface forcing site, with perturbations propagating downstream smoothly, demonstrating physically consistent behavior devoid of numerical noise. In contrast, the Ideal_Test simulation (Fig. 7b) shows that additional 2Δx grid-scale waves emerge near the local surface friction anomaly, superimposed on the expected frictional wind response. These 2Δx waves persist throughout the integration period with temporally quasi-constant amplitude, consistent with their behavior in CMA-GFS 3.0. Figure 7c and d demonstrates that the influence of surface friction is confined to the boundary layer (below the 50th model level). Comparing Fig. 7d with Fig. 7c, the staggered coupling generates 2Δx oscillation within the boundary layer at the point of non-uniform surface friction and these spurious noise patterns propagate to adjacent grid points, primarily near the surface. The 6 h simulated lowest level u-wind (Fig. 7e) confirms that these 2Δx waves represent spurious numerical oscillations induced by dispersive energy propagation, which are most pronounced adjacent to the inhomogeneous surface forcing region. The Ideal_Test simulation results in Fig. 7 successfully reproduce both the distinctive 2Δx noise and their characteristic spatial distribution observed in CMA-GFS 3.0.

Figure 7 Idealized experiment results (domain (grid point number): x-grid = 41–59; z-grid = 1–60): (a, b) Temporal evolution of winds at the lowest model level in (a) Ideal_Ctrl and (b) Ideal_Test; (c, d) Vertical structure at 6 h for (c) Ideal_Ctrl and (d) Ideal_Test; (e) 6 h wind at the lowest model level.

The idealized experiments reveal that the noise generation involves two coupled processes: (1) dispersive effects inherent to the advection-diffusion discretization (grid staggering), and (2) selective amplification and phase organization of small-scale components in surface forcing through these dispersive mechanisms. As demonstrated in Sect. 2.2, the stripe patterns in wind fields of CMA-GFS 3.0 predominantly occur over landmasses or near islands. This strongly indicates that wind stripe patterns in CMA-GFS 3.0 stem from inconsistent grid staggering in its physics-dynamics coupling. When combined with non-uniform surface forcing, this grid mismatch triggers dispersive wave propagation – ultimately producing the observed 2Δx stripe patterns.

3.4 Analysis of one-dimensional linear wave solution

For the one-dimensional linear damped wave equation:

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\partial u}{\partial t}}+c{\displaystyle \frac{\partial u}{\partial x}}=-\mathit{\alpha }u,\mathit{\alpha }>\mathrm{0}\end{array}$$

where c is the phase speed, −αu the damping term and α the damping coefficient.

In Eq. (6), the term $c\frac{\partial u}{\partial x}$ represents the advection term, while the damping term −αu acts as physical forcing. Based on the idealized experimental configuration outlined in the preceding section but with explicit damping time discretization, the right-hand side term of Eq. (6) can be discretized through two distinct numerical treatments. The first approach computes the damping term directly at the advection scheme's collocated grid points (x_j), giving

$$\begin{array}{}\text{(7)}& \mathrm{\Delta }{u}_j^{\mathrm{P}}=-\mathrm{\Delta }t{\mathit{\alpha }}_j{u}_j\end{array}$$

where the subscript j represents the value at grid point x_j. Equation (7) corresponds to the Ideal_Ctrl configuration in the idealized experiment. This lead to the solution of Eq. (6) as:

$$\begin{array}{}\text{(8)}& u_j=u_j^{\mathrm{D}}{\mathrm{e}}^{-\mathit{\alpha }t}\end{array}$$

where $u_j^{\mathrm{D}}=u_{\mathrm{0}}{e}^{ik(x_j-c^{\mathrm{D}}t)}$ is the wave solution of the advection term at x_j after spatiotemporal differencing. The specific form of phase velocity c^D depends on the finite difference scheme employed (Durran, 2010).

The alternative approach, analogous to the treatment in Eqs. (1) and (2), involves performing two averaging operations during the physics-dynamics coupling, similar to the Ideal_Test configuration. Consequently, the damping term in Eq. (6) leads to the following change in u at grid point x_j:

$$\begin{array}{}\text{(9)}& \mathrm{\Delta }{u}_j^{\mathrm{P}}=-{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}\left({\mathit{\alpha }}_{j-\frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{u_{j-\mathrm{1}}+u_j}{\mathrm{2}}}+{\mathit{\alpha }}_{j+\frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{u_j+u_{j+\mathrm{1}}}{\mathrm{2}}}\right).\end{array}$$

For numerically resolved scales ($\left|k\mathrm{\Delta }x\right|\ll \mathrm{1}$) where Taylor expansions remain valid, the variable u at neighboring grid points can be approximated as:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{u}_{j\pm \mathrm{1}}& =u_j\pm \mathrm{\Delta }x{\left.{\displaystyle \frac{\partial u}{\partial x}}\right|}_j+{\displaystyle \frac{\mathrm{\Delta }{x}^{\mathrm{2}}}{\mathrm{2}\mathrm{!}}}{\left.{\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial x^{\mathrm{2}}}}\right|}_j\pm {\displaystyle \frac{\mathrm{\Delta }{x}^{\mathrm{3}}}{\mathrm{3}\mathrm{!}}}{\left.{\displaystyle \frac{{\partial }^{\mathrm{3}}u}{\partial x^{\mathrm{3}}}}\right|}_j\\ & +\mathcal{O}\left(\mathrm{\Delta }{x}^{\mathrm{4}}\right),\end{array}\end{array}$$

where the notation ${\left.\frac{\partial u}{\partial x}\right|}_j$ denotes the partial derivatives of u at location x_j and $\mathrm{\Delta }x=x_{j+\mathrm{1}}-x_j$ is the Equidistant grid space. Substituting Eq. (10) into Eq. (9) and systematically neglecting fourth-order terms (𝒪(Δx⁴)), the individual wavelike solution to Eq. (6) becomes:

$$\begin{array}{}\text{(11)}& u_j=u_j^{\mathrm{D}}{\mathrm{e}}^{-{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{j}t}{\mathrm{e}}^{-i\frac{k\mathrm{\Delta }{\mathit{\alpha }}_j\mathrm{\Delta }x}{\mathrm{4}}t}{\mathrm{e}}^{\frac{k^{\mathrm{2}}{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_j\mathrm{\Delta }{x}^{\mathrm{2}}}{\mathrm{4}}t}{\mathrm{e}}^{i\frac{k^{\mathrm{3}}\mathrm{\Delta }{\mathit{\alpha }}_j\mathrm{\Delta }{x}^{\mathrm{3}}}{\mathrm{24}}t},\end{array}$$

where ${\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_j=\frac{{\mathit{\alpha }}_{j-\frac{\mathrm{1}}{\mathrm{2}}}+{\mathit{\alpha }}_{j+\frac{\mathrm{1}}{\mathrm{2}}}}{\mathrm{2}}$, $\mathrm{\Delta }{\mathit{\alpha }}_j={\mathit{\alpha }}_{j+\frac{\mathrm{1}}{\mathrm{2}}}-{\mathit{\alpha }}_{j-\frac{\mathrm{1}}{\mathrm{2}}}$. Compared to Eq. (8), Eq. (11) exhibits three additional terms in its wave solutions at x_j:

1. The term ${\mathrm{e}}^{-i\frac{k\mathrm{\Delta }{\mathit{\alpha }}_j\mathrm{\Delta }x}{\mathrm{4}}t}$ introduces a non-dispersive phase velocity modification.

2. The term ${\mathrm{e}}^{\frac{k^{\mathrm{2}}{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_j\mathrm{\Delta }{x}^{\mathrm{2}}}{\mathrm{4}}t}$ enhances the wave amplitude (equivalently reducing damping).

3. The term ${\mathrm{e}}^{i\frac{k^{\mathrm{3}}\mathrm{\Delta }{\mathit{\alpha }}_j\mathrm{\Delta }{x}^{\mathrm{3}}}{\mathrm{24}}t}$ generates wavenumber-dependent phase modifications, a characteristic signature of dispersive wave propagation.

The phase velocity in Eq. (11) combines both non-dispersive and dispersive modifications to the phase velocity of advection c^D as $c^{\prime }=c^{\mathrm{D}}+\frac{\mathrm{\Delta }{\mathit{\alpha }}_j\mathrm{\Delta }x}{\mathrm{4}}-\frac{k^{\mathrm{2}}\mathrm{\Delta }{\mathit{\alpha }}_j\mathrm{\Delta }{x}^{\mathrm{3}}}{\mathrm{24}}$.

When α exhibits no spatial variation (hence Δα_j=0), Eq. (11) reduces to:

$$\begin{array}{}\text{(12)}& u_j=u_j^{\mathrm{D}}{\mathrm{e}}^{-{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{j}t}{\mathrm{e}}^{\frac{k^{\mathrm{2}}{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_j\mathrm{\Delta }{x}^{\mathrm{2}}}{\mathrm{4}}t},\end{array}$$

where only the amplitude-modifying term persists, producing exclusively damping reduction without any phase velocity modifications (neither non-dispersive shifts nor dispersive effects).

Comparing Eqs. (8), (11), and (12) leads to the following conclusions:

1. The staggered-grid discretization scheme, where advection and damping terms are computed at offset grid points, fundamentally modifies wave solutions when compared to collocated approaches.

2. For spatially varying damping coefficients (α), this numerical framework introduces coupled amplitude-phase distortions, including wavenumber-dependent propagation speeds characteristic of dispersive systems.

3. In contrast, uniform α fields restrict the influence of staggered-grid effects solely to amplitude modulation, preserving the original non-dispersive wave kinematics.

The reformulated Eq. (9) can be expressed as:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}\mathrm{\Delta }{u}_j^{\mathrm{P}}& =-\mathrm{\Delta }t{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_j{u}_j-{\displaystyle \frac{{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_j\mathrm{\Delta }t}{\mathrm{4}}}(u_{j-\mathrm{1}}-\mathrm{2}{u}_j+u_{j+\mathrm{1}})\\ & -{\displaystyle \frac{\mathrm{\Delta }{\mathit{\alpha }}_j\mathrm{\Delta }t}{\mathrm{8}}}(u_{j+\mathrm{1}}-u_{j-\mathrm{1}}).\end{array}\end{array}$$

The right-hand side of Eq. (13) contains three distinct terms. The first term is identical to Eq. (7). The second term exhibits the form of second-order horizontal diffusion – note its negative coefficient, which actually suppresses small-scale fluctuations in the damping process itself and clearly cannot generate 2Δx oscillations. This term corresponds to the amplitude modification term in Eq. (11), reducing the damping influence. The third term represents a second-order centered difference scheme. As discussed by Durran (2010), second-order centered differencing introduces dispersion in solutions to one-dimensional linear wave equations, producing upstream-propagating noise near spikes – precisely as demonstrated in our idealized experiments. Evidently, the noise source we examine originates specifically from the third term of the right-hand side of Eq. (13).

Idealized experiments and linear wave analysis identify wind stripe noise in CMA-GFS 3.0 as requiring the coexistence of two independent factors: Physical forcing inhomogeneity (e.g., non-uniform surface friction) and physics-dynamics coupling on staggered grids. Critically, neither factor alone can generate dispersion noise – their combined presence is strictly necessary. This explains why stripes emerge preferentially over landmasses or near islands where sharp forcing gradients intersect with the model's inherent grid architecture (Fig. 2).

4 Confirming the origin of wind stripe patterns in CMA_GFS 3.0: sensitivity experiments

The preceding analysis demonstrates that low-level wind stripe noise over terrestrial region and near islands in CMA-GFS 3.0 stems from horizontal grid mismatch in physics-dynamics coupling. To validate this mechanism, we conduct a targeted sensitivity experiment (Exp_Test) by implementing conformal grid alignment for coupled processes.

In Fig. 1, the physical-point u-component can be expanded as $u_{\mathrm{p}}(x+\frac{\mathrm{\Delta }x}{\mathrm{2}})=u\left(x\right)+\mathcal{O}\left(\mathrm{\Delta }x\right)$. In Exp_Test, we retain solely the zeroth-order approximation, enforcing $u_{\mathrm{p}}(x+\frac{\mathrm{\Delta }x}{\mathrm{2}})=u\left(x\right)$. The same treatment is applied to the v-component ($v_{\mathrm{p}}(y+\frac{\mathrm{\Delta }y}{\mathrm{2}})=v\left(y\right)$), ensuring consistent physics-dynamics coupling for both horizontal velocity components. In correspondence with Eqs. (1) and (2), we now present the modified formulations:

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{\displaystyle }P\left(u_{i,j}\right)& {\displaystyle }=P\left(u_{i-\frac{\mathrm{1}}{\mathrm{2}},j}\right)\\ {\displaystyle }P\left(v_{i,j}\right)& {\displaystyle }=P\left(v_{i,j-\frac{\mathrm{1}}{\mathrm{2}}}\right)\end{array}\end{array}$$

When mapping from physics point back to the dynamics point, we have:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{\displaystyle }\mathrm{\Delta }{u}_{i-\frac{\mathrm{1}}{\mathrm{2}},j}^{\mathrm{P}}& {\displaystyle }=P\left(u_{i,j}\right)\mathrm{\Delta }t\\ {\displaystyle }\mathrm{\Delta }{v}_{i,j-\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{P}}& {\displaystyle }=P\left(v_{i,j}\right)\mathrm{\Delta }t\end{array}\end{array}$$

This conformal alignment ensures dynamical and physical processes are computed at coincident grid points, eliminating the need for interpolation-based coupling shown in Eqs. (1) and (2). The experimental configuration replicates the Exp_Ctrl setup described in Sect. 2.2, maintaining identical initial conditions and physical parameterizations.

Compared to the results from Exp_Ctrl (Fig. 2), Exp_Test maintains small-scale variability in terrestrial wind fields while completely eliminating 2Δx stripe patterns (Fig. 8). The residual inhomogeneity over land, as evidenced in Fig. 8, reflects authentic physical responses to subgrid surface forcing inhomogeneity (e.g., topographic roughness, land-use variations), distinguishing it from numerical dispersion noise.

Figure 8 Same as Fig. 2, but for Exp_Test.

Spectral analysis serves as a rigorous tool for quantifying energy distribution across frequencies. To diagnose the impacts of conformal physics-dynamics coupling on wind stripe patterns we conduct one-dimensional energy spectrum analysis separately for the u-component (zonal direction) and v-component (meridional direction). We selected the East Asian domain (70–145° E, 10–65° N) as a representative non-uniform surface region and the tropical Pacific (160° E–120° W, 20° S–20° N) as a typical uniform surface region for comparative analysis. The u-component was processed as a one-dimensional east-west oriented sequence for Fast Fourier Transform (FFT) analysis, while the v-component underwent identical treatment along the north-south axis. Through spectral analysis, the power spectral density (PSD) was estimated by computing the squared modulus of the FFT coefficients (Fig. 9).

Figure 9 PSD of (a) u- and (b) v-components at 18 h over tropical Pacific (TP: 20° S–20° N, 160° E–120° W) and East Asia (EA: 10–65° N, 70–145° E) domains for Exp_Ctrl and Exp_Test.

The key results of PSD investigation can be summarized as follows:

1. Over land, both u- and v-components in Exp_Ctrl demonstrate considerably strong energy at small scales (high wavenumbers), with no energy decay observed from 4Δx to 2Δx. In contrast, Exp_Test achieves significant suppression of small-scale energy, effectively mitigating high-frequency computational noise in terrestrial regions.

2. Over oceanic regions, the homogeneous ocean surfaces produce broadly similar spectral characteristics in both experiments. Although this general spectral agreement suggests that the horizontally homogeneous sea surface lacks the strong physical forcing inhomogeneity required to generate prominent wind stripe patterns, minor differences can still be observed at the smallest resolved scale (2Δx). While most of the central Pacific study region features homogeneous oceanic surfaces, resolvable islands (e.g., Hawaii and numerous small islands east of Australia) introduce observable heterogeneity. As demonstrated in Fig. 10 (depicting near-surface winds over the eastern Australian waters), the Exp_Ctrl exhibits clear stripe patterns in low-level winds over some of the islands and their upwind regions, such as parts of Fiji (16–19° S, 177–180° E) and Vanuatu (14–20° S, 166–171° E), whereas Exp_Test shows no such artifacts. This contrast highlights the role of surface inhomogeneity in noise generation.

3. Terrestrial spectra exhibit distinct energy distributions between Exp_Ctrl and Exp_Test across scales ranging from 2Δx to 8Δx. This spectral difference indicates that when inhomogeneous surface forcing combined with staggered-grid physics-dynamics coupling, dominant 2Δx-scale perturbations emerge consistently while statistically significant oscillations develop at discrete larger scales (e.g., 4Δx).

Figure 10 18 h forecasts of (a) u- and (b) v-component for Exp_Ctrl, and (c) u- and (d) v-component for Exp_Test, at the lowest model level from CMA-GFS 3.0.

Our sensitivity analysis confirms that the noise in CMA-GFS 3.0 originate specifically from the staggered-grid coupling between dynamic and physical processes over inhomogeneous surfaces. However, the noise-free Exp_Test configuration artificially decouples the wind field from mass fields (e.g., temperature and moisture) spatially in physical processes, violating fundamental physical constraints. We emphasize that Exp_Test serves solely as a mechanistic diagnostic to isolate the noise source rather than a viable approach for operational model improvement.

5 Conclusion and discussion

The removal of fourth-order horizontal diffusion in CMA-GFS 3.0 (Shen et al., 2023), while enhancing small-scale system simulations, has inadvertently exposed underlying 2Δx stripe patterns in low-level wind fields that were previously suppressed by the horizontal diffusion scheme. The results of a case study demonstrate that these stripe patterns exhibits one-dimensional oriented distribution characteristics, with meridional alignment in u-component and zonal alignment in v-component. These stripes maintain continuous presence throughout the 120 h period, with spatially predominant occurrence over land and near islands, while exhibiting diurnal amplitude modulation. We systematically elucidate the generation mechanism of these stripe patterns through a combination of numerical experiments and theoretical analysis.

Through controlled experiments with the dynamic core alone (no physics) and examination of associated surface static fields, we conclusively exclude the possibility that either the dynamic numerics or physical parameterizations independently generate the observed 2Δx noise. This leaves physics-dynamics coupling as the sole plausible origin, consistent with our previous findings (Chen et al., 2020) that 2Δx computational modes may emerge specifically when dynamics and physics interact through staggered-grid coupling. Two-dimensional idealized experiments successfully reproduced 2Δx waves generated at inhomogeneous surface friction points when combining advection and friction on staggered grids. Linear wave theory analysis further revealed that this dispersion effect originates from the coupling between dynamic advection and physical forcing on staggered grids under non-uniform physical forcing. These mechanistic analyses indicate that the grid mismatch in dynamic-physical coupling is the fundamental cause of noise generation over inhomogeneous surfaces.

Sensitivity experiments confirm that when dynamic-physical grid consistency is maintained and interpolation coupling is avoided, the wind field noise can be completely eliminated. Spectral analyses reveal land-ocean contrasts: Over oceans, staggered-grid coupling shows negligible impacts due to surface homogeneity. Over land, grid mismatch in physics-dynamics coupling under inhomogeneous surface elevates small-scale energy (2–8Δx), generating broadband noise, while unstaggered physics-dynamics coupling restores the energy spectrum to a reasonable decay curve with decreasing scale.

Our results demonstrate that the stripe noise in low-level wind field of CMA-GFS 3.0 stems from grid mismatch in physics-dynamics coupling under non-uniform forcing conditions. These findings offer critical insights for other NWP models employing similar grid configurations, revealing potential numerical dispersion risks in physics-dynamics coupling on traditional staggered grids.

Building on these findings, we systematically evaluate potential solutions to mitigate wind-field noise in CMA-GFS 3.0. Since surface inhomogeneity is inherently present in NWP models, one essential approach to noise reduction involves resolving the grid mismatch between physics parameterizations and dynamical core. In latitude-longitude grid configurations, the Arakawa C-grid exhibits superior dispersion characteristics compared to unstaggered A-grid arrangements, particularly in maintaining proper phase relationships for geostrophic adjustment processes. However, adopting an A-grid configuration may introduce computational challenges in the dynamical core that could degrade numerical accuracy.

It must be emphasized that grid configuration represents a fundamental architectural feature of numerical modeling systems. As the cornerstone of model design, grid setup is determined at the earliest development stage. Redesigning grids equates to building an entirely new model system, making modifications to CMA-GFS's existing Arakawa C-grid configuration unfeasible. An alternative approach would involve computing wind fields directly on the dynamics grid within physical parameterizations. However, this strategy encounters substantial technical obstacles owing to the inherently tight coupling between wind and temperature/moisture fields in physical processes, especially when managing their interactions on staggered grids.

A compromise solution involves employing higher-order horizontal diffusion or filtering which is widely implemented in NWP models. Methods like high-order horizontal diffusion can effectively eliminate 2Δx high-frequency noise while retaining most resolvable-scale energy. Specifically, as demonstrated by Xue (2000), 4th-order horizontal diffusion preserves about 80 % of 4Δx-scale spectral energy, and 6th-order diffusion achieves ∼90 % retention at this scale. However, one-dimensional linear wave analysis shows that staggered grid coupling under non-uniform forcing not only produces dispersion effects but also alters wave amplitude and phase speed, and these systematic deviations cannot be corrected by simple noise-removal tools.

Higher-order horizontal diffusion can serve as a practical remedy for unexplained noise. However, this study has definitively pinpointed the specific sources of wind-field noise in CMA-GFS. We therefore propose the following targeted recommendations:

1. While the piecewise-constant sampling method (Eqs. 14 and 15) effectively suppresses numerical noise, it may introduce directional biases. A more physically consistent approach would be upwind sampling, where:

   $$\begin{array}{}\text{(16)}& \begin{array}{rlrl}{\displaystyle }{u}_{\mathrm{p}}\left(x+{\displaystyle \frac{\mathrm{\Delta }x}{\mathrm{2}}}\right)& {\displaystyle }=u\left(x\right)& {\displaystyle }& {\displaystyle }\text{if}u\ge \mathrm{0}\\ {\displaystyle }{u}_{\mathrm{p}}\left(x+{\displaystyle \frac{\mathrm{\Delta }x}{\mathrm{2}}}\right)& {\displaystyle }=u(x+\mathrm{\Delta }x)& {\displaystyle }& {\displaystyle }\text{otherwise}\end{array}\end{array}$$

   This straightforward approach offers immediate operational feasibility and can be rapidly implemented to address the issue. We recommend trial implementation in the operational system followed by comprehensive impact assessment.

2. The strong connection between wind and heat transfer in the boundary layer turbulent diffusion makes it difficult to compute momentum diffusion directly at wind grid points. As demonstrated in Chen et al. (2020), interpolating the diffusivity (rather than prognostic variables) effectively eliminated vertical grid-scale noise in thermodynamic fields – a numerical artifact originally induced by staggered-grid coupling between dynamic and physical processes. Following this approach, we recommend averaging the diffusion coefficient back to the wind points and performing the vertical diffusion on the wind points, thereby avoiding interpolation of prognostic wind variables. Given its demonstrated efficacy in addressing similar discretization challenges in our previous studies (Chen et al., 2020), this approach merits implementation and systematic evaluation.

3. In developing next-generation model frameworks, the coordination of physics-dynamics coupling should be a core design consideration to fundamentally prevent such issues.