Yttrium plays an important role as a radiochemical dosimetry detector for determining high energy 14 MeV neutron fluences, through measurement of (n,2n) activation products. The total (n,2n) cross section is known rather well from extensive activation measurements on the stable 89Y isotope, and from a previous activation measurement on the unstable 88Y ground state. However, until now the branching ratios to the ground state and excited isomers in 88Y via the 89Y (n,2n) reaction were not well known, and furthermore, uncertainty estimates were not available for these cross sections and branching ratios. This paper describes how gamma-ray transitions between states in (n,2n) and (n,n') reactions measured using the GEANIE detector at Los Alamos' LANSCE facility, together with theory predictions using the GNASH code, enable us to determine these quantities for the ENDF/B-VII evaluation. A previous measurement by Dietrich, at Livermore, provided important complementary information to the GEANIE analysis. We describe an uncertainty quantification analysis that uses the GNASH-KALMAN approach to evaluate cross sections for the 89Y (n,2n) population of the 88Y ground state and two meta-stable isomers (m1 and m2), along with their uncertainties. Our new results agree with Arthur's historic Los Alamos evaluated cross sections within a few percent below 15 MeV (with larger differences above 15 MeV). The (n,2n) cross sections to the 88Y ground state and m1, m2 isomers impact the average 88Y(n,2n) 87Y cross section at leading-order; we determine this 14.1 MeV average cross section 88Y(n,2n) 87Y = 1107 mb (± 4%), which agrees with the value obtained from Arthur's evaluation to 0.7%. An alternative method to predict cross sections, uncertainties, and covariance data, is described that uses the European TALYS reaction modeling code and a Backward-Forward Monte-Carlo uncertainty quantification technique. This approach uses a microscopic optical model, together with Hauser-Feshbach and preequilibrium reaction mechanisms, and the underlying model parameters and their uncertainties and correlations are determined through a Monte-Carlo filtering method based on comparisons with measured cross section data. We compare the results obtained using this approach with the GNASH-KALMAN method. The evaluated cross sections are rather similar in the two approaches. We show how the uncertainty information, as embodied in the resulting covariance matrices, is also qualitatively similar in both approaches.