During the Late Pleistocene and Holocene retreat of paleo-ice sheets in North America and Europe, vast proglacial lakes existed along the land terminating margins. These proglacial lakes impacted ice sheet dynamics by imposing boundary conditions analogous to a marine terminating margin. Such lacustrine boundary conditions cause changes in the ice sheet geometry, stress balance and frontal ablation and therefore affect the mass balance of the entire ice sheet. Despite this, dynamically evolving proglacial lakes have rarely been considered in detail in ice sheet modeling endeavors. In this study, we describe the implementation of an adaptive lake boundary in the Parallel Ice Sheet Model (PISM), which we call PISM-LakeCC. We test our model with a simplified glacial retreat setup of the Laurentide Ice Sheet (LIS). By comparing the experiments with lakes to control runs with no lakes, we show that the presence of proglacial lakes locally enhances the ice flow, which leads to a lowering of the ice sheet surface. In some cases, this also results in an advance of the ice margin and the emergence of ice lobes. In the warming climate, increased melting on the lowered ice surface drives the glacial retreat. For the LIS, the presence of lakes triggers a process similar to marine ice sheet instability, which caused the collapse of the ice saddle over Hudson Bay. In the control experiments without lakes, Hudson Bay is still glaciated when the climate reaches present-day (PD) conditions. The results of our study demonstrate that glacio-lacustrine interactions play a significant role in the retreat of land terminating ice sheet margins.