Engineers often use seismic site response analyses to evaluate the propagation of seismic waves through a soil deposit. These analyses have traditionally relied on dynamic properties (e.g., the maximum shear modulus, Gmax, and shear modulus reduction curves G/Gmax) measured on homogeneous soil specimens. However, soil properties are not uniform, and it is common to see variabilities at different scales. This study uses physical models with controlled variability to evaluate the Gmax and G/Gmax in heterogeneous soil samples to reflect the effects of soil variability. We used a 3D soil printer to create heterogeneous soil models of 20x20x20 cm that match target random field simulations generated by spatially varying the liquid limit on clayey-type soils. The models were used to sample cylindrical specimens, then tested in a torsional shear rheometer to measure Gmax and G/Gmax. The results show that Gmax and the G/Gmax curves are affected by the internal specimen variability, particularly by the coefficient of variation (COV) in the underlain random fields. Using the experimental results, we also show that the Hill  approach for estimating properties of composite materials provides a good estimation for the experimentally measured G/Gmax curves in a heterogeneous soil media, especially in specimens with low to medium variability. Finally, we evaluated the performance of the homogenization bounding theory of Hashin and Shtrikman , also formulated for composite materials, on estimating Gmax. We find that the Gmax of specimens with low and medium variability can be reasonably estimated using the upper bound. In contrast, samples with high variability are better estimated using the lower bound.