The equations for a vertical cut in cohesive material are given in several publications including Heyman (1973): The stability of a vertical cut (see attachment Equations1.jpeg) and Verruijt (2001) Soil mechanics (see attachment Equations2.jpeg). I would like to try to replicate this analytical example as closely as possible using FEM. I am using Plaxis2D and the results for vertical and horizontal stress are given in Result1.jpeg and Result2.jpeg respectively. The material defined is Mohr-Coulomb with Undrained C material type. The parameters are as follows:

gamma_sat = 18 kN/m3
E_undr = 5000 kN/m2
nu_undr = 0.495
su (undrained shear strength) = 1000 kN/m2

The geometry is 20 meters high (total) and 40 meters long (total).

The vertical stresses (Result1.jpeg) look promising, since at the bottom of the left part (zone 3 according to Equations2.jpeg) we have a value nearing 180 kN/m3 and at the right part (zone 2 according to Equations2.jpeg) we have a value nearing 360 kN/m3. The horizontal stresses (Result1.jpeg) however look completely wrong. I was expecting the right top part (zone 1 according to Equations2.jpeg) to show a value of 0 kN/m3 throughout and the bottom part (zones 2 and 3 according to Equations 2.jpeg) to show a minimum of 180 kN/m3 throughout. The bottom right part (zone 2 according to Equations2.jpeg) shows decent values which correspond to the minimum result of 360 kN/m3 but the bottom left part (zone 3 according to Equations2.jpeg) does not.

Why are the horizontal stresses deviating from the analytical result? How can I replicate the distributions predicted by the publications on an FEM program?