We consider the interaction-round-a-face version of the isotropic six-vertex model. The associated spin chain is made of two coupled Heisenberg spin chains with different boundary twists. The phase diagram of the model and the long distance correlations were studied in [Nucl. Phys. B, 995 (2023) 116333]. Here, we compute the short-distance correlation functions of the model in the ground state for finite system sizes via non-linear integral equations and in the thermodynamic limit. This was possible since the model satisfies the face version of the discrete quantum Knizhnik-Zamolodchikov (qKZ) equation. A suitable ansatz for the density matrix is proposed in the form of a direct sum of two Heisenberg density matrices, which allows us to obtain the discrete functional equation for the two-site function ω(λ1,λ2). Thanks to the known results on the factorization of correlation functions of the Heisenberg chain, we are able to compute the density matrix of the IRF model for up to four sites and its associated spin chain for up to three sites.