One of the most critical challenges in today's CMOS VLSI design is the lack of predictability in chip performance at design stage. One of the process variabilities comes from the voltage drop variations in on-chip power distribution networks. In this paper, we present a novel analysis approach for computing voltage drops of large power grid networks under process variations. The new algorithm is very efficient and scalable for huge networks with a large number of variational variables. This approach, called variational extended truncated balanced realization ( varETBR), is based on model order reduction techniques to reduce the circuit matrices before the variational simulation. It performs the parameterized reduction on the original system using variation-bearing subspaces. After the reduction, Monte Carlo based statistical simulation is performed on the reduced system and the statistical responses of the original system are obtained thereafter. varETBR calculates variational response Grammians by Monte Carlo based numerical integration considering both system and input source variations in generating the projection subspace. varETBR is very scalable for the number of variables and flexible for different variational distributions and ranges as demonstrated in experimental results. Experimental results, on a number of IBM benchmark circuits up to 1.6 million nodes, show that the varETBR can be 1900 X faster than the Monte Carlo method and is much more scalable than one of the recently proposed approaches.