We describe an Enriched Algebraic Multiscale Solver (EAMS) that overcomes the deficiency of existing multiscale methods for flow in heterogeneous media with large coherent correlation structures and high contrasts. For a given multiscale method, EAMS enriches the coarse space with local basis-functions specifically aimed at the largest error components in the solution space. For this purpose, the discrete error equation is used to identify the solution modes that are missing from the multiscale operator. The identified error modes, which are complex combinations of a spectrum of wave numbers, are then localized (truncated) and added to the prolongation operator. The enrichment process is repeated iteratively until the desired convergence rate is reached. The identification and enrichment processes are algebraic, and they are performed adaptively during the iterative solution process. Using challenging test cases from the literature, we show that EAMS leads to great improvements in the robustness and efficiency of existing state-of-the-art multiscale linear solvers. In most settings, the convergence rate of AMS is improved significantly by supplementing the standard basis functions with a few basis functions guided by the error equation. Since the enrichment is adaptive and algebraic, it can be integrated into any existing multiscale linear-solver implementation. EAMS is expected to be most useful in modeling evolution multiphase problems in heterogeneous reservoirs, whereby the changes in the character of the linear system - across Newton iterations and time steps - are relatively mild.