Rare earth elements (REE) are critical elements found in monazite, xenotime, and hydrated REE phosphates which typically form in hydrothermal mineral deposits. Accurate predictions of the solubility of these REE phosphates and the speciation of REE in aqueous fluids are both key to understanding the controls on the transport, fractionation, and deposition of REE in natural systems. Previous monazite and xenotime solubility experiments indicate the presence of large discrepancies between experimentally derived solubility constants versus calculated solubilities by combining different data sources for the thermodynamic properties of minerals and aqueous species at hydrothermal conditions. In this study, these discrepancies were resolved by using the program GEMSFITS to optimize the standard partial molal Gibbs energy of formation (Δ_fG^°₂₉₈) of REE aqueous species (REE³⁺ and REE hydroxyl complexes) at 298.15 K and 1 bar while keeping the thermodynamic properties fixed for the REE phosphates. A comprehensive experimental database was compiled using solubility data available between 25 and 300 °C. The latter permits conducting thermodynamic parameter optimization of Δ_fG^°₂₉₈ for REE aqueous species. Optimal matching of the rhabdophane solubility data between 25 and 100 °C requires modifying the Δ_fG^°₂₉₈ values of REE³⁺ by 1-6 kJ/mol, whereas matching of the monazite solubility data between 100 and 300 °C requires modifying the Δ_fG^°₂₉₈ values of both REE³⁺ and REEOH²⁺ by ∼ 2-10 kJ/mol and ∼ 15-31 kJ/mol, respectively. For xenotime, adjustments of Δ_fG^°₂₉₈ values by 1-26 kJ/mol are only necessary for the REE³⁺ species. The optimizations indicate that the solubility of monazite in acidic solutions is controlled by the light (L)REE³⁺ species at <150 °C and the LREEOH²⁺ species at >150 °C, whereas the solubility of xenotime is controlled by the heavy (H)REE³⁺ species between 25 and 300 °C. Based on the optimization results, we conclude that the revised Helgeson-Kirkham-Flowers equation of state does not reliably predict the thermodynamic properties of REE³⁺, REEOH²⁺, and likely other REE hydroxyl species at hydrothermal conditions. We therefore provide an experimental database (ThermoExp_REE) as a basic framework for future updates, extensions with other ligands, and optimizations as new experimental REE data become available. The optimized thermodynamic properties of aqueous species and minerals are available open access to accurately predict the solubility of REE phosphates in fluid-rock systems.