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.

@DanielLichtblau, perfect! That's the answer I was looking for. My goal is indeed to find all real solutions of underdetermined polynomial systems. Since NSolve was giving some results for certain cases, I got confused. But since you confirmed that NSolve is not going to do that for me, that's all I wanted to know. A suggestion though, it may be good to mention this in the Mathematica documentation so that people don't get confused.

Thanks again.

@MichaelRogers, Thanks for your response. No, I did not get that error, though I used to get that error in earlier versions of Mathematica (now I have 11.2). Regardless of that error, I think technically these algorithms shouldn't give real isolated solutions for underdetermined systems. The Groebner basis method is defined over complex variables and gives complex solutions. If the original system has only isolated solutions, the Groebner Basis method gives all isolated complex solutions which would include isolated real solutions as well. So all is good.

However, if the original system has infinitely many solutions (underdetermined), and we add as many random hyperplanes (or linear conditions, as you say) as required to completely determine the system (i.e., only isolated solutions exist), the Groebner basis method will then give us complex isolated solutions for the new and bigger system. That's because you are adding complex hyperplanes (remember, the GB mehod is strictly defined over complex variables).

The fact that the NSolve command gives us real isolated solutions back is suggesting that it may be using some additional method/trick to extract the real isolated solutions from complex solution curves of the original system. So, my question is much more specific than Frank's response. In any case, Mathematica documentation should at least tell us what method it is using, rather than its user figuring out things. Otherwise, I can't use the results generated by it into a scientific publication. Again, I don't need any proprietary information here obviously. I just need to know if the method makes sense and kosher.

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