Upon reviewing the three main methods of estimating vertical motion in prior papers (kinematic, thermodynamic, and QG-omega), I've decided that implementing the kinematic method might be the most straightforward method to understand and also seems to be best suited to my needs. I initially tried solving for the traditional QG-Omega equation, but I am really struggling with "inverting" the laplacian on the left-hand side of the equation, and so have now opted to approach this problem using the kinematic method. However, I'm running into some trouble working out the logic of using this method in my head.

Its derivation seems easy enough. Starting with the continuity equation in pressure coordinates: $$ \frac{dw}{dp} + \frac{du}{dx} + \frac{dv}{dy} = 0 $$

we can rearrange and integrate to yield:

$$ w_{p} = w_{p+\Delta p} + \int_{p}^{p+\Delta p} ( \frac{du}{dx} + \frac{dv}{dy})dp $$

where w is vertical velocity, p is some reference pressure, u is zonal wind, and v is meridional wind.

Now the problem I'm running into is that much of the literature recommends against using the kinematic method for acquiring a realistic estimate of vertical motion. This seems mainly due to the fact that the geostrophic wind is non-divergent, and errors in the ageostrophic wind are large enough to significantly affect horizontal divergence. This makes sense to me for instances where one is using observations, since the goal is to match reality as closely as possible, and observational readings are prone to error both due to instrumentation and error that likely arises from the interpolation process from discrete observations onto a grid.

However, for model data, would this drawback still apply? There are no "errors" since the fields are continuous and the model is simulating everything based on its own assumptions. In that sense, I'm not looking to match reality, but rather match whatever the model's reality is. It seems to me like this is how the WRF model already calculates omega (indicated within the comments of the calc_ww_cp subroutine in the WRF code, although my Fortran is a bit too rusty to actually make sense of it). Basically, I'm looking for the most efficient way of obtaining a realistic quantity representing vertical motion within the model output.

In that sense, is the integration of the continuity equation the best way of accomplishing this for numerical model output? Or is it worthwhile to pursue a different method?

Thanks in advance, and I look forward to any fruitful discussions that can come out of this!

The major drawback of the kinematic method (errors in the ageostrophic wind propagating into estimates of $\omega$) does not apply to model data where the accuracy in horizontal wind is approximated to several orders of digits. In observational data, where this order of accuracy does not exist, it's recommended to avoid the kinematic method. In some personal tests, I found that integrating divergence on isobaric levels in reanalysis data (ERA5) successfully reproduces the $\omega$ field within the reanalysis itself. Digging through the ECMWF technical notes confirms that continuity is used for their purpose of estimating $\omega$ in gridded data (specifically, they refer to Simmons & Burridge (1980) - An Energy and Angular-Momentum Conserving Vertical Finite-Difference Scheme and Hybrid Vertical Coordinates [Specifically eqn (2.5). This is the approach ECMWF uses to solve for omega on hybrid vertical levels]).

Speaking subjectively, this method seems to be the easiest to implement and the most intuitive. In instances where $\omega$ is not available for your data and the phenomena you wish to study satisfies the hydrostatic approximation (such as in most GCMs or studies examining synoptic-scale features), the kinematic approach is justifiable.