## 3 Numerical Discrete Element Method (DEM)  
The numerical simulations have been conducted on the Discrete Element environment offered by 3 DEC [37]. This section summarises the general features of 3 DEC related to the specificities of the present numerical simulations. The considered specimens are described using the so-called simplified micro-modelling strategy. Each block is modelled in 3D by a separate entity, while zero-thickness interfaces model the dry joints. In 3 DEC, blocks can be either deformable or rigid. In the present paper, given the low stresses acting on them, blocks are considered fully rigid and characterised by only six degrees of freedom and a unique mechanical parameter: their density (\(\rho=1800\) kg/m\({}^{3}\)). Furthermore, a single numerical simulation using deformable blocks with very low elastic rigidity (E = 1 \(\times\) 10\({}^{8}\) Pa) confirmed the negligible influence of the block deformation itself. At the interface between blocks (either for bed or head joints), contact forces are modelled through equivalent springs (Figure 2). One normal spring relates the elastic block interpenetration in the normal direction of the joint (u\({}_{\text{n}}\)) with the normal stress (\(\sigma\)). Similarly, two tangential springs relate the shear elastic relative displacement (u\({}_{\text{s}}\)) with the shear stresses (\(\tau\)) in the two tangential directions. Equation (1) describes the link between the increments in stress and displacement, in the process of moving from step t to step t + \(\Delta\)t. Only one shear equation is shown for brevity, instead of the two implemented in 3 DEC.  
\[ \begin{split}\Delta\sigma=\text{k}_{\text{n}}\cdot\Delta\text{u }_{\text{n}}&=\text{k}_{\text{n}}\cdot[\text{u}_{\text{n}}(\text {t}+\Delta\text{t})-\text{u}_{\text{n}}(\text{t})]\\ \Delta\tau=\text{k}_{\text{s}}\cdot\Delta\text{u}_{\text{s}}& =\text{k}\text{s}\cdot[\text{u}_{\text{s}}(\text{t}+\Delta\text{t})-\text{u}_ {\text{s}}(\text{t})].\end{split}  \tag{1}\]  
Here, k\({}_{\text{n}}\) and k\({}_{\text{s}}\) are the joint's normal and tangential elastic stiffness and are sufficiently large to avoid unrealistic interpenetration of blocks. The experimental characterisation of joint stiffness parameters is not provided in [22]. Therefore, their values are varied through a parametric analysis presented in the next section. When two blocks are automatically determined to be in contact, a few contact points are generated, and the contact forces are transmitted discretely at these locations. By default, 3 DEC creates contact points at each vertex of blocks, but it has been found that more contact points were required to model the  
Figure 2: WD 6 after vertical gravity stabilisation (k\({}_{\text{n}}=\) k\({}_{\text{s}}=\) 3 \(\times\) 10\({}^{7}\)). Boundary conditions and reference axis are shown. The monitored point is also shown (see Section 5). actual stress distribution accurately [7]. This issue is further discussed in the next section concerning the modelled experiments. The Discrete Element algorithm uses an explicit algorithm to solve the equations of motion. At a given step t, all forces (including contact, body and other applied forces) and position of blocks are known. For each block, the equations of motion are integrated to compute the new acceleration of blocks (Equation (2)). Afterwards, a central difference algorithm calculates the block velocities, which enable finding the new block positions at step t + \(\Delta\)t.  
\[ \begin{split}\text{m}_{\text{i}}\text{a}_{\text{i}}(\text{t})= \Sigma\text{F}_{\text{i}}(\text{t})\\ \text{v}_{\text{i}}(\text{t}+\Delta\text{t}/2)=\text{v}_{\text{i} }(\text{t}-\Delta\text{t}/2)+\text{a}_{\text{i}}(\text{t})\cdot\Delta\text{t} \\ \text{u}_{\text{i}}(\text{t}+\Delta\text{t})=\text{u}_{\text{i}}( \text{t})+\text{v}_{\text{i}}(\text{t}+\Delta\text{t}/2)\cdot\Delta\text{t}. \end{split}