In this work, we present results from molecular dynamics simulations of aqueous solutions of furfural and 5-hydroxymethylfurfural, which are important intermediates in the hydrothermal carbonization processes of biomass conversion. The computations were performed both at ambient and hydrothermal conditions using a two-level factorial design varying concentration, temperature, and pressure. A number of equilibrium and dynamic properties have been computed including enthalpies and free energies of vaporization, free energies of solvation, diffusion coefficients, and rotational/reorientational correlation times. Structural properties of solutions were analyzed using radial and spatial distribution functions. It was shown that the formation of hydrogen bonds among 5-hydroxymethylfurfural molecules is preferred compared to hydrogen bonding between 5-hydroxymethylfurfural and water. In addition, our results suggest that the oxygen atoms in the furan rings of furfural and 5-hydroxymethylfurfural do not participate in hydrogen bonding to the same extent as the oxygen atoms in the hydroxyl and carbonyl groups. It is also observed that furfural molecules aggregate under certain conditions, and we show how this is affected by changes in temperature, pressure, and concentration in agreement with experimental solubility data. The analysis of the computational results provides useful insight into the structure and dynamics of the considered molecules at conditions of hydrothermal carbonization, as well as at ambient conditions.