In the present study, we investigated the solubility of celecoxib the in the aqueous solution of two deep eutectic solvents (DESs) at various DES weight fractions within the temperature range of 283.2-323.2 K. Choline chloride (ChCl) as hydrogen bond acceptor (HBA) and ethylene glycol (EG) and 1,2-propanediol (PD) as hydrogen bond donors (HBDs) were used for the preparation of DESs. Experimental results showed that the solubility of celecoxib is significantly enhanced by approximately one million times in DES-contained PD compared to pure water. Also, celecoxib's solubility was enhanced by increasing temperature and the weight fraction of DESs. Besides, DES2 (ChCl + PD) showed a higher solubility for celecoxib compared to DES1 (ChCl + EG) in all experiments. Moreover, two temperature-dependent Buchowski-Ksiazaczak λh and modified Apelblat equations and the local composition models such as Wilson, non-random two-liquid theory (NRTL), and universal quasi-chemical (UNIQUAC) were applied to correlate the experimental solubility data. The correlated solubility was in good agreement with the experimental data. In addition, two-segment-based models, including non-random two-liquid-segment activity coefficient (NRTL-SAC) and universal quasi-chemical segment activity coefficient (UNIQUAC-SAC) models provided relatively good performance in solubility prediction. Furthermore, the calculation of the apparent standard thermodynamic properties showed that enthalpy is the primary contribution term to the apparent standard Gibbs free energy in the dissolution process of celecoxib. Finally, we applied molecular dynamics simulations to describe the drug solubility from the molecular point of view. The radial distribution function (RDF) analysis revealed that the experimental solubility is consistent with the MD simulation results, so the attractive interactions between celecoxib-DES1 and celecoxib-DES2 are much higher than the water-celecoxib.