The nonlinear gravitational-wave memory effect is a prediction of general relativity in which test masses are permanently displaced by gravitational radiation. We implement a method for calculating the expected memory waveform from an oscillatory gravitational-wave time series. We use this method to explore the phenomenology of gravitational-wave memory using a numerical relativity surrogate model. Previous methods of calculating the memory have considered only the dominant oscillatory (=2, m=|2|) mode in the spherical harmonic decomposition or the post-Newtonian expansion. We explore the contribution of higher-order modes and reveal a richer phenomenology than is apparent with =|m|=2 modes alone. We also consider the "memory of the memory" in which the memory is, itself, a source of memory, which leads to a small, O(10-4), correction to the memory waveform. The method is implemented in the python package GWMemory, which is made publicly available.