Digital rock physics provides a useful way to compute effective petrophysical properties (i.e., elastic modulus and permeability) of porous media from digital images. With increase in sample size to capture multimodal and multiscale microstructures, conventional approaches based on direct numerical simulation (DNS) are becoming very computationally intensive or even infeasible. To address this computational challenge, we propose a hierarchical homogenization method (HHM) to estimate effective petrophysical properties of digital rocks with large sample sizes up to billions of voxels.  The hierarchical approach divides large digital rocks into small sub-volumes and separately predicts their effective properties using DNS or data-driven methods. The effective petrophysical properties of the full digital rock is then predicted from the upscaled model at the intermediate scale. Compared with the classic homogenization method directly from the micro-scale to macro-scale, this hierarchical homogenization approach considerably reduces the computational costs and makes it applicable to large digital rocks.