What is the purpose of a 1x1 convolutional layer?
First of all, Conv2D(16, kernel_size=(1,1)) applied to 28x28x1 will produce 28x28x16 as it doesn't change the number of dimensions. Secondly, I took a look at the paper you refer to, and haven't found that they apply 1x1 filters to an input with just one channel as in the example you provided. In case of just one channel in the input your reasoning is correct, but it cannot be transferred to input with several channels.

Conv2D by definition, assign weights between all channels of a selected area (in case of 1x1 kernel to a whole column of nodes) and corresponding column of nodes in the output grid. Then it applies these weights subsequently to all valid areas according to stride and padding parameters. Thus applying 1x1 kernel on 1-channel input doesn't make sense as all output channels will be proportional to each other and to the input channel as well. In case of multiple channels, as per definition Conv2D each node of the output will be a weighted (non-trivial) sum of a corresponding input column and won't be proportional to either channel of the input in most cases.