These few-site entropies tend to converge fairly fast with finite size, because they ultimately depend on very local physics. If you switch to DMRG you should easily see the values converge fully. You can also analytically obtain the single- and two-site reduced density matrices of transverse-field Ising and XY chains in the thermodynamic limit, see Osborne and Nielsen, Phys. Rev. A 66, 032110 (2002). In their notation the Ising chain reads [Eq. (1) with $\gamma=1$]
$$
H = -\sum_j \left( \lambda \sigma^x_j \sigma^x_{j+1} + \sigma_j^z \right),
$$
so the case $\lambda=1$ corresponds to your $h=1$.
The one-site reduced density matrix is given by [Eq. (25)]
$$
\rho_1 = \frac{I+\langle \sigma^z\rangle \sigma^z}{2},
$$
where in the ground state [Eq. (35)]
$$
\langle \sigma_z \rangle = \frac{1}{\pi} \int_0^\pi \mathrm{d}\phi \frac{1+\lambda\cos \phi }{\sqrt{1+\lambda^2+2\lambda \cos\phi}}.
$$
Evaluating the entropy at $\lambda=1$, where $\langle \sigma^z\rangle = 2/\pi$, one finds
$$
\begin{align}
S_1 &= -Tr \left[ \rho_1 \log \rho_1 \right] / \log{2}\\
&= \frac{\pi \left(\log (4)+2 \log (\pi )-\log \left(\pi ^2-4\right)\right)-4 \tanh ^{-1}\left(\frac{2}{\pi }\right)}{\pi \log (4)} \approx 0.68376.
\end{align}
$$
Your largest system size gets quite close!
The two-site reduced density matrix for two sites separated by a distance $r$ is given by [Eq. (26)]
$$
\rho_{0r} =\frac{I_{0r}+\left<\sigma^z\right>\left(\sigma_0^z+\sigma_r^z\right)+\sum_{k=1}^3\left<\sigma_0^k\sigma_r^k\right>\sigma_0^k\sigma_r^k}{4}.
$$
The expression involves spin-spin correlations $\langle \sigma_0^k \sigma_r^k\rangle$, which are known at criticalty ($\lambda=1$), see Pfeuty, Ann. Phys. (N.Y.) 57, 79-90 (1970) [Eqs. (3.3)-(3.5)].
$$
\begin{align}
\langle \sigma_0^x\sigma_r^x\rangle_c &=\left(\frac{2}{\pi}\right)^r 2^{2r(r-1)}\frac{H(r)^4}{H(2r)}\\
\langle\sigma_0^y\sigma_r^y\rangle_c &=-\frac{\langle \sigma_0^x\sigma_r^x\rangle }{4r^2-1}\\
\langle\sigma_0^z\sigma_r^z\rangle_c&=\frac{4}{\pi^2}\frac{1}{4r^2-1}
\end{align}
$$
where $H(r)=1^{r-1}2^{r-2}\dots (r-1)$. (Osborne and Nielsen reproduce these as Eqs. (49)-(51), but lack a factor of $1/\pi$ for $\langle \sigma^z_0 \sigma^z_r\rangle_c$.) The subscript ${}_c$ denotes that these are connected correlation functions, i.e.
$$
\langle \sigma_0^k \sigma_r^k \rangle_c = \langle \sigma_0^k \sigma_r^k \rangle - \langle \sigma_0^k\rangle \langle \sigma_r^k\rangle,
$$
but note that the two-site reduced density matrix uses the full expression $\langle \sigma_0^k \sigma_r^k \rangle$, see the answer to this related question. At criticality only $\langle \sigma^z\rangle$ is non-zero, so the disconnected part only matters for the $zz$-correlation.
Approximate values for the nearest-neighbor correlations are $\langle \sigma_0^x\sigma_1^x\rangle_c\approx 0.63662$, $\langle \sigma_0^y\sigma_1^y\rangle_c\approx -0.212207$, and $\langle \sigma_0^z\sigma_1^z\rangle_c\approx 0.135095 \quad \Rightarrow \quad \langle \sigma_0^z\sigma_1^z\rangle\approx 0.54038$.
Evaluation of the nearest-neighbors two-site entropy with normalization as in the question then yields
$$
S_{12} = -Tr \left[ \rho_{01} \log \rho_{01} \right]/(2\log 2) \approx 0.428054,
$$
which is also consistent with the numerical result.
