Investigating molecular interactions between cellulose surfaces and hemicelluloses is crucial for understanding plant cell wall architecture. Here, we employed molecular dynamics (MD) simulations to study how strongly xylan, the most common type of hemicelluloses, interact with cellulose surfaces. We systematically investigated effects from xylan substitution pattern, degree of polymerization, cellulose surface structure and the explicit water models. We studied four xylan models: a bare xylan without substitutions, an acetylated xylan where every other xylose unit is acetylated, and two acetylated xylan models with a Methyl-D-glucuronic unit (either in the middle of the chain or at the end). We then calculated the free energy of adsorption to see how strongly each of these xylans adhere to the hydrophobic (100) and the hydrophilic (110) cellulose surfaces. Our results show that the acetylated xylan adheres more strongly to the cellulose surfaces, compared to the unsubstituted xylan. This is more pronounced on the (100) cellulose surface and is consistent for all the common water models including SPC, SPC/E, TIP3P and TIP4P_2005. The presence of Methyl-D-glucuronic unit also enhanced the adhesion to the (100) surface, but it reduces the adhesion to the (110) surface. We also observed a linear relationship between the xylan degree of polymerization and the predicted adsorption energies, which can be used for predicting the adhesion of xylan polymers to the cellulose surfaces. Such molecular level understanding of xylan-cellulose adhesion is provide insights into plant cell wall biosynthesis, and is importnat for better understanding of the plant cell wall architecture.