A new many-valley theoretical formulation of the strain-dependent shallow donor polarizability is developed which should yield better agreement with experiment for the case when the $1s\ensuremath{-}{A}_{1}\ensuremath{-}1s\ensuremath{-}{E}_{a}$ valley-orbit splitting $6\ensuremath{\Delta}(0)$ is much larger than the strain-induced splitting of the valleys. The ground-state wave function ${\ensuremath{\psi}}_{\mathrm{GS}}(\stackrel{\ensuremath{\rightarrow}}{\mathrm{r}},x)$, $x$ the reduced valley strain, is taken as ${\ensuremath{\psi}}_{\mathrm{GS}}(\stackrel{\ensuremath{\rightarrow}}{\mathrm{r}},x)=\ensuremath{\gamma}(x){\ensuremath{\psi}}_{A1}(\stackrel{\ensuremath{\rightarrow}}{\mathrm{r}},{a}_{A}(x))+\ensuremath{\beta}(x){\ensuremath{\psi}}_{{E}_{a}}(\stackrel{\ensuremath{\rightarrow}}{\mathrm{r}},{a}_{E}(x))$ where $\ensuremath{\gamma}(x)$ and $\ensuremath{\beta}(x)$ are determined by the valley-repopulation model. ${a}_{A}(x)$ and ${a}_{E}(x)$ are strain-dependent Bohr radii, which differ by 16% for the P donor in Si for $x=0$ and this difference grows with $x$. Unlike the previous approach each valley wave function component consists of two parts with different strain-dependent Bohr radii. ${\ensuremath{\psi}}_{\mathrm{GS}}(\stackrel{\ensuremath{\rightarrow}}{\mathrm{r}},x)$ is employed, using the Hass\'e method, to calculate the strain-dependent donor polarizability for P and Sb donors in Si. The results are compared to the experimental data for P and Sb donors in Si. For the isocoric P donor the new calculated results are in better agreement with the experimental results. However, one must still incorporate a strain-dependent valley-orbit splitting parameter $\ensuremath{\Delta}(x)$ to fit the experimental data. The results of the force-fit $\ensuremath{\Delta}(x)$ case show ${a}_{A}(x)$ increasing with $x$ and the ground-state energy ${E}_{\mathrm{GS}}(x)$ decreasing with $x$---both results are the opposite of those expected for the case $\ensuremath{\Delta}(x)=\ensuremath{\Delta}(0)=\mathrm{constant}$. The donor piezohyperfine data of Wilson and Feher are reanalyzed utilizing independently measured values of $6\ensuremath{\Delta}(0)$ and the shear deformation potential ${\ensuremath{\Xi}}_{u}$. The results show that ${a}_{A}(x)$ must slowly increase with $x$, which is in qualitative agreement with the ${a}_{A}(x)$ increase inferred from the polarizability results.