This pattern could be replicated in any language, but is implemented here in Python. It is a tabular manipulation, requiring no GIS software, but relies on the id and idDown attributes to define routing connectivity.
The hydrographic network is known to the function as a array, networkArr, where each record has two attributes: 1) the flowline identifier, id, and 2) the identifier of the downstream flowline, idDown. The array is "typed", meaning that the columns have names (strings).
The function will add the starting currComid to the list of all upstream feature identifiers, comidList. If there are any flowlines whose idDown points to the currComid, then the function will call itself once for each inflowing neighbor using that neighbor's id as the currComid argument. If there is a flowline draining into one of the inflowing neighbors, then the function will call itself using that neighbor's id. This process will climb upslope according to the networkArr connectivity, until the networkArr indicates that no records have an idDown value that eqauls the current currComid (i.e., that the current flowline is a headwater). The final product of the function is the comidList. Note that the client should provide an empty list to the function on the first iteration.
Below is an example of some arcpy code for extracting the relevant information from a featureclass.