A major challenge in molecular simulations is to describe denaturant-dependent folding of proteins order to make direct comparisons with {\it in vitro} experiments. We use the molecular transfer model, which is currently the only method that accomplishes this goal albeit phenomenologically, to quantitatively describe urea-dependent folding of PDZ domain, which plays a significant role in molecular recognition and signaling. Experiments show that urea-dependent unfolding rates of the PDZ2 domain exhibit a downward curvature at high urea concentrations, which has been interpreted by invoking the presence of a sparsely populated high energy intermediate. Simulations using the MTM and a coarse-grained model of PDZ2 are used to show that the intermediate, which has some native-like character, is present in equilibrium both in the presence and absence of urea. The free energy profiles show that there are two barriers separating the folded and unfolded states. Structures of the transition state ensembles, (\(TSE1\) separating the unfolded and \(I_{EQ}\) and \(TSE2\) separating \(I_{EQ}\) and the native state), determined using the \(P_{fold}\) method, show that \(TSE1\) is expanded; \(TSE2\) and native-like. Folding trajectories reveal that PDZ2 folds by parallel routes. In one pathway folding occurs exclusively through \(I_1\), which resembles \(I_{EQ}\). In a fraction of trajectories, constituting the second pathway, folding occurs through a combination of \(I_{1}\) and a kinetic intermediate. The radius of gyration (\(R_g^{U}\)) of the unfolded state is more compact (by \(\sim\) 9\%) under native conditions. Decrease in \(R_g^{U}\) occurs on the time scale on the order of utmost \(\sim\) 20 \(\mu s\). The modest decrease in \(R_g^{U}\) and the rapid collapse suggest that high spatial and temporal resolution are needed to detect compaction in finite-sized proteins.