penetration is positive, Hertzian contact theory can be used to determine the resulting force $f_i$ and contact radius $a_i$ as: $$f_i = c_i\langle g_i\rangle^{3/2}, \quad a_i = \sqrt{r_i  \langle g_i\rangle}, \quad\mbox{ with } \quad c_i = \frac43 E^* \sqrt{r_i}   \label{eq:force}$$ where $\langle g_i\rangle = \max\{0,g_i\}$ denotes Macaulay brackets, therefore the forces are induced only if the penetration is positive. We assume that the half-space’s surface deforms under action of contact forces following Boussinesq solution, thus the vertical displacement due to a single concentrated force $f$ located at $x_f,y_f$ is given by: $$\Delta z_f = \frac{f}{\pi E^*\sqrt{(x-x_f)^2+(y-y_f)^2}}.
  \label{eq:dz}$$ This coupling between additional displacement and forces makes the system of equations - strongly coupled. Then, following the superposition principle, at a given asperity $i$, its total vertical shift is given by: $$\Delta z_i = \frac{1}{\pi E^*}\sum\limits_{j \in \mathcal I\setminus i} \frac{f_j}{d_{ij}} =  {\boldsymbol{B}} {\boldsymbol{F}},$$ where $d_{ij} = \sqrt{(x_i-x_j)^2+(y_i-y_j)^2}$ and the symmetric interaction matrix ${\boldsymbol{B}}$ (containing the inverse distances) and the vector of forces ${\boldsymbol{F}}$ are given by: $${\boldsymbol{B}} = 
  \frac{1}{\pi E^*}\left[\begin{array}{ccccc}\displaystyle        0 & d^{-1}_{12} & d_{13}^{-1} & \quad\dots & d_{1N}^{-1}\\[3pt]
        d^{-1}_{12} & 0 & d_{23}^{-1} & \quad\dots & d_{2N}^{-1}\\[3pt]
        d^{-1}_{13} &  d_{23}^{-1} & 0 & \quad\dots & d_{2N}^{-1}\\[3pt]
        \vdots & \vdots  & 