This is workflow for single glacier:

arcpy.PolygonToLine_management("Glacier", "./OUTLINE.shp")
arcpy.InterpolateShape_3d("DEM", "OUTLINE", "./outline3d.shp")
arcpy.FeatureVerticesToPoints_management("outline3d", "./points3d.shp", "ALL")
arcpy.AddZInformation_3d("points3d", "Z")
arcpy.Sort_management("points3d", "./sorted.shp", "Z DESCENDING")
# select extreme points
arcpy.SelectLayerByAttribute_management("sorted", "NEW_SELECTION", '"FID" in ( 0, 772)')
arcpy.SplitLineAtPoint_management("OUTLINE", "sorted", "./threeOrMore.shp", "0.5 Meters")
# set environment extent to glacier
arcpy.gp.EucAllocation_sa("threeOrMore", "./alloc", "", "", "5", "FID", "", "")
arcpy.gp.FocalStatistics_sa("ALLOC", "./fstat", "Rectangle 3 3 CELL", "RANGE", "DATA")
arcpy.gp.Thin_sa("fStat", "./thinned", "ZERO", "NO_FILTER", "ROUND", "50")
arcpy.RasterToPolyline_conversion("thinned", "./Skeleton.shp", "ZERO", "0", "SIMPLIFY", "VALUE")
arcpy.Clip_analysis("Skeleton", "Glacier", "./clipped.shp")
arcpy.SelectLayerByLocation_management("clipped", "INTERSECT", "sorted", "10 Meters")
arcpy.Dissolve_management("clipped", "./c_line.shp", "", "", "SINGLE_PART")

It's enough to create a loop through individual features to achieve your goal. You still need to decide what to do with holes in polygons caused by nunataks.

This is workflow for single glacier:

arcpy.PolygonToLine_management("Glacier", "./OUTLINE.shp")
arcpy.InterpolateShape_3d("DEM", "OUTLINE", "./outline3d.shp")
arcpy.FeatureVerticesToPoints_management("outline3d", "./points3d.shp", "ALL")
arcpy.AddZInformation_3d("points3d", "Z")
arcpy.Sort_management("points3d", "./sorted.shp", "Z DESCENDING")
# select extreme points
arcpy.SelectLayerByAttribute_management("sorted", "NEW_SELECTION", '"FID" in ( 0, 772)')
arcpy.SplitLineAtPoint_management("OUTLINE", "sorted", "./threeOrMore.shp", "0.5 Meters")
# set environment extent to glacier
arcpy.gp.EucAllocation_sa("threeOrMore", "./alloc", "", "", "5", "FID", "", "")
arcpy.gp.FocalStatistics_sa("ALLOC", "./fstat", "Rectangle 3 3 CELL", "RANGE", "DATA")
arcpy.gp.Thin_sa("fStat", "./thinned", "ZERO", "NO_FILTER", "ROUND", "50")
arcpy.RasterToPolyline_conversion("thinned", "./Skeleton.shp", "ZERO", "0", "SIMPLIFY", "VALUE")
arcpy.Clip_analysis("Skeleton", "Glacier", "./clipped.shp")
arcpy.SelectLayerByLocation_management("clipped", "INTERSECT", "sorted", "10 Meters")
arcpy.Dissolve_management("clipped", "./c_line.shp", "", "", "SINGLE_PART")

It's enough to create a loop through individual features to achieve your goal. You still need to decide what to do with holes in polygons caused by nunataks.