def METHOD_NAME(bathymetry_2d, elev_init, store_station_time_series=True, **model_options):
    """
    Construct a *linear* shallow water equation solver for tsunami
    propagation modelling.
    """
    mesh2d = elev_init.function_space().mesh()
    t_end = 2 * 3600.0
    u_mag = Constant(5.0)
    t_export = 60.0
    dt = 60.0
    if os.getenv("THETIS_REGRESSION_TEST") is not None:
        t_end = 5 * t_export
    # Bathymetry
    bathymetry_2d.interpolate(bathymetry_2d - elev_init)
    # Create solver
    solver_obj = solver2d.FlowSolver2d(mesh2d, bathymetry_2d)
    options = solver_obj.options
    options.use_nonlinear_equations = False
    options.element_family = "dg-dg"
    options.simulation_export_time = t_export
    options.fields_to_export = ["elev_2d"]
    options.simulation_end_time = t_end
    options.horizontal_velocity_scale = u_mag
    options.swe_timestepper_type = "CrankNicolson"
    if not hasattr(options.swe_timestepper_options, "use_automatic_timestep"):
        options.timestep = dt
    options.swe_timestepper_options.solver_parameters = {
        "snes_type": "newtonls",
        "ksp_type": "gmres",
        "pc_type": "bjacobi",
        "sub_pc_type": "ilu",
    }
    options.update(model_options)
    solver_obj.create_equations()
    # Set up gauges
    if store_station_time_series:
        for name, data in read_station_data().items():
            sta_x, sta_y = coord_system.to_xy(data["lon"], data["lat"])
            cb = TimeSeriesCallback2D(
                solver_obj,
                ["elev_2d"],
                sta_x,
                sta_y,
                name,
                append_to_log=False,
                start_time=data["start"],
                end_time=data["end"],
            )
            solver_obj.add_callback(cb)
    # Set boundary conditions
    solver_obj.bnd_functions["shallow_water"] = {
        100: {"un": Constant(0.0), "elev": Constant(0.0)},
    }
    # Set initial condition
    solver_obj.assign_initial_conditions(elev=elev_init)
    return solver_obj