A Review of  “Three-dimensional variational assimilation with a multivariate background error covariance for the Model for Prediction Across Scales – Atmosphere with the Joint Effort of Data Assimilation Integration (JEDI-MPAS 2.0.0-beta)” submitted to Geoscientific Model Development by Jung et al. (2023)

Review Decision: Major Revisions

Manuscript type:  Development and technical papers

General Comments:

This manuscript is a part of a series of JEDI-MPAS development papers that include the first one introducing the first JEDI-MPAS implementation by Liu et al. (2022) and the second one demonstrating the EDA of 3DEnVar with JEDI-MPAS by Guerrette et al. (2023, under review). Compared to the two companion papers mentioned above, this manuscript focuses on documenting 3DVar and the associated static background error covariance (BEC) of JEDI-MPAS using BUMP. While the paper is well-written, easy to follow, and would be useful to the NWP community interested in using JEDI-MPAS, my biggest concern is regarding the experimental design and the presentation of the results from using the static BEC generated for JEDI-MPAS. Comparing 3DVar with Pure 3DEnVar and hybrid 3DEnVar, all of which used the same static BEC generated via BUMP, really does not say much about the novelty of using BUMP. It is more of an assurance that the static BEC produced using BUMP works as expected. I suspect if one produces static BEC following Wu et al. (2002) and uses recursive filters, they are likely to reach similar conclusions that hybrid 3DEnVar > Pure 3DEnVar > 3DVar. That said, I think it is more important to show the novelty of BUMP compared to other approaches in generating static BEC for JEDI-MPAS (this is also suggested by the authors in their future work), given that this study is the first evaluation of BUMP for use in atmospheric DA. It is also important to demonstrate the efficiency, usefulness, and flexibility of BUMP. Overall, I think this manuscript is of great importance to both research and operational NWP communities and its topic also aligns with the aims and scope of GMD, however, some work, which is not minor, is needed to reformulate the presentation of the manuscript. With that, I recommend major revisions with a few comments and questions listed below.

1. Page 1, Line 14: “hybrid covariance that combines that used …” should be “hybrid covariance, which combines that used…”
2. Page 1, Line 19: missing “centers” after numerical weather prediction
3. Page 2, Lines 44-46: What prevents the products of vectors with the univariate spatial correlation matrices from being computed on the native mesh? It is not clear here as to why it is required to compute in a thinned subset of mesh and interpolated to back to full-resolution mesh. Is this related to BUMP NICAS mentioned later?
4. Page3, Line 56: “heigh-based” should be “height-based”
5. Page 3, Line 67: “the United Forward Operator” should be “the Unified Forward Operator”
6. Page 3, Line 79: “distinguishes” should be “distinguished”
7. In section 2.3, It is a little difficult to understand the new approach (relative to the previous one used in Liu et al. 2022) as described in the second and third paragraphs. The second paragraph talks about analysis variable change, not increments, however, the third paragraph is all about increments even though the variables being described in the two paragraphs are the same. The most unclear part is from line 92 to line 95. How do p, rho_d, and theta_d change relative to background forecast even if no observations were assimilated. Why would background forecast change in this case?  And how does this new approach contribute to temperature bias reduction in the stratosphere? Assuming the stratospheric temperature bias was an issue in Liu et al. (2022), please provide more context and consider re-write these two paragraphs for more consistency.
8. Page 4, Lines 100-101: what is the purpose of mentioning “without a halo region” here? Was it meant to indicate that there are no operations of derivatives or interpolations required for the state(x) and increment (delta x) in JEDI-MPAS as such no halo region is needed? Please explain.
9. Page 4, Line 106: “of” is missing after “independent”
10. Page 4, Lines 112-113: Please explain in more details on how BUMP models the univariate spatial correlations differently from the recursive filters, which was used in GSI? Perhaps try to make efforts to connect this paragraph with those on the last paragraph (Lines 165-170) of Page 6. 
11. Section 3.2: Another concern I have is regarding the parameters of the multivariate background error covariance being diagnosed from NCEP GFS forecasts, which is a different model from MPAS. The authors mentioned in their future work that they plan to use an ensemble from JEDI-MPAS to train the covariance model, but why not use the MPAS forecasts to diagnose these parameters in the first place? Isn’t the MAPS forecast a natural and more straightforward choice that can realistically represent the multivariate structure of the forecast errors of MPAS compared to NCEP GFS? In addition, knowing that the covariance is meant to be used for 6-hourly cycling DA (which is typical), why would the author follow the traditional approach and train the covariance model using data from 24 hour forecast differences and then apply a re-scaling factor to address the gap? This begs the question whether the re-scaling could be avoided if the covariance was trained from 6-h differences of MPAS forecasts. In addition, what is the spatial resolution of the NCEP GFS samples?
12. Page 8, Lines 211-219: Does this modification suggest that the assumed correlation function (the fifth-order, compactly supported function from Gaspari and Cohn 1999) is not optimal for stream function and velocity potential, but okay for other variables such as temperature, specific humidity and surface pressure?
13. The authors didn’t show analysis increments of specific humidity (q) from the two single observation tests. Is it because there is no correlation between q and other analysis variables (meaning q is univariate) so there is zero q increment? In addition, there is no single observation test for observation of q. Is it due to the same reason?
14. Page 9, Lines 253-257: As stated here, a one-month “cycling” experiments were performed, but at “each” cycle, a 20-member ensemble of 6-hour MPAS forecast was performed using IC from GEFS. Are the experiments fully cycled or partial cycled? The latter one means each cycle always cold-start from GEFS. I would think the experiments are fully cycled, but it is not clear to me which approach was really taken by the authors.
15. As stated in the conclusions, the formulation of the JEDI-MPAS static B generally follows Wu et al. (2002), but with the novel use of BUMP for multiple elements of the covariance model. Although the BUMP package including VBAL and VAR are introduced to perform variable transforms, these transformations generally follow Wu et al. (2002). The most novel parts of BUMP are perhaps the NICAS and HDIAG drivers that are used to model the univariate correlation in place of recursive filters. The NICAS and HDIAG drivers allow one to specify spatial correlation functions and compute convolutions on (semi-) native mesh, which should be considered an enhancement over recursive filters used in GSI. However, the subsequent modification which halves the diagnosed horizontal correlation length for stream function and velocity potential makes this enhancement less promising. If the specified correlation function (i.e., the fifth-order compactly supported function from GC1999 chosen in this study) leads to results that are far from the sample statistics, then why not make efforts to find a more appropriate correlation function that may be more suitable for certain variables? It gives the impression that the novelty that BUMP brought about was not fully exploited.
16. Figure 1 Caption: regression coefficients of (a) delta T, (b) delta phi, and (c) delta kai are the nonzero elements at this mesh cell of the submatrices should correspond to M, L, and N, respectively…, not L, M and N, according to Equation (4).
17. Figure 6.: “using the length scale the gives the best fit” should be “using the length scale that gives the best fit”
18. Figure 7: Please draw the reduced horizontal correlation length in a different color to show the effect of the additional modifications made to the raw statistics.
19. Figures 7-9: Please include variable units.

Review of Three-dimensional variational assimilation with a multivariate background error covariance for the Model for Prediction Across Scales-Atmosphere with the Joint Effort for data Assimilation Integration (JEDI-MPAS 2.0.0-beta) by Jung et al.

This manuscript describes the 3DVar system within JEDI-MPAS and the developments for the static multivariant covariances. The developed static covariances include variable changes for winds from sf/vp to u/v, balance operators representing wind-mass balance, and error standard deviations. Examinations on single-observation tests provide similar results to previous global 3DVar systems. Cycling DA experiments over a month-long period show that hybrid with the benefits of both static and ensemble covariances yields better performance than 3DVar and pure EnVar. I think this manuscript introduces a part of JEDI developments and can contribute much to the research community. However, the authors did not provide details of the novel developments. This manuscript can be perfected if those details can be added. I suggest a major revision before accepting this manuscript for publication.

1. Line 56: “heigh” -> “height”
2. Line 67: “United” -> “Unified”
3. Line 132-133: Can you clarify the purpose of using the same level only for calculating regression coefficients for δψ and δχ? Do you think their vertical cross-correlations are weak/negligible?
4. Line 138: I think this manuscript can contribute much to the research community. If the authors can include how BUMP implements these operators from the algorithm perspective in detail, this manuscript can be at a higher level.
5. Line 147: You directly used GFS forecasts to calculate the static error statistics. Are GFS forecasts appropriate to be used to represent MPAS model errors? I doubt it.
6. Line 152: Using NCL first seems to make the procedure complicated. Is it an essential step, or the alternative strategies exist?
7. Line 166-171: This is one novelty part relative to the other utilities (e.g., gen_be_v2 in GSI), right? I suggest the authors give more details about HDIAG.
8. Figure 8: Is it the same observation location used as Fig.7, but for the zonal wind? Can you clearly state the location of the single zonal wind observation and mark it in Fig. 8?
9. Figure 9: Same as the comment for Fig. 8. Please mark the location of the assimilated observations.