One must remark that derivatives in Sobolev spaces are usually taken in the sense of distributions: given $k\in\mathbb{N}_0=\{0,1,2,\ldots\}$, $H^k(\mathbb{R}^n)$ is the space of tempered distributions $u$ on $\mathbb{R}^n$ such that all partial derivatives of order $\leq k$ belong to $L^2(\mathbb{R}^n)$, where these derivatives are in the distribution (i.e. weak) sense: for all $f\in\mathscr{S}(\mathbb{R}^n)$ we have that $$u(f)=\int_{\mathbb{R}^n}\overline{u(x)}f(x)dx$$ and $$u_\alpha(f)=(-1)^{|\alpha|}u(\partial^\alpha f)=(-1)^{|\alpha|}\int_{\mathbb{R}^n}\overline{u(x)}\partial^\alpha f(x)dx=\int_{\mathbb{R}^n}\overline{u_\alpha(x)}f(x)dx$$ define $u(x),u_\alpha(x)\in L^2(\mathbb{R}^n)$ for each $$\begin{split}\alpha&=(\alpha_0,\ldots,\alpha_{n-1})\in\mathbb{N}_0^n\ ,\,\partial^\alpha=\partial_0^{\alpha_0}\cdots\partial_{n-1}^{\alpha_{n-1}}\ ,\,\\1&\leq|\alpha|=\|\alpha\|_{l^1}=\alpha_0+\cdots+\alpha_{n-1}\leq k\ .\end{split}$$ Just think of the above formulae as in integration by parts. This means that, as a rule, solutions $u$ of the Klein-Gordon equation in $\mathbb{R}^{n+1}$ that belong to $H^k(\mathbb{R}^{n+1})$ are weak solutions: for all $f\in\mathscr{S}(\mathbb{R}^{n+1})$ we have that $$u((\Box+m^2)f)=0\ .$$ In order to convert weak derivatives to classical derivatives so as to conclude that $u$ is a classical solution of the Klein-Gordon equation, one needs $u$ to be at least $\mathscr{C}^2$, which is the case if $u\in H^k(\mathbb{R}^{n+1})$ for $k>\frac{n+1}{2}+2$ thanks to the Sobolev inequality. That being said, the well-posedness of the Cauchy problem for the Klein-Gordon equation stems from the so-called energy estimates: for all $k\in\mathbb{N}$, $T>0$ we have $C_k>0$ such that for all $u\in\mathscr{S}(\mathbb{R}^{n+1})$ $$\begin{split}\|u(t,\cdot)\|_{H^k(\mathbb{R}^n)}&+\|\partial_0 u(t,\cdot)\|_{H^{k-1}(\mathbb{R}^n)}\leq C_k\left(\|u(0,\cdot)\|_{H^k(\mathbb{R}^n)}+\|\partial_0 u(0,\cdot)\|_{H^{k-1}(\mathbb{R}^n)}\phantom{\int^T_0}\right.\\ &\left.+\int^T_0\|(\Box+m^2)u(s,\cdot)\|_{H^{k-1}(\mathbb{R}^n)}ds\right)\ ,\quad |t|\in[0,T]\ .\end{split}$$ $H^k$ norms of higher-order time derivatives of $u$ for $k\geq 0$ also satisfy similar estimates once we appeal to the identity $\partial_0^2=(\Box+m^2)+\sum^n_{j=1}\partial_j^2-m^2$ and use the bound $$\begin{split}\|\partial_0^2 u(0,\cdot)\|_{H^k(\mathbb{R}^n)}&\leq\|(\Box+m^2)u(0,\cdot)\|_{H^k(\mathbb{R}^n)}+C\|u(0,\cdot)\|_{H^{k+2}(\mathbb{R}^n)}\\ &\leq C\|u(0,\cdot)\|_{H^{k+2}(\mathbb{R}^n)}+D\int^T_0\left(\|(\Box+m^2)u(s,\cdot)\|_{H^k(\mathbb{R}^n)}\right.\\ &\phantom{\leq}\left.+\|(\Box+m^2)\partial_0 u(s,\cdot)\|_{H^k(\mathbb{R}^n)}\right)ds\end{split}$$ for suitable constants $C,D>0$ (the last estimate comes from the one-dimensional Sobolev inequality on $[0,T]$, see e.g. Section 8.2 of the book Functional Analysis, Sobolev Spaces and Partial Differential Equations by H. Brézis (Springer-Verlag, 2011)). Unfortunately, the control of higher-order time derivatives by means of energy estimates is omitted or at best only glossed over in the vast majority of textbooks ("details are left as an exercise to the reader"). Putting all these together allows us to get the following, slightly weaker form for the energy estimates: for all $T>0$, $k\in\mathbb{N}$ we have constants $C'_k>0$ such that for all $u\in\mathscr{S}(\mathbb{R}^{n+1})$ $$\begin{split}\|u\|_{H^k([0,T]\times\mathbb{R}^n)}&\leq C'_k\left(\|u(0,\cdot)\|_{H^k(\mathbb{R}^n)}+\|\partial_0 u(0,\cdot)\|_{H^{k-1}(\mathbb{R}^n)}\right.\\&\phantom{\leq}\left.+\|(\Box+m^2)u\|_{H^{k-1}([0,T]\times\mathbb{R}^n)}\right)\ .\end{split}$$ These estimates by themselves already guarantee (by density of $\mathscr{S}$ in $H^k$ for all $k\in\mathbb{Z}$) uniqueness of $H^k$ weak solutions of the Klein-Gordon equation in $[0,T]\times\mathbb{R}^n$. They can also be extended to $k\in-\mathbb{N}_0$ by appealing to the Green operator ( = fundamental solution) $(1-\Delta)^{-1}$, where $\Delta$ is the ($n$- or $(n+1)$-dimensional) Laplacian, observing that such an operator commutes with all partial derivatives and therefore with the Klein-Gordon operator as well.
The existence of weak solutions in $H^k([0,T]\times\mathbb{R}^n)$ for initial data in $H^k,H^{k-1}$ can be achieved by appealing to the duality between positive-order and negative-order Sobolev spaces (up to an appropriate prescription of initial conditions at $t=0$ and $t=T$). Since the Klein-Gordon operator is formally self-adjoint, the $H^k$ energy estimate allows us to define a weak solution as a continuous linear functional in the image of (the formal adjoint of) $\Box+m^2$ in $H^{-k+1}$ for all $k\in\mathbb{Z}$, suitably extended to the whole codomain $H^{-k-1}$ by means of the Hahn-Banach theorem (uniqueness of the extension is guaranteed by the energy estimates in $H^{k+1}$). The importance of extending the energy estimates to $k\in-\mathbb{N}_0$ is that it allows us to obtain weak solutions in $H^k$ for positive $k$ by this duality technique. As mentioned above, if $k>\frac{n+1}{2}+2$ we conclude that the weak solution $u$ is $\mathscr{C}^2$ and therefore is a classical solution of the Klein-Gordon operator, thanks to the Sobolev inequality.
(Edit: on the other hand, also due to the Sobolev inequality, one obtains a unique distributional (weak) solution $u\in\mathscr{D}'(\mathbb{R}^{n+1}$) for given distributional sources and initial data, since the energy estimates for (the formal adjoint of) $\Box+m^2$ allow one to control a complete family of seminorms for test functions supported in any given compact subset. Uniqueness of distributional (weak) solutions then follows from the energy estimates by density of $\mathscr{D}$ in $\mathscr{D}'$. Moreover, one can also localize the energy estimates in space - more precisely, within slices of past- or future-directed lightcones -, which then also yield the finite speed of propagation property of solutions that characterizes hyperbolic PDE's such as the Klein-Gordon equation)
This duality method for proving well-posedness of the Cauchy problem for the Klein-Gordon equation (as well as other linear PDE's) goes back to Garding and Lax. One can find a terse exposition of this strategy for the (possibly variable-coefficient) wave equation in Chapter I of the book Lectures on Non-Linear Wave Equations by C. D. Sogge (2nd. edition, International Press, 2008).
