Why are the following outputs different?
N[(x^2 + y^2) /. Thread[{x, y} -> {1., 1.}], 10]
gives 2..
N[(1^2 + 1^2), 10]
gives 2.000000000.
How can I use replacement but still retain numerical precision?
Once you've lost precision, you cannot get it back. Use
N[(x^2 + y^2) /. Thread[{x, y} -> {1, 1}], 10]
or
(x^2 + y^2) /. Thread[{x, y} -> {1`10, 1`10}]