There is considerable potential for data-driven modelling to describe pathdependent soil response. However, the complexity of soil behaviour imposes significant challenges on the training efficiency and the ability to generalise. This study proposes a novel physics-constrained hierarchical (PCH) training strategy to deal with existing challenges in using data-driven models to capture soil behaviour. Different from previous strategies, the proposed hierarchical training involves 'low-level' and 'high-level' neural networks, and linear regression, in which the loss function of the neural network is constructed using physical laws to constrain the solution domain. Feedforward and long short-term memory (LSTM) neural networks are adopted as baseline algorithms to further enhance the present method. The data-driven model is then trained on random strain loading paths and subsequently integrated within a custom finite element (FE) analysis to solve boundary value problems by way of validation. The results indicate that the proposed PCH-LSTM approach improves prediction accuracy, requires much less training data and has a lower performance sensitivity to the exact network architecture compared to traditional LSTM. When coupled with FE analysis, the PCH-LSTM model is also shown to be a reliable means of modelling soil behaviour under loading-unloading-reloading and proportional strain loading paths. In addition, strain localisation and instability failure mechanisms can be accurately identified. PCH-LSTM modelling overcomes the need for adhoc network architectures thereby facilitating fast and robust model development to capture complex soil behaviours in engineering practice with less experimental and computational cost.