Parabolic equation (PE) methods have been widely applied to the modeling of wireless propagation in tunnel environments. However, the relevant literature does not include concrete guidelines for the choice of the parameters of these methods and the tradeoffs involved. This paper provides a comprehensive analysis of the two sources of error that arise when PE methods are employed for the modeling of radio-wave propagation scenarios: the well-known numerical dispersion error stemming from the finite-difference solvers for PE and the approximation error stemming from the use of PE for the solution of wave propagation problems that are subject to Maxwell's equations. The analysis is performed for four methods, three of which have been already used in PE-based propagation studies, namely, the Crank-Nicolson (CN) scheme, the alternative-direction-implicit (ADI) method, and its locally one-dimensional (LOD-ADI) version. The fourth method is the Mitchell-Fairweather (MF)-ADI scheme that has been recently shown to be a promising alternative technique for tunnel propagation modeling. The proposed method leads to robust criteria for the choice of spatial discretization in realistic propagation scenarios, as shown via numerical examples.