This manuscript presents and compares several methods to smooth terrain data to be used in weather models with terrain-following coordinates. The contents of the paper are useful for modeling and the paper's structure is coherent. However, the presentation is not always clear. The manuscript therefore requires minor revisions before it can be published. The issues that need to be addressed are:

line 21: the 35 degree as threshold value needs some more context, since it is heavily relied upon throughout the paper. What mechanism causes the instability at 35 degrees, and why can some models tolerate steeper gradients? 

line 32: what is a 'closest approach'?

line 38: 'thorough' - I do not disagree, but adjectives like this do not belong in scientific text.

line 71: The grid spacing, probably not the resolution, increases with height.

Section 2.1.2: I find the inner and extended FastEddy domains a bit confusing: it seems that only results from the extended domain are presented, so what's the use of the inner domain in this paper? Also, the naming implies that these are two (nested) simulations (like for WRF), which is not the case, if I understand it correctly.

Figure 1 caption: Typo (crush instead of crash).

line 139: How exactly is the terrain upscaled? Since the paper is about terrain data processing, this is an important detail.

line 149: it should be more explicitly stated what convergence entails.

Figure 6: What is the purpose of the metric relative elevation differences? A certain elevation difference caused by the smoothing will have a relative elevation difference that depends on the location in the domain. Why not present absolute elevation difference or slope difference?

Figure 7:  In the slope density distribution, it is not clear what each color means.

Figure 7:  All three resolutions have the same maximum wavenumber: 2 10^-2 m^-1, which corresponds to 50 m. Since the three panels use different resolutions, it is not clear what is meant with the wavenumber. Also, the y axes lack units.

line 265: 'the developed local smoothing method ensures numeral stability', this conclusion is a bit too strong (with one case study).

Summary

This manuscript addresses a practical challenge in high-resolution atmospheric modeling: numerical instabilities arising from steep terrain slopes in terrain-following coordinate systems. The authors develop and evaluate a localized terrain smoothing approach that applies a Gaussian filter with adaptive standard deviation over a 3×3 stencil, using a blending factor to preserve original terrain at the edges of the smoothing window.  Fifteen different smoothing configurations across various parameters are compared for convergence speed, computational cost, and terrain distortion. The selected method (simultaneous 3×3, BF = 0.2) is implemented in both WRF(v4.6.0) and FastEddy (integrated into v3.0), and demonstrated on a real-terrain case study over the city of Murcia, Spain, where simulations are able to complete successfully after applying the local smoothing.

Overall, the study focuses on creating a practical terrain filtering approach designed to minimize the number of topography grid points that are adjusted. Major questions/concerns that should be addressed:

1. The Nyquist cutoff limits information represented on a grid to 2*dx. If topography data from a higher resolution dataset is sampled onto a dx resolution grid, there will be aliasing of higher wavenumbers onto the grid. These can create noise that is generally removed by applying terrain smoothing. If the local filtering approach is applied, how is aliasing addressed across the rest of the domain?
2. Explain why global smoothing is a bad thing, given that aliasing is an issue when sampling from a higher resolution dataset. There could be an argument that there are many uncertainties in the topography and land cover datasets, but it seems odd to adjust the grid points only in the steep terrain sections of the domain.
3. There is no evaluation of the meteorological forecast differences with global vs local smoothing. Does the local smoothing actually give better simulation results? Why bother with a more expensive local method if the global method works quickly, and the raw topo data needs filtering anyway to remove aliasing? Can the benefits be highlighted more?

Additional specific comments are given below. Major revisions are recommended. 

Specific comments

1. Lack of meteorological validation of the smoothed-terrain simulations: the paper demonstrates the locally smoothed terrain allows previously failing simulations to complete without CFL errors with real cases. With that said the paper shows no comparison of any comparison with observation. Even a qualitative comparison (10m-Wind, T2, vertical wind components…) would substantially strengthen the paper. 
2. Line 32 - "introduced a closest approach to localized" Since the idea of localized smoothing has not yet been introduced, perhaps rephrase this.

1. Lines 35-45 - The last paragraph with the 'table of contents' could be merged with the previous paragraph to make it seem less repetitive.
2. Mention which namelist parameters specify the typical terrain smoothing in WRF and FastEddy.
3. Line 55-56: The vertical grid mentions 45 levels with increased resolution near the surface with no further detail. What is the lowest model level height? Further details on vertical levels seem to be relevant information because the slope angle that causes instability depends on the grid aspect ratio. The stretching factor and other grid parameters should be given.
4. Line 68-70: Same applies with FastEddy domain using 80 vertical levels with stretching to achieve ~10m resolution near the surface. What’s the near-surface vertical spacing exactly? Any detailed information needed to assess the effective resolution and the relevance of the 35° slope threshold.
5. Figure 1 caption - 'model crush' → crash
6. Line 69 - 2700 m is quite a low model top for complex terrain!
7. Line 72 - explain why you need an extended domain?
8. Line 75 - 'once resampled' --> the 2m resolution topography should be smoothed to 10 m -- it shouldn't be resampled because that retains wavelengths smaller than the resolution of the grid.
9. Line 98 - "the use of a localized method ensures that the majority of the grid points are not being modified in any way, as only the steep-slope points (and their immediate surroundings) are smoothed out" - again this means you have aliasing noise. 
10. Line 101 - "the maximum number of iterations is set to 1.5 times the number of grid points with steep slopes." Why? 
11. Line 139 - Why not just start with the 25 m resolution domain and put that in Figure 1 instead? It seems unnecessary to add it here and add an extra figure (Fig 3) for this.
12. Line 166 - "The global smoothing method is by far the fastest, followed by the simultaneous methods" and Line 170 - "Based on computational time alone, any of the simultaneous methods would be suitable candidates for selection." But 859 seconds is pretty long if you can do it in 0.39 s.
13. Can you prove that local smoothing gives you better simulation results? Why bother with a more expensive local method if the global method works quickly, and the raw topo data needs filtering anyway to remove aliasing?
14. Figure 4 - Why does the max slope increase with # of iterations? (bottom right panel)
15. Figure 7 - Explain why there is an increase in light blue near the cutoff.
16. Line 190 "increase in power is observed for wavenumbers 190 greater than 10−2 m−1, likely reflecting small-scale noise introduced during the blending process" - this extra noise could be introducing error into the simulation as well. Why is this noise ok to include but the error from global smoothing is not ok?
17. Figure 8 - this is presented as though the original terrain is the correct one. But this is higher resolution data sampled on the coarser domain and thus contains incorrect frequencies/values.
18. Line 199 - "The subsequent steps of the WPS workflow remain unchanged" - doesn't this mean that WPS applies its own smoothing on top of this?
19. Line 216 - "(epssm ≥ 0.9)." This is really high! There are other ways to make the simulation run, such as coarsening dz, because the ratio of dz/dx and terrain slope is what affects stability. 
20. Line 217 - "may result in unphysical behavior" Why unphysical? It's just reducing the order of accuracy of the time advancement for the vertical implicit scheme from 2nd to 1st order?
21. Table B1 and B2  - this information doesn't seem relevant to this paper unless you add actual simulation results to this paper. (Typo: Table B2: Deciduoud broadleaf forest -> Deciduous broadleaf forest)