ABSTRACT With increasing modelling capabilities in geomechanics incorporating non-linear material behaviour, an increasing complexity of models representing excavation of tunnels can be analysed. Often, analysis is using finite elements, where problems size can be a controlling factor due to limited computational resources. An approach to analysing these problems is to optimise the discretisation of the differential equations by concentrating elements where solution accuracy counts the most. This paper describes a priori p-refinement method for mesh improvement for stress analysis of prismatic tunnels. The improvement results in a mesh with higher-order elements near regions of interest and lower-order elements elsewhere. The focus is on the automated insertion of transitional elements at the interface of the two regions, based on the zone of disturbance in the geomaterial. The method was incorporated into a finite element code, and its performance was tested on practical problems, where the global stiffness matrix size was reduced by 70%, and the calculation times were significantly shortened. The average difference compared to complex models without mesh improvement was less than 3%. It is anticipated that the presented method enables analysis of problems that would otherwise be beyond the current computing capabilities leading to safer and economical designs.