Eumelanin, a member of the melanin family, is a black-brown insoluble pigment. It possesses a broad range of properties such as antioxidation, free radical scavenging, photoprotection, and charge carrier transportation. Surprisingly, the exact molecular structure of eumelanin remains undefined. It is, however, generally considered to consist of two main building blocks, 5,6-dihydroxyindole (DHI) and 5,6- dihydroxyindole carboxylic acid (DHICA). We focus on DHI and report, for the first time, a computational investigation of the structural properties of DHI-eumelanin aggregates in aqueous solutions. First, multimicrosecond molecular dynamics (MD) simulations at different concentrations were performed to investigate the aggregation and ordering of tetrameric DHI-eumelanin protomolecules. This was followed by umbrella sampling (US) and density functional theory (DFT) calculations to study the physical mechanisms of stacking. Aggregation occurs through formation of nanoscale stacks and was observed in all systems. Further analyses showed that aggregation and coarsening of the domains is due to a decrease in hydrogen bonds between the eumelanins and water; while domains exist, there is no long-range order. The results show noncovalent stacks with the interlayer distance between eumelanin protomolecules being less than 3.5 Å. This is in good agreement with transmission electron microscopy data. Both free energy calculations and DFT revealed strong stacking interactions. The electrostatic potential map provides an explanation and a rationale for the slightly sheared relative orientations and, consequently, for the curved shapes of the nanoscale domains.