We study the resistance force of cylindrical objects penetrating quasi-statically into granular media experimentally and numerically. Simulations are validated against experiments. In contrast to previous studies, we find in both experiments and simulations that the force-depth relation consists of three regimes, rather than just two: transient and steady-state. The three regimes are driven by different dynamics: an initial matter compression, a developing stagnant zone, and an increase in steady-state force with a fully developed stagnant zone. By simulations, we explored the effects of a wide range of parameters on the penetration dynamics. We find that the initial packing fraction, the inter-granular sliding friction coefficient, and the grain shape (aspect ratio) have a significant effect on the gradient Kφ of the force-depth relation in the steady-state regime, while the rolling friction coefficient noticeably affects only the initial compression regime. Conversely, Kφ is not sensitive to the following grain properties: size, size distribution, shear modulus, density, and coefficient of restitution. From the stress fields observed in the simulations, we determine the internal friction angles φ, using the Mohr-Coulomb yield criterion, and use these results to test the recently-proposed modified Archimedes' law theory. We find excellent agreement, with the results of all the simulations falling very close to the predicted curve of φ vs. Kφ. We also examine the extreme case of frictionless spheres and find that, although no stagnant zone develops during penetration into such media, the value of their internal friction angle, φ = 9° ± 1°, also falls squarely on the theoretical curve. Finally, we use the modified Archimedes' law theory and an expression for the time-dependent growth of the stagnant zone to propose an explicit constitutive relation that fits excellently the force-depth curve throughout the entire penetration process.