In this paper we derive analytical solutions to statistical moments for transient saturated flow in two-dimensional, bounded, randomly heterogeneous porous media. By perturbation expansions, we first derive partial differential equations governing the zeroth-order head h(0) and the first-order head term h(1), where orders are in terms of the standard deviation of the log transmissivity. We then solve h(0) and h(1) analytically, both of which are expressed as infinite series. The head perturbation h(1) is then used to derive expressions for autocovariance of the hydraulic head and the cross covariance between the log transmissivity and head. The expressions for the mean flux and flux covariance tensor are formulated from the head moments based on Darcy's law. Using numerical examples, we demonstrate the convergence of these solutions. We also examine the accuracy of these first-order solutions by comparing them with solutions from both Monte Carlo simulations and the numerical moment equation method.