Try

Trace[
 NSolve[x^2 + y^2 - 1 == 0, {x, y}, Reals],
 _GroebnerBasis,
 TraceInternal -> True
 ]

It does use GroebnerBasis. I do not see why NSolve cannot use complex methods and then verify "all variables, parameters, constants, and function values are restricted to be real."

@DanielLichtblau, Thanks for your message! For this particular example, it would be straightforward, sure. That's why I had made a comment in my original question "While this example may be easy enough, Mathematica also finds real solutions for larger systems." How does Mathematica find isolated real solutions out of underdetermined system of polynomial equations when the number of equations/variables are larger? The GB method assumes complex variables. For underdetermined systems, it is very tricky to extract real isolated solutions which would now be embedded in complex positive dimensional components. I can send you a more complicated problem if you want. Or I think you would already know such examples anyway.

I get this message when I run your code:

NSolve::infsolns: Infinite solution set has dimension at least 1. Returning intersection of solutions with
-((92291 x)/87992)-(121001 y)/175984 == 1.

It is omitted from your quoted result. Did you not get it? Obviously, NSolve adds linear conditions until the system is not underdetermined. (I say, obviously, because that's what the message says.) With that in mind, Frank's response seems to be exactly what is done to find the solutions (at least the Gröbner basis part).

For the following equation, it adds two hyperplanes:

NSolve[x^2 + y^2 + z^2 - 1 == 0, {x, y, z}, Reals]

By the way, you can turn off that method and get parametrized solutions (at least for many algebraic equations):

SetSystemOptions["NSolveOptions" -> "UseSlicingHyperplanes" -> False]
NSolve[x^2 + y^2 + z^2 - 1 == 0, {x, y, z}, Reals]

NSolve::svars: Equations may not give solutions for all "solve" variables.

Out[] = {
 {z ->  ConditionalExpression[-1.` Sqrt[1.` - 1.` x^2 - 1.` y^2],
   -1.` < x < 1.` && -1.` Sqrt[1.` - 1.` x^2] <= y <= Sqrt[1.` - 1.` x^2]]},
 {z -> ConditionalExpression[Sqrt[1.` - 1.` x^2 - 1.` y^2],
   -1.` < x < 1.` && -1.` Sqrt[1.` - 1.` x^2] <= y <= Sqrt[1.` - 1.` x^2]]},
 {x -> -1.`, y -> 0, z -> 0},
 {x -> 1.`, y -> 0, z -> 0}}

For sparse linear systems, Solve and NSolve use several efficient numerical methods, mostly based on Gauss factoring with Markowitz products (approximately 250 pages of code). For systems of algebraic equations, NSolve computes a numerical Gröbner basis using an efficient monomial ordering, then uses eigensystem methods to extract numerical roots.

http://reference.wolfram.com/language/tutorial/SomeNotesOnInternalImplementation.html