We calculate Faraday rotation in global axisymmetric magnetohydrodynamic simulations of geometrically thick accretion flows. These calculations are motivated by the measured rotation measure (RM) of -6 * 105 rad m-2 from Sgr A* in the Galactic center, which appears to have been stable over the past 7 years. In our numerical simulations, the quasi-steady state structure of the accretion flow, and the RM it produces, depends on the initial magnetic field threading the accreting material. In spite of this dependence, we can draw several robust conclusions about Faraday rotation produced by geometrically thick accretion disks: i) the time averaged RM does not depend that sensitively on the viewing angle through the accretion flow, but the stability of the RM can. Equatorial viewing angles show significant variability in RM (including sign reversals), while polar viewing angles are relatively stable if there is a large scale magnetic field threading the disk at large radii. ii) Most of the RM is produced at small radii for polar viewing angles while all radii contribute significantly near the midplane of the disk. Our simulations confirm previous analytic arguments that the accretion rate onto Sgr A* must satisfy \dot M_in \ll \dot M_Bondi 10-5 \mpy in order to not over-produce the measured RM. We argue that the steady RM -6 * 105 rad m-2 from Sgr A* has two plausible explanations: 1) it is produced at 100 Schwarzschild radii, requires \dotM_in 3 * 10-8 Mo yr-1, and we view the flow at an angle of 30^o relative to the rotation axis of the disk; in our simulations, the variation in RM across a finite-sized source is sufficient to depolarize the emission below 100 GHz, consistent with observations. 2) Alternatively, the RM may be produced in the relatively spherical inflowing plasma near the circularization radius at 103-104 Schwarzschild radii, with the magnetic field perhaps amplified by the magnetothermal instability. Time variability studies of the RM can distinguish between these two possibilities.