I am very new to python coding and have been struggling with this for a while I cannot see the reason for this error! 

Here is the code that is causing this error:

    import arcpy
    import numpy as np
    import math
    from arcpy.sa import *

    # Get the input arguements

    # Input DEM
    input_raster = arcpy.GetParameterAsText(0)

    # Output 
    output_raster = arcpy.GetParameterAsText(1)

    # Input Discharge
    mean_discharge = float(arcpy.GetParameterAsText(2))

    # Input number of timesteps
    input_timestep = int(arcpy.GetParameterAsText(3))

    #Other stuff to put outside while loop

    # The below text takes the input raster and calculates the bottom left corner
    extent_xmin_result = arcpy.GetRasterProperties_management(input_raster, "LEFT")
    extent_xmin = float(extent_xmin_result.getOutput(0))
    extent_ymin_result = arcpy.GetRasterProperties_management(input_raster, "BOTTOM")
    extent_ymin = float(extent_ymin_result.getOutput(0))
    
    # Turns the corner into a point
    bottom_left_corner = arcpy.Point(extent_xmin, extent_ymin)

    arcpy.overwriteOutput = True

    arcpy.env.workspace = r"C:\SedMod\SedMod\Workspace.gdb"
    arcpy.env.scratchWorkspace = r"C:\SedMod\SedMod\Workspace.gdb"

    actual_timestep = 0

while actual_timestep <= input_timestep:
    # Slope calculation

    # Set local variables
    out_measurement = "PERCENT_RISE" 

    # Check out the ArcGIS Spatial Analyst extension license
    arcpy.CheckOutExtension("Spatial")

    #Execute Slope
    out_slope = Slope(input_raster, out_measurement)

    # Convert to gradient
    slope_percentage = out_slope / 100

    #Save slope
    save_location = r"C:/SedMod/SedMod/Workspace.gdb/test" + str(actual_timestep)
    slope_percentage.save(save_location)

    #Execute Flow Direction

    output_flowdirection = FlowDirection(input_raster)

    #Execute Drainage Area

    # Method for describing the raster and then using the cellsize to determine the Drainage factor for calculating Drainage area
    desc_raster = arcpy.Describe(input_raster)
    raster_cell_width = desc_raster.MeanCellWidth # Cell width
    raster_cell_height = desc_raster.MeanCellHeight # Cell height
    spatial_reference = desc_raster.spatialReference # Spatial reference
    multi_factor = (raster_cell_height / 1000) * (raster_cell_width / 1000)
                
    # Run the flow direction
    output_flowaccumulation = FlowAccumulation(output_flowdirection)
    output_drainage_area = output_flowaccumulation * multi_factor

    #Execute Discharge Calculation
    
    mean_discharge_float = float(mean_discharge)

    #Calculating the discharge based on the average taken from the Avon catchment to bath (1605km2) and a mean flow of 20.762 m3/s giving a discharge per km2 of 0.01294 m3/s
    discharge_raster = output_drainage_area * mean_discharge_float

    #Execute StreamPower

    # Set variables for calculating stream power    
    density_of_water = 1000
    acceleration_gravity = 9.81
        
    # Stream power calculation = stream power, g is the density of water (1000 kg/m3), g is acceleration due to gravity (9.8 m/s2), discharge (m3/s) and channel slope.
    unit_width_stream_power = (density_of_water * acceleration_gravity * discharge_raster * slope_percentage) / raster_cell_width
    
    #Execute calculating Sediment Transport        
            
    # Set variables for calculating dimensionless stream power    
    pw = 1000 # Density of water
    ps = 2650 # Density of sediment
    g = 9.81 # Acceleration due to gravity 
    Di = 0.01 # Sediment size
                                                          
    # Dimensionless unit width stream power
    dimensionless_stream_power = unit_width_stream_power / ((16186.5) * math.sqrt((1.65 * 9.81 * (Di * Di *Di ))))
    
    # Additional variables to calculating sediment transport
    timestep_hour = 3600
            
    # Calculates dimensionless sediment transport       
    dimensionless_sediment_transport = Con(dimensionless_stream_power < 0.25, ((dimensionless_stream_power ** 6) * 100), (0.2 * (dimensionless_stream_power ** 1.5)))
                       
    # Calculates bed transport rate
    bed_material_transport = dimensionless_sediment_transport * ((16186.5) * math.sqrt((1.65 * 9.81 * (Di * Di *Di ))))            
           
    # Calculates it for the cell size in kg per second        
    bed_material_transport_kg = (bed_material_transport * raster_cell_width)  / 9.81
    
    # Converts it into volume m3
    bed_material_transport_volume = (bed_material_transport_kg / 2650) * timestep_hour
          

    
    # Convert the sediment transport and the flow direction rasters into Numpy arrays
    sediment_transport_np = arcpy.RasterToNumPyArray(bed_material_transport_volume, '#', '#', '#', -9999)
    flow_direction_np = arcpy.RasterToNumPyArray(output_flowdirection, '#', '#', '#', -9999)
    
    [rows,cols] = flow_direction_np.shape
    elevation_gain = np.zeros((rows,cols), np.float)
    elevation_gain_32 = np.zeros((rows,cols), np.float)
    elevation_gain_64 = np.zeros((rows,cols), np.float)
    elevation_gain_128 = np.zeros((rows,cols), np.float)
    elevation_gain_16 = np.zeros((rows,cols), np.float)
    elevation_gain_1 = np.zeros((rows,cols), np.float)
    elevation_gain_2 = np.zeros((rows,cols), np.float)
    elevation_gain_4 = np.zeros((rows,cols), np.float)
    elevation_gain_8 = np.zeros((rows,cols), np.float)
            
    # Main body for calculating elevation change        
    
    # Attempt 4
    
    # Move the sediment transport in the direction for 32 direction        
    for [i, j], flow in np.ndenumerate(flow_direction_np):
        try:
            if flow == 32:
                elevation_gain_32[i - 1, j - 1]  = sediment_transport_np[i, j]
        except IndexError:
                elevation_gain_32[i - 1, j - 1] = 0
                
    # Move the sediment transport in the direction for 64 direction
    for [i, j], flow in np.ndenumerate(flow_direction_np):
        try:
            if flow == 64:
                elevation_gain_64[i - 1, j]  = sediment_transport_np[i, j]
        except IndexError:
            elevation_gain_64[i, j] = 0
    
    # Move the sediment transport in the direction for 128 direction                  
    for [i, j], flow in np.ndenumerate(flow_direction_np):
        try:
            if flow == 128:
                elevation_gain_128[i - 1, j + 1]  = sediment_transport_np[i, j]
        except IndexError:
            elevation_gain_128[i, j] = 0                
    
    # Move the sediment transport in the direction for 16 direction
    for [i, j], flow in np.ndenumerate(flow_direction_np):
        try:
            if flow == 16:
                elevation_gain_16[i, j - 1]  = sediment_transport_np[i, j]
        except IndexError:
            elevation_gain_16[i, j] = 0
    
    # Move the sediment transport in the direction for 1 direction
    for [i, j], flow in np.ndenumerate(flow_direction_np):
        try:
            if flow == 1:
                elevation_gain_1[i, j + 1]  = sediment_transport_np[i, j]
        except IndexError:
            elevation_gain_1[i, j] = 0
    
    # Move the sediment transport in the direction for 2 direction
    for [i, j], flow in np.ndenumerate(flow_direction_np):
        try:
            if flow == 2:
                elevation_gain_2[i + 1, j + 1]  = sediment_transport_np[i, j]
        except IndexError:
            elevation_gain_2[i, j] = 0
    
    # Move the sediment transport in the direction for 4 direction
    for [i, j], flow in np.ndenumerate(flow_direction_np):
        try:
            if flow == 4:
                elevation_gain_4[i + 1, j]  = sediment_transport_np[i, j]
        except IndexError:
            elevation_gain_4[i, j] = 0
    
    # Move the sediment transport in the direction for 8 direction
    for [i, j], flow in np.ndenumerate(flow_direction_np):
        try:
            if flow == 8:
                elevation_gain_8[i + 1, j - 1]  = sediment_transport_np[i, j]
        except IndexError:
            elevation_gain_8[i, j - 1] = 0                                   
    
    # Covert numpy arrays back into ArcGIS raster format as I do not know how to successfully add together multidimensional numpy arrays                  
    elevation_gain_32_ras = arcpy.NumPyArrayToRaster(elevation_gain_32, bottom_left_corner, raster_cell_width, raster_cell_height, -9999)
    elevation_gain_64_ras = arcpy.NumPyArrayToRaster(elevation_gain_64, bottom_left_corner, raster_cell_width, raster_cell_height, -9999)
    elevation_gain_128_ras = arcpy.NumPyArrayToRaster(elevation_gain_128, bottom_left_corner, raster_cell_width, raster_cell_height, -9999)
    elevation_gain_16_ras = arcpy.NumPyArrayToRaster(elevation_gain_16, bottom_left_corner, raster_cell_width, raster_cell_height, -9999)
    elevation_gain_1_ras = arcpy.NumPyArrayToRaster(elevation_gain_1, bottom_left_corner, raster_cell_width, raster_cell_height, -9999)
    elevation_gain_2_ras = arcpy.NumPyArrayToRaster(elevation_gain_2, bottom_left_corner, raster_cell_width, raster_cell_height, -9999)
    elevation_gain_4_ras = arcpy.NumPyArrayToRaster(elevation_gain_4, bottom_left_corner, raster_cell_width, raster_cell_height, -9999)
    elevation_gain_8_ras = arcpy.NumPyArrayToRaster(elevation_gain_8, bottom_left_corner, raster_cell_width, raster_cell_height, -9999)
    
    # Calculate  volume in raster for whole catchment        
    volume_in = elevation_gain_32_ras + elevation_gain_64_ras + elevation_gain_128_ras + elevation_gain_16_ras + elevation_gain_1_ras + elevation_gain_2_ras + elevation_gain_4_ras + elevation_gain_8_ras
    
    #State what volume out is and the elevation           
        
    net_volume = volume_in - bed_material_transport_volume
    
    elevation_change = net_volume / (raster_cell_width * raster_cell_height)
    
    new_elevation = input_raster + elevation_change
    
    actual_timestep = actual_timestep + 1
    
    new_elevation.save(output_raster)   

The error I am getting is:

    Traceback (most recent call last):
      File "C:\SedMod\SedMod\Slope1.py", line 47, in <module>
        out_slope = Slope(input_raster, out_measurement, )
      File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\sa\Functions.py", line 5523, in Slope
        z_factor)
      File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\sa\Utils.py", line 47, in swapper
        result = wrapper(*args, **kwargs)
      File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\sa\Functions.py", line 5518, in wrapper
        z_factor)
      File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\geoprocessing\_base.py", line 498, in <lambda>
        return lambda *args: val(*gp_fixargs(args, True))
    RuntimeError: ERROR 010240: Could not save raster dataset to C:\SedMod\SedMod\Workspace.gdb\numpy_ras6 with output format FGDBR.
    
    Failed to execute (SedModv1).