In this study, a vertical velocity (w) assimilation scheme is developed in the 3DVar method. In this w assimilation scheme, the Richardson equation is employed as the operator for w, enabling the 3DVar system to update dynamic and mass variables from w assimilation. The w assimilation scheme looks reasonable to me. However, I have concerns about the experimental configurations for both the single case study and batch experiments in the manuscript. Therefore, I recommend a major revision.

Major comments:

1. The authors highlight the use of the Richardson equation as the operator for w to reduce imbalance, but there is a lack of corresponding discussion and figures to support this assertion. I strongly recommend the authors incorporate related discussions and provide figures that illustrate the effectiveness of employing the Richardson equation for w in reducing imbalance.
2. In the single observation test, only horizontal increments are presented, while the w observation captures vertical motions. It would be more informative to display vertical cross-sections of the increments, revealing the vertical spread of assimilating w single observations. I am particularly interested in understanding the vertical range of assimilating w Does it extend from the bottom level to the top level? Are there any constraints limiting the impacts of w assimilation? Additionally, based on Fig. 1, the range of significant increments appears to be approximately 8 degrees. Have you considered tuning the decorrelation scales of BEC for w assimilation in real case experiments? I believe adjusting these scales could be crucial for constraining the impact radii of w observations in convective-scale data assimilation.
3. The assimilation configurations in both the heavy rainfall case study and batch experiments appear unconventional. In the case study, CNTL performs no assimilation, while DA-W assimilates w hourly for only three hours. This leads to DA-W being initialized four times additionally, potentially introducing significant impacts. It's challenging to attribute these impacts solely to the assimilation of w In the batch experiments, both cases seem to involve single-time assimilation each day, likely contributing to minimal impacts from w assimilation. I recommend the authors rerun both batch and single-case experiments, considering an increased number of assimilation cycles per day. For consistency, the single case experiments could from one day of the batch experiments.

Minor comments:

Line 42: Change “vertical velocity” to “w”. Check it throughout the manuscript.

Line 54: Change “computing” to “computational”.

Line 54: The physical constraint is implicitly considered in the ensemble BEC in the EnKF-based methods.

Line 58: the observation operator of w is not the problem. The issue is the effectiveness of assimilating w in the forecasts.

Line 86: “wind” should be “winds”.

Line 91: It is not accuracy. H links the model state variables to the observed variables, not the control variables.

Line 96-97: This sentence is not accuracy enough. The observation operator combines the dynamic and mass fields, but it cannot adjust these variables. Maybe it could be modified to “enabling the 3DVar method to adjust ….”.

Line 100: add “height” after “top”.

Line 116-119: Make it simple. Maybe the (2) and (3) can be combined.

Line 122: It is common to use bold H to represent the tangent linear observation operator and with the subscript T to represent its adjoint operator.

Eq. (9): The cost function J in Eq. (1) corresponds to analysis state x, while in Eq. (9) it is for the control variable. Modify it to make it consistent.

Fig. 3: How is the bias score calculated? Is it the frequency bias? For me, the values of bias are too large, indicating a significant overestimation of both experiments.

Fig. 4: The flowchart is a little bit confusing to me. Is it continuous cycling or partial cycling? Were the assimilations done only at 06Z every day? The description of the flowchart is not clear enough to me.

Fig. 5 and 6: The authors may enhance clarity by presenting the average forecast skills of ETS, FSS, and BIAS over a ten-day cycling period instead of displaying results for each case. By utilizing ten days' samples, the inclusion of error bars in the figures can provide a more comprehensive representation of variability and uncertainty.