When dealing with angular momentum of a combined system you have two possible basis. Let $\mathscr{H}_1$ be the Hilbert space of particle 1 and $\mathscr{H}_2$ be the Hilbert space of particle 2.
On $\mathscr{H}_1$ you have the basis $|j_1 m_1\rangle$ where $j_1 = 1$ and $m_1 = -1,0,1$. On $\mathscr{H}_2$ you have the basis $|j_2 m_2\rangle$ where $j_2 = 1$ and $m_2 = -1,0,1$. In that sense both $\mathscr{H}_1$ and $\mathscr{H}_2$ are three-dimensional Hilbert spaces.
On the other hand, the combined system is $\mathscr{H}_1\otimes \mathscr{H}_2$. One obvious basis for this space is that associated to the complete set of commuting observables $J_1^2,J_2^2,J_{1z}J_{2z}$. This is the basis
$$|j_1,j_2;m_1,m_2\rangle=|j_1,m_1\rangle\otimes |j_2,m_2\rangle.$$
It is clear that $j_1,j_2=1$ and $m_1,m_2=-1,0,1$. You thus have $9$ basis states. But all basis of a finite dimensional Hilbert space have the same number of element which is its dimension, whatever basis set you use, it will have $9$ states.
For completeness, another natural basis is that of the total angular momentum. You define
$$\mathbf{J}=\mathbf{J}_1\otimes \mathbf{1}+\mathbf{1}\otimes \mathbf{J}_2$$
to get operators $J_z$ and $J^2$. They commute with $J_1^2$ and $J_2^2$ so that you get a basis
$$|j_1,j_2;j,m\rangle$$
This is the basis of total angular momentum. It is a result then, that you can consult for example in Cohen's book Volume 2 on the "Addition of Angular Momentum" chapter, that the possible values for $j$ the eigenvalues of $J^2$ are
$$j=|j_1+j_2|,|j_1+j_2-1|,\dots, |j_1-j_2|$$
Here $j_1=j_2=1$ hence the possible values for $j$ are $$j=2,1,0.$$
Now for $j = 0$ you have just $m = 0$ (one state), for $j = 1$ you have $m = -1,0,1$ (three states) and for $j = 2$ you have $ m = -2,-1,0,1,2$ (five states). This gives a total of $9$ states on the basis of total angular momentum as anticipated.
If
$$
\vert j_1m_1;j_2 m_2\rangle =\sum_{J(M)} C^{JM}_{j_1m_1;j_2m_2} \vert JM\rangle \tag{1}
$$
where the Clebsch-Gordan coefficients are the inner products
$$
C^{JM}_{j_1m_1;j_2m_2} =\langle JM\vert j_1m_1;j_2m_2\rangle =
\langle j_1m_1;j_2m_2\vert JM\rangle
$$
since the CG's are real. Next, start with $\vert JM\rangle$ and expand it in the uncoupled basis using the identity operator
$$
\hat{\mathbb{1}}=\sum_{j_1m_1;j_2m_2}\vert j_1m_1;j_2m_2\rangle
\langle j_1m_1;j_2m_2 \vert
$$
to get
\begin{align}
\vert JM\rangle &= \sum_{j_1m_1;j_2m_2}\vert j_1m_1;j_2m_2\rangle
\langle j_1m_1;j_2m_2 \vert JM\rangle \\
&=\sum_{j_1m_1;j_2m_2}\vert j_1m_1;j_2m_2\rangle
C^{JM}_{j_1m_1;j_2m_2}\, . \tag{2}
\end{align}
Thus, the CGs allow you to go in both direction: from the uncoupled to the coupled basis as in (1), and from the coupled to the uncoupled basis as in (2).
Best Answer
The dimension of the Hilbert space in the stated setup is four, and thus all matrix representations of the operators should be of dimension four. This is independent of basis representation.
To wit: let's start with the spin-addition basis $|J, M\rangle$ where we have the $4$ states $|J=1, M=1\rangle$, $|J=1, M=0\rangle$, $|J=1, M=-1\rangle$ and $|J=0, M=0\rangle$, and we associate them with the basis column vectors in that order (that is $|J=1, M=1\rangle = (1, 0, 0, 0)^T$ etc.). Then $$\vec{S}^2 = \hbar^2\begin{pmatrix} 2 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 \\ 0 & 0 & 2 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix} $$ and $$S_z^2 = \hbar^2\begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix} $$
Now we can work in the separated basis with the four basis states $|m_1=\hbar/2, m_2=\hbar/2\rangle$, $|m_1=\hbar/2, m_2=-\hbar/2\rangle$, $|m_1=-\hbar/2, m_2=\hbar/2\rangle$, $|m_1=\hbar/2, m_2=-\hbar/2\rangle$ and we associate them again with the basis column vectors in that order. The matrix representation of $S_z^2$ is still diagonal $$S_z^2 = \hbar^2\begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} $$ but we get non-diagonal terms for the representation of $\vec{S}^2$. However calculating $\langle m_1, m_2 | \vec{S}^2 | m_1', m_2' \rangle$ is quite straight-forward, using CG. We have two states that can only have non-zero matrix elements with themselves (as they are also states in the momentum-added basis): $$ \langle m_1=\uparrow, m_2=\uparrow | \vec{S}^2 | m_1', m_2' \rangle = 2\hbar^2\delta_{m_1', \uparrow}\delta_{m_2',\uparrow}$$ $$ \langle m_1=\downarrow, m_2=\downarrow | \vec{S}^2 | m_1', m_2' \rangle = 2\hbar^2\delta_{m_1', \downarrow}\delta_{m_2',\downarrow}$$ and have a non-diagonal matrix elements when we look at the subspace associated with total $S_z = 0$: $$ \langle m_1=\uparrow, m_2=\downarrow | \vec{S}^2 | m_1=\uparrow, m_2=\downarrow \rangle = \langle m_1=\downarrow, m_2=\uparrow | \vec{S}^2 | m_1=\downarrow, m_2=\uparrow \rangle = \hbar^2$$ $$ \langle m_1=\hbar/2, m_2=\downarrow | \vec{S}^2 | m_1=\downarrow, m_2=\uparrow \rangle = \hbar^2$$ so we end up getting $$\vec{S}^2 = \hbar^2\begin{pmatrix} 2 & 0 & 0 & 0 \\ 0 & 1 & 1 & 0 \\ 0 & 1 & 1 & 0 \\ 0 & 0 & 0 & 2 \end{pmatrix} $$
You can now convince yourself that the two representations are related to one another by basis transformations, which is given exactly by the CG coefficients (which is not surprising - this is how we got the matrices in the first place).