Let $A=\bigl[a_{ik}\bigr]$ be the unknown $(4\times4)$-matrix in question. Given that the $10$ described sums all have the same value $s$ we have a system ${\cal S}$ of $10$ linear equations in the $16$ unknown variables $a_{ik}$. The $(10\times16)$-matrix $M$ (all entries are $0$ or $1$) of this system has only rank $9$. In fact the four rows of $M$ corresponding to the four rows of $A$ sum to $(1,1,\ldots,1)$, as do the next four rows of $M$ corresponding to the four columns of $A$. The rank of $M$ could be even smaller because we do not really understand the influence of the two rows of $M$ corresponding to the main diagonals of $A$.
When believe in ${\rm rank}(M)=9$ then we can argue as follows: Setting up for ${\cal S}$ a Gaussian reduction process for the $16$ unknowns $a_{ik}$ will in the end present $9$ of these unknowns as functions of the $7$ others. A non resolvable case cannot result when all row and column sums have the equal value $s$.
According to the algorithm of the calculation we therefore can choose the values of $7$ suitable $a_{ik}$, and we have to make all these choices in order to determine the $16$ $a_{ik}$ completely. I solved the system with Mathematica. Here is the system:

And these are the solutions:

You see that $9$ variables $a_{ik}$ appear as functions of the seven others and $s$. Which variables are the "free" ones depends on the algorithm of the solver. Here the variables $a_{11}$, $a_{12}$, $a_{13}$, $a_{21}$, $a_{22}$, $a_{23}$, $a_{31}$ are "free".