Simulations of ice sheet evolution over glacial cycles requires integration of observational constraints using ensemble studies with fast ice sheet models. These include physical parameterisations with uncertainties, for example, relating to grounding line migration. Ice dynamically more complete models are slow and have thus far only be applied for <1,000 years, leaving many model parameters unconstrained. Here we apply a 3D thermomechanically coupled full-Stokes ice sheet model to the Ekström Ice Shelf embayment, East Antarctica, over a full glacial cycle (40,000 years). We test the model response to differing ocean bed properties that provide an envelope of potential ocean substrates seawards of today’s grounding line. The end member scenarios include a hard, high friction ocean bed and a soft, low friction ocean bed. We find that predicted ice volumes differ by >50% under almost equal forcing. Grounding line positions differ by up to 49 km, show significant hysteresis, and migrate non-steadily in both scenarios with long quiescent phases disrupted by leaps of rapid migration. The simulations quantify evolution of two different ice sheet geometries (namely thick and slow vs. thin and fast), triggered by the variable grounding line migration over the differing ocean beds. Our study extends the timescales of 3D full-Stokes by an order of magnitude to previous studies with the help of parallelisation. The extended time frame for full-Stokes models is a first step towards better understanding other processes such as erosion and sediment redistribution in the ice shelf cavity impacting the entire catchment geometry.