Harmonic Earth tide components in well water levels have been used to estimate hydraulic and geomechanical subsurface properties. However, the validity of various methods based on analytical solutions has not been established. First, we review the theory and examine the latest analytical solution used to relate well water levels to Earth tides. Second, we develop and verify a novel numerical model coupling hydraulics and geomechanics to Earth tide strains. Third, we assess subsurface conditions over depth for a range of realistic properties. Fourth, we simulate the well water level response to Earth tide strains within a 2D poroelastic layered aquifer system confined by a 100 m thick aquitard. We find that the analytical solution matches two observations (amplitudes and phases) to multiple unknown parameters leading to non-unique results. We reveal that undrained and confined conditions are necessary for the validity of the analytical solution. This occurs for the dominant M_2 frequency at depths >50 m and requires specific storage at constant strain of S_ε ≥ 10^-6 m^-1, in combination with hydraulic conductivity of the aquitard k_l ≥ 5*10^-5 ms^-1 and aquifer ka ≥ 10^-4 ms^-1. We further illustrate that the analytical solution is valid in unconsolidated systems, whereas consolidated systems require additional consideration of the Biot modulus. Overall, a priori knowledge of the subsurface system supports interpretation of the groundwater response. Our results improve understanding of the effect of Earth tides on groundwater systems and its interpretation for subsurface properties.