Optimal vertical coordinate for complex topography
This is a big question, but there are benefits and difficulties with all representations, and the efficacy of any coordinate will also depend upon other numerical choices such as the calculation of horizontal pressure gradients [Zängl 2012].  There are broadly two categories: terrain following layers and Cartesian systems such as step terrain and cut cell grids.

Terrain following layers are typically implemented using a coordinate transformation such that the computational domain is rectangular.  There are essentially two choices to make here.  First, the vertical coordinate may be geometric height, pressure, or entropy (potential temperature), or some other monotonic function.  Second, the decay function controls how the influence of the terrain on computational surfaces decays with height: this might be a linear decay (called Basic Terrain Following (BTF) or 'sigma' coordinates [Gal-Chen &amp; Somerville 1975]), or something more elaborate such as Smooth LEvel VErtical (SLEVE) [Schär et al 2002] or Smoothed Terrain Following (STF) [Klemp 2011].

One of the motivations for smoothing coordinate surfaces is to make horizontal pressure gradient calculations more accurate: in areas of steep terrain metric terms in the coordinate transformation can increase errors [Dempsey &amp; Davis 1998].  This is also one of the motivations for step terrain and cut cell grids, since the grid is orthogonal everywhere except in the cells next to the ground.

Another issue with representing terrain is that cells can become very thin (over steep slopes with terrain following layers) or very small (in a cut cell or partial-step grid) which constrains the timestep for explicit numeric methods due to the CFL criteria [Klein 2009].  There are several methods that can overcome this 'small cell problem' [Jebens et al 2011, Yamazaki &amp; Satomura 2010, Steppeler et al. 2002] 

This is by no means an exhaustive answer but I hope it gives you some understanding of the issues involved with vertical coordinates.  Additionally, I gave a presentation in 2015 on this topic, the slides are available online.  Another useful review article on nonhydrostatic models, including terrain representation and vertical coordinates, is presented in Steppeler et al 2003.

References


Dempsey, D. and C. Davis, 1998: Error analyses and test of pressure gradient force schemes in nonhydrostatic, mesoscale model. 12th Conf. on Numerical Weather Prediction, Phoenix, AZ. Amer. Meteor. Soc., 236–239
Gal-Chen, T. and R. C. Somerville, 1975: On the Use of a Coordinate Transformation for the Solution of the Navier-Stokes Equations. J. Comp. Phys., 17, 209–228
Jebens, S., O. Knoth and R. Weiner, 2011: Partially implicit peer methods for the compressible Euler equations. J. Comp. Phys., 230, 4955–4974.
Schär, C., D. Leuenberger, O. Fuhrer, D. Lüthi and C. Girard, 2002: A new terrain-following vertical coordinate formulation for atmospheric prediction models. Mon. Wea. Rev., 130, 2459–2480
Steppeler, J., H.-W. Bitzer, M. Minotte and L. Bonaventura, 2002: Nonhydrostatic Atmospheric Modeling using a z-Coordinate Representation. Mon. Wea. Rev., 130, 2143–2149
Steppeler, J., R. Hess, U. Schättler and L. Bonaventura, 2003: Review of numerical methods for nonhydrostatic weather prediction models. Meteorology and Atmospheric Physics, 82, 287–301
Yamazaki, H. and T. Satomura, 2010: Nonhydrostatic Atmospheric Modeling Using a Combined Cartesian Grid. Mon. Wea. Rev., 138, 3932–3945
Zängl, G., 2012: Extending the numerical stability limit of terrain-following coordinate models over steep slopes. Mon. Wea. Rev., 140, 3722–3733