<p>The stress integration of critical soil model is usually based on implicit Euler algorithm, where the stress predictor is corrected by employing a return mapping algorithm. In the case of a large load step, the solution of the local nonlinear system to compute the plastic multiplier may not be attained. To overcome this problem, a sub-stepping scheme is generally used to improve the convergence of the local nonlinear system solution strategy. Nevertheless, the complexity of the tangent operator of the sub-stepping scheme is high. This complicates the use of Newton–Raphson algorithm to obtain global quadratic convergence. In this paper, a formulation for consistent tangent operator is developed for implicit sub-stepping integration for the modified Cam-Clay model and unified Clay and Sand model. This formulation is highly efficient and can be used with problem involving arbitrary large load step, such as tunnel simulation.</p>