Molecular dynamics simulations have been performed to investigate the atomistic deformation mechanisms of hierarchically nanotwinned Cu under nanoindentation. When the grain size (d) and the spacing of primary twins (lambda (1)) are fixed, the hardness is observed to increase with decreasing spacing of secondary twins (lambda (2)) until a critical lambda (2), and then decrease with further decreasing lambda (2). The size effect of lambda (2) on the plastic area beneath the indenter is observed to be exactly opposite to the trend of the size effect on the hardness. There exist two plastic zones beneath the indenter: the severe plastic zone and the moderate plastic zone. In the severe plastic zone, high density of dislocation networks are observed and the deformation mechanisms are independent of lambda (2). The deformation mechanisms in the moderate plastic zone are highly dependent on the lambda (2), which is the origin of the size effect on the hardness. Below the critical lambda (2), the deformation mechanisms are dominated by the softening mechanisms with decreasing lambda (2): (i) detwinning of secondary twins and (ii) nucleation and propagation of partial dislocations with a small angle to the boundaries of secondary twins. Above the critical lambda (2), the deformation mechanisms are dominated by the strengthening mechanisms with decreasing lambda (2): partial dislocations are blocked by the boundaries of primary twins or secondary twins.